Simulation Extrapolation Method for Cox Regression Model with a Mixture of Berkson and Classical Errors in the Covariates using Calibration Data

1 Introduction

Many epidemiological or biomedical studies often aim to investigate the association between the time to an event of interest and some covariates. The examination of the association is usually performed under the popular Cox proportional hazards model [1] that stipulates proportional covariate effects on the hazard function for the time-to-event. A commonly encountered problem in such an examination is that some covariates potentially involve measurement error. This may be due to imprecision of the data measuring instrument, high costs to obtain precise measurements or the nature of the covariates themselves. For example, CD4 counts in HIV studies [2], dietary protein intake in nutrition studies and radiation dose in radiation epidemiology studies [3] are well-known to be subject to measurement error. Naively omitting the measurement error and simply using the error-affected covariate data in the survival analysis could lead to invalid inference [4].

There are two different types of measurement errors, namely the classical type and the Berkson type [3]. The main difference between these two types of errors lies in their relationship with the unobserved covariates. The classical error is assumed to be independent of the unobserved covariates whereas the Berkson error is positively correlated with the unobserved covariates. Self-reported data, collected through questionnaire typically involve recall error that is usually assumed to be of the classical type [5]. Also, grouped-exposure data obtained by assigning the same value of the unobserved exposure variable (e. g. air pollution or radiation dose) to individuals sharing some characteristics are often assumed to be contaminated with Berkson error [6]. In many other situations, the covariate data may be affected by a mixture of errors of both classical and Berkson types. For example, DS02 radiation dose estimates for atomic-bomb survivors, who were followed up by the Radiation Effects Research Foundation (RERF) are contaminated with both errors [7]. The classical error is thought to originate from the survivors’ individual recollections of location and shielding at the time of the bombings and the Berkson error is believed to be associated with averaging radiation dose estimates among survivors with common location and shielding type. Taken individually and if omitted, measurement error of each of these two types could affect the estimation of the survival model parameters. The impact of the classical error in the covariates of a Cox regression model has been well-documented and includes the attenuation of the covariate effect [8, 9]. Moreover, a few references have reported that the Berkson error may also have a similar effect [6, 10]. However, the effect of a mixture of classical and Berkson errors in the covariates of a proportional hazards model still remains to be well investigated. The presence of errors of both types in the covariates of the regression model is expected to worsen the attenuation effect as will become clear in the simulation study section. It is therefore important to simultaneously account for both types of errors when there are present in the covariates of a Cox proportional hazards model.

Several methods have been developed to deal with covariate measurement error in a Cox proportional hazards model and extensive attention has been given to the classical error. Popular methods adjusting for classical error in survival analysis are the regression calibration (RC) method [11, 12], simulation extrapolation (SIMEX) method [3, 13, 14], parametric corrected score [15], non-parametric corrected score [8, 16], conditional score [17, 18], maximum likelihood method [19, 20] and non-parametric maximum likelihood method [21]. On the other hand, methods addressing the issue of Berkson error in the covariates in Cox regression model are very scarce and this problem still has to be well-discussed in the literature. Among the very limited number of references, [22] discussed an extension to a Berkson error of a corrected profile likelihood method for Cox model with classical error in the covariates. It is important to note that methods that simply deal with classical error or Berkson error exclusively may not work well when errors of both types are involved in the covariates. Moreover, the existing methods accounting for a mixture of errors of both types are developed for generalized linear models [23, 24] and may not be directly applicable in the context of survival analysis. Therefore, new statistical methods that jointly adjust for both classical and Berkson errors need to be developed for a reliable inference in this situation.

In this paper, we address the problem of parameter estimation in Cox regression analysis when some covariates are subject to a combination of both Berkson and classical errors. We are not aware of an existing reference in the statistical literature that has properly dealt with this problem. We assume that two replicates of the error-prone covariates are available for all individuals. Also, we assume that data on an instrumental variable for the mis-measured covariate are available only for some individuals in a subset of the study cohort that is here referred to as calibration subsample. For illustration, we consider the AIDS Clinical Trials Group (ACTG) 175 study, a randomized double-blind clinical trial that was conducted to compare four antiretroviral treatments with either a single nucleoside or two nucleosides in (HIV-1) infected adult participants [25]. In addition to sociodemographic data, the study also collected data on survival time and virologic variables such as CD4 cell count and CD8 cell count. An important interest is in assessing the effect of the true baseline CD4 cell count on time to AIDS disease development or death under a Cox proportional hazards model adjusting for the treatment effect. A complication related to the use of standard survival analysis methods for this assessment is the possible contamination of the CD4 cell count data with measurement error, which has been assumed to be of the classical type in previous applications [8, 16]. Here we allow the measurement error to incorporate features of both classical and Berkson error types. We use a mixture of classical and Berkson errors model, which is more general than the classical model. Moreover, we treat the two CD4 cell counts taken before the start of the treatment and within three weeks of randomization as replicates. Also, the CD8 cell count measured at randomization is likely correlated with the baseline CD4 cell count and may be used as instrumental variable for the true baseline CD4 count. There were 2441 participants with two replicates for the true baseline CD4 cell count and only 1459 of them had CD8 cell count measured at randomization. We propose a SIMEX-based new approach to adjusting for both classical error and Berkson error under these settings. It places no assumption on the mixture proportion of the error variances.

In the next section, we provide the model formulations for the survival data and the mixture of classical and Berkson errors. In Section 3, we present the naive and RC methods and develop the SIMEX method for this problem in Section 4. We present a simulation study in Section 5 and apply the developed method to the data from the ACTG 175 study in Section 6. Finally, we conclude with some discussions in Section 7.

2 Model formulations

Suppose that there are n study subjects, which are followed-up over a study period of length τ. For subject i, let Ti0 denote the unobserved survival time, which is assumed to be right-censored and Ci be the censoring time. The observed survival data are (Ti,δi), where Ti=min(Ti0,Ci) is the observed failure time and δi=I(Ti0≤Ci) is the non-censoring indicator, i = 1,…, n. Also, let Xi denote a vector of time-independent covariates that are subject to measurement error, and Zi represent a vector of other time-independent covariates that can be measured accurately. It is assumed that replicates Wij, j=1,…,ki(≥2) of Xi are available for all subjects and data on an instrumental variable Qi for Xi are available only for some subjects in a calibration subsample. Letting ηi indicate whether Qi is observed (ηi=1) or not (ηi=0) and Wi=(Wi1′,…,Wiki′)′, the observed data are (Ti,δi,Wi,Zi,ηiQi), i = 1,…, n, which are assumed to be independent and identically distributed random vectors. In the ACTG 175 data example, δi is the indicator of the event of interest, which is AIDS disease progression or death during the study follow-up and Ti is the time to the event if δi=1 or end of follow-up (due to lost to follow up or end of the study observation period) if δi=0. Also, Xi is a single covariate representing the log-transformation of baseline CD4 counts, Zi is a vector of indicator variables of randomized treatment groups, Wij (j = 1,2) are the log-transformation of two CD4 count measurements taken before treatment initiation and Qi is the log-transformation of baseline CD8.

We consider the following mixture of Berkson and classical errors model for the relation between Wij and Xi:

where Li is a latent variable with mean μl and variance Σl, and Ub,i is a zero-mean error with variance Σb. Also, Uc,ij, j=1,…,ki, are independent and identically normally distributed with mean 0 and variance Σc. Moreover, Uc,ij is independent of Ub,i, j=1,…,ki, (Uc,i1,…,Uc,iki) and Ub,i are independent of (Ti0,Ci,Li,Zi), i = 1,…, n. Model (1), which was also studied by [23, 26] and [27] in different contexts encompasses features of both classical and Berkson error models. It reduces to a Berkson measurement error model when Σc=0 and a classical error model when Σb=0. Furthermore, Qi is assumed to be associated with Xi according to the following model:

where α=(α0,α1′)′ is a vector of unknown parameters and ϵi is a zero mean normally distributed random variable with variance σϵ2 and independent of (Ti0,Ci,Li,Zi,Ub,i,Uc,i1,…,Uc,iki), i = 1,…, n. The time-to-event is associated with Xi and Zi via the Cox proportional hazards model, which is a commonly used tool for survival data analysis. Assuming that Ti0 is conditionally independent of Ci, Wi and Qi given Xi and Zi, the model formulates the hazard function as follows.

where λ0(.) is an unspecified baseline hazards function and β=(βx′,βz′)′ is a vector of unknown parameters. Our interest lies in the estimation of β based on the observed data. In what follows, we let Yi(u)≡I(Ti≥u) denote the at-risk process and define the counting process Ni(u)≡I(Ti≤u,δi=1), where I(.) is the indicator function.

3 Naive and RC methods

The naive and RC estimations are two simple imputation-based methods that could be used to estimate the parameter β although, they might not handle correctly the measurement errors, especially when the error magnitudes are not small.

4 SIMEX approach

The SIMEX method is a two-step procedure consisting of a simulation step and a subsequent extrapolation step. It was originally developed by [13] to deal with classical error in the covariates of a linear or nonlinear regression model and later studied by [3, 14, 31] and [32]. In addition, applications of this method to a mixture of errors of both types in the covariates of a generalized linear models are discussed in [3] and [24]. Here we provide an extension of the method to the Cox proportional hazards model when some covariates are subject to both classical and Berkson errors. We describe the implementation of the two steps of the procedure for our problem in what follows. The simulation step requires an estimating function that would yield a consistent estimator of β if there were no classical error (Σc=0) and Xi were substituted with Wˉi in that estimating function. The partial likelihood score in (4) would have been a suitable choice for this purpose if the error were simply of classical type. However, its use in our context is complicated by the presence of Berkson-type error in addition to classical error. We instead use estimation equations based on a maximum likelihood function.

It should be noted that if Xi were observed and λ0(.) known, the maximum likelihood estimator of β would solve the following equation:

Also, it is useful to note that Li=Wˉi if the error were simply of Berkson type (Σc=0) under Model (1). Assuming that λ0(.) and the distribution of (Li,Zi) do not involve β, a consistent estimator for β in the presence of simply Berkson error (Σc=0) would be obtained by solving

Also, a Breslow type estimator for the baseline hazards function λ0(.) could be expressed as λˆ0(u)=∑i=1ndNi(u)/∑i=1nYi(u)Eexp(βx′Xi+βz′Zi)|Ni,Yi,Li,Zi. However, Li cannot be observed due to the presence of classical error (in addition to the Berkson error) and simply replacing Li with Wˉi in (5) will potentially lead to a biased estimator for β in our situation. In the simulation step of the proposed SIMEX-based method, a number of R (with R>1) naive estimates of β are obtained by solving (5), in which Li is imputed with Wζ,r,i=Wˉi+ζUr,i, where ζ is a non-negative scalar and Ur,i∼N(0,ki−1Σc), i = 1,…, n, r=1,…,R. Let βˆr(ζj) denote the naive estimator for  β using Wζj,r,i in place of Li in (5) and denote  βˆ(ζj)=R−1∑r=1Rβˆr(ζj) the pseudo-estimate of β using ζj, where 0=ζ0<ζ1<…<ζJ, and J>1 is the number of pseudo-estimates. The extrapolation step consists in fitting a regression model of each component of  βˆ(ζj) on the ζj’s, j = 1,…, J, using the ordinary least squares estimation method. The SIMEX estimate of each component of β is then obtained by extrapolating back to the case when ζ=−1, which represents the situation of no classical error. Common extrapolation functions are polynomial in ζ of the form P(ζ,a)=a0+a1ζ+…+aqζq, where q≥1 is the degree of the polynomial. For a particular component of β, an estimator aˆ=(aˆ0,…,aˆq)′ for the parameter a=(a0,…,aq)′ is obtained via the regression of the corresponding components of the βˆ(ζj)′s on the ζj’s, j = 1,…, J. For example, letting βˆx(ζj) represent the component of βˆ(ζj) that corresponds to βx, and βˆx(ζ)={βˆx(ζ1),…,βˆx(ζJ)}′, an estimator for a would be aˆ={D(ζ)′D(ζ)}−1D(ζ)′βˆx(ζ), where D(ζ) is the J × q-matrix whose jth row is (1,ζj,…,ζjq)′. The SIMEX estimate of that particular component of β is then obtained as P(−1,aˆ)=aˆ0−aˆ1+…+aˆq(−1)q. Here R, J and q are assumed known. The standard errors can be obtained based on the bootstraps method. The SIMEX method leads to a consistent estimator for β if the true expression of the extrapolation function is used. In practice, approximations to the form of the extrapolation function are necessary as its closed form expression is unknown. It is important to note that the performance of the method could be affected by the choice of the extrapolation function.

A conditional expectation of the form E{h(Xi)|Ni,Yi,Li,Zi} involved in (5) is given as follows.

where

The nuisance parameters μ, α, Σl, Σb, Σc and σϵ2 are needed for the evaluation of the conditional expectations involved in the implementation of the RC and SIMEX methods. They may be estimated from the data if unknown. For the single covariate case for example, the vector of nuisance parameters is ν=(μ,α0,α1,σl2,σb2,σc2,σϵ2), where σl2=Σl, σb2=Σb and σc2=Σc. All the components of ν can be identified using the observations (Wi1,…,Wiki,ηiQi,Ti), i = 1,…, n. Moreover, a consistent estimator of ν can be obtained based on the method of moments, which leads to the following estimating equations:

where Tˉ=n−1∑i=1nTi.

5 Simulation study

We performed a simulation study to examine the finite-sample performance of the proposed SIMEX method and compare it to the naive and RC methods for the case of a single covariate with a mixture of Berkson and classical measurement errors. The latent variable Li was generated from N(μl,σl2) with μl=0 and σl2=1. The unobserved covariate Xi was simulated as Xi=Li+Ub,i, where the Berkson error, Ub,i followed a zero-mean normal distribution with variance σb2. Two replicates, Wij, j=1,2=ki, were generated from the classical error model Wij=Li+Uc,ij with Uc,ij generated from N(0,σc2). We set σb2=0.3 or 0.5, and σc2=0.3 or 0.5 to study the separate and combined effects of the Berkson error and classical error on the performances of the various estimators. The survival times were simulated following the Cox proportional hazards model with unit baseline hazard λ(u|Xi,Zi)=exp(βxXi+βzZi), where Zi followed a Bernoulli trial with a success probability of 0.5, βx=log(2) and βz=0. A common censoring was applied to subjects inducing a censoring rate of 50%. The instrumental variable was simulated following the model Qi=α0+α1Xi+ϵi, where ϵi was normal with mean 0 and variance σϵ2. We took α0=1, α1=0.8 and ϵi=0.5. The variable ηi, indicating whether Qi is available or not was generated from the Bernoulli distribution with probability of success P(ηi=1)=c. We set c = 0.5 or 0.7 to explore how the proportion of the calibration data influences the results.

A total of 500 Monte Carlo samples of size n = 500 or 1000 were generated in the simulations. We estimated the parameter of interest β=(βx,βz)′ based on the naive method (Naive), which replaces Xi by Wˉi=(Wi1+Wi2)/2, RC method (RC1) that uses E(Xi|Wˉi) in place of Xi, RC method (RC2) that substitutes Xi by E(Xi|Wˉi,Qi) if ηi=1 and E(Xi|Wˉi) otherwise and the SIMEX method with quadratic polynomial (SIM1) or polynomial of degree 4 (SIM2) as extrapolation function in extrapolation step. For the SIMEX method, we created R=50 additional data sets of measurement error at each point ζ∈0,0.25,0.5,0.75,1 in the simulation step. In addition, the vector of nuisance parameters ν=(μl,α0,α1,σl2,σb2,σc2,σϵ2)′ was estimated by the method of moments. Integrals involved in the evaluation of condition expectations for the implementation of RC and SIMEX methods were computed using Gauss-Hermite integration techniques with 20 quadrature points. The estimators were evaluated with regard to their biases (Bias), sample standard deviation of the estimates (SD), average of the estimated standard erros (ASE) of the estimators, mean squared errors (MSE=Bias2+SD2) and coverage probabilities (CP) of their 95% Wald confidence intervals. The standard errors for RC1, RC2 estimators were computed using the sandwich method, whereas the estimation of the standard error for the SIMEX estimators was based on the bootstrap method, re-sampling 30 times.

Table 1 shows the results of the simulations pertaining to the estimation of the parameter β for 500. For βx, the naive estimator noticeably shows the worst performance among all the methods with respect to biases and coverage probabilities. Its biases and MSE’s are substantially large and increase as the variance of the Berkson error or that of the classical error becomes large. The attenuation effect seems to be more serious with the classical error than the Berkson error. The bias problem is severely accentuated by the presence of errors of both types. Also, it can be seen that its coverage probabilities are well below the nominal level of 95%. Hence, it is important to correct for both types of errors when they are present in the covariates of a Cox regression model. RC methods appear to have reduced the bias from the naive estimator although not completely. RC2, which uses the data in the calibration subsample generally exhibits smaller biases (except when σb2=0.3 and P(η = 1) = 0.5) and larger coverage probabilities than RC1. Both SIM1 and SIM2 display smaller biases than the RC-based estimators. SIM1 works well in correcting the bias from the Berkson error and performs acceptably when the magnitude of the classical error is small. However, its bias becomes noticeable when the variance of the classical error is large. This highlights the importance of the choice of the extrapolation function for SIMEX method, which may perform well in terms of bias reduction when the choice is appropriate. Meanwhile, SIM2 appears to show very small biases and good coverage probabilities, although its MSE is larger than the RC and SIM1 counterparts. It can also be noted that the biases and ASE’s of RC2, SIM1 and SIM2 decrease to some extent as the proportion of the calibration data increases. Moreover, all methods, including the naive one perform equally well with respect to biases for the estimation of the parameter βz. The results of the simulations for the sample size n = 1000 are shown in Table 2. The performances of the methods are similar to those observed in Table 1. Also, the biases and ASE’s are smaller than those presented in Table 1 for all the methods.

The simulation results for the estimation of the nuisance parameters by the method of moments are presented in Table 3 for the setting with σb2=0.5 and σc2=0.5. The estimators generally show small bias and acceptable coverage probabilities. Also, it can be noted that the standard errors for the nuisance parameters involved in the distribution of the instrumental variable get smaller as the size of the calibration subsample increases. The biases and ASE’s of the estimators decrease as the sample size increases.

In addition, we conducted some simulations to investigate the performances of the methods when Xi is not normally distributed. The simulation settings were similar to those in Table 1 with the exception that Li was simulated from gamma distribution G(1,1), Ub,i followed the uniform U(−σb3,σb3), σb2=0.5, σc2=0.5 and P(ηi=1)=0.5. We estimated the parameter β by all the methods under the correctly specified or misspecified models for Li and Ub,i to gain insight into potential effects of models misspecification on the performances of the methods. The misspecified model for Ub,i assumed normality for Ub,i even though this assumption did not hold. The vector of nuisance parameters was estimated using the method of moments. The results of these simulations are reported in Table 4. The performances of the methods under the correctly specified models are similar to those observed in Table 1 and Table 2 with the SIMEX methods showing smaller biases and better coverage probabilities than the RC methods. Also, the SIMEX estimators show similar performances under both correct and misspecified models situations, suggesting that the method may not be very affected by a misspecification of the distribution for the latent variable or the measurement error.