Semiparametric Estimation for the Transformation Model with Length-Biased Data and Covariate Measurement Error
Li-Pang Chen








Semiparametric Estimation for the Transformation Model with Length-Biased Data and Covariate Measurement Error
    
Li-Pang ChenaCorrespondence author: Li-Pang Chen. 
Email: L358CHEN@uwaterloo.ca
aDepartment of Statistics and Actuarial Science, 
University of Waterloo 

200 University Ave W, Waterloo, ON N2L 3G1

    Received; accepted
===========================================================================================================================================================================================================




Analysis of survival data with biased samples caused by left-truncation or length-biased sampling has received extensive interest. Many inference methods have been developed for various survival models. These methods, however, break down when survival data are typically error-contaminated. Although error-prone survival data commonly arise in practice, little work has been available in the literature for handling length-biased data with measurement error. In survival analysis, the transformation model is one of the frequently used models. However, methods of analyzing the transformation model with those complex features have not been fully explored. In this paper, we study this important problem and develop a valid inference method under the transformation model. We establish asymptotic results for the proposed estimators. The proposed method enjoys appealing features in that there is no need to specify the distribution of the covariates and the increasing function in the transformation model. Numerical studies are reported to assess the performance of the proposed method.



Length-based sampling; measurement error; prevalent cohort; survival analysis; transformation model.






§ INTRODUCTION
 

In the study of a disease history, the time from the onset of an initiating event to the focused disease event (or failure) is usually the interest in the epidemiological and biomedical researches. One of the most attractive data comes from the prevalent sampling design, in which individuals only experience the initiating event but not the failure event before the recruiting time. Under this sampling scheme, individuals might not be observed because they experience the failure event before the recruiting time. Such a phenomenon caused by the delayed entry is called left-truncation and tends to produce a biased sample. Meanwhile, individuals who are recruited in the study may drop out or may not experience the failure event at the end of the study. It is called right-censoring in the dataset. To be specific, let ξ be the calendar time of the recruitment and let u and r denote the calendar time of the initiating event (or the disease incidence) and the failure event, respectively, where u<r, and u < ξ < r.  Let T^* = r-u be the failure time, and let A^* = ξ -u denote the truncation time. Assume that A^∗ and T^∗ are independent. Let X^* be the associated covariate of dimension p × 1.  Let f(t) and S(t) be the density function and the survivor function of the failure time T^*, respectively. Let A,T, and X represent, respectively, the truncation time, the failure time, and the covariates for those subjects who are recruited in the study. That is, (A,T,X) has the same joint distribution as (A^∗, T^∗, X^∗) given T^∗≥ A^∗. If T^∗ < A^∗, then such an individual is not included in the study to contribute any information.
 
 Specifically, in the case of stable disease, the occurrence of disease onset follows the stationary Poisson process. To see this, we first impose two assumptions (e.g., Huang et al. 2012), including (a) the variables (T^∗, X^∗) are independent of the disease incidence u, and (b) disease incidence occurs over calendar time at a constant rate. Based on these two assumptions, the joint density function of (A^∗, T^∗) given X^∗ = x is given by (Lancaster 1990, Chapter 3)
 
    f(t|x)/μ(x) I(t > a > 0),

 where μ(x) = ∫_0^∞α f(α | x) d α. Moreover, by (<ref>), the failure time T^∗ given X^∗ = x has a length-biased density function tf(t|x)/μ(x).
  It implies that the truncation time A^∗ follows the uniform distribution, and the survival time in the prevalent cohort has a length-biased sampling distribution since the probability of the survival time is proportional to the length of survival time. 
 
  In addition, we let C denote the censoring time for a recruited subject. Let Y = min{ T, A+C } be the observed survival time and let δ = I(T ≤ A+C) be the indicator of a failure event. Figure <ref> gives an illustration of the relationship among those variables.


In the past literature, several models for the analysis of  biased samples were developed. To name a few, Huang and Qin (2013) and Chen (2018+) developed the estimation procedures for the additive hazards model for general left-truncated data where the distribution of the truncation time is unspecified. Qin and Shen (2010) and Huang et al. (2012) studied the Cox model for length-biased sampling. In addition, the transformation model is also a widely used model in survival analysis, which is formulated by

    H( T^∗) = -X^∗^⊤β + ϵ,

where H(·) is an unknown increasing function, ϵ is a random variable with a known distribution, and β is the p × 1 vector of parameters of main interest. The transformation model gives a broad class of some frequently used models in survival analysis. Be more specific, when ϵ has an extreme value distribution, then T^∗ follows the proportional hazards (PH) model; whereas when ϵ has a logistic distribution, then T^∗ follows the proportional odds (PO) model. When biased samples occur, Shen et al. (2009) and Wang and Wang (2015) proposed valid methods to estimate the parameter of main interest based on length-biased sampling. In this paper, we mainly focus on the development of estimation method for the transformation model.

On the other hand, the other important feature in survival analysis is measurement error in covariates. Assume that X is an unobserved covariate and W is a surrogate version of X. The classical additive model is frequently used, which is given by

    W = X + η,

where η follows the normal distribution with mean zero and covariance matrix Σ_η. Assume that X and η are independent and Σ_η is known for ease of discussion.

When the absence of left-truncation or length-biased sampling, there are many contributions on correcting the mismeasurement. For example, Nakamura (1992) proposed an approximate corrected partial likelihood method for the Cox model. Yi and Lawless (2007) developed the conditional expectation approach to correct the error effect for the Cox model. More discussions can be found in Yi (2017, Chapter 3).

However, little work has been available when those two features happen simultaneously. In the past literature, Chen and Yi (2018) proposed the corrected pseudo-likelihood estimation to estimate the parameter for the Cox model subject to left-truncation and right-censoring data and covariate measurement error. However, other survival models, such as the transformation model, have not been fully explored when those two complex features occur in the dataset. In this paper, our main purpose is to develop a valid estimation procedure to estimate the parameter for the transformation model with length-biased sampling and measurement error in covariate.

Our work is partially motivated by the Worcester Heart Attack Study (WHAS500) data (Hosmer et al. 2008) which involves biased and incomplete data. Basically, three types of time are recorded:  time of the hospital admission, time of the hospital discharge, and time of the last follow-up (which is either death or censoring time). The total follow-up length is defined as the time gap between the hospital admission and the last follow-up, and the hospital stay time is defined as the time length between the hospital admission and the hospital discharge. Data can only be collected for those individuals whose total follow-up length is larger than the hospital stay time, creating biased sample (e.g., Kalbfleisch and Prentice 2002, Section 1.3; Lawless 2003, Section 2.4).  It is interesting to study how the risk factors are associated with the survival times after the patients are discharged from the hospital. In addition, `blood pressure' and `body mass index' are variables in the data set, and they may be measured with error. To conduct sensible analyses, it is imperative to account for possible measurement error effects that are induced from error-prone covariates.  

The remainder is organized as follows. In Section <ref>, we first present the proposed method to correct the error effect and derive the estimator. After that, we develop the theoretical result for the proposed method. Numerical results, including simulation study and real data analysis, are provided in Sections <ref> and <ref>, respectively. Finally, we conclude the paper with discussions in Section <ref>.







§ MAIN RESULTS
 
In this section, we propose a method to deal with the measurement error in covariates, and then derive the estimators of β and H(·). The theoretical results of the estimators are also established in this section.



 §.§ Construction of The Estimating Function
 
In this section, we review a method proposed by Wang and Wang (2015) and assume that the covariate X_i is collected without mismeasurement. 

 Suppose we have a sample of n subjects where for i=1, ⋯, n, (Y_i, A_i, δ_i, X_i, W_i) has the same distribution as ( Y,A,δ,X ,W). Let R_i(t) = I( A_i ≤ t ≤ Y_i ) and N_i(t) = I ( Y_i ≤ t, δ_i = 1 ) for i = 1,⋯,n. Define 

    M_i(t) = N_i(t) - ∫_0^t R_i(u) δ_i w(u)/w(Y_i)λ(u|X_i)

for i = 1,⋯,n, where λ(y|X_i) is the conditional hazard function of Y_i given X_i, w(t) = ∫_0^t S_C(u)du, and S_C(t) is the survivor function of the variable C. It can be shown that (<ref>) is a square-integrable martingale (e.g., Anderson et al. 1993). Define r(t,Y_i,δ_i) = δ_i w(t)/w(Y_i), and let r(t,Y_i,δ_i) = δ_i w(t)/w(Y_i) denote the estimator of r(t,Y_i,δ_i), where w(t) = ∫_0^t S_C(u)du and S_C(·) is the Kaplan-Meier estimator. Therefore, based on (<ref>) and the counting process theory (e.g., Anderson et al. 1993), we can construct the following estimating equations

    ∑_i=1^n dM_i(t) = 0, 
           ∑_i=1^n X_i dM_i(t) = 0,

or equivalently,

    ∑_i=1^n [ dN_i(t) -  R_i(t) r(t,Y_i,δ_i) dΛ{X_i^⊤β + H(t)}] = 0,

and

    ∑_i=1^n X_i [ dN_i(t) -  R_i(t) r(t,Y_i,δ_i) dΛ{X_i^⊤β + H(t)}] = 0,

where Λ(·) is the cumulative hazard function of ϵ. When β is fixed, let H_X(t;β) denote the unique solution of (<ref>). Next, replacing H(·) in (<ref>) by H_X(t;β) gives

    U_X(β) ≜∑_i=1^n X_i [ dN_i(t) -  R_i(t) r(t,Y_i,δ_i) dΛ{X_i^⊤β + H_X(t;β)}] = 0.

Consequently, let β_true denote the solution of U_X(β)=0 and it is the estimator of β based on the covariate X_i. Moreover, H_X(t;β_true) is the estimator of H(t).

Based on the true covariates X_i, we have the following proposition which is given in Theorem 1 of Wang and Wang (2015).

Let β_0 be the true parameter. Under regularity conditions in  Wang and Wang (2015), β_true is a consistent estimator and √(n)( β_true - β_0 ) converges weakly to a normal distribution with mean zero and variance-covariance matrix Σ_∗^-1Σ^∗ (Σ_∗^-1)^⊤, where Σ_∗ and Σ^∗ can be found in  Wang and Wang (2015).


However, as described in Section <ref>, we only observe W_i instead of X_i. Therefore, replacing X_i by W_i in (<ref>) gives the naive estimating function

    U_naive(β) ≜∑_i=1^n W_i [ dN_i(t) -  R_i(t) r(t,Y_i,δ_i) dΛ{W_i^⊤β + H_W(t;β)}],

where H_W(t;β) is the solution of (<ref>) replacing X_i by W_i.
We let β_naive denote the naive estimator which is a solution of U_naive(β) = 0.




 §.§ Corrected Estimating Function
 
It is well-known that the naive estimator β_naive incurs the tremendous bias (e.g., Carroll et al. 2006; Yi 2017). To correct the mismeasurement and reduce the bias of the estimator, we employ the simulation-extrapolation (SIMEX) method (Cook and Stefaski 1994). The proposed procedure is the following three stages: 


  Stage 1Simulation

Let B be a given positive integer and let 𝒵 = {ζ_0,ζ_1,⋯,ζ_M } be a sequence of pre-specified values with 0 = ζ_0<ζ_1<⋯<ζ_M. where M is a positive integer, and ζ_M is pre-specified positive number such as ζ_M = 1.

For a given subject i with i=1,⋯,n and b = 1,⋯,B, we generate η_b,i from N(0,Σ_η). Then for observed vector of covariates W_i, we define W_i(b,ζ) as

    W_i(b,ζ) = W_i + √(ζ)η_b,i

for every ζ∈𝒵. Therefore, the conditional distribution of W_i(b,ζ) given X_i is N( X_i, (1+ζ) Σ_η).


  Stage 2 Estimation

In this stage, replacing X_i in (<ref>) and (<ref>) by W_i(b,ζ) defined in  (<ref>) gives new estimating equations

    ∑_i=1^n [ dN_i(t) -  R_i(t) r(t,Y_i,δ_i) dΛ{W_i^⊤(b,ζ) β + H(t)}] = 0,

and

    ∑_i=1^n W_i(b,ζ) [ dN_i(t) -  R_i(t) r(t,Y_i,δ_i) dΛ{W_i^⊤(b,ζ) β + H(t)}] = 0.

Similar to the procedure in Section <ref>, we first fix β and let H(t;b,ζ,β) denote the solution of (<ref>). After that, replacing H(·) in (<ref>) by H(t;b,ζ,β) gives the estimating equation 

    U_SIMEX(β) 
       ≜   1/n∑_i=1^n W_i(b,ζ) [ dN_i(t) -  R_i(t) r(t,Y_i,δ_i) dΛ{W_i^⊤(b,ζ) β + H(t;b,ζ,β)}] 
       =    0.

for every b = 1,⋯,B and ζ∈𝒵. Let β(b,ζ) denote the solution of U_SIMEX(β) = 0. Moreover, we define

    β(ζ) = 1/B∑_b=1^B β(b,ζ).



  Stage 3 Extrapolation

By (<ref>), we have a sequence {( ζ, β(ζ)) : ζ∈𝒵}. Then we fit a regression model to the sequence

    β(𝒵) = φ(𝒵, Γ) + ϱ,

where φ(·) is the user-specific regression function, Γ is the associated parameter, and ϱ is the noise term. The parameter Γ can be estimated by the least square method, and we let Γ denote the resulting estimate of Γ.

Finally, we calculate the predicted value

    β_SIMEX = φ(-1, Γ)

and take β_SIMEX as the SIMEX estimator of β.


  Stage 4 Estimation of H(·) 

Furthermore, we can also derive the estimator of the unknown function H(·). To do this, we first replace β by β_SIMEX in H(t;b,ζ,β), which gives H(t;b,ζ,β_SIMEX). For every t and ζ, taking average with respect to b gives H(t;ζ,β_SIMEX) = 1/B∑_b=1^B H(t;b,ζ,β_SIMEX). Finally, similar to Stage 3 above, fitting a regression model and taking ζ = -1 as a predicted value yields a final estimator H(t;β_SIMEX), also denoted as H_SIMEX(t).






 §.§ Theoretical Results
 
In this section, we present the consistency and the asymptotic distribution of the proposed estimators β_SIMEX and H_SIMEX(·), and the proofs are available in Appendix <ref>.

 
Let β_0 be the true parameter and let H_0(t) denote the true increasing function. Under regularity conditions in Appendix <ref>, estimators β_SIMEX and H_SIMEX(·) have the following properties:


  (1) β_SIMEXp⟶β_0  as  n →∞;

  (2) H_SIMEX(t) p⟶ H_0(t)  as  n →∞;

  (3) √(n)(β_SIMEX - β_0 ) d⟶ N(0, {∂φ/∂Γ(-1, Γ) }𝒬{∂φ/∂Γ(-1, Γ) }^⊤)  as  n →∞;

  (4) √(n){H_SIMEX(t) - H_0(t) } converges to the Gaussian process with mean zero and covariance function E{ℋ_i(t) ℋ_i(s) },

where the exact formulations of 𝒬 and ℋ_i(t) are placed in Appendix <ref>.



§ NUMERICAL STUDY
 



 §.§ Simulation Setup

 We examine the setting where ϵ is generated from the extreme value distribution and the logistic distribution, and the truncation time A^∗ is generated from the uniform distribution U(0,1). Let β = (β_1, β_2)^⊤ denote a two-dimensional vector of parameters, and let β_0 = (β_10, β_20)^⊤ be the vector of true parameters where we set β_0 = ( 1,1 )^⊤. We consider a scenario where X^∗ = (X_1^∗, X_2^∗)^⊤ are generated from a bivariate normal distribution with mean zero and variance-covariance matrix Σ, which is set as ( [   4 0.7; 0.7   3 ]). 
 Given ϵ, X^∗ and β_0, the failure time T^∗ is generated from the model: 

    log T^∗ = - ( X_1^∗β_10 + X_2^∗β_20) + ϵ.

Based on our two settings of ϵ, the failure time T^∗ follows the PH model and the PO model, respectively.
  Therefore, the observed data (A,T,X) is collected from (A^∗,T^∗,X^∗) by conditioning on that T^∗≥ A^∗. We repeatedly generate data these steps we obtain a sample of a required size n=200. For the measurement error process, we consider model (<ref>)  with error η∼ N ( 0, Σ_η), where  Σ_η is a diagonal matrix of σ_η which is taken as 0.01, 0.5, and 0.75, respectively.
  
We consider two censoring rates, say 25% and 50%, and let the  censoring time C be generated from the uniform distribution U(0,c), where c is determined by a given censoring rate.
Consequently, Y and δ are determined by Y = min{ T, A+C } and δ = I ( T ≤ A + C ). In implementing the proposed method, we set B = 500 and partition the interval [0,2] into subintervals with width 0.25, and let the resulting cutpoints be the values of ζ. We take the regression function φ(·) to be the quadratic polynomial function, which is a widely used function in many cases (e.g., Cook and Stefaski 1994; Carroll et al. 2006). Finally, 1000 simulations are run for each parameter setting. 



 §.§ Simulation Results


We mainly examine three estimators which are discussed in Sections <ref> and <ref>, and they are labeled by Wang-Wang (β_true), Naive (β_naive), and Chen (β_SIMEX). We report the biases of estimates (Bias), the empirical variances (Var), the mean squared errors (MSE), and the coverage probabilities (CP) of those three estimators under the measurement error model (<ref>). The results are reported in Table <ref>. 

First, the censoring rate and measurement degree have noticeable impact on each estimation methods.
As expected, biases and variance estimates increase as the censoring rate increases. 
When the measurement degree increases, biases of both β_naive and β_SIMEX are increasing, and the impact of the measurement error degrees seems more obvious on the naive estimator β_naive.

Within a setting with a given censoring rate and a measurement error degree, the naive method and the proposed method perform differently. When measurement error occurs, the performance of the proposed method is better than the naive method. The naive method produces considerable finite sample biases with coverage rates of 95% confidence intervals significantly departing from the nominal level. The proposed method outputs satisfactory estimate with small finite sample biases and reasonable coverage rates of 95% confidence intervals. Compared to the variance estimates produced by the naive approach, the proposed method which accounts for measurement error effects yield larger variance estimates, and this is the price paid to remove biases in point estimators. This phenomenon is typical in the literature of measurement error models. However, mean squared errors produced by the proposed method tends to be a lot smaller than those obtained from the naive method. Finally, we also present the numerical results of the estimator β_true with the true covariate X. In general, β_true gives the smallest bias and it is an efficient estimator with small variance. This result also verifies the validity of method proposed by Wang and Wang (2015).


 




 


§ DATA ANALYSIS
 
 
In this section, we apply the proposed method to analyze the data arising from the Worcester Heart Attack Study (WHAS500), which are described in Section <ref>. 

Discussed by Hosmer et al. (2008), a survival time was defined as the time since a subject  was admitted to the hospital. We are interested in studying survival times of patients who were discharged alive from the hospital. Hence, a selection criterion was imposed that only those subjects who were discharged alive were eligible to be included in the analysis. That is, individuals were not enrolled in the analysis if they died before discharging from the hospital, hence biased sample occurs. With such a criterion, a sample of size 461 was available.
In this data set, the censoring rate is 61.8%. 
To be more specific, the total length of follow-up (lenfol) is the survival time (i.e., Y_i = min(T_i,C_i )), the length of hospital stay (los) is the truncation time (i.e., A_i), and the vital status at last follow-up (fstat) is δ_i. In our analysis, the covariates include the body mass index (BMI) and the blood pressure (BP) of a patient. Suppose that BMI and BP are subject to measurement error, we let W = ( W_1, W_2 )^⊤ where W_1 and W_2 denote BMI and BP, respectively, and consider the measurement error model (<ref>). Let β = ( β_1,β_2 )^⊤, where β_k is the associated parameter with W_k for k=1,2. We fit the regression model (<ref>) with ϵ being the extreme value distribution and the logistic distribution, respectively.

In this data set, there is no additional data source, such as a validation subsample or replicated measurements which is often required to describe the measurement error process (e.g., Carroll et al. 2006; Yi 2017). To get around this and understand the impact of measurement error on estimation, we carry out sensitivity analyses.  Specifically, let Σ be the sample covariance matrix of W, and for sensitivity analyses we consider Σ_η = Σ + Σ_e to be the covariance matrix for the measurement error model (<ref>), where Σ_e is the diagonal matrix with diagonal elements being a common value σ_e ∈ [0,1].



In Table <ref>, we report the point estimates (EST), the standard errors and p-values for the estimators β_SIMEX for the cases with σ_e = 0.15, 0.50 and 0.75, respectively, corresponding to minor, moderate and large measurement error, and report the results of β_naive by directly using W in the analysis. All the point estimates produced by our proposed method are fairly close as observed, and the results are fairly stable regardless of the degree of measurement error. Both our proposed method and the naive estimator suggest that two covariates are significant regardless of the specification of the distribution of ϵ. 


  
 


§ DISCUSSION
 

In this article, we focus the discussion on the transformation model with length-biased sampling and develop a valid method to correct the covariate measurement error and derive an efficient estimator. In this article, we also establish the large sample properties, and the numerical results guarantee that our proposed method outperforms. Although we only focus on the simple structure of the measurement error model, our method can easily be extended to complex measurement error models or additional information, such as repeated measurement or validation data. In addition, there are still many challenges in this topic, such as analysis without specifying the distribution of the truncation time A^∗ or the discussion of time-dependent covariates with mismeasurement. These topics are also our researches in the future.








§ REGULARITY CONDITIONS
 


  (C1) Θ is a compact set, and the true parameter value β_0 is an interior point of Θ.

  (C2) Let τ be the finite maximum support of the failure time.

  (C3) The { A_i, Y_i,X_i } are independent and identically distributed for i=1,⋯,n. 

  (C4) The covariate X_i is bounded.

  (C5) Conditional on X_i^∗, ( T_i^∗,  X_i^∗) are independent of A_i^∗.

  (C6) Censoring time C_i is non-informative. That is, the failure time T_i and the censoring time C_i are independent, given the covariates X_i.

  (C7) The regression function φ(·) is true, and its first order derivative exists.

Condition (C1) is a basic condition that is used to derive the maximizer of the target function. (C2) to (C6) are standard conditions for survival analysis, which allow us to obtain the sum of i.i.d. random variables and hence to derive the asymptotic properties of the estimators. Condition (C7) is a common assumption in SIMEX method.







§ PROOF OF THEOREM <REF>
 
Before presenting the proof, we first define some notation. For any s,t ∈ [0,τ], b = 1,⋯,B, and ζ∈𝒵, let

    B_1(t;b,ζ)    =    E [ R_i(t) r(t,Y_i,δ_i) λ'{ W_i^⊤(b,ζ) β(b,ζ) + H_0(t) }], 
    
    B_2(t;b,ζ)    =    E [ R_i(t) r(t,Y_i,δ_i) λ{ W_i^⊤(b,ζ) β(b,ζ)  + H_0(t) }], 
    
    B_1^W(t;b,ζ)    =    E [ W_i(b,ζ) R_i(t) r(t,Y_i,δ_i) λ'{ W_i^⊤(b,ζ) β(b,ζ)  + H_0(t) }], 
    
    B_2^W(t;b,ζ)    =    E [ W_i(b,ζ) R_i(t) r(t,Y_i,δ_i) λ{ W_i^⊤(b,ζ) β(b,ζ)  + H_0(t) }], 
    
    B(t,s;b,ζ)    =   exp{∫_s^t B_2^-1(u;b,ζ) B_1(u;b,ζ) dH_0(u)  },

where λ'(·) is the derivative of λ(·). We further define

    𝐰(t)    =    B_2^-1(t;b,ζ) [ B_2^W(t;b,ζ) . 
            . + ∫_t^τ{ B_1^W(s;b,ζ) - B_2^-1(s;b,ζ) B_2^W(s;b,ζ) B_1(s;b,ζ) } B(t,s;b,ζ) dH_0(s) ], 
    
    a(t)    =   1/π(t)∫_0^τ E [ W_i(b,ζ) δ_i R_i(s)/w(Y_i)( I(t ≤ s) ∫_t^s S_C(u) du . . 
            . . - w(s)/w(Y_i) I(t ≤ Y_i ) ∫_t^Y_i S_C(u) du) dΛ{ W_i^⊤ (b,ζ) β(b,ζ) + H_0(s) }],

where π(t) = P( Y_i - A_i ≥ t ).




Proof of Theorem <ref> (1):

The Kaplan-Meier estimator S_C(t) over [0,τ] is uniformly consistency to S_C(t) in the sense that sup_t ∈ [0,τ]| S_C(t) - S_C(t) | p⟶ 0 as n →∞ (Pollard 1990; van der Vaart 1998). It implies that as n →∞,

    sup_t ∈ [0,τ]| w(t) - w(t) | p⟶ 0

and

    sup_t ∈ [0,τ]| r(t,Y_i,δ_i) - r(t,Y_i,δ_i) | p⟶ 0

since sup_t ∈ [0,τ]| ∫_0^t S_C(u)du - ∫_0^t S_C(u)du | ≤sup_t ∈ [0,τ]∫_0^t| S_C(u)- S_C(u)|du  p⟶ 0 as n →∞. Furthermore, by the similar derivation in Step 1 of Wang and Wang (2015), we have

    .∂H(t;b,ζ,β)/∂β|_β = β_0   =    - ∫_0^t B(s,t;b,ζ)/B_2(s;b,ζ) B_1^W (s;b,ζ) dH_0(s) + o_p(1) 
       ≜    A(t) + o_p(1).



Let β(b,ζ) be a solution of E{U_SIMEX(β)} = 0. Since β(b,ζ) is a solution of U_SIMEX(β) = 0. By (<ref>), (<ref>), (<ref>), and the Uniformly Law of Large Numbers (van der Vaart 1998), we have that 1/n U_SIMEX(β) converges uniformly to E{U_SIMEX(β)}. Then we have that as n →∞,

    β(b,ζ) p⟶β(b,ζ).

By definition (<ref>), taking averaging with respect to b on both sides of (<ref>) gives that as n →∞,

    β(ζ) p⟶β(ζ)

for every ζ∈𝒵.
By (<ref>), we can show that as n →∞,

    Γp⟶Γ.

Since β_SIMEX = φ(-1,Γ), therefore, by the continuous mapping theorem, we have that as n →∞,

    β_SIMEXp⟶β_0.

 

Proof of Theorem <ref> (2):

By (<ref>), we have H(t;b,ζ,β_SIMEX) - H(t;b,ζ,β_0) = o_p(1) for every t ∈ [0,τ], b, and ζ. Taking average with respect to b gives H(t;ζ,β_SIMEX) - H(t;ζ,β_0) = o_p(1). On the other hand, by the Uniformly Law of Large Numbers and similar derivations in Wang and Wang (2015) with ζ→ -1, we have that as n →∞, H(t;β_0) - H_0(t) p⟶ 0 for all t ∈ [0,τ]. Therefore, we conclude that as n →∞ and ζ→ -1, H_SIMEX(t) - H_0(t) p⟶ 0 by the fact that H(t;-1,β_SIMEX) - H_0(t) = H(t;-1,β_SIMEX) - H(t;-1,β_0) + H(t;-1,β_0) - H_0(t).




Proof of Theorem <ref> (3):

For b = 1,⋯,B and ζ∈𝒵, applying the Taylor series expansion on (<ref>) around β(b,ζ) gives

    0    =    U_SIMEX( β(b,ζ) ) 
       =    U_SIMEX( β(b,ζ) ) + ∂ U_SIMEX( β(b,ζ) )/∂β{β(b,ζ) - β(b,ζ) } + o_p( 1/√(n)),

or equivalently,

    √(n){β(b,ζ) - β(b,ζ) }   =   (-∂ U_SIMEX( β(b,ζ) )/∂β)^-1√(n) U_SIMEX( β(b,ζ) ) 
             + o_p( 1 ).


By (<ref>), (<ref>), (<ref>), and the Uniformly Law of Large Numbers, we have that as n →∞, 
    (-∂ U_SIMEX( β(b,ζ) )/∂β) p⟶𝒜(b,ζ),

where 

    𝒜(b,ζ)    =    - E[ ∫_0^τ{ W_i(b,ζ) - 𝐰(t) } W_i^⊤(b,ζ) R_i(t) . 
            ×. r(t,Y_i,δ_i) λ' { W_i^⊤ (b,ζ) β (b,ζ) + H_0(t) } dH_0(t) ].



On the other hand, by the similar derivations in Wang and Wang (2015) and standard algebra, we can derive U_SIMEX(β(b,ζ)) as a sum of i.i.d. random functions. The exact formulation is given by

    √(n) U_SIMEX(β(b,ζ)) = 1/√(n)∑_i=1^n Ψ_i(b,ζ) + o_p(1),

where

    Ψ_i(b,ζ)    =   ∫_0^τ[ {W_i(b,ζ) - 𝐰(t)} dM_i(t) a(t) dM_C(t) ],

M_C = I(Y_i ≤ t, δ_i = 0 ) - ∫_0^t I(Y_i ≥ u ) dΛ_C(u), and Λ_C(·) is the cumulative hazards function of C.

Combining (<ref>) and (<ref>) with (<ref>) yields

    √(n){β(b,ζ) - β(b,ζ) } = 1/√(n)∑_i=1^n 𝒜^-1(b,ζ) Ψ_i(b,ζ) + o_p( 1 ).

By (<ref>), taking average with respect to b on both sides of (<ref>) gives

    √(n){β(ζ) - β(ζ) } = 1/√(n)∑_i=1^n Φ_i(ζ) + o_p( 1 )

for ζ∈𝒵, where Φ_i(ζ) =  1/B∑_b=1^B 𝒜^-1(b,ζ) Ψ_i(b,ζ).


Let β(𝒵) = {β(ζ) : ζ∈𝒵} denote the vectorization of estimator β(ζ) with every ζ∈𝒵. By the Central Limit Theorem on (<ref>), we have that as n →∞,

    √(n){β(𝒵) - β(𝒵) }d⟶ N(0, Ω( 𝒵) ),
 
where Ω( 𝒵) = cov{Φ_i(𝒵) }. By the Taylor series expansion on φ( 𝒵, Γ) with respect to Γ, we have

    φ( 𝒵, Γ) - φ( 𝒵, Γ) ≈∂φ( 𝒵, Γ)/∂Γ(Γ - Γ).

Let 𝒞 = ∂φ( 𝒵, Γ)/∂Γ and 𝒟 = {∂φ( 𝒵, Γ)/∂Γ}^⊤∂φ( 𝒵, Γ)/∂Γ. Combining (<ref>) and (<ref>) gives that as n →∞,

    √(n)(Γ - Γ) d⟶ N(0, 𝒟^-1𝒞Ω(𝒵) 𝒞^⊤𝒟^-1).

 Finally, let 𝒬 = 𝒟^-1𝒞Ω(𝒵) 𝒞^⊤𝒟^-1. Since the SIMEX estimator is defined by β_SIMEX = φ(-1, Γ), then combining (<ref>) and (<ref>) with ζ→ -1 and applying the delta method give that as n →∞,
 
    √(n)(β_SIMEX - β_0 ) d⟶ N(0, {∂φ/∂Γ(-1, Γ) }𝒬{∂φ/∂Γ(-1, Γ) }^⊤).

  
 
 Proof of Theorem <ref> (4):

 We first consider the expression of √(n){H(t;b,ζ,β_SIMEX) - H(t;b,ζ,β_0) }. By the Taylor series expansion with respect to β, we have
 
    √(n){H(t;b,ζ,β_SIMEX) - H(t;b,ζ,β_0) }
       =   ∂H(t;b,ζ,β_0)/∂β√(n)( β_SIMEX - β_0 ) 
       =    A(t) √(n)( β_SIMEX - β_0 ) + o_p(1) 
       =    A(t) 1/√(n)∑_i=1^n {∂φ/∂Γ(-1, Γ) }𝒟^-1𝒞Φ_i(𝒵) + o_p(1),

where the second equality is due to (<ref>) and the third equality is due to (<ref>).

By similar derivations in Step 3 of Wang and Wang (2015), we have

    √(n){H(t;b,ζ,β_0) - H_0(t) } = 1/√(n)∑_i=1^n ∫_0^t B(s,t;b,ζ)/B_2(s;b,ζ) dM_i(s) + o_p(1).

Then combining (<ref>) and (<ref>) gives

    √(n){H(t;b,ζ,β_SIMEX) - H_0(t) } = 1/√(n)∑_i=1^n 𝒯_i(t;b,ζ) + o_p(1),

where 𝒯_i(t;b,ζ) = A(t) {∂φ/∂Γ(-1, Γ) }𝒟^-1𝒞Φ_i(𝒵) + ∫_0^t B(s,t;b,ζ)/B_2(s;b,ζ) dM_i(s). Taking average on both sides of (<ref>) with respect to b yields

    √(n){H(t;ζ,β_SIMEX) - H_0(t) } = 1/√(n)∑_i=1^n 𝒯_i(t;ζ) + o_p(1),

where 𝒯_i(t;ζ) = 1/B∑_b=1^B 𝒯_i(t;b,ζ).

Suppose that φ_H(ζ,Γ_H(t)) is a function with the same conditions in (C7), and Γ_H(t) is the associated parameter depending on time t. For t ∈ [0,τ] and ζ∈𝒵, we fit a regression model on H(t;ζ,β_SIMEX) and φ_H(ζ,Γ_H(t)), and derive the estimator of Γ_H(t) which is denoted by Γ_H(t). Furthermore, similar to the derivations in (<ref>), we have  

    φ_H( 𝒵, Γ_H(t) ) - φ_H( 𝒵, Γ_H(t) ) ≈∂φ_H( 𝒵, Γ_H(t) )/∂Γ_H(t){Γ_H(t) - Γ_H(t) }.

Let 𝒰(t) = ∂φ_H( 𝒵, Γ_H(t) )/∂Γ_H(t) and 𝒱(t) = {∂φ_H( 𝒵, Γ_H(t) )/∂Γ_H(t)}^⊤∂φ_H( 𝒵, Γ_H(t) )/∂Γ_H(t). Combining (<ref>) and (<ref>) yields

    √(n){Γ_H(t) - Γ_H(t) } = 1/√(n)∑_i=1^n 𝒱^-1(t) 𝒰(t) 𝒯_i(t;𝒵) + o_p(1),

and since the estimator H_SIMEX(t) is a predicted value of φ_H(ζ,Γ_H(t)) by taking ζ→ -1, then by (<ref>) and the delta method with ζ→ -1, we obtain

    √(n){H_SIMEX(t) - H_0(t) }
       =   1/√(n)∑_i=1^n {∂φ_H(-1,Γ_H(t))/∂Γ_H(t)}𝒱^-1(t) 𝒰(t) 𝒯_i(t;𝒵) + o_p(1) 
       ≜   1/√(n)∑_i=1^n ℋ_i(t) + o_p(1).

where ℋ_i(t) = {∂φ_H(-1,Γ_H(t))/∂Γ_H(t)}𝒱^-1(t) 𝒰(t) 𝒯_i(t;𝒵).
Finally, by the Central Limit Theorem and similar derivations in Wang and Wang (2015), we conclude that √(n){H_SIMEX(t) - H_0(t) } converges to the Gaussian process with mean zero and covariance function E{ℋ_i(t) ℋ_i(s) }. 

□



=22.5pt

16

[Anderson et al. (1993)]R0 
Anderson, P.K., Borgan, O., Gill, R.D. and Keiding, N. (1993) Statistical Models Based on Counting Processes. Springer-Verlag New York.


[Carroll et al. (1996)]R1
Carroll, R. J., Kchenho, H., Lombard, F. and Stefanski, L. A. (1996) Asymptotics for
the SIMEX estimator in nonlinear measurement error model. Journal of the American
Statistical Association, 91, 242 - 250.

[Carroll et al. (2006)]R2
Carroll, R. J., Ruppert, D., Stefanski, L. A., and Crainiceanu, C. M. (2006) Measurement Error in Nonlinear Model. Chapman & Hall/CRC, New York.


[Chen (2019)]R3
Chen, L.-P. (2019). Pseudo likelihood estimation for the additive hazards model with data subject to left-truncation and right-censoring.  Statistics and Its Interface, 12, 135-148.

[Chen and Yi (2018)]R4
Chen, L.-P., Yi, G. Y. (2018). Semiparametric methods for left-truncated and right-censored survival data  with covariate measurement error. Submitted.

[Cook and Stefaski (1994)]R5
Cook, J. R. and Stefaski, L. A. (1994). Simulation-extrapolation estimation in parametric
measurement error models. Journal of the American Statistical Association, 89, 1314-1328.

[Hosmer et al. (2008)]Hosmer2008
Hosmer, D. W., Lemeshow, S. and May, S. (2008) Applied Survival Analysis : Regression Modeling of Time to Event Data. Wiley, New York.


[Huang et al. (2012)]R6
Huang, C. Y., Qin J., and Follmann, D. A. (2012) A maximum pseudo-profile likelihood estimator for the Cox model under length-biased sampling. Biometrika, 99, 199-210.

[Huang and Qin  (2013)]R7
Huang, C. Y. and Qin J. (2013) Semiparametric estimation for the additive hazards model with left-truncated and right-censored data. Biometrika, 100, 877-888.

[Kalbfleisch and Prentice (2002)]Kalbfleisch2002
Kalbfleisch, J. D. and Prentice, R. L. (2002) The Statistical Analysis of Failure Time Data. Wiley, New York.

[Lancaster (1990)]Lancaster1990
Lancaster, T. (1990). The Econometric Analysis of Transition Data. Cambridge University Press.

[Nakamura (1992)]R8
 Nakamura, T. (1992). Proportional hazards model with covariates subject to measurement error. Biometrics, 48, 829-838.

[Ning et al. (2014)]R9
Ning, J., Qin, J., and Shen, Y. (2014) Semiparametric accelerated failure time model for length-biased data with application to dementia study. Statistica Sinica, 24, 313–333.

[Pollard (1990)]Pollard1990
Pollard, D. (1990) Empirical Process: Theory and Applications. NSF-CBMS Regional Conference Series in Probability and Statistics, Vol. 2. IMS, Hayward, CA.



[Qin and Shen (2010)]R10
Qin, J. and Shen, Y. (2010) Statistical methods for analyzing right-censored length-biased data under Cox model. Biometrics, 66, 382-392.


[Qin et al. (2011)]R11
Qin, J., Ning, J., Liu, H. and Shen, Y. (2011). Maximum likelihood estimations and EM algorithms with length-biased data. Journal of the American Statistical Association, 106, 1434–1449.

[Shen et al. (2009)]R12
Shen, Y.,  Ning, J., and Qin, J. (2009). Analyzing length-biased data with semiparametric transformation and accelerated failure time models. Journal of the American Statistical Association, 104, 1192–1202.
[Su and Wang (2012)]R13
Su, Y. and Wang, J. (2012) Modeling left-truncated and right-censored survival data with longitudinal covariates. The Annals of Statistics, 40, 1465 - 1488.

[Yi and Lawless (2007)]R14
Yi, G. Y. and Lawless, J. F. (2007)  A corrected likelihood method for the proportional hazards model with covariates subject to measurement error. Journal of Statistical Planning and Inference, 137, 1816-1828.

[van der Vaart (1998).]R-van
van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press, New York.

[Yi and Reid (2010)]R17
Yi, G. Y. and Reid, N. (2010). A note on misspecified estimating functions. Statistica Sinica, 20, 1749–1769.


[Yi (2017)]R15
Yi, G. Y. (2017).  Statistical Analysis with Measurement Error and Misclassication: Strategy, Method and Application. Springer, New York.

[Wang and Wang(2015)]R16
Wang, X. and Wang, Q. (2015) Estimation for semiparametric transformation models with length-biased sampling. Journal of Statistical Planning and Inference, 156, 80–89.










