Cox Regression Model Under Dependent Truncation
Lior Rennert Sharon X. Xie

Entanglement spectrum of mixed states
    Oded Zilberberg
    December 30, 2023
=====================================




Truncation is a statistical phenomenon that occurs in many time to event studies. For example, autopsy-confirmed studies of neurodegenerative diseases are subject to an inherent left and right truncation, also known as double truncation. When the goal is to study the effect of risk factors on survival, the standard Cox regression model cannot be used when the data is subject to truncation. Existing methods which adjust for both left and right truncation in the Cox regression model require independence between the survival times and truncation times, which may not be a reasonable assumption in practice. We propose an expectation-maximization algorithm to relax the independence assumption in the Cox regression model under left, right, or double truncation, to an assumption of conditional independence. The resulting regression coefficient estimators are consistent and asymptotically normal. We demonstrate through extensive simulations that the proposed estimators have little bias and, in most practical situations, have a lower mean-squared error compared to existing estimators. We implement our approach to assess the effect of occupation on survival in subjects with autopsy-confirmed Alzheimer's disease.




§ INTRODUCTION

Truncation is a statistical phenomenon that has been shown to occur in a wide range of applications, including survival analysis, epidemiology, economics, and astronomy. Individuals who are subject to truncation provide no information to the investigator. Left truncation occurs when data is only recorded for individuals whose event time exceeds a random time (i.e. left truncation time). Under left truncation, individuals with smaller event times are less likely to be observed, resulting in a study sample that is biased towards larger event times and risk factors associated with larger event times. Right truncation occurs when data is only recorded for individuals whose event time proceeds a random time (i.e. right truncation time). Under right truncation, individuals with larger event times are less likely to be observed, resulting in a study sample that is biased towards smaller event times and risk factors associated with smaller event times. When both left and right truncation are present, this is known as double truncation.

Double truncation is inherent in autopsy-confirmed studies of neurodegenerative diseases <cit.>. Left truncation occurs because individuals enter the study after the onset of the disease, and therefore those who succumb to the disease before they enter the study are unobserved. The right truncation occurs because individuals who live past the end of the study date do not receive a pathological diagnosis of the disease. Since these subjects cannot be definitively diagnosed with a particular disease, they are excluded from the autopsy-confirmed study sample and therefore provide no information to the investigator. This is contrary to censored individuals, who provide partial information about their survival time. We note, however, that right censoring is not possible in autopsy-confirmed studies, since any individual who has an autopsy performed will also have a known survival time. This truncation scheme is illustrated in Figure 1, where only individuals whose time of death falls between the study entry time and end of study time are observed. 

The aim of our data analysis is to get accurate estimates of the effect of risk factors on survival from disease symptom onset in subjects with autopsy-confirmed Alzheimer's disease (AD), the most common neurodegenerative disease. Because individuals with shorter survival times are less likely to enter the study, left truncation leads to a study sample that is biased towards larger survival times and risk factors associated with larger survival times. Similarly, individuals with longer survival times are more likely to live past the end of the study, and therefore right truncation leads to a study sample that is biased towards smaller survival times and risk factors associated with smaller survival times. If double truncation is not accounted for, then the regression coefficient estimators from the Cox regression model <cit.> will be biased.

Methods to handle double truncation have recently started gaining traction in the literature. In 2017, three methods were published to adjust the Cox model under double truncation <cit.>. The estimation procedure for all three methods rely on estimating the joint distribution of the left and right truncation times, which is used to compute the probability that a subject is observed (i.e. not truncated). These probabilities are then used as weights or offsets in the Cox model. However, the estimation of the truncation distribution relies on the assumption of independence between the observed survival and truncation times, which may not a reasonable assumption in practice. For example, according to the Alzheimer's association and discussions with our clinical investigators, factors such as lower age of symptom onset, depression, and stress are associated with delayed study entry. Since these factors are associated with survival, this induces a dependence between the left truncation times and survival times. As shown in the simulation studies in Section 3, the regression coefficient estimators from <cit.> are sensitive to violations of this independence assumption. Therefore, the existing literature is unable to address the unique challenges present in our study.  

In this paper, we propose a novel method to relax the assumption of independence between the observed survival and truncation times in the Cox proportional hazards model under left, right, and double truncation. Specifically, by conditioning on the observed truncation times, our method relaxes the independence assumption to an assumption of conditional independence. Treating the truncated survival times as missing, we introduce an expectation-maximization (EM) algorithm to estimate the regression coefficients and baseline hazard rates. This approach, which completely avoids the estimation of the truncation distribution, yields consistent and asymptotically normal regression coefficient estimators. We show through extensive simulation studies that our proposed estimators have little bias in small samples, while the estimators based on the methods introduced in <cit.>, and the standard model which ignores double truncation, can be heavily biased under violations of the independence assumption. We show that even if the independence assumption is satisfied, our proposed method performs as well as the existing approaches. We illustrate our method by analyzing the effect of cognitive reserve on survival in an autopsy-confirmed AD cohort.

Cognitive reserve (CR) is a widely used hypothetical construct intended to account for individual differences in cognitive decline and clinical manifestations of dementia among individuals with AD <cit.>. CR hypothesizes that individuals develop cognitive strategies and neural connections throughout their life times through experience such as occupation, education, and other forms of mental engagement <cit.>. This may modulate the effects of AD because of compensatory strategies obtained from a higher level of professional performance or a good education <cit.>. For example, CR may have a protective role in the brain and therefore lengthen survival from disease symptom onset <cit.>. 

Occupation, often used as a proxy for CR, has been shown to modulate survival in healthy aging and AD <cit.>. Several studies in healthy aging have examined the possible protective influence of higher occupational attainment on survival <cit.>. However other studies have shown that for individuals with AD, those with a higher occupational attainment had a higher mortality rate than those with a lower occupation attainment <cit.>. The caveat to previous studies assessing the effect of occupation on survival is that most consisted of populations with clinically diagnosed AD subjects, which can be unreliable <cit.>. Due to the inaccuracy of clinical diagnosis of AD, autopsy-confirmation is used for a definitive diagnosis <cit.>. Without an accurate diagnosis of AD, any estimates of factors affecting survival are not reliable. In this paper, we aim to get improved estimates of the effect of occupation on survival from an autopsy-confirmed AD sample, adjusting for both truncation and dependence.

The remainder of this paper is organized as follows. In Section 2 we introduce the proposed EM method, including the estimation procedure and the large sample properties of the resulting estimators. In Section 3, we conduct a simulation study to assess the finite sample performance of the proposed estimators under dependent truncation. In Section 4, we apply the proposed method to estimate the effect of occupation on survival in individuals with autopsy-confirmed AD. Discussion and concluding remarks are given in Section 5. Proofs of the large sample results are outlined in the Appendix.



§ METHODS
  
We first introduce notation and assumptions. Let T denote the survival time of interest (e.g. survival time from disease symptom onset), L denote the left truncation time (e.g. time from disease symptom onset to entry into the study), R denote the right truncation time (e.g. time from disease symptom onset to the end of study date), and  denote a p× 1 vector of covariates. Let N denote the size of the target population – the population that would have been observed had there been no truncation present in the study. Due to double truncation, we only observe (T_i,L_i,R_i,_i) for i=1,...,n≤ N individuals who live long enough to enter the study (i.e. T ≥ L) and do not live past the end of the study (i.e. T≤ R). Here we have denoted the population random variables from the target population without subscripts, and the sampling random variables from the observed sample with subscripts.

The proportional hazards model <cit.> is considered the standard regression model for analyzing traditional right-censored survival data. The model assumes that the covariate-specific hazard function is given by λ_(t) = λ(t)exp(β'Z), where  is a p× 1 regression parameter vector, and λ(t) is the baseline hazard function and is unspecified. When the survival data are subject to selection bias, Cox's partial likelihood approach <cit.> cannot be directly applied. This is because the observed data are not a representative sample of the target population, and therefore the observed, biased data do not follow the model that is assumed for the unbiased data from the target population. When the data is biased due to double truncation, the distribution of the observed survival time T_i is given by:

    P(T_i≤ t|_i) = P(T≤ t|_i,L≤ T≤ R) = P(T≤ t, L≤ T≤ R|_i)/P(L≤ T≤ R|_i)≠ P(T≤ t|_i),

which differs from the distribution of the survival time T from the target population. Therefore the resulting estimates of the regression coefficients based on data from the observed sample will be biased estimators of the regression coefficients from the target population. 





Under the assumption of independence between the survival and truncation times, <cit.> adjust for double truncation by estimating the probability that a subject with survival time T_i is observed, defined by π_i = P(L≤ T≤ R|T=T_i), i=1,...,n. These probabilities are then used as weights or offsets in the Cox regression model. For example, under double truncation and independence between the survival times and truncation times, Rennert and Xie (2017) consistently estimate the true p× 1 regression coefficient vector β_0 by , the solution to 

    U_w(β,π)    = ∑_i=1^n∫_0^τπ_i^-1{Z_i(t) - ∑_j=1^n π_j^-1Y_j(t)exp{β'Z_j(t)}Z_j(t)/∑_j=1^n π_j^-1Y_j(t)exp{β'Z_j(t)}}dN_i(t) = 0,

where π = (π_1,...,π_n), Y_i(t) = I(T_i ≥ t), N_i(t) = I(T_i ≤ t), and I is the indicator function. Here τ is the maximum of the observed event times. The standard Cox regression estimator <cit.> which ignores double truncation, , is the solution to U_w(β,1) = 0, where U_w(β,1) is the score equation from the standard Cox model. 

The caveat of the approaches which adjust for double truncation is that they require estimating the distribution of the truncation times, which is needed to obtain the estimator of the selection probabilities π. Existing methods to estimate the truncation distribution require independence between the survival and truncation times. When this independence assumption is violated, the estimator of the truncation distribution will be biased, and therefore the estimator of the selection probabilities π will be biased. Because the methods in <cit.> depend on π, the resulting regression coefficient estimators will also be biased. The severity of this bias is demonstrated by the simulation studies in the next section. 

When the survival times are conditionally independent of the truncation times given the covariate , the likelihood of the observed survival times, conditional on the truncation times and covariates, is given by

    L_n(,Λ) = ∏_i=1^nλ(T_i)exp(β'Z_i)exp{-Λ(T_i)exp(β'Z_i)}/α_i(,λ),
  
where α_i(,λ) = exp{-Λ(L_i)exp(β'Z_i)} - exp{-Λ(R_i)exp(β'Z_i)} and Λ(t) = ∫_0^tλ(u)du. That is, α_i(,λ) = P(L≤ T≤ R|_i,L_i,R_i;,λ) is the probability of observing a random subject from the target sample with covariate _i and truncation times L_i and R_i. Conditioning on the truncation times allows us to utilize the information in the covariates  to relax the assumption of independence to an assumption of conditional independence. Furthermore, this conditioning completely avoids the need to estimate the distribution of the truncation times.

The log-likelihood function, log L_n(,Λ), can be expressed as

    l_n(,Λ) = n^-1∑_i=1^n[∫_0^τ{logλ(t)+β'Z_i-Λ(t)exp(β'Z_i)}dN_i(t)-logα_i(,λ)].
  
Due to the difficulties of maximizing the log-likelihood (<ref>) over all absolutely continuous cumulative hazard functions, we allow the estimator of λ to be discrete. Because the maximum likelihood estimation (MLE) of β and λ may be computationally intractable if directly solving the score equations for (<ref>), we estimate _0 and λ using an EM algorithm. This has the advantage that its maximization step (M-step) only involves the complete-data likelihood. Based on the EM algorithm, we provide a convenient estimation approach to obtain estimators of the regression coefficients and baseline hazard function under left, right, or double truncation. This approach allows the survival and truncation times to be dependent through the covariate vector . Furthermore, it does not require the estimation of the truncation time distribution. The estimation approach given here can easily be implemented using standard software for the Cox regression model. 



 §.§ Proposed EM Algorithm

Motivated by the approach in <cit.>, who proposed EM algorithms for length-biased and right-censored data, Shen and Liu <cit.> proposed an EM algorithm to obtain pseudo MLEs of the regression coefficients from the Cox model under independent left and right truncation. They referred to their MLEs as pseudo because their proposed likelihood included the plug-in value of the estimator of the selection probabilities π. However, as the authors point out, the estimated selection probabilities will be biased if the truncation times depend on the covariates Z. Hence, the resulting pseudo MLEs of the regression coefficients from the Cox model will also be biased.

We propose an EM algorithm for obtaining the MLE of (,) based on (<ref>). This allows us to relax the assumption of independence required by the methods in <cit.> to an assumption of conditional independence by avoiding the estimation of the truncation distribution (and corresponding selection probabilities). Similar to the approaches of <cit.> and <cit.>, we let t_1 <  ... < t_d denote the ordered, distinct failure times for {T_1,...,T_n}. We develop the EM algorithm based on the discrete version of Λ, which we redefine as a step function only taking jumps at t_1,...,t_d. Specifically, we set Λ(t) = ∑_t_j≤ tλ_j, were λ_j is the positive jump at time t_j for j=1,...,n. 

Our observed data consists of O = {_1,...,_n}, where O_i ≡ (T_i,L_i,R_i,Z_i) for i=1,...,n. Let O^* = {T_ir^*;i=1,...,n,r=1,...,m_i} denote the truncated latent data, where T_ir^* is the missing survival time for a subject with truncation times (L_i,R_i) and covariate vector Z_i for i=1,...,n and r=1,...,m_i. For notational convenience, we set θ = (β,λ) and define the density of T at time t_j, given Z_i, as f_i(t_j;) = λ_jexp(β'Z_i)exp{-∑_s=1^jλ_sexp(β'Z_i)}, where = (λ_1,...,λ_d). Assuming the latent survival times T_ir^* take their values in {t_1,...,t_d}, the complete data log-likelihood is given by

    l_full(;,^*) = ∑_j=1^d∑_i=1^n[I(T_i=t_j) + ∑_r=1^m_iI(T_ir^* = t_j)]log f_i(t_j;β,λ)



To estimate the parameter , the EM algorithm begins by choosing an initial value for , say ^(0). In our setting, we can choose ^(0) = (_s,_s), which are the estimates from the standard Cox model. For k=0,1,2,..., the expectation step (E-step) consists of calculating the expected value of the complete data log likelihood function l_full(;,^*) with respect to the missing data T_ir^*, i=1,...,n,r=1,...,m_i, conditional on the observed data (T_i,L_i,R_i,_i), i=1,...,n, under the current estimate ^(k). That is, we compute:

    Q(θ;θ^(k)) = E_θ^(k)[l_full(;,^*)|]
 
In the maximization step (M-step), we choose ^(k+1) to maximize Q(θ;θ^(k)). That is, we set

    ^(k+1) = arg max_Q(θ;θ^(k))

The E- and M-steps are carried out again, but this time with ^(k) replaced by ^(k+1). The E- and M-steps are then alternated repeatedly until ^(k+1)-^(k)<ϵ for some prespecified error ϵ>0.

The EM algorithm described here is under the double truncation setting. If only left truncation is present, the algorithm is easily adjusted by setting R_i = ∞ for i=1,...,n. When only right truncation is present, we set L_i = -∞ for i=1,...,n. Note that when only left truncation is present, the standard Cox regression estimator can account for dependent left truncation by adjusting the risk set at a given time point to include all individuals who are alive and in the study at that time <cit.>. We denote this estimator by _s,l. We show through simulations in the next section that when only left truncation is present, our proposed estimator and _s,l yield nearly identical results. 



  §.§.§ E-step

At the k^th iteration, define θ^(k) =  (β^(k),λ^(k)). Then,

    Q(θ;θ^(k))     =  E_θ^(k)[∑_j=1^d∑_i=1^n[I(T_i=t_j) + ∑_r=1^m_iI(T_ir^* = t_j)]log f_i(t_j;)|O] 
        = ∑_j=1^d∑_i=1^n {I(T_i=t_j) + ∑_j=1^d∑_i=1^nE_θ^(k)[∑_r=1^m_iI(T_ir^* = t_j)|O_i]}log f_i(t_j;) 
        = ∑_j=1^d∑_i=1^n {I(T_i=t_j) + ∑_j=1^d∑_i=1^nE_m_i[m_i× E_[I(T_ir^* = t_j)|O_i]]}log f_i(t_j;),

where

    E_[I(T_ir^* = t_j)|O_i]     =  P_(T_ir^*=t_j|L_i,R_i,_i) = P_(T=t_j|L_i,R_i,_i,{T<L}∪{T>R}) 
        = P_(T=t_j,{T<L}∪{T>R}|L_i,R_i,_i)/P_({T<L}∪{T>R}|L_i,R_i,_i) = f_i(t_j;)×[I(t_j<L_i)+I(t_j>R_i)]/1-α_i().

Since m_i is the number of missing/truncated subjects with covariate values _i and truncation times L_i and R_i, m_i follows a geometric distribution with success rate α_i(). Therefore when =, E[m_i] = 1-α_i()/α_i(). 

The complete data log likelihood is then given by 

    Q(;) = {I(T_i=t_j)+I(t_j<L_i)+I(t_j>R_i)/α_i()f_i(t_j;)}log f_i(t_j;)






  §.§.§ M-step

Let =I(T_i=t_j)+I(t_j<L_i)+I(t_j>R_i)/α_i()f_i(t_j;). The complete data log likelihood can be written as 

    Q(;) = log f_i(t_j) = λ_j + '_i + ∑_s=1^jexp('_i)λ_s,

where =∑_i=1^n and =∑_j=1^d.

Treating  as constant, we set ∂ Q(;)/∂λ_j = 0 to get a closed form solution to λ_j as a function of :

    λ_j = /∑_s=j^d w_is^(k)exp('_i),    j=1,...,d.

Differentiating Q(;) with respect to  yields  

    ∂ Q(;)/∂ = _i + ∑_s=1^j_iexp('_i)λ_s.

Setting the equation above equal to 0 and inserting the equation for λ_j yields

    _i - ∑_s=1^dw_+s^(k){∑_j=s^d_iexp('_i)/∑_j=s^dexp('_i)} = 0.

The estimating equation (<ref>) can be solved by specifying the “weights” option in the “coxph” function in R. First, a weight vector of length nd must be created: w_nd^(k) = (w_11^(k),...,w_1d^(k),...,w_n1^(k),...,w_nd^(k)). The corresponding failure time data and covariate vectors are also created with length nd as follows: T_nd = (t_1,...,t_d,...,t_1,...,t_d) and _nd = (_1,...,_1,...,_n,...,_n). Letting Δ_nd be the identity vector of length nd, the solution to (<ref>), which we denote by ^(k+1), can be obtained with the following command: 

    coxph(Surv(T_nd,Δ_nd)∼Z_nd, weights = w_nd^(k), subset = which(w_nd^(k)>0)).

Plugging ^(k+1) into (<ref>) yields an updated estimator for , ^(k+1). We then set θ^(k+1) =  (β^(k+1),λ^(k+1)), and repeat the E- and M-steps. We continue to alternate between the E- and-M steps until θ^(k+1)  - θ^(k) < ϵ, for some prespecified error ϵ > 0. The MLE of the hazard ratio is then given by = β^(k+1). We denote the corresponding baseline hazard by _em = λ^(k+1), and the cumulative baseline hazard function by (t)=∑_t_j≤ tλ_j^(k+1). 

The EM algorithm presented here falls into the general scheme of the ECM algorithm, and therefore its convergence to the local maximizer is guaranteed by the same conditions required for convergence of the ECM algorithm <cit.>. The uniqueness of the resulting estimators are guaranteed by the regularity conditions in Appendix A.1. 




 §.§ Asymptotic Properties

In this section, we establish the strong consistency and asymptotic normality of the proposed EM estimators. Here we denote the proposed estimators by θ = (,), and denote the true regression coefficients and cumulative baseline hazard function _0 = (_0,Λ_0). The asymptotic properties of the proposed estimators refer to the situation when the total number of observed (non-truncated) subjects n→∞. The following theorems assume that the regularity assumptions in Appendix A.1 hold.

Theorem 1: Under the regularity assumptions given in the Appendix,  is consistent: As n→∞,  converges to _0, and (t) converges to Λ_0(t) almost surely and uniformly in t for t∈ [0,τ]. 

The existence and uniqueness of the MLE can be proved based on the log-likelihood function 

    l_n(,) = n^-1∑_i=1^n[   ∫_0^τ'_idN_i(t) + ∑_s=1^dlogλ_s× I(T_i=t_s) - ∑_t_s≤ t_iλ_sexp('_i) 
    
    -    log{exp(-exp('_i)∑_t_s<L_iλ_s) - exp(-exp('_i)∑_t_s≤ R_iλ_s)}].
  
Theorem 1 can then be proved by applying the classical Kullback-Leibler information approach as in <cit.>.

Theorem 2: Under the regularity assumptions given in the Appendix, √(n)(-_0) converges weakly to a tight mean-zero Gaussian process.

Theorem 2 is proved using the Z-theorem for infinite dimensional equations <cit.>. The proofs of Theorem 1 and Theorem 2 are outlined in Appendix A.2 and A.3, respectively. 

To obtain an estimate of the standard deviation of , we apply the simple bootstrap technique. In our setting, the bootstrap sample is obtained by drawing n independent vectors (T_j^b,L_j^b,R_j^b,_j^b), j=1,...,n, from the observed data vectors (T_i,L_i,R_i,_i), i=1,...,n, with replacement. These data vectors are then used to obtain an estimate of regression coefficients, denoted by ^(b). This process is repeated B times to obtain the B estimators ^(1),....,^(B). The estimate of the standard deviation of  is computed by taking the standard deviation of the ^(b), b=1,...,B. We denote this estimate by σ_. We show through simulation studies in the next section that the standard deviation of  is accurately estimated by σ_.



§ SIMULATIONS

In this section we examine the performance of the proposed estimator under dependent truncation. We compare our proposed estimator to the weighted estimator which adjust for double truncation but assumes independence between the survival and truncation times. We also compare the proposed estimator to the estimator from the standard Cox regression model. In all simulations, the survival times were generated from a proportional hazards model with hazard function λ(t)exp(β_1Z_1+β_2Z_2), and follow a Weibull distribution with scale parameter ν = 0.001 and shape parameter κ = 5. We set β_1=β_2 = 1, and generated the risk factors Z_1 and Z_2 from independent Unif[0,5] distributions. The truncation times were also simulated from Weibull distributions with scale parameter ν = 0.001 and shape parameter κ = 5. The left truncation times were generated from a proportional hazards model with hazard function λ_L(l)exp(β_L1Z_1+β_L2X), and the right truncation times were simulated from a proportional hazards model with hazard function λ_R(r)exp(β_R1Z_1+β_R2Y). Here X and Y were generated from independent Unif[0,5] distributions, with β_L2 = β_R2 = 1. To adjust the proportion of missing data due to left and right truncation, the truncation times were multiplied by constants c_l and c_r, respectively. A higher value of c_l induced a higher proportion of missing data due to left truncation, while a lower value of c_r induced a higher proportion of missing data due to right truncation. Because the survival, left, and right truncation times are all functions of Z_1 for β_1≠ 0, β_L1≠ 0, and β_R1≠ 0, they are dependent. However, the survival and truncation times are conditionally independent given Z_1. To adjust the degree of dependence between T and L, and T and R, we varied the regression coefficients β_L1 and β_R1, respectively.  

We conducted 1000 simulation repetitions with sample sizes of n = 100 and 250. To obtain n observations after truncation, we simulated N=n/1-q observations, where q is the proportion of truncated data. For each simulation, we estimated =(β_1,β_2), using the proposed EM estimator =(,), the weighted Cox regression estimator =(,), and the standard Cox regression estimator =(,). Of the estimators which adjust for double truncation under the independence assumption, we only focus on  from <cit.>, as previous simulations (not shown here) have concluded that this estimator and that in <cit.> are nearly identical, and both outperform the estimators in <cit.>. For each estimator, we calculated the estimated bias, observed sample standard deviations (SD), estimated standard errors (SE), and the average empirical coverage probability of the 95% confidence intervals (Cov). To compare the efficiency of the estimators which adjust for double truncation to the efficiency of the standard estimator, we calculated the relative mean-squared error (MSE) of β_j to , j=1,2. That is, we computed rMSE(β_j) = MSE(β_j)/MSE() for β_j = and β_j =. We used 200 bootstrap resamples to estimate the standard error of  and , j=1,2.

Table 1 shows the results of the simulations described above. In the first model, we set β_L1 = -1 to induce a negative dependence between the survival times and left truncation times, which resulted in a correlation of -0.35. In the second model, we set β_L1 = 0 to induce independence between the survival and left truncation times. In the third model, we set β_L1 = 1 to induce a positive dependence between the survival times and left truncation times, which resulted in a correlation of 0.35. Here c_l and c_r were chosen such that 25% of the survival times were left truncated and 25% of the survival times were right truncated, which resulted in q≈ 0.50. The parameter β_R1 was set to 1 in all models, which resulted in a correlation of 0.35 between the survival times and right truncation times.    

In all models, the proposed EM estimators  and  had little bias, while the standard estimators  and  were biased. The weighted estimator  was heavily biased in all models, while  was biased in the first set of models (ρ_LT = -0.35). The observed sample standard deviations of the proposed estimators were accurately estimated by the bootstrap technique, and the coverage probabilities of the proposed estimators were all close to the nominal level of 0.95. The coverage probabilities of  and  were well below the nominal level, as were the coverage probabilities for . Furthermore, the mean-squared errors of the proposed estimators were lower than those of the weighted and standard estimators in almost all settings, indicating that the proposed EM method is more efficient.

We further explored the bias and MSE of these estimators as a function of left and right truncation proportion (Figure 2). We set β_L1 = β_R1 = 1 and n=250, which corresponded to the setting of the last model in Table 1, inducing a positive dependency between the survival times and both left and right truncation times. The proposed estimators had little bias, regardless of truncation proportion. Even under mild truncation, the weighted estimator  of the regression coefficient corresponding to Z_1, which is correlated with the truncation times, was biased. This bias increased drastically as the proportion of right truncation increased. The bias was relatively small for both the proposed and weighted estimator of the regression coefficient corresponding to Z_2, which is uncorrelated with the truncation times. Both standard estimators  and  were heavily biased in this setting. The MSE of  was significantly lower than the MSE of  for j=1,2. Furthermore, in most cases the MSE of  was lower than the MSE of the standard estimator , i.e. rMSE() < 1 for j=1,2. 

In Figure 3, we compared the bias and MSE of these estimators under varying truncation proportions, when the assumption of independence holds (i.e. β_L1 = β_R1 = 0). The proposed EM estimators  and weighted estimators  had little bias, while the standard estimators  were biased for j=1,2. We also compared the rMSE of  and  to , j=1,2. As indicated by the bottom row of Figure 3, the proposed EM estimators had similar efficiency to the weighted estimators when the independence assumption holds. When the proportion of missing data due to left and right truncation were approximately equal, the standard estimator was more efficient than the proposed EM estimator and the weighted estimator. This is because the bias due to left truncation canceled out with the bias due to right truncation when these proportions were equal, which yielded a lower MSE. 

The standard Cox regression model can accommodate left truncation when the left truncation time is conditionally independent of the survival times given the observed risk factors. We compare the estimator from this model to our proposed estimator under dependent left truncation only in Appendix A.4. As shown in Figure 4, the results are nearly identical, and both estimators outperform the weighted estimators in this setting.




§ APPLICATION TO ALZHEIMER'S DISEASE

We illustrate our method by considering an autopsy-confirmed AD study conducted by the Center for Neurodegenerative Disease Research at the University of Pennsylvania. The target population for the research purposes of this study consists of all subjects with AD symptom onset before 2012 that met the study criteria and therefore would have been eligible to enter the center. Our observed sample contains all subjects who entered the center between 1995 and 2012, and had an autopsy performed before July 1, 2012. Thus one criterion for a subject to be included in our sample is that they did not succumb to AD before they entered the study, yielding left truncated data. In addition, our sample only contains subjects who had an autopsy-confirmed diagnosis of AD, and therefore we have no knowledge of subjects who live past the end of the study. Thus our data is also right truncated. Our data consists of n=91 subjects, all of whom have event times. The event time of interest is the survival time (T) from AD symptom onset. The left truncation time (L) is the time between the onset of AD symptoms and entry into the study (i.e. initial clinic visit). The right truncation time (R) is the time between the onset of AD symptoms and the end of the study, which is taken to be July 1, 2012. Due to double truncation, we only observe subjects with L≤ T≤ R.

We are interested in assessing the effect of occupation on survival in AD. Occupation is often used as a proxy for cognitive reserve (CR), which hypothesizes that individuals develop cognitive strategies and neural connections throughout their life times through experience such as occupation, education, and other forms of mental engagement <cit.>. A common hypothesis in the literature is that CR has protective role in the brain and modulates the effects of AD because of compensatory strategies obtained from a higher level of professional performance and therefore lengthens survival during the course of the disease <cit.>.  

However some studies have shown a higher mortality rate in AD individuals with higher occupational attainment <cit.>. This supports an alternative theory of CR; individuals with higher CR tolerate more pathology which delays the onset of the disease. Because higher age of AD symptom onset is associated with an increased risk of mortality, this would support the hypothesis that those with higher CR would have an increased risk of mortality. There are two caveats to the studies described above. The first is that these studies consisted of populations with clinically diagnosed AD subjects, which can be unreliable. The second caveat is that the statistical analyses were subject to confounding, since age of AD symptom onset was not recorded nor adjusted for. 

Here we are interested in obtaining improved estimates of the effect of occupation on survival from an autopsy-confirmed cohort of individuals with AD who have a known age of disease symptom onset. We use the highest occupational attainment for a given subject as a proxy for their CR. Primary occupation was classified and ranked based on the US census categories. In the following analyses, subjects who were classified as manager, business/government, and professional/technical workers were labeled as having high occupational attainment in our study. Subjects classified as unskilled/semiskilled, skilled trade or craft, and clerical/office workers were classified as having low occupational attainment. This classification is consistent with previous studies <cit.>. Age at AD symptom onset was estimated based on a family report at first contact with the individual. 

We first check the assumption of independence between the observed survival and truncation times using the conditional Kendall's tau proposed by Martin and Betensky <cit.>. The resulting p-value is 0.038, and therefore we reject this independence assumption. The corresponding Kendall's tau statistic is τ_K = (0.20,0.16), indicating positive dependence between the survival times and truncation times. The positive dependence between the left truncation times and survival times is clinically plausible because doctors often attribution the symptoms of early onset AD (onset of AD before 65 years of age) to other causes such as depression and stress, hence delaying the study entry time. Since younger age at onset is also associated with higher survival, this induces a positive dependence between the left truncation times and survival times.


Due to the dependence between the survival and truncation times, we apply the proposed method to estimate the effect of occupation on survival, adjusting for age at AD symptom onset and sex. Table 2 displays the results from the Cox regression model using the proposed EM estimators, weighted estimators, and the standard estimators. Using the proposed method, the estimated log hazard ratio for age at AD symptom onset is 0.029 (p-value = 0.016), indicating that AD individuals who have symptom onset one year later are roughly 3% more likely to die than subjects who have symptom onset a year earlier (e^0.029=1.03). The estimated effect of female is -0.636 (p-value = 0.023), indicating that males are almost twice as likely to die than females (e^0.636=1.89). These effects are nearly doubled using the weighted method which assumes independence, however the effects are not statistically significant (p-values = 0.117 and 0.088, respectively). 

High occupational attainment is associated with increased survival in all models. Under the proposed method, the effect of high occupational attainment on survival is -0.673 (p-value = 0.009), indicating that those with a low occupational attainment are approximately twice as likely to die than those with a high occupational attainment (e^0.673=1.96). This effect is attenuated under the weighted and standard methods, and neither method yielded statistically significant estimates (p-values are 0.186 and 0.158, respectively).  



§ DISCUSSION

We proposed a novel method which relaxes the independence assumption between the observed survival and truncation times in the Cox model under left, right, or double truncation to an assumption of conditional independence between the observed survival and truncation times. We obtained consistent and asymptotically normal estimators of the regression coefficients and baseline hazard function by maximizing the conditional likelihood of the observed survival times using an EM algorithm. The simulation studies confirmed that the proposed estimators had little bias in small samples, while the naïve estimators from the Cox models which ignore truncation or assume independence were biased. The existing methods which adjust for truncation but assume independence resulted in heavily biased estimators of the regression coefficients for risk factors of survival that were also correlated with the truncation times. Furthermore, the proposed estimators were more efficient than the naïve estimators in most of the simulation settings. 

We applied our proposed method to an autopsy-confirmed sample of individuals with Alzheimer's disease (AD). AD is a major neurodegenerative disease which currently affects 5.3 million people in the United States according to the Alzheimer's Association. In 2017 alone, AD and other dementias will have cost the nation an estimated $259 billion. Autopsy-confirmation is needed for a definitive diagnosis of AD, and a definitive diagnosis is necessary to accurately estimate the effect of potential risk factors associated with a given neurodegenerative disease. However, autopsy-confirmed samples of neurodegenerative diseases are subject to an inherent selection bias due to double truncation. Existing methods which adjust the Cox model in the presence of double truncation assume that the observed survival and truncation times are independent. This assumption may not be reasonable for studies of neurodegenerative diseases. In our data example, this independence assumption was rejected. Therefore, previous methods are not appropriate for our setting.

Given the severity of Alzheimer's disease on patients, their caregivers, and society, accurate estimation of the effects of risk factors on survival is crucial. One such factor, cognitive reserve (CR), is hypothesized to lengthen survival during the course of the disease. Using occupation as a proxy for CR, we estimated the effect of CR on survival in an autopsy-confirmed AD sample. Using our proposed method to adjust for both left and right truncation and dependence between the survival and truncation times, we found that a low occupational attainment was associated with shortened survival. Compared to existing methods, the estimated hazard ratios for occupation on survival were larger under our proposed method. This is consistent with many studies concluding that an individual's occupation may provide a protective effect and lengthen survival in AD. These findings suggest the importance of incorporating occupation in treatment trials and prognostic considerations in individuals with AD. 

A limitation of our proposed method is that it in its current form, it cannot properly handle time-varying covariates measured after study entry, such as cognitive test scores. This is a consequence of the estimation procedure, which uses an expectation-maximization algorithm to estimate the latent survival times conditional on the observed truncation times and risk factors. This leads to predicting survival times based on risk factors measured after death for those missing subjects whose survival time is less than their left truncation time, which may yield biased regression coefficient estimators. 

The proposed method has useful implications for observational studies. Double truncation has been shown to be present in a variety of studies, such as studies of clinically diagnosed Parkinson's disease <cit.>, childhood cancer <cit.>, astronomy data <cit.>, and studies based on registry data <cit.>. In fact, any data pulled from a disease registry will be subject to inherent right truncation, since data is only recorded for subjects who have the disease and are entered in the registry by the time the data is extracted <cit.>. In certain cases, the data will also be subject to left truncation <cit.>. In a similar fashion, studies which only include data from individuals whose event times fall within the time course of the study are subject to double truncation <cit.>. Therefore careful consideration of the study design must be taken into account when fitting the Cox proportional hazards model. Furthermore, the assumption of independence should always be tested, given the high sensitivity of existing methods to this assumption. For example, a quick application of a Kendall's conditional Tau test <cit.> revealed this independence assumption is violated in the AIDS data used in <cit.>. We therefore recommend using the proposed estimators in most practical settings, since they have little bias, and in most situations, have a lower mean-squared error compared to existing estimators under left, right, or double truncation, under a wide range of dependence structures.




§ APPENDIX

We first define some notation. Let _n= and Λ_n = Λ_em(·), and set =(,). The true parameter is denoted by _0 = (_0,Λ_0).

A.1 Regularity Assumptions


  * The true hazard function λ_0(·) is continuously differentiable, Λ_0(0) = 0, and Λ_0(τ) < ∞.

  * The true parameter vector _0 lies in a compact set 𝔹. The set 𝔸 contains all nondecreasing functions Λ satisfying regularity assumption 1.

  * EZ and Eexp(') are bounded, where ≡√(z_1^2+...+z_p^2).

  * The information matrix -∂^2E[l_n(,())]/∂^2|_=_0 is positive definite. Here () = _em is used to emphasize the dependence on .

  * If P(b' = c) = 1 for some constant c, then b=0. 

Assumptions 1 and 2 are required for stochastic approximation. Assumptions 3 and 4 are needed to establish the asymptotic properties of the regression parameter estimates from the Cox model <cit.>. Assumption 5 implies no covariate colinearity and thus ensures that the model is identifiable. 


A.2 Consistency: Proof of Theorem 1

Since each function of  in l_n(,) is concave or strictly concave, and the summation of concave functions is concave, the log-likelihood function l_n(,) is strictly concave in . Therefore we can find a unique maximizer (·,) of l_n(,) for each  in a compact set 𝔹. The existence of the NPMLE for (,) follows by compactness of 𝔹 for the likelihood l_n(,(·,)), which is continues in . Uniqueness is guaranteed by Assumption for in A.1 for large samples.  

Here we show that if  converges, it must converge to _0. As _n maximizes the log-likelihood given in (<ref>), , the empirical Kullback-Leibler distance - must be nonnegative. Suppose  converges to some ^* = (^*,Λ^*). Then by the strong law of large numbers (SLLN), - must converge to the negative Kullback-Leibler distance between P_^* and P__0. Here P_ is the probability measure under the parameter . Since the Kullback-Leibler distance and - are nonnegative, the Kullback-Leibler distance between P_^* and P__0 must be zero. Therefore P_^*=P__0 almost surely, and it then follows from model identifiability that ^* = _0. Therefore if  converges, it must converge to _0.

The technical details to show that  indeed converges are similar to those in <cit.>. The idea is to find a further convergent subsequence for any subsequence of , and then apply Helly's selection theorem. Here we provide only a sketch of the proof. The first step is to show that  stays bounded. By regularity assumption 3,  is in a compact set and is therefore bounded. To show  is bounded, we make use of the fact that the empirical Kullback-Leibler distance l_n() - l_n() is always non-negative for each  in the parameter set. Using the approach of Murphy (1994), it can be shown that if  does indeed diverge to ±∞, then it is possible to construct some sequence _n such that l_n() - l_n() eventually becomes negative infinity, which contradicts the nonnegativity of the empirical Kullback-Leibler distance. 

Since  stays bounded, we can apply Helly's selection principal to find a further convergent subsequence _n_k = (_n_k,Λ_n_k) for any subsequence of  indexed by {1,...,n}. By the classical Kullback-Leibler information approach, and the SLLN, _n_k must converge to _0. It then follows from Helly's selection theorem that the entire sequence (,(t)) must converge to (_0,Λ_0(t)) for every t∈[0,τ], where τ=t_d is the maximum of the observed event times. Since Λ_0(·) is assumed to be monotone and continuous, the convergence of (t) is uniformly in t∈[0,τ]. Because the proof is carried out for a fixed ω in the underlying probability space Ω, where the SLLN is applied countably many times, the convergence here is also almost surely a true convergence. 


A.3 Asymptotic Normality: Proof of Theorem 2

Here we outline the proof for the weak convergence of _n, which follows the proof for weak convergence in <cit.>. The proof consists of the application of empirical process theory and the Z-theorem for infinite dimensional estimating equations <cit.>. 

Denote the score equation for  by U_1n() = ∂/∂. To obtain the score equation Λ(·), we define the submodel dΛ_ϵ = (1+ϵ h)dΛ, where h is a bounded and integrable function. Setting h(·)=I(·≤ t), the score equation for Λ is given by
U_2n(t,) = ∂ l_n(,Λ_ϵ)/∂ϵ|_ϵ=0. 

We denote the vector of the score functions by U_n(·,) = [U_1n() , U_2n(t,)]. The expectation E_0 under the true value _0 is given by U_0(·,) = [U_10() , U_20(t,)], where U_10() = E_0[U_1n()] and U_20(t,) = E_0[U_2n(·,)].

By the definition of the MLE, U_n(·,_n)=0. Since U_0(·,_0)=0, we can show that |√(n){U_0(·,_n)-U_n(·,_n)} -√(n){U_n(·,_0)-U_0(·,_0)}|=o_p(1). The estimating equation evaluated at _0, √(n)U_n(·,_0) = √(n){U_n(·,_0)-U_0(·,_0)}, is a sum of iid terms. We can therefore use empirical process theory to show that √(n)U_n(·,_0) converges weakly to 𝕎=(𝕎_1,𝕎_2), where 𝕎_1 is a Gaussian random vector and 𝕎_2 is a tight Gaussian process. The covariance matrix for 𝕎_1 is given by Σ_11 = E_0{U_1n(_0)^⊗ 2}, and the covariance between 𝕎_2(s) and 𝕎_2(t) is given Σ_22(s,t) = E_0{U_2n(s,_0)U_2n(t,_0)}.   

Applying the Z-theorem for the infinite dimensional estimating equations (theorem 3.3.1 in <cit.>), we have that under the regularity conditions in A.1, √(n)(_n-_0) converges weakly to a tight mean-zero Gaussian process -U̇__0^-1(𝕎). 

Here U̇__0 is the Fréchet derivative of the map   U_0(·,) evaluated at _0. Using arguments similar to Appendix A.5 in <cit.>, we can show U_0(·,) is Fréchet differentiable and the Fréchet derivative, U̇__0, is continuously invertible. By definition of the Fréchet derivative, we have that U̇__0{√(n)(_n-_0)} = -√(n){U_n(·,_0)-U_0(·,_0)} + o_p(1). This completes the proof. 


A.4 Simulations under dependent left truncation

Under dependent left truncation only, we compared our proposed method to the weighted method and the standard method which accounts for dependent left truncation. To adjust the correlation between the left truncation times and survival times, we varied the parameter β_L1 between -1 and 1. In this setting, a value of β_L1=0, which yields a correlation of 0, indicates independence between the left truncation times and survival times. We denote the standard regression coefficient estimator which adjusts for dependent left truncation as _s,l = (β_s,l,1,β_s,l,2). As shown in Figure 4,  and β_s,l,j had little bias, while the weighted estimators  were biased for j=1,2. As indicated by the bottom row of Figure 4, the proposed EM estimators had similar efficiency to the standard estimators which accounted for dependent left truncation, and both estimators were more efficient than the weighted estimators. 



§ ACKNOWLEDGEMENTS

We would like to thank Dr. Murray Grossman for his contribution to the clinical aspect
of this paper. Mr. Rennert received support from NIH National Institute of Mental Health
grant T32MH065218 and Dr. Xie from NIH grant R01-NS102324, AG10124, AG17586, and
NS053488.








 

 


 












