Semiparametric Regression Estimation for Recurrent Event Data with Errors in Covariates under Informative Censoring

1 Introduction

In many longitudinal follow-up studies, recurrent event data are collected when subjects experience an event multiple times. For example, patients with superficial bladder cancer may experience tumor recurrence many times; patients with cystic fibrosis may experience repeated lung exacerbations; patients with chronic granulomatous disease may experience repeated pyogenic infections [1, 2]. Models for recurrent event data can be categorized into two different classes: time-to-event or gap time models. In time-to-event models, interest focuses on the occurrence rate of an event over time [3–5]. In gap time models, interest lies in the gap time between two consecutive events [6].

In this study, we focus on the time-to-event models. The time-to-event models may be constructed on the basis of an intensity function [7] or a rate function [3, 8]. The intensity function uniquely determines the probability structure of the recurrent event process. However, it needs to specify the occurrence of an event given the prior event history correctly. On the other hand, the rate function allows for arbitrary dependence among the recurrent events and provides a direct interpretation on the occurrence rate without conditioning on the prior event history. Our primary focus is to assess the average effects of treatments or risk factors, that is, we are mainly interested in the inference of the rate function. Lawless and Nadeau [9] estimated the cumulative rate function nonparametrically, and applied their approach to industrial warranty data. In addition, Hu and Lagakos [8] proposed a nonparametric method to study the rate function of viral load changing process for HIV infected patients. Nevertheless, all of the above approaches need to assume non-informative censoring or the observation mechanism is independent of the recurrent process. In practice, the assumption is usually violated; for example, when the recurrent event process is interrupted by some terminal events that are related to the recurrent events. A potential remedy is to consider a frailty model which allows dependence between the recurrent event process and the informative drop-out through a non-negative frailty variable. In general, the distribution of the frailty variable is assumed to be known [10] and thus the likelihood-based approach [11] is preferred. More recently, Kalbfleisch et al. [12] proposed a weighted estimating equation approach with the weight specified by a gamma frailty distribution. However, in general it is not easy to verify the frailty distribution due to invisibility of the frailty variables. To avoid specification of the frailty distribution, Wang et al. [13] and Wang and Huang [14] considered a conditional likelihood approach, where the unobserved frailty variables are “conditioned away” in their proposed estimating equations.

The aforementioned approaches, nevertheless, require that the covariates are correctly measured. In many epidemiologic or medical studies, the covariates may suffer from measurement errors. For example, baseline plasma selenium level is an important predictor for the occurrence of skin cancers in the Nutritional Prevention of Cancer (NPC) trial study [15]. However, the true value of plasma selenium level can never be measured because of intrinsic biological variability or limited instrumental precision. Instead, the values we observed are contaminated with measurement errors. The most convenient approach is to treat the observed covariates as the true covariates in the regular estimating procedure, which is also referred to as the naive approach. However, the naive estimator obtained from this approach is generally known to be inconsistent ([16], Chapter 3). In survival and longitudinal data analysis, intensive research has been done to deal with measurement error problems. For Cox regression, Prentice [17] proposed likelihood approaches with normal measurement error and rare disease assumptions. Wang et al. [18] applied regression calibration to the partial score function, and investigated the performance of the regression calibration estimator through simulation studies; whereas, Nakamura [19] constructed unbiased estimating equations based on the concept of corrected scores. For nonlinear mixed models, Wu [20], Liu and Wu [21] and Wu et al. [22] proposed estimating approaches for longitudinal response data when the covariates are measured with errors, which can also handle censoring in the response and missing data. In recurrent event analysis, nonetheless, little has been addressed for measurement error problems. Under a normal measurement error assumption, Jiang et al. [23] proposed a moment corrected method to adjust for the bias of a naive estimator under a semi-parametric model. However, their approach not only requires the assumption of non-informative censoring but also assumes that the censoring distribution is independent of the covariates.

The present study is motivated by the NPC trial study, which aimed to assess the efficacy of oral supplement of plasma selenium in preventing the development of skin cancers such as squamous cell carcinoma (SCC). This clinical trial began in 1983 and had included approximately 1,300 patients with dermatologic cancer histories. Nearly half of the patients in the NPC trial were randomly assigned to the placebo or treatment group respectively. Patients in the treatment arm were supposed to take 200 μg of plasma selenium supplement per day. In the study period, the patients in the trial might experience SCC events repeatedly. Each incidence of a new SCC was diagnosed and recorded by certified doctors. The medical records were reviewed by the clinical coordinators at the semi-annual visit, the annual contact or by self-report to ensure the completeness of the data. At the time of randomization, many prognostic risk factors of SCC were recorded including the baseline plasma selenium level. As we mentioned, the plasma selenium level may be measured with errors. In the original study, Clark et al. [15] did not take measurement error into account and found a nonsignificant negative plasma selenium effect on developing SCC. The result contradicted the evidence of the previous studies which showed high correlation between plasma selenium level and several kinds of cancer. Later, many studies focused on the effect of plasma selenium level on the recurrences of SCC by assuming an independent censoring assumption, some of which also took measurement error into account [23]. However, we found a significant negative relationship between the censoring time and the SCC occurrence rate. This implies that the independent censoring assumption is not satisfied. Therefore, the existing methods are not appropriate for the NPC trial data.

This paper is organized as follows. In Section 2, statistical models for recurrent events and measurement errors are given. In Section 3, we propose a regression calibration method and a moment corrected method to correct the measurement errors in the presence of informative censoring. The simulation results are given in Section 4 to investigate the finite sample performance of the proposed methods. Then, we applied the proposed methods to the NPC trial data to evaluate the effect of selenium on the recurrence of SCC in Section 5, and concluded with a discussion in Section 6. The regularity conditions and technical proofs are provided in the Appendix and the Supplementary Information.

2 Model illustration

3 Correction for errors-in-variable

Assume the observed data {(Ci,(Ti1,…,Timi),{Wi1,…,Wiki},Zi),i=1,…,n} are independent and identically distributed (iid), where Tij denotes the observed event times for j=1,…,mi, and mi denotes the number of recurrent events occurred before Ci for subject i. As we mentioned in Section 2.1, C is conditionally independent of the recurrent event process N(t) given (ν,X,Z). Then, by eq. (2) we have

If Λ0(t) and X are known, the estimating equations ∑i=1n(1,Xi′,Zi′){miΛ0−1(Ci)−ea0+βX′Xi+βZ′Zi}=0 for (βX,βZ) are unbiased. In practice, they cannot be implemented since Xi is unobserved and Λ0(t) is unknown. To deal with the unknown function Λ0(t), we start with the conditional likelihood function of (Ti1,…,Timi) given (Ci,νi,mi,Xi,Zi). Under the Poisson process assumption, such a conditional likelihood can be constructed from a set of iid random variables with truncated density ∏j=1miλ0(Tij)/Λ0(Ci)I(0≤Tij≤Ci). Define a rescaled baseline function ϕ(t)≡λ0(t)/Λ0(τ) and Φ(t)=∫0tϕ(u)du=Λ0(t)/Λ0(τ) for t∈[0,τ], where Φ(τ)=1. The conditional likelihood is given by ∏i=1nP(Ti1,…,Timi|Ci,νi,mi,Xi,Zi) which is proportional to ∏i=1n∏j=1miϕ(Tij)/Φ(Ci). As pointed out by Wang et al. [13], the conditional likelihood shares the same form as the nonparametric likelihood for right-truncated data. Thus, Φ(t) can be consistently estimated by the product limit estimator

where {T(l)} are the ordered and distinct values of {Tij}i=1,…,n;j=1,…,mi, n(l) is the number of events occurred at T(l), and N(l) is the number of events which satisfy Tij≤T(l)≤Ci. Note that the non-parametric estimation of Φ does not require any information from the covariates and the unobserved frailty variable. Hence, Φˆ(t) is a consistent estimator even if X is measured with errors or the frailty distribution is unspecified.

For the issue of identifiability, let μν=1 without loss of generality. The expectation of the event number divided by the rescaled baseline function before time C is

where β0=log(Λ0(τ)). With the above equation, we can construct the unbiased estimating equations by Φ(t) instead of the unknown Λ0(t). After replacing the unknown X with the average of the replicates W‾i=∑j=1kiWij/ki, we can obtain the naive estimating equations

Then, the naive estimator βˆN=(βˆN,0,βˆN,X′,βˆN,Z′)′ is obtained by solving eq. (3) and Λ0(t) can be estimated by Λˆ0N(t)=Φˆ(t)exp(βˆN,0). Due to the measurement error in W‾, it can be shown that βˆN does not converge to the true parameter β, where β=(β0,βX′,βZ′)′. Based on eq. (3), we develop a regression calibration method and a moment corrected method to adjust for covariate measurement errors in the following subsections.

4 Simulation study

In this section, we evaluate the performance of the RC and MC methods with the naive approach under the semi-parametric model via the simulation studies. Additionally, the corrected partial likelihood (CPL) approach proposed by Jiang et al. [23] is also listed for comparison. The CPL estimator takes measurement error into account but assumes non-informative and covariate-independent censoring.

We consider a regression model with a continuous covariate X and a discrete covariate Z. Let X∼N(0,σX2=1/3) be the error-prone covariate which is unobserved, while Z∼Bin(0.5) be a random treatment assignment and is precisely obtained. For subject i, we generate ki repeated surrogates Wij=Xi+Uij for Xi where ki is generated from a discrete uniform distribution ranging from 1 to 4 and Uij∼N(0,σU2). With the repeated surrogates, we estimate the nuisance parameter γ by solving ∑i=1nΨi(γ)/n=0, where Ψ is shown in the Appendix. We conduct the simulations with reliability ratio (RR) σX2/(σX2+σU2)=0.8 and 0.5. The reliability ratio is used to represent the magnitude of the error contamination, and lower reliability ratio indicates higher error contamination. We generate νi∗ from a mixture model of which νi∗ follows a uniform distribution ranging from 0.5 to 1.5 when Zi=0, and follows a uniform distribution ranging from 1.5 to 4 otherwise. Then, the frailty variable is νi=exp(−Zlog(2.75))νi∗. When (νi,Xi,Zi) is given, the recurrent event process {Ni(t)} is generated with intensity function λ(t|νi,Xi,Zi)=νiλ0(t)exp(βXXi+βZZi) in which λ0(t)=(t−6)3/360+0.6,t∈[0,τ], τ=10 for i=1,…,n. We consider two distinct coefficient parameters (βX,βZ)=(log(1.5),log(1.5)) and (βX,βZ)=(log(3),log(1.5)). To show the robustness of the proposed estimators, the first two scenarios are conducted under different censoring time settings. In Scenario 1, we let the censoring time C depend on W. When Wi1>0, Ci is generated from an exponential distribution with mean 10νi−1 and is truncated after τ=10; otherwise, Ci is generated from an exponential distribution with mean 0.5νi−1 and is truncated after τ=10. In Scenario 2, we let the censoring time C depend on X. We generate Ci from the mixed exponential distribution in the same way as in Scenario 1 with Wi replaced by Xi. In addition, we conduct two cases to investigate the sensitivity of the conditional normal assumption imposed on the covariate X. In Scenario 3, X is uniformly distributed over the interval −3σX2,3σX2 and Z is allowed to be correlated with X. Let Z∗=X+ε where ε∼N(0,σX2), and Z=1 if Z∗≤0 and Z=0 otherwise. The other variables are generated the same as those in Scenario 2. Further, a non-normal measurement error case is considered in Scenario 4. We generate measurement error U from a skew normal distribution with mean 0, variance σU2 and skewness parameter α=−2, and X from N(0,σX2=1/3). The remaining variables are generated the same as those in Scenario 3. A total of 200 replicates with sample sizes n=300 and n=600 are generated in each simulation configuration. In the tables, BIAS denotes the average bias, ASE denotes the average standard error estimation, ESD denotes the empirical sample standard deviation, and CP and CL denote the coverage probability and average interval length of the 95 % confidence interval based on 200 runs. The standard errors of the proposed estimators are obtained by taking the square roots of the diagonal elements from the sandwich variance estimators given in Appendices B and C.

The results of Scenarios 1 to 4 are demonstrated in Tables 1 to 4. In general, the naive estimator for the error-prone covariate X has large biases and disastrous coverage probabilities as shown in all tables. This phenomenon is due to the common attenuation effect. The degree of bias becomes critical when the error-prone covariate effect is large and the reliability ratio is low. In Scenarios 1 and 2, the naive estimation of the effect of Z is not affected by the measurement errors since X and Z are generated to be mutually independent. In Scenario 3 in which X and Z are correlated, the naive estimator for βZ also has low coverage probabilities which is shown in Tables 3 and 4. Further, the numerical equivalence of the RC and MC estimators is also seen in the simulation results.