Rheumatoid arthritis is a chronic inflammatory disease that affects the joints. Disease activity leads to pain, swelling, and joint damage. The cause of rheumatoid arthritis is unknown but it is known to be an autoimmune disease. Disease activity score based on 28 joints (DAS28) is a composite score measuring disease activity that integrates measures of physical examination (number of tender and swollen joints), erythrocyte sedimentation rate (a hematology test that is a non-specific measure of inflammation), and patient reported outcome (patient self reported global health on a visual analogue scale).
DAS28 scores can range from 0.49 to 9.07. A score greater than 5.1 is considered high disease activity, 3.2-5.1 moderate disease activity, 2.6-3.2 low disease activity, and a score less than 2.6 is considered to be clinical remission. The goal of treatment is to achieve the lowest possible level of disease activity and clinical remission if possible. Figure 1 displays the 28 joints considered in the calculation of DAS28 scores.
Clinical trials are strategically designed to collect the disease activity score of each patient over multiple clinical visits, meanwhile a patient may drop out before their intended completion due to various reasons. The dropout terminates the longitudinal data collection on the patient’s disease activity score. In the presence of informative dropout (Rubin, 1976; Little and Rubin, 1987), that is, the dropout depends on latent variables from the longitudinal process, simply applying mixed effects models to analyze longitudinal outcomes may lead to biased results because the assumption of random dropout is violated (Laird and Ware, 1982). It would be more appropriate to jointly model the DAS28 scores and the dropout to account for the dependency between the two processes.
Several methods have been proposed in the literature for jointly modeling longitudinal and time-to-event data, as reviewed by Hogan and Laird (1997); Tsiatis and Davidian (2004); Yu et al. (2004); Ibrahim et al. (2005); Hanson et al. (2011). Joint models have been widely used in HIV/AIDS clinical trials, cancer vaccine (immunotherepy) trials and quality of life studies. References include Schluchter (1992); Pawitan and Self (1993); De Gruttola and Tu (1994); Tsiatis et al. (1995); Faucett and Thomas (1996); Wulfsohn and Tsiatis (1997); Faucett et al. (1998); Henderson et al. (2000); Wang and Taylor (2001); Xu and Zeger (2001); Law et al. (2002); Song et al. (2002); Chen et al. (2002, 2004); R Brown and G Ibrahim (2003); Brown and Ibrahim (2003); Brown et al. (2005); Chi and Ibrahim (2006, 2007); Zhu et al. (2012), and many others.
Our approach combines different methodologies to deal with the complexity of the data. The joint model has two parts: the first is a longitudinal random change point model for repeated measures DAS28 scores; the second part is a time to event hazard risk log normal model that models two competing risk events. Both parts are linked through a patient specific random effect. Similar approaches are described next. Jacqmin-Gadda et al. (2006) propose a joint model for repeated measurements (cognitive test score) and time to event (dementia). The part of their model for repeated measures, as ours, assumes a linear trend before a random change point and, in contrast to ours, assumes a polynomial trend after it. They assume a lognormal model for the time to event as we do for our competing risk model. Faucett et al. (2002) apply a similar random effects model for the AIDS trial along with a time-dependent proportional hazards model for the time to AIDS. In their random effects longitudinal model they assume, as we do, two linear trends before and after a random change point. Although, as ours, their model contains a subject-specific random intercept, in contrast to ours, all subjects have a common transition point. Pauler and Finkelstein (2002) adopts a similar random effect for the longitudinal PSA trajectories and Cox proportional hazard model for the progression times. By comparing the data pattern of people who relapse and who do not, they add a strong constraint to the random effects model with a new indicator – the slope after the transition point has to be greater than the slope before the transition point. For event-time data, Cox proportional hazards model with a piecewise exponential model is applied to increase flexibility.
Currently, the most commonly used method for jointly modeling lontidudinal and time-to-event data are likelihood based joint models. Likelihood based joint models are classified as either selection models or pattern mixture models(Little, 1995)
. Selection models factor the joint distribution through the marginal distribution model of the longitudinal outcome and the conditional distribution of time-to-event given longitudinal latent variables. Pattern-mixture models(Glynn et al., 1986; Little, 1993) factor the joint distribution differently by modeling the longitudinal data stratified by time-to-event (e.g., time-to-dropout) patterns. Our interest in this paper is towards marginal inference on the longitudinal profiles. This is, we are interested in assessing the treatment effects while adjusting for the imbalanced drop out. For this reason we consider selection models over pattern-mixture models. We also consider a competing risk survival model for informative dropout because, in practice, a patient may drop out before their intended completion date due to various competing risks. Joint analysis of longitudinal measurements and competing risks failure time data has been studied recently both in terms of frequentist approaches (Elashoff et al., 2007, 2008) and Bayesian approaches (Huang et al., 2010)
. Such approaches of joint modeling can make the estimates for the longitudinal marginal profile less biased. Therefore we can draw more reliable conclusions from clinical trials with longitudinal outcomes.
The motivation comes from a randomized double-blinded clinical Trial of Etanercept and Methotrexate with radiographic Patient Outcomes (TEMPO, Keystone et al., 2009). The primary objective of this trial is to evaluate the efficacy and safety of a combination therapy of methotrexate and etanercept for the treatment of rhuematoid arthritis versus methotrexate alone and etanercept alone. Methotrexate, a disease-modifying antirheumatic drug (DMARD), indirectly inhibits the enzyme adenosine deaminase, resulting in anti-inflammatory effects. Etanercept, a tumor necrosis factor (TNF) inhibitor, is a large molecule biologic therapy that acts like a sponge to remove TNF- molecules from the joints and blood, thereby reducing disease activity and slowing joint destruction. The two drugs have different mechanisms of action, therefore it is hypothesized that the combination therapy will be more effective in the treatment of rheumatoid arthritis. In the TEMPO study, a total of 682 subjects were randomized at a 1:1:1 ratio of methotrexate alone, etanercept alone or the combo-therapy. At the end of the three year study, approximately 52% of the subjects dropped out from the study due to the following reasons: inefficacious treatment, adverse events, administrative decisions, and others (Table 1). Dropout due to inefficacious treatment was highest in the methotrexate group (18.4%) compared to etanercept (17.5%) and the combination treatment (5.6%). Dropout due to adverse events was lowest in the combination treatment (20%) compared to methotrexate (25.4%) and etanercept (20.2%).
|Completed study||88 (38.6%)||107 (48.0%)||132 (57.1%)||327 (47.9%)|
|Discontinued study||140 (61.4%)||116 (52.0%)||99 (42.9%)||355 (52.1%)|
|adverse event||58 (25.4%)||45 (20.2%)||46 (20%)||149 (21.8%)|
|inefficacious treatment||42 (18.4%)||39 (17.5%)||13 (5.6%)||94 (13.8%)|
|administrative decision||9 (4.0%)||13 (5.8%)||19 (8.2%)||41 (6%)|
|other reasons||31 (13.6%)||19 (8.5%)||21 (9.1%)||71 (10.4%)|
The remainder of this paper is organized as follows. Section 2 presents our proposed Bayesian joint model. In Section 3, we apply the model to the TEMPO study and concluding remarks are provided in Section 4.
2 Bayesian Joint Model
2.1 Model for Longitudinal DAS28 Scores
Let denote the DAS28 scores for subject where is the number of DAS28 scores collected over time for subject in the study. Let
denote the corresponding vector of observation times anddenote the covariate vector for subject .
In the subject level data, we observed a change-point pattern in the DAS28 scores with a steep decrease in DAS28 scores after initial treatment followed by stability or a more gradual change. For this reasons we consider the following random change-point model
where is the trajectory function on the log scale for subject and is the associated random residual error. We modeled on the log scale to improve linearity and to remove the constraint of positive values. The trajectory function on the log scale corresponds to two straight lines that meet at a random change-point . The parameters , , , and represent the subject level intercept, slopes before and after the random change-point respectively. Here we use the notation and.
We use normal prior distributions for the parameters , , , and .
We let the means of the parameters vary for each treatment group by letting
be dummy variables for treatment. Theparameters are then coefficients that will represent treatment effect. In our application, with three treatment groups,
. We complete the model by placing normal hyper-priors with means of zero and large variances on the hyperparameters, this is . We also choose vague priors for the variances:
. We note that the above model is equivalent to a mixed-effect model where the ”random-effect parameters” are re-parameterized and centered around the ”fixed-effect parameters”. This hierarchical centering reparameterization technique reduces autocorrelation between consecutive Markov chain Monte Carlo (MCMC) samples(Gelfand et al., 1995).
2.2 Model for Competing Risk Dropout
Competing risk is defined by Gooley et al. (1999)
as the situation where one type of event precludes the occurrence or alters the probability of occurrence of another event(Pintilie, 2006). We consider the following competing risk model for informative dropout. For the motiviating example, dropout due to adverse event and dropout due to inefficacious treatment is considered to be informative. Dropout due to administrative decisions and other are considered uninformative, that is, dropout is random and does not depend on disease activity outcome.
Let denote the time-to-dropout for subject due to risk where for adverse event and for inefficacious treatment. Let denote the corresponding censoring indicators for subject where if subject dropped out of the study due to risk , 0 otherwise. Because dropout tends to occur in the earlier part of the study, we use a lognormal survival regression model for to allow for a humped-shaped hazard function,
Here are the cause-specific effects on dropout, are the dummy variables for treatment. The last term, , is the component for the joint analysis where is the latent variable (or vector) that, in order to link the longitudinal and dropout processes, is a function of the parameters in the RHS of (1), and is the coefficient (or vector of coefficients). See next section for more detail about . We use normal hyper-priors for and with means of zero and large variances, say 1000. We also assume a vague prior for .
3 Application to the TEMPO Study
We coded the joint model described in Section 2 in WinBUGS (Lunn et al., 2000, see Appendix for code). The covariates in the model are and representing the two dummy variables required for three treatment groups. For example, if subject received methotrexate alone, 0 otherwise, and if subject received etanercept alone, 0 otherwise. Here, the combo-therapy is the reference group.
For , we consider different latent variables described in the literature. In the works of Faucett and Thomas (1996); Wang and Taylor (2001); Ibrahim et al. (2004), and R Brown and G Ibrahim (2003), missingness (or dropout) is considered to be associated with outcome. This is known as an outcome dependent selection model where the estimated trajectory at the time of dropout is used as the latent variable. Wu and Carroll (1988); Wulfsohn and Tsiatis (1997); Follmann and Wu (1995); Henderson et al. (2000); Guo and Carlin (2004) considered an alternative approach known as the shared parameter selection model where the missingness (or dropout) is directly related to a trend over time. In the shared parameter selection models, the parameters from the trajectory function are used as the latent variables. We considered the different latent variable(s) in the literature and used a selection model criterion, Deviance Information Criterion (DIC), to determine which latent variable(s) to use in the final model. The DIC provides a measure of goodness of fit penalized by the effective number of parameters. The model with the smallest DIC is preferred (Spiegelhalter et al., 2002).
For each model we ran 10,000 MCMC iterations with two chains, each with different initial values, retaining every fifth sample. We discarded the first 5,000 iterations as burn-in. We assessed convergence using standard diagnostic tools, that is, converging and mixing of the two chains. We also assessed autocorrelation as a function of time-lag. Model 1 sets , representing separate analysis where the longitudinal model and the dropout process are not linked. Model 1 assumes that dropout is not informative of outcome. Model 2 considers dropout is associated with the outcome at the time of dropout, denoted as , . Model 3 considers dropout is associated with the disease activity at baseline, . In Model 4, dropout is associated with the initial change in disease activity prior to the change-point, . In Model 5, dropout is associated with the change in disease activity following the change-point, . Model 6 shares all parameters from the trajectory function where is now a vector with . Convergence and acceptable autocorrelation was met by all six models. Model 6 is preferred with the smallest DIC (Table 2).
|Model: Type of Analysis||DIC|
|1: Separate analysis||0||781.43|
|Joint analysis where dropout is related to:|
|2: disease activity at time of dropout||488.98|
|3: baseline disease activity||790.63|
|4: initial change in disease activity after treatment||574.10|
|5: the long-term change in disease activity||406.45|
|6: sharing all parameters from the trajectory||391.04|
Visually, Model 6 fits well to the subject-level data (Figure 2). Figure 3 shows the population-level curves for each treatment group using separate analysis (Model 1) versus joint analysis (Model 6). One would conclude that subjects exposed to the combo-therapy are more likely to reach clinical remission earlier (DAS28 score ). However, Model 1 tends to overestimate treatment effectiveness in the methotrexate group (MTX) where dropout due to inefficacy was the highest 18.4% compared to 17.5% in the etanercept group (ETAN) and 5.6% in the combo-therapy (MTXETAN).
|(a) Patient completed study||(b) Patient dropped out due to AE|
|(c) Patient dropped out due to||(d) Patient dropped out due to|
|inefficacious treatment||administrative reasons|
Clinical trials are prone to incomplete outcome data due to patient drop. In the presence of informative dropout and when drop out is differential across the treatment groups, separate analysis of the longitudinal data may lead to biased results. In the example TEMPO study, separate analysis tends to overestimate the treatment effect in the methotrexate group (MTX) where dropout due to inefficacy was the highest 18.4% compared to 17.5% in the etanercept group (ETAN) and 5.6% in the combo-therapy (MTXETAN). The proposed joint model provides a way to perform longitudinal data analysis while adjusting for informative dropout with competing risk. Such approach of joint modeling can make the effect estimates less biased and therefore we can draw more reliable conclusions from the clinical trial. In the case where the proportion of dropout is small and there is no association between dropout and patient outcome, the joint model analysis should give similar results to the separate model analysis.
In this paper, we also considered a linear regression random-change point model for the longitudinal DAS28 scores to account for the observed change-point pattern in the DAS28 scores with a steep decrease in DAS28 scores after initial treatment followed by stability or a more gradual change. The complexity of the model is easily handled under a Bayesian framework. The time to run the model is within minutes and convergence was observed for all parameters. Our proposed Bayesian joint model can also be extended to safety endpoints where safety endpoints are modeled as a survival process.
- Brown and Ibrahim (2003) Brown, E. R. and J. G. Ibrahim (2003). Bayesian approaches to joint cure-rate and longitudinal models with applications to cancer vaccine trials. Biometrics 59(3), 686–693.
- Brown et al. (2005) Brown, E. R., J. G. Ibrahim, and V. DeGruttola (2005). A flexible b-spline model for multiple longitudinal biomarkers and survival. Biometrics 61(1), 64–73.
et al. (2002)
Chen, M.-H., J. G. Ibrahim, and D. Sinha (2002).
Bayesian inference for multivariate survival data with a cure
Journal of Multivariate Analysis80(1), 101–126.
- Chen et al. (2004) Chen, M.-H., J. G. Ibrahim, and D. Sinha (2004). A new joint model for longitudinal and survival data with a cure fraction. Journal of Multivariate Analysis 91(1), 18–34.
- Chi and Ibrahim (2007) Chi, Y. and J. Ibrahim (2007). A new class of joint models for longitudinal and survival data accomodating zero and non-zero cure fractions: A case study of an international breast cancer study group trial. Stat Sin 17, 445–462.
- Chi and Ibrahim (2006) Chi, Y.-Y. and J. G. Ibrahim (2006). Joint models for multivariate longitudinal and multivariate survival data. Biometrics 62(2), 432–445.
- De Gruttola and Tu (1994) De Gruttola, V. and X. M. Tu (1994). Modelling progression of cd4-lymphocyte count and its relationship to survival time. Biometrics, 1003–1014.
- Elashoff et al. (2007) Elashoff, R. M., G. Li, and N. Li (2007). An approach to joint analysis of longitudinal measurements and competing risks failure time data. Statistics in medicine 26(14), 2813–2835.
- Elashoff et al. (2008) Elashoff, R. M., G. Li, and N. Li (2008). A joint model for longitudinal measurements and survival data in the presence of multiple failure types. Biometrics 64(3), 762–771.
- Faucett et al. (1998) Faucett, C. L., N. Schenker, and R. M. Elashoff (1998). Analysis of censored survival data with intermittently observed time-dependent binary covariates. Journal of the American Statistical Association 93(442), 427–437.
et al. (2002)
Faucett, C. L., N. Schenker, and J. M. Taylor (2002).
Survival analysis using auxiliary variables via multiple imputation, with application to aids clinical trial data.Biometrics 58(1), 37–47.
- Faucett and Thomas (1996) Faucett, C. L. and D. C. Thomas (1996). Simultaneously modelling censored survival data and repeatedly measured covariates: a gibbs sampling approach. Statistics in medicine 15(15), 1663–1685.
- Follmann and Wu (1995) Follmann, D. and M. Wu (1995). An approximate generalized linear model with random effects for informative missing data. Biometrics, 151–168.
et al. (1995)
Gelfand, A. E., S. K. Sahu, and B. P. Carlin (1995).
Efficient parametrisations for normal linear mixed models.Biometrika, 479–488.
- Glynn et al. (1986) Glynn, R. J., N. M. Laird, and D. B. Rubin (1986). Selection modeling versus mixture modeling with nonignorable nonresponse. In Drawing inferences from self-selected samples, pp. 115–142. Springer.
- Gooley et al. (1999) Gooley, T. A., W. Leisenring, J. Crowley, and B. E. Storer (1999). Estimation of failure probabilities in the presence of competing risks: new representations of old estimators. Statistics in medicine 18(6), 695–706.
- Guo and Carlin (2004) Guo, X. and B. P. Carlin (2004). Separate and joint modeling of longitudinal and event time data using standard computer packages. The american statistician 58(1), 16–24.
- Hanson et al. (2011) Hanson, T. E., A. J. Branscum, and W. O. Johnson (2011). Predictive comparison of joint longitudinal-survival modeling: a case study illustrating competing approaches. Lifetime data analysis 17(1), 3–28.
- Henderson et al. (2000) Henderson, R., P. Diggle, and A. Dobson (2000). Joint modelling of longitudinal measurements and event time data. Biostatistics 1(4), 465–480.
- Hogan and Laird (1997) Hogan, J. W. and N. M. Laird (1997). Mixture models for the joint distribution of repeated measures and event times. Statistics in medicine 16(3), 239–257.
- Huang et al. (2010) Huang, X., G. Li, and R. M. Elashoff (2010). A joint model of longitudinal and competing risks survival data with heterogeneous random effects and outlying longitudinal measurements. Statistics and its interface 3(2), 185.
- Ibrahim et al. (2004) Ibrahim, J. G., M.-H. Chen, and D. Sinha (2004). Bayesian methods for joint modeling of longitudinal and survival data with applications to cancer vaccine trials. Statistica Sinica, 863–883.
- Ibrahim et al. (2005) Ibrahim, J. G., M.-H. Chen, and D. Sinha (2005). Bayesian survival analysis. Wiley Online Library.
- Jacqmin-Gadda et al. (2006) Jacqmin-Gadda, H., D. Commenges, and J.-F. Dartigues (2006). Random changepoint model for joint modeling of cognitive decline and dementia. Biometrics 62(1), 254–260.
- Keystone et al. (2009) Keystone, E., B. Freundlich, M. Schiff, J. Li, and M. Hooper (2009). Patients with moderate rheumatoid arthritis (ra) achieve better disease activity states with etanercept treatment than patients with severe ra. The Journal of rheumatology 36(3), 522–531.
- Laird and Ware (1982) Laird, N. M. and J. H. Ware (1982). Random-effects models for longitudinal data. Biometrics, 963–974.
- Law et al. (2002) Law, N. J., J. M. Taylor, and H. Sandler (2002). The joint modeling of a longitudinal disease progression marker and the failure time process in the presence of cure. Biostatistics 3(4), 547–563.
- Little (1993) Little, R. J. (1993). Pattern-mixture models for multivariate incomplete data. Journal of the American Statistical Association 88(421), 125–134.
- Little (1995) Little, R. J. (1995). Modeling the drop-out mechanism in repeated-measures studies. Journal of the American Statistical Association 90(431), 1112–1121.
- Little and Rubin (1987) Little, R. J. and D. B. Rubin (1987). Statistical analysis withmissing data.
- Lunn et al. (2000) Lunn, D. J., A. Thomas, N. Best, and D. Spiegelhalter (2000). Winbugs-a bayesian modelling framework: concepts, structure, and extensibility. Statistics and computing 10(4), 325–337.
- Pauler and Finkelstein (2002) Pauler, D. K. and D. M. Finkelstein (2002). Predicting time to prostate cancer recurrence based on joint models for non-linear longitudinal biomarkers and event time outcomes. Statistics in medicine 21(24), 3897–3911.
- Pawitan and Self (1993) Pawitan, Y. and S. Self (1993). Modeling disease marker processes in aids. Journal of the American Statistical Association 88(423), 719–726.
- Pintilie (2006) Pintilie, M. (2006). Competing risks: a practical perspective, Volume 58. John Wiley & Sons.
- R Brown and G Ibrahim (2003) R Brown, E. and J. G Ibrahim (2003). A bayesian semiparametric joint hierarchical model for longitudinal and survival data. Biometrics 59(2), 221–228.
- Rubin (1976) Rubin, D. B. (1976). Inference and missing data. Biometrika, 581–592.
- Schluchter (1992) Schluchter, M. D. (1992). Methods for the analysis of informatively censored longitudinal data. Statistics in medicine 11(14-15), 1861–1870.
- Song et al. (2002) Song, X., M. Davidian, and A. A. Tsiatis (2002). A semiparametric likelihood approach to joint modeling of longitudinal and time-to-event data. Biometrics 58(4), 742–753.
- Spiegelhalter et al. (2002) Spiegelhalter, D. J., N. G. Best, B. P. Carlin, and A. Van Der Linde (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 64(4), 583–639.
- Tsiatis et al. (1995) Tsiatis, A., V. Degruttola, and M. Wulfsohn (1995). Modeling the relationship of survival to longitudinal data measured with error. applications to survival and cd4 counts in patients with aids. Journal of the American Statistical Association 90(429), 27–37.
- Tsiatis and Davidian (2004) Tsiatis, A. A. and M. Davidian (2004). Joint modeling of longitudinal and time-to-event data: an overview. Statistica Sinica, 809–834.
- Wang and Taylor (2001) Wang, Y. and J. M. G. Taylor (2001). Jointly modeling longitudinal and event time data with application to acquired immunodeficiency syndrome. Journal of the American Statistical Association 96(455), 895–905.
- Wu and Carroll (1988) Wu, M. C. and R. J. Carroll (1988). Estimation and comparison of changes in the presence of informative right censoring by modeling the censoring process. Biometrics, 175–188.
- Wulfsohn and Tsiatis (1997) Wulfsohn, M. S. and A. A. Tsiatis (1997). A joint model for survival and longitudinal data measured with error. Biometrics, 330–339.
- Xu and Zeger (2001) Xu, J. and S. L. Zeger (2001). The evaluation of multiple surrogate endpoints. Biometrics 57(1), 81–87.
- Yu et al. (2004) Yu, M., N. J. Law, J. M. Taylor, and H. M. Sandler (2004). Joint longitudinal-survival-cure models and their application to prostate cancer. Statistica Sinica, 835–862.
- Zhu et al. (2012) Zhu, H., J. G. Ibrahim, Y.-Y. Chi, and N. Tang (2012). Bayesian influence measures for joint models for longitudinal and survival data. Biometrics 68(3), 954–964.