 Research article
 Open Access
 Published:
A multiple imputation method based on weighted quantile regression models for longitudinal censored biomarker data with missing values at early visits
BMC Medical Research Methodology volume 18, Article number: 8 (2018)
Abstract
Background
In patientbased studies, biomarker data are often subject to left censoring due to the detection limits, or to incomplete sample or data collection. In the context of longitudinal regression analysis, inappropriate handling of these issues could lead to biased parameter estimates. We developed a specific multiple imputation (MI) strategy based on weighted censored quantile regression (CQR) that not only accounts for censoring, but also missing data at early visits when longitudinal biomarker data are modeled as a covariate.
Methods
We assessed through simulation studies the performances of developed imputation approach by considering various scenarios of covariance structures of longitudinal data and levels of censoring. We also illustrated the application of the proposed method to the Prospective Study of Outcomes in Ankylosing spondylitis (AS) (PSOAS) data to address the issues of censored or missing Creactive protein (CRP) level at early visits for a group of patients.
Results
Our findings from simulation studies indicated that the proposed method performs better than other MI methods by having a higher relative efficiency. We also found that our approach is not sensitive to the choice of covariance structure as compared to other methods that assume normality of biomarker data. The analysis results of PSOAS data from the imputed CRP levels based on our method suggested that higher CRP is significantly associated with radiographic damage, while those from other methods did not result in a significant association.
Conclusion
The MI based on weighted CQR offers a more valid statistical approach to evaluate a biomarker of disease in the presence of both issues with censoring and missing data in early visits.
Background
With advances in biotechnology, biological markers (biomarkers) continue to play an important role in an increasing number of biomedical studies. Biomarkers have led to a better understanding of the natural history and development of acute and chronic diseases, providing insights of the mechanism of treatment effects to identify and classify patients into different risk categories, and potential biological pathways that can be used to guide the therapeutic strategies for future treatment targets. Biomarker data are often measured over a period of time in biomedical studies to determine if the temporal changes differ between the patients who develop disease and those who do not. For example, Creactive protein (CRP) is one of the primary biomarkers known to reflect the degree of inflammation in the body and it has been widely used for studies of Ankylosing spondylitis (AS) to monitor disease activity, assess response to treatment and predict radiographic progression. However, in a longitudinal study, biomarker data may not be collected in certain time points for some patients. Furthermore, the biomarker data may be subject to censoring due to limits of detection (LOD). For example, in the Prospective Study of Outcomes in Ankylosing Spondylitis (PSOAS) [7], CRP data not only were censored due to limits of detection but also were incompletely collected at early study visits because blood sample collection was not a part of the original study design.
Compared to single imputation, it has been shown that modelbased imputation techniques such as multiple imputation (MI) methods can account for the uncertainty about the prediction of the unknown missing values [21] and provide more valid statistical inference (Lubin [15]). Likelihoodbased MI approaches have been proposed to address censoring issues due to limits of detection when the biomarker data are considered as covariates in a model. For example, Lee et al. [8] proposed MI for the multiplecensored correlated covariates based on the Gibbs sampling method. However, these methods focus on estimation of mean of biomarkers and assume normality of the distribution of biomarker data, which may not be valid as most biomarker data are highly skewed even after logtransformation. These limitations prompted development of alternative methods for nonnormal biomarker data. For example, Powell [19, 20] proposed using quantile regression models for censored data (i.e., censored quantile regression) with detection limit and it has been extended for longitudinal data using improved computational methods (Wang and Fygenson [27]; Lee and Kong [9]; Sun et al. [23]). As an important alternative to the mean regression models, quantile regression models are increasingly used in longitudinal study due to its robustness to nonnormality or heteroscedasticity and minimal assumption imposed on the quantiles of data. Estimating different quantiles should be of more practical interest especially in the presence of censoring issue, providing a broader picture of data distribution; specifically, it is common that some quantiles of biomarker data show significant effects that are not significant in other quantile. MI approaches that are based on censored quantile regression to handle censored covariates have been also proposed. For example, Wang and Feng [28] developed a multiple imputation approach for Mregression models. However, these approaches cannot handle multiple imputation in the presence of both censored and missing covariates in longitudinal data setting.
Lee and Kong [10] proposed an estimation approach based on censored quantile regression using the inverse probability weighting technique to handle longitudinal response variable with both censoring and monotone missingness [10] which is mainly caused by dropout; the basic concept of this method is that an individual’s contribution to the estimating equations is incorporated by the inverse probability weights for dealing with missing data at a dropout time. Since in PSOAS, CRP data for some patients were not completely collected at early study visits due to study design, this necessitated development of a new approach to handle missing data while controlling for censoring issue simultaneously.
In this paper, we propose a specific multiple imputation strategy that not only account for censoring, but also missing data at early visits when longitudinal biomarker data are modeled as a covariate. Assuming monotone missing pattern holds, we adopt the idea of Lee and Kong’s estimation method to establish weighted censored quantile models which are incorporated into our developed multiple imputation process. The focus here is assessing through simulation studies the performances of our multiple imputation approach by comparing relative efficiency of our method with that of complete case analysis and other traditional multiple imputation methods. We also illustrate application of the method to real life data from PSOAS to achieve realistic situations while specifically evaluating the association between CRP and radiographic damage.
Methods
The prospective study of outcomes in Ankylosing Spondylitis (PSOAS)
Ankylosing spondylitis is a chronic inflammatory disease characterized by inflammatory spinal pain usually beginning in the second to fourth decades of life that can result in chronic pain, that can result in functional impairment and diminished quality of life, and, in some patients, complete spinal fusion. In PSOAS, participants meeting the modified New York (mNY) Classification Criteria for AS [26] were enrolled from one of the five study sites (CedarsSinai Medical Center in Los Angeles, California, the University of Texas Health Science Center at Houston (UTH), the NIH Clinical Center, the University of California at San Francisco (UCSF), and the Princess Alexandra Hospital in Brisbane, Australia (PAH)^{Footnote 1} and were followed for up to 13 years (through two cycles of NIH funding: 2002–2006 and 2007–2016). At each study visit, spaced 6 months apart, the patients underwent comprehensive clinical evaluation for disease activity and functional impairment. Selfreported outcomes were measured at 6month intervals and radiographic data, including AP pelvis Xrays, AP and lateral lumbsacral spine films and lateral cervical spinal films were collected every 2 years in order to assess longitudinal radiographic damage which was defined by scoring the Bath Ankylosing Spondylitis Radiology Index (BASRI) [17] and the modified Stoke Ankylosing Spondylitis Spine Score (mSASSS) [3]. Creactive protein (CRP) levels and erythrocyte sedimentation rate (ESR) as well as medication usage were also determined at each clinical visit.
Analysis cohort
One of the objectives of PSOAS was to evaluate factors associated with longitudinal radiographic severity and rate of progression in AS patients. Specifically, we focused on evaluating longitudinal association between CRP level and radiographic damage which is assessed by mSASSS values (range 0–72) at each Xray visit. We considered analysis cohort of 295 patients who were confirmed AS by mNY criteria and had at least 4 years of radiologic follow up data (as of August 2016). However, we faced with two challenges in analyzing CRP data in relation to mSASSS. First, we found that 13.3% of CRP values were leftcensored due to being below the limit of detection. Another issue is related to unobserved CRP data during early visits for some patients which was by study design. Specifically, of the 295 patients with at least 4 years of radiologic follow up, 37% have been followed since first cycle of funding, i.e., Study I (enrolled from 2002–2006)^{Footnote 2} and then consented to Study II (enrolled from 2007–2016) for their continued participation. The inclusion of these patients increases the study power by increasing the number of patients who have been followed for 10+ years. However, CRP levels for 111 patients among our study cohort were not collected for up to first two consecutive Xray visits because initially blood sample collection (e.g., CRP and ESR) was not a part of the protocol for some patients in Study I. Different scenarios in which missing CRP data are generated are presented with detailed descriptions in Fig. 1. An important feature of CRP data was that it has a monotonic missing pattern; i.e., if a value was missing at visit t then the values for all previous visits (i.e., k,1≤k≤t−1) were also missing. This monotonic pattern was found in 92% of 111 patients. We also discovered that the number of patients who had CRP data missing at their first visit only (73.9%) was higher than the number of those who had missing data for both first and second visits (26.1%). There were only 5 patients who had CRP levels missing at all first three visits. Also we noted that patients who were enrolled earlier had a higher number of visits with unobserved CRP data.
Statistical approach
Our approach for assessing the longitudinal association between CRP and radiographic damage (i.e., mSASSS) includes four steps: Step (1) modeling missing data processes, Step (2) applying the inverse weighting techniques to censored quantile regression (CQR) using the probabilities of missing early visits that are estimated from the modeled missing data process, Step (3) employing multiple imputation process for both censored and missing CRP data at early visits, based on quantile estimates from established weighted CQR in Step (2), and Step (4) conducting longitudinal regression analyses where imputed CRP data are treated as an independent variable and mSASSS as a response variable. Natural logtransformed CPR was used in our analyses to reduce its highly rightskewed distributions.

Step (1)
Since the true probability of missing data is unknown, we modeled missing data process for each visit separately, through binary logistic regression with a response variable which indicates whether CRP data are observed or not (i.e., 1 = observed; 0 = missing) and independent variables that include observed CRP data at later times (i.e., the first available CRP) and other variables that were associated with missingness of CRP. Using the predicted probabilities (=η_{ it } for the ith patient and the tth visit) that were estimated from these models, we derived the patientlevel probability weights (=π_{ i }) under monotone missing data mechanism. Details regarding derivation of probability weights are shown in Appendix A.2. Although other reasons may lead to informative missing, we predicted the probability of missing CRP at early visits conditional on the observed CRP data at the rest of follow up visits (if censored, imputed by half of detection limit (DL), i.e., DL/2), and a set of covariates including study sites and the study group that patients were first enrolled in.

Step (2)
Censored quantile estimating equation incorporated by inverse probability weights was defined as a function of the variables that were significantly associated with CRP. Details regarding specific models and estimation procedure are presented in Appendices A.1 and A.2. In practice, it is important to use all available information to build the best imputation model [18, 21, 22]; we conducted the weighted CQR model based on the variables that include the covariates and the outcome of the potential analysis models even if they have limited predictive power. Once weighted CQR model was established, parameter estimates of different quantile levels were obtained by implementing function ‘crq’ in the R package quantreg for the existing optimization algorithms [4, 21]. Specifically, we used the option ‘Powell’ for method and ‘left’ for censoring type (i.e., ctype). It is well known that even if the missing data depend on the observed data, the weighted estimating equations provide unbiased estimation, when the missing data process is modeled with correctly specified probability [10, 12].

Step (3)
Missing CRP data were imputed by the uth conditional quantiles based on quantilespecific parameter estimates of aforementioned weighted CQR, where random variable U was generated from a uniform distribution between 0 and 1. We can estimate quantilespecific parameters using R function ‘crq’ with a function argument called ‘tau’, the quantile level at which the model is to be estimated. For censored CRP data, we first estimated the conditional probability of censoring, denoted by ω, using longitudinal logistic regression model with adjustment of potential predictors of censoring, such as study sites (because censoring rates varied over study sites from 2% to 29%), ESR levels, functional outcomes [2], disease activity [6] and medications usage. Then we used ω to randomly generate values of v from a uniform distribution between 0 and ω, which were used to impute the censored value by the vth conditional quantile. Since the conditional probability of censoring was estimated from logistic regression model and the imputations are obtained from a separated CQR, in few samples, the imputed values may not be less than a desired detection limit. When this situation occurs, we discarded the corresponding cases and used another v value, sampled from a uniform distribution between 0 and ω.

Step (4)
After Step (3) was repeated M=5 times to generate five imputed datasets, we conducted longitudinal regression analyses to evaluate association between CRP and dependent variable, mSASSS for each of imputed datasets. Details of analysis model are described in “Analysis of PSOAS data” section. To obtain the parameter estimates of interest, we defined the combined MI estimator as a mean of five estimates. Variance of MI estimators was determined based on 500 bootstrap samples, by resampling the observations with replacement and pvalues were calculated assuming the normality of estimated parameters. Additional details related to our imputation procedures are discussed in Appendix A.3.
Simulation studies
We conducted simulation studies to investigate the performance of our developed MI methods through different scenarios. We considered the following four different scenarios of longitudinal data structures, as well as different levels of censoring for generating biomarker data.

Scenario 1
Multivariate normal (MVN) distribution ; exchangeable covariance structure

Scenario 2
Multivariate normal (MVN) distribution ; unstructured covariance

Scenario 3
Multivariate exponential (MVE) distribution; exchangeable covariance

Scenario 4
Multivariate exponential (MVE) distribution; heteroscedastic covariance structure (i.e., covariance depends on a set of covariates)
In order to generate a longitudinal outcome variable that mimics the distribution of mSASSS (=y_{ it }) in PSOAS, we used the following regression model
for the ith patient and the tth visit, where α_{0}=5, α_{1}=−4, α_{2}=−6, w_{ it } represents the longitudinal structure variable time (t=1, ⋯,4) and \(z^{*}_{it}\) denotes complete biomarker data which have been generated based on the aforementioned four scenarios. An error term ε_{ it } was generated from multivariable normal distribution based on exchangeable covariance structure with a correlation coefficient ρ of 0.3. We then produced missing and censored values for biomarker data z_{ it } that mimic the missing data pattern of CRP levels in PSOAS. We used a logistic regression to model probability of observed biomarker data η_{ it }, based on variables w_{ it }, y_{ it } and the first observed biomarker data after time t (i.e., \(z_{it'}, t<t'\leq 4\phantom {\dot {i}\!}\)). Based on probability η_{ it }, we calculated the patientlevel probability weights (=π_{ i }) under monotone missing data mechanism through a specific function of η_{ ij }, as shown in Appendix B. For generating censored data, we chose the detection limit c, as the (100×r)th percentile of the simulated biomarker data, where r is the censoring rate (i.e., r=0.1, 0.15, 0.2, 0.3). We simulated data for 75% of patients who had missing data for up to first 3 visits (i.e., missing at first visit only, first and second visits, or all first three visits), and 25% of patients who had complete measurements up to visit 4. For each scenario, three hundred simulation datasets with sample size of 250 were generated. Details regarding parameters used for data generation and covariance structures for each scenario are described in Appendix B.
For multiple imputation, we fitted weighted CQR (wCQR) models using both known probability weights π (MIwCQR _{1}) and \(\hat {\pi }\) that was estimated through the aforementioned logistic regression model. Since the observed biomarker data z_{i,t+j} in this logistic regression model had also censored values, we considered fitting two separate wCQR models, one based on imputed censored data by DL/2 (MIwCQR _{2}) and the other using uncensored data only (MIwCQR _{3}), in order to see the impact of these two approaches on parameter estimates. We also considered unweighted CQR method (MICQR) which accounts for censoring but ignores the missing data mechanism. Other imputation methods were further applied, that included Markov Chain Monte Carlo (MCMC)based MI methods [11, 14] through on Bayesian frame work for a monotone missing data, which were implemented in the function ‘mice’ of the R package mice [25], imputing only missing values (MIMCMC _{1}), as well as imputing both censored and missing values (MIMCMC _{2}). MCMCMI algorithm obtains the posterior distribution of parameters by sampling iteratively from conditional distributions based on Gibbs sampling method.
Using imputed biomarker data generated from different MI methods described above, we conducted longitudinal regression analysis using model 1, for each scenario. For assessing the performance of each estimator, we calculated bias and ratio of the mean squared error (MSE) of the omniscient estimator (OMNI), the gold standard, which is based on the complete data without censoring, to that of each estimator. Throughout we refer to this ratio of MSEs as relative efficiency (RE), which will be used for comparing the performance of the aforementioned methods. Moreover, we also conducted complete case analysis using only observed biomarker data where censored data were imputed by a single value of DL/2 (CCDL/2).
Analysis of PSOAS data
Censored or missing CRP data at early visits in PSOAS were imputed based on six different approaches (CCDL/2, MIMCMC _{1}, MIMCMC _{2}, MICQR, MIwCQR _{2}, MIwCQR _{3}) that are described in “Simulation studies” section. For each imputed dataset, we assessed the longitudinal association between natural logtransformed CRP levels and mSASSS while controlling for potential confounding factors, that included Bath Ankylosing Spondylitis Disease Activity Index (BASDAI), medication usages of Tumor Necrosis Factors inhibitors (TNFi), and Nonsteroidal AntiInflammatory drugs (NSAIDs) as well as demographic information such as sex, race, disease duration, comorbidity, education and smoking status. Multivariable mixed effect Poisson regression models were conducted for each imputed dataset, to account for the correlations of repeated measures within a patient.
Results
Simulation study results
Table 1 presents the results of simulation study across different censoring rates based on aforementioned Scenario 1; Bias and relative efficiency (100 × RE) are presented for each parameter, α_{0}, α_{1}, α_{2}. Overall, our proposed three MIwCQR approaches (i.e., MIwCQR _{1}, MIwCQR _{2}, MIwCQR _{3}) produced more efficient estimators (i.e., higher REs) than other MI methods that were used for comparison. Specifically, with 10% censored data, RE of our MI methods for α_{1}, the coefficient of biomarker, ranged from 53.4% to 53.8%, while for CCDL/2 RE was 3.2% and for that of MIMCMC _{1}, MIMCMC _{2} and MICQR were 5.0, 39.5, and 49.6%, receptively. Although as the censoring rate increased the magnitude of RE for our methods decreased slightly, we still observed higher REs, ranging from 40.1% to 44.9% for α_{1}, compared to other methods with REs ranging from 8% for CCDL/2 to 39.98% for MICQR, when 30% of data were censored. We also observed similar patterns in REs for other two coefficients, α_{0} and α_{2}.
Similar findings were observed under Scenario 2Scenario 4, as shown in Table 2. Our MI methods provided higher REs compared to the other methods across all these three scenarios in the presence of 30% censoring. It also demonstrates that our weighted CQR methods is not sensitive to the choice of covariance structure, as compared to MCMCbased methods that assume normality of biomarker data. For example, RE for MIMCMC _{2} under Scenario 3 (MVE) was about 49% lower than that of Scenario 2 (MVN) (i.e., from 21.52 for Scenario 2 to 11.07 for Scenario 3), while our methods provided consistent REs (< 0.5% change) over all three scenarios.
Figure 2 displays distribution of biomarker data from one of simulated datasets, distinguishing the observed data from imputed data by our MIwCQR _{2} method; the distribution of data after imputation was very similar to the complete data distribution that were originally simulated. Similar findings were observed for our other methods, MIwCQR _{1} and MIwCQR _{3} (Figures are not shown).
Results of applying the proposed methods to PSOAS data
General characteristics of patients in PSOAS
Among 295 AS patients who had at least 4 years of radiologic follow up, the mean follow up time was 6.49 years (standard deviation (SD) = 2.37) with maximum years of 13.5 and mean number of Xray visits was 3.6 (SD = 1.2). The cohort was 76.3% male, 81% white, 8.8% Hispanic, with a mean age 42.6 years (SD = 13.1) and a mean disease duration 18.0 years (SD = 12.7) at baseline. Of the 295 patients, 54 (18.3%) were from UCSF, 83 (28.1%) from UTH, 106 (35.9%) from Cedars Sinai, 42 (14.2%) from the NIH Clinical Center and 10 (3.4%) were from PAH. At baseline visit, 71.5% of patients had at least one comorbidity, 41.5% were eversmokers and 10.7% were current smokers. A median mSASSS at baseline visit was 5 (interquartile range (IQR) = [0, 24]), the first observed CRP level with censored values imputed by DL/2 had a median of 0.4 (IQR = [0.2, 0.8]) and patients with mSASSS ≥ 4 had higher median CRP level compared to those with mSASSS < 4 (0.43 vs. 0.31).
Analysis results
Table 3 shows the adjusted rate ratios (RR) and pvalues of complete case analysis using censored CRP data imputed by DL/2 (CCDL/2), and those from imputed CRP levels by three other methods: MIMCMC2, MICQR and MIwCQR _{2}. The results from MIMCMC _{2} and MIwCQR _{2} were very similar to those from MIMCMC _{1} and MIwCQR _{3} respectively (data not shown). This may be because the censoring rate of CRP in PSOAS data is not high enough to cause differences between these methods. However, there were noticeable differences in the estimates and the corresponding pvalues for CRP across these four methods. The results from our method (MIwCQR _{2}) suggest that higher CRP is significantly associated with radiographic damage (adjusted RR = 1.018; 95% confidence interval (CI) = [1.004, 1.031]; p=0.0095), while the other methods did not result in a statistically significant association (p=0.99 for CCDL/2; p=0.56 for MIMCMC _{2}; p=0.18 for MICQR).
Discussion and conclusion
Biomarker data are often subject to left censoring due to inability to obtain complete data when the measurements are below the limit of detection. In longitudinal studies, it is also possible that biomarker data are not completely collected during early study visits which introduces a monotonic missing data pattern. Both likelihood based joint modeling techniques ([5, 16, 24]) and quantile regression approaches ([10, 12, 29]) have been used to deal with monotonic missing data. However, most of these methodological developments have dealt with monotone missing data caused by termination from a trial or study (e.g., dropout), to our knowledge there are no published studies that have developed imputation approaches that specifically accommodate both censoring and missing data at early follow up visits. In this article we have developed the use of multiple imputation procedure that is based on weighted censored quantile regression model to account for both leftcensored and monotone missing biomarker data during early visits. Specifically, we applied inverse probability weighting techniques to incorporate missing data in early visits through a multiple imputation based on censored quantile regression.
Our findings from the simulation study indicate that our proposed method performs better than other MI methods as assessed by higher RE. Further, our approach is not sensitive to the choice of covariance structure as compared to other methods that assume normality of biomarker data. The results of our method MIwCQR _{2}, where missing data probability weights were estimated based on the imputed censored data by DL/2, were similar to those of MIwCQR _{1}, where the true probability weights were used. This is reassuring to use estimated probabilities for missing data based on the logistic regression model for missing data process, as in real life applications we may not have information about the true probabilities for missing data. However, when the uncensored data were only used for the missing data model (MIwCQR _{3}), the results were not as good as the ones from MIwCQR _{1} and MIwCQR _{2} (i.e., lower RE). Since missing data model is defined by linear covariate effects, uncensored data can be used to obtain a consistent estimate of probability of missing, but it may introduce some bias when there is a strong underlying nonlinear relationship or the censoring rate is considered high (e.g., > 30%). This is consistent with literature that indicate accurate estimates of probability is critical when the inverse weighted probability based approach is applied [27, 29].
Moreover, we demonstrated application of our methods to real data from the PSOAS cohort by examining the longitudinal association between CRP levels and radiographic damage in a situation where CRP levels for some patients were either not collected in the early visits or leftcensored due to the detection limit. The results from our method indicated that higher CRP is significantly associated with radiographic damage, while the other methods did not result in a significant association. This finding is also consistent with clinical expectation that CRP is associated with radiographic severity in patients with AS [1].
Though developed method could be implemented with standard software package quantreg in R that fits quantile regression, we are currently developing R package for users to easily implement the proposed MI approaches. Censored quantile regression has been extended to data censored at both lower and upper thresholds [4], therefore our method can be also directly extended to doubly censored biomarker data. There is a growing interest in developing MI methods that impute missing data across multiple medications while accounting for the correlations among them, which can be also extended by our proposed method.
Based on our earlier work (Lee and Kong [10]), we expect our method to provide consistent estimates for the parameters in the weighted quantile regression model, assuming the missing data model is correctly specified [10, 12]. Despite aforementioned advantages for our method shown earlier, we acknowledge that the misspecification of the model for missing data process may introduce bias in the estimation of parameters. Therefore, it is important to identify the model carefully and interpret the analysis results cautiously [8, 10].
Appendix A: Model formulation/ Estimation procedure
A.1 Weighted censored quantile regression model accounting for missing early visits
Let \(z^{*}_{it}\) be the biomarker measurement for the ith subject at time t assuming all subjects are to be observed at the same time. Suppose we define the linear regression model \(z^{*}_{it}= \boldsymbol x_{it}^{T}\boldsymbol \beta +e_{it},\quad i=1,\cdots, n; \quad t=1,\cdots, m \label {latent}\), where x_{ it } is a p×1 vector of covariates that can include the time of measurement, β is an unknown p×1 vector of regression parameters and the random errors e_{ it } are correlated within the subject to reflect the serial correlations of repeated measurements within each individual. If the τth conditional quantile of e_{ it } given \(\boldsymbol x_{it}^{T}\) is assumed to be zero, a quantile regression model related to the τth quantile of response variable, \(q_{\tau }\left (z^{*}_{it}\right)\), conditional on x_{ it } has the form \(q_{\tau }\left (z^{*}_{it}\right)=\boldsymbol x_{it}^{T}\boldsymbol \beta _{\tau }, \quad 0<\tau <1\), where β_{ τ } is a vector of quantile specific regression parameters corresponding to the coefficient β in the linear regression model above. When there exists a lower detection limit, say c, \(z^{*}_{it}\) is a latent variable and we cannot observe the biomarker measurement if it has a value below c and we only observe \(z_{it}=z^{*}_{it}\), if \(z^{*}_{it}>c\). This leads to the longitudinal censored quantile regression (CQR) model defined as \(z_{it}= max\left (c, \boldsymbol x_{it}^{T}\boldsymbol \beta +e_{it}\right)\). We can define the objective function for longitudinal censored data as
The loss function ρ_{ τ }(u)=u{τ−I(u≤0)}, with I(·) being an indicator function, represents the contribution by residuals. The estimates resulting from (2) are equivalent to the solution of estimating equation
To apply the weighting techniques to the censored quantile regression model for handling missing data at early visits, let O_{ i } be a random variable indicating the time point when the data collection was started for the ith subject. O_{ i } can take the values between 1 and m. If the subject has completed 1 ∼m followup visits then O_{ i }=1, and if the subject had missing data from visit 1 to m1 then O_{ i }=m. We denote \(\boldsymbol z_{i}^{o}\) as the observed response history since the data collection was started, and X_{ i }={x_{i1},⋯,x_{ im }}^{T} as a set of covariates that were observed from the complete study visits 1 ∼m. When the biomarker measurements are MAR, the conditional probability of missing early visit from the baseline to the o_{ i }−1 occasion for the ith subject is \(\pi _{io_{i}}=Pr\{O_{i}=o_{i}\boldsymbol z_{i}^{o}\), X_{ i },γ} (o_{ i }=1,⋯,m), where \(\pi _{io_{i}}>0\) and γ is a parameter vector of the regression model. Now the weighted estimating equations for censored quantile regression model can be defined as
The basic idea of weighted estimating equations is to weight each subject’s contribution by the inverse probability of missing early visits to a given occasion. After we define \(\boldsymbol x_{it}^{w}=\frac {1}{\pi _{io_{i}}}\boldsymbol x_{it}, z_{it}^{w}=\frac {1}{\pi _{io_{i}}}z_{it}\) and \( {c_{i}}^{w}=\pi _{io_{i}}^{1} {c}\), Eq. (4) can be written in the same form as the unweighted estimating Eq. (3) as follows:
Thus, the corresponding objective function is in the form of
Now the traditional censored quantile regression estimation algorithm ([4, 20]) can be straightly applied to minimize this objective function. Details related to estimation procedures, inference, and asymptotic properties of parameter estimators were discussed in Lee and Kong [10].
A.2 Missing data process
If the missing early visit data arise from the MAR mechanism, estimation of probability of missing early visits is straightforward. To illustrate the missing data process, denote R_{ it } as the missing status of response variable z_{ it }, i.e., R_{ it }=1 if z_{ it } is observed and 0 otherwise. Then R_{ ij }=1 implies that \(\phantom {\dot {i}\!}R_{ij'}=1\) for all j^{′}>j given the monotone missing pattern. To indicate when the data collection is started, we define a random variable O_{ i } as \(O_{i}=1+(m\sum _{t=1}^{m} R_{it})\). The probability of missing early visit \(\pi _{io_{i}}\) from baseline to occasion o_{ i }−1 can be given by
where \(z_{io_{i}'}^{o}\) is the first observed z value after time o_{ i } (i.e., o_{ i }<oi′≤m). When we define the probability of being observed at time t for the ith subject as \(\eta _{it}=Pr\left (R_{it}=1 R_{i,t+1}= \cdots = R_{im}=1, z_{it'}^{o}, \boldsymbol X_{i}, \boldsymbol \gamma \right)\), t<t^{′}≤m, we can carry out the probability of missing early visits in terms of η_{ it } as follows:
where γ is the parameter vector of the regression model for η_{ it }. Then appropriate regression models such as logistic regression model can be used to estimate η_{ it }, and then we can calculate \(\pi _{io_{i}}\) based on the equation above.
A.3 Multiple imputation process based on weighted censored quantile regression
Multiple imputation techniques [21] have been widely used for the general handling of missing data. However, the censoring and monotone missing mechanisms should be incorporated in the imputation models to deal with the complexity of missingness. Based on the weighted censored regression of quantiles that was introduced in the previous section, we propose the multiple imputation procedure to fill in the data that are leftcensored or missing at early visits. The conditional censoring probability of z_{ it }, ω(X_{ it })=Pr(z_{ it }<cX_{ it }) can be estimated using a logistic regression model \(logit[\omega (\boldsymbol X_{it}) ]=\boldsymbol X_{it}^{T}\boldsymbol \delta \), where δ is unknown parameter vector and then v is sampled from uniform distribution UNIF (0,ω(X_{ it })) in order to impute the censored value z_{ it } by its conditional quantile of \({z}_{it}^{*}=\boldsymbol X_{it}^{T} \hat {\boldsymbol \beta }_{v}\) which is estimated through fitting weighted censored quantile regression model where z_{ it } is treated as the dependent variable as described in the previous section. We also draw u from UNIF(0,1) to fill in the missing value by \({z}_{it}^{*}=\boldsymbol X_{it}^{T} \hat {\boldsymbol {\beta }}_{u}\), that is the uth conditional quantile of z_{ it } given X_{ it }. Once the imputed datasets are generated, any analysis designed for the complete dataset can be applied to each of M imputed datasets. To obtain the parameter estimators of interest in the regression model \(y_{it}= \alpha _{0}+\alpha _{1} z^{*}_{it}+\alpha _{2} w_{it}+\epsilon _{it}\), we define the combined MI estimator as \(\hat {{\boldsymbol {\alpha }}}_{MI}=M^{1}\sum _{k=1}^{M} \hat {{\boldsymbol {\alpha }}}_{k}\). Wang and Feng [28] discussed the asymptotic properties of their proposed multiple imputation procedure based on the conditional quantile function and suggested bootstrapping because the asymptotic variance of MI estimators takes complex forms and it is difficult to estimate directly. We adopted a bootstrap method by resampling the paired observations with replacement based on 500 bootstrap samples to obtain the standard errors for estimated parameters and pvalues that were calculated by using the normality of estimated parameter \(\hat {{\boldsymbol {\alpha }}}_{MI}\).
Appendix B: Simulation study design
Suppose the subjects are to be observed at the same m time points. The latent longitudinal biomarker data are generated from the model
where the covariates include a variable x_{ it } with Poisson(20) distribution and w_{ it } representing the tth assessment time which is set equal to t. Given the covariates, random error vectors, e_{ i }=(e_{i1},⋯,e_{ im })^{T} for i=1,⋯,n, are assumed to be mutually independent and have conditional τth quantile equal to zero. Let us consider an error term \(e_{it}F^{1}_{e_{it}}(\tau)\), where F(·) denotes the cumulative distribution function and \(F^{1}_{e_{it}}(\tau)\) is the τth quantile of e_{ it } given x_{ i } and t. We simulated the random variable e_{ it } from each of the following distributions and calculated \(F^{1}_{e_{it}}(\tau)\) with τ=0.5 for each scenario.

Scenario 1
Multivariate normal distribution (MVN); exchangeable covariance structure:
e_{ i }={e_{i1},⋯,e_{i4}}^{T}∼MVN(0,σ^{2}R), where σ^{2}=1 and correlation matrix R is exchangeable with ρ=0.3
\( \boldsymbol R=\left (\begin {array}{cccc} 1.0 & 0.3 & 0.3 & 0.3\\ & 1.0 & 0.3 & 0.3 \\ & & 1.0 & 0.3 \\ & & & 1.0 \\ \end {array}\right).\)

Scenario 2
Multivariate normal (MVN) distribution ; unstructured covariance:
e_{ i }={e_{i1},⋯,e_{i4}}^{T}∼MVN(0,σ^{2}R), where σ^{2}=1 and correlation matrix
\( \boldsymbol R=\left (\begin {array}{cccc} 1.00 & 0.75 & 0.44 & 0.54\\ & 1.00 & 0.37 & 0.46 \\ & & 1.00 & 0.08 \\ & & & 1.00 \\ \end {array}\right).\)

Scenario 3
Multivariate exponential distribution; exchangeable covariance:
e_{ it }=exp(ξ_{ it })−1 and ξ_{ i }={ξ_{i1},⋯,ξ_{i4}}^{T}∼MVN(0,σ^{2}R), where σ^{2}=1 and R is exchangeable with ρ=0.3.

Scenario 4
Multivariate exponential distribution; heteroscedastic covariance structure (i.e., covariance depends on a set of covariates X_{ i }:
e_{ it }=exp(ξ_{ it })−1 and ξ_{ i }={ξ_{i1},⋯,ξ_{i4}}^{T}∼MVN(0,1/(1+x_{i1})R), where R is exchangeable with ρ=0.3. Note that in this case, \(F^{1}_{e_{it}}(\tau)\) varies with x_{i1}.
We set β=(β_{0},β_{1},β_{2})^{T}=(2.3,−0.25,−0.1)^{T}, m=4 and overall censoring percentage r=10,15,20 or 30%. We chose the detection limit c as the (100 ×r)th sample percentile of the simulated biomarker data \(z^{*}_{it}\). Using latent variable \(z^{*}_{it}\), we finally generated \(y_{it}= \alpha _{0}+\alpha _{1} z^{*}_{it}+\alpha _{2} w_{it}+\epsilon _{it}\), (α_{0}=5,α_{1}=−4,α_{2}=−6). As for the missing data process, the logistic regression model below was postulated,
where the parameter vector is α=(γ_{0},γ_{1},γ_{2})^{T}=(1,−8.5,0.5)^{T}, η_{ it } is the conditional probability of being observed at time t, and \(z^{*}_{i,t+1}\) is the observed biomarker data at the time point t+1. Under this setting, we assumed the subjects with higher level of marker z^{∗} are more likely to have missing data.
Notes
 1.
Participants from PAH have been enrolled since 2007.
 2.
Study IA: 5year funded study (enrolled from 2002–2006) for AS patients with disease symptom duration of > 20 years. Patients were initially enrolled for one visit but the protocol was amended to include the second follow up visit about 2 to 3 years after their initial enrollment; Study IB: 2year longitudinal study (enrolled from 2003–2006) for AS patients with disease symptom duration of < 20 years.
Abbreviations
 AS:

Ankylosing spondylitis
 BASRI:

Bath Ankylosing Spondylitis radiology index
 CC:

Complete case
 CQR:

Censored quantile regression
 CRP:

Creactive protein
 DL:

Detection limit
 ESR:

Erythrocyte sedimentation rate
 LOD:

Limits of detection
 MCMC:

Markov Chain Monte Carlo
 MI:

Multiple imputation
 mNY:

Modified New York
 mSASSS:

Modified Stoke Ankylosing Spondylitis spine score
 MSE:

Mean squared error
 MVE:

Multivariate exponential
 MVN:

Multivariate normal
 NSAIDs:

Nonsteroidal AntiInfammatory drugs
 OMNI:

Omniscient
 PAH:

Princess Alexandra Hospital in Brisbane, Australia
 PSOAS:

Prospective study of outcomes in Ankylosing spondylitis
 RE:

Relative efficiency
 RR:

Rate ratio
 TNFi:

Tumor Necrosis factors inhibitors
 UCSF:

University of California at San Francisco
 UTH:

University of Texas health science center at Houston
 wCQR:

Weighted CQR
References
 1
Braun J, Baraliakos X, Hermann KA, Xu S, Hsu B. Serum Creactive Protein Levels Demonstrate Predictive Value for Radiographic and Magnetic Resonance Imaging Outcomes in Patients with Active Ankylosing Spondylitis Treated with Golimumab. J Rheumatol. 2016. doi:https://doi.org/10.3899/jrheum.160003.
 2
Calin A, Farrett S, Whitelock H, Kennedy LG, OHea J, Mallorie P, Jenkinson T. A new approach to defining functional ability in ankylosing spondylitis: the development of the Bath Ankylosing Spondyitis Functional Index. J Rheumatol. 1994; 21:2281–5.
 3
Creemers MC, Franssen MJ, vant Hof MA, Gribnau FW, van de Putte LB, van Riel PL. Assessment of outcome in ankylosing spondylitis: an extended radiographic scoring system. Ann Rheum Dis. 2005; 64:1279.
 4
Fitzenberger B. A guide to censored quantile regressions In: Maddala GS, Rao CR, editors. Handbook of Statistics, Volume 15 (Robust Inference). Amsterdam: Elsevier Science; 1997. p. 405–37.
 5
Gao S, Thiebaut R. Mixedeffect Models for Truncated Longitudinal Outcomes with Nonignorable Missing Data. J Data Sci. 2009; 7(1):13–25.
 6
Garrett S, Jenkinson T, Kennedy LG, Whitelock H, Gaisford P, Calin A. A new approach to defining disease status in ankylosing spondylitis: the Bath Ankylosing Spondylitis Disease Activity Index. J Rheumatol. 1994; 21:2286–91.
 7
Gensler L, Ward MM, Reveille JD, Weisman MH, Davis Jr JC. Clinical, radiographic and functional differences between juvenileonset and adultonset ankylosing spondylitis: results from the PSOAS cohort. Ann Rheum Dis. 2008; 67(2):233–7.
 8
Lee M, Kong L, Weissfeld L. Multiple Imputation For LeftCensored Biomarker Data Based On Gibbs Sampling Method. Stat Med. 2012; 31(17):1838–48.
 9
Lee M, Kong L. Median Regression for Longitudinal Leftcensored Biomarker Data subject to Detection Limit. J Stat Biopharm Res. 2011; 3(2):363–71.
 10
Lee M, Kong L. Quantile Regression For Longitudinal Biomarker Data Subject to Left Censoring and Dropouts. Commun Stat Theory Methods. 2014; 43(21):4628–41.
 11
Li KH. Imputation Using Markov Chains. J Stat Comput Simul. 1988; 30:5779.
 12
Lipsitz SR, Fitzmaurice GM, Molenberghs G, Zhao LP. Quantile regression methods for longitudinal data with dropouts: Application to CD4 cell counts of patients infected with the human immunodeficiency virus. J R Stat Soc: Ser C: Appl Stat. 1997; 46(4):463–476.
 13
Little RJA, Rubin DB. Statistical analysis with missing data, 2nd edition. New York: Wiley; 2002.
 14
Liu M, Wei L, Zhang J. Review of guidelines and literature for handling missing data in longitudinal clinical trials with a case study. Pharm Stat. 2006; 5:718.
 15
Lubin JH, Colt JS, Camann D, Davis S, Cerhan JR, Severson RK, Bernstein L, Hartge P. Epidemiologic evaluation of measurement data in the presence of detection limits. Environ Health Perspect. 2004; 112(17):1691–6.
 16
Lyles RH, Lyles CM, Taylor DJ. Random regression models for human immunodeficiency virus ribonucleic acid data subject to left censoring and informative dropouts. J R Stat Soc: Ser C. 2000; 49(4):485–97.
 17
Mackay K, Mack C, Brophy S, Calin A. The Bath Ankylosing Spondylitis Radiology Index (BASRI): a new, validated approach to disease assessment. Arthritis Rheum. 1998; 41:2263–70.
 18
Meng X. Multipleimputation inferences with uncongenial sources of input. Stat Sci. 1994; 9(4):538–73.
 19
Powell JL. Least Absolute Deviations Estimation for the Censored Regression Model. J Econ. 1984; 25(3):303–25.
 20
Powell JL. Censored regression quantiles. J Econ. 1986; 32(1):143–55.
 21
Rubin DB. Multiple Imputation for Nonresponse in Surveys. New York: Wiley; 1987.
 22
Rubin DB. Multiple imputation after 18 + years. J Am Stat Assoc. 1996; 91(434):473–89.
 23
Sun X, Peng L, Manatunga A, Marcus M. Quantile regression analysis of censored longitudinal data with irregular outcomedependent followup. Biometrics. 2016; 72:64–73.
 24
Thiebaut R, JacqminGadda H, Babiker A, Commenges D. Joint modelling of bivariate longitudinal data with informative dropout and leftcensoring, with application to the evolution of CD4+ cell count and HIV RNA viral load in response to treatment of HIV infection. Stat Med. 2005; 24(1):65–82.
 25
van Buuren S, GroothuisOudshoorn K. mice: Multivariate Imputation by Chained Equations in R. J Stat Softw. 2011; 45(3):1–67.
 26
van der Linden S, Valkenburg H, Cats A. Evaluation of diagnostic criteria for ankylosing spondylitis. A proposal for modification of the New York criteria. Arthritis Rheum. 1984; 27:361–8.
 27
Wang JH, Fygenson M. Inference for Censored Quantile Regression Models in Longitudinal studies. Ann Stat. 2009; 37(2):756–81.
 28
Wang JH, Feng X. Multiple Imputation for Mregression with Censored Covariates. J Am Stat Assoc. 2012; 107(497):194–204.
 29
Yi GY, He W. Median Regression Models for Longitudinal Data with Dropouts. Biometrics. 2009; 65(2):618–25.
Acknowledgments
The authors are grateful to the editors and referees for their thoughtful comments and constructive suggestions that have led to a considerable improvement of the earlier version. We acknowledge the support provided by the Biostatistics/ Epidemiology/ Research Design (BERD) component of the Center for Clinical and Translational Sciences (CCTS) for this project. CCTS is mainly funded by the NIH Centers for Translational Science Award (UL1 TR000371) by the National Center for Advancing Translational Sciences (NCATS). The authors also would like to recognize that this study was presented at 2016 Joint Statistical Meetings (JSM) held on July 30  August 4, 2016 in Chicago, IL.
Funding
This work was supported by grants from the United States Department of Health and Human Services, National Institutes of Health (NIH), National Institute of Arthritis and Musculoskeletal and Skin Diseases (NIAMS) (P01AR05291506) and the Spondylitis Association of America.
Availability of data and materials
The datasets generated and analyzed during the current study are not publicly available due to ongoing process of PSOAS data collection and quality assurance. However, the datasets are available from the corresponding author on reasonable request.
Author information
Affiliations
Contributions
ML made substantial contributions to conception and design, implemented the simulation study, performed statistical analyses, and was a major contributor in writing the manuscript. MHR commented and provided edits on the manuscript at all stages. JDR and MB edited the manuscript. LG and MW reviewed the findings of PSOAS data analysis. LD contributed to the manuscript drafting process, specifically for illustrations of PSOAS cohort. All authors read and approved the final manuscript.
Corresponding authors
Correspondence to MinJae Lee or Mohammad H. Rahbar.
Ethics declarations
Ethics approval and consent to participate
This study received approval from the institutional review boards of each institution and the respective institutional ethics boards approved the study [CedarsSinai IRB; Committee on Human Research at UCSF; Princess Alexandra Hospital (PAH) Human Research Ethics Committee; Committee for the Protection of Human Subjects at the University of Texas at Houston; National Institutes of Health (NIH)]. All study patients provided written informed consent.
Consent for publication
Not applicable (no individual person’s data contained in the manuscript).
Competing interests
No potential conflict of interest was reported by the authors.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Lee, M., Rahbar, M.H., Brown, M. et al. A multiple imputation method based on weighted quantile regression models for longitudinal censored biomarker data with missing values at early visits. BMC Med Res Methodol 18, 8 (2018). https://doi.org/10.1186/s1287401704639
Received:
Accepted:
Published:
Keywords
 Limit of detection
 Leftcensoring
 Missing early visits
 Quantile regression
 Multiple imputation