A Bayesian Joint model for Longitudinal DAS28 Scores and Competing Risk Informative Drop Out in a Rheumatoid Arthritis Clinical Trial
Violeta G. Hennessey Luis G. Leon-Novelo Juan Li Li Zhu Eric Chi Joseph G. Ibrahim
stat.AP
















A Bayesian Joint model for Longitudinal DAS28 Scores and Competing Risk Informative Drop Out in a Rheumatoid Arthritis Clinical Trial
    Violeta Hennessey
Warner Home Video,  Los Angeles, CA, USA. E-mail: v.g.hennessey@gmail.comLuis León NoveloUniversity of Texas Health Science Center at Houston- School of Public Health, Houston, TX, USA
Juan Li Department of Biostatistics, Amgen Inc., Thousand Oaks, CA, USA
Li Zhu [3]Xin Huang [3]
Eric ChiDepartment of Biostatistics, Amgen Inc.,  South San Francisco, CA, USA
 Joseph G. Ibrahim Department of Biostatistics, The University of Noth Carolina at Chapel Hill, Chapel Hill, NC, USA
    
==================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================


Rheumatoid arthritis 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 patients activity score.   In the presence of informative dropout, that is, the dropout depends on latent variables from the longitudinal process, simply applying a model to analyze the longitudinal outcomes may lead to biased results because the assumption of random dropout is violated. In this paper we develop a data driven Bayesian joint model for modeling DAS28 scores and competing risk informative drop out. The motivating example is a clinical trial of Etanercept and Methotrexate with radiographic Patient Outcomes  <cit.>.


Keywords:
joint selection models;
longitudinal data;
patient reported outcomes; 
random change-point.



  (-27mm,-1.5cm)[0pt][r]Hennessey et.al.







§ INTRODUCTION

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= 0.56√(Tender28)+0.28√(Swollen28)+0.70log(ESR);                                                        +0.014(Patient Global Health on VAS) ]

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.



r0.4
 

 
    < g r a p h i c s >

 

  
 []Twenty-eight 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 
<cit.>, 
 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
 <cit.>.  
 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 
<cit.>. 
 Joint models have been widely used in HIV/AIDS clinical trials, cancer vaccine (immunotherepy) trials and quality of life studies. References include 
<cit.>, 
 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.
<cit.> 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.
<cit.> 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. 
<cit.> 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 <cit.>. 
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 
<cit.> 
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 
 <cit.> 
 and Bayesian approaches 
 <cit.>. 
 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  <cit.>. 
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 <ref>). 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%).


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.


§ BAYESIAN JOINT MODEL




 §.§ Model for Longitudinal DAS28 Scores

Let 𝐲_i=(y_i1,…,y_in_i) denote the DAS28 scores for subject i where n_i is the number of DAS28 scores collected over time for subject i in the study. Let 𝐭_i=(t_i1,…,t_in_i) denote the corresponding vector of observation times and 𝐱_i=(x_i1,…,x_im) denote the covariate vector for subject i.

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

    log(y_ij)=ψ_ij+ϵ_ij



    ψ_ij=α_i+ β_1iκ_i+β_1i(t_ij-κ_i)_- +β_2i(t_ij-κ_i)_+

where ψ_ij=ψ_i(t_ij) is the trajectory function on the log scale for subject i and ϵ_ij is the associated random residual error. We modeled y_ij 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 κ_i.
The parameters α_i, β_1i, β_2i, and κ_i represent the subject level intercept, slopes before and after the random change-point respectively. Here we use the notation (t_ij-κ_i)_- = min{t_ij-κ_i, 0} and
(t_ij-κ_i)_+ = max{t_ij-κ_i, 0}. The random residual errors are assumed to be 
independent and
normally distributed with mean zero and variance σ^2.

We use normal prior distributions for the parameters α_i, β_1i, β_2i, and κ_i.


    [                        α_i∼            N(μ_α_i, σ^2_α),                      μ_α_i=   γ_α0+γ_α1x_i1+…+γ_α mx_im;                       β_1i∼         N(μ_β_1i, σ^2_β_1),                     μ_β_1i= γ_β10+γ_β11x_i1+…+γ_β1mx_im;                       β_2i∼         N(μ_β_2i, σ^2_β_2),                     μ_β_2i= γ_β20+γ_β21x_i1+…+γ_β2mx_im;                        κ_i∼            N(μ_κ_i, σ^2_κ),                      μ_κ_i=   γ_κ0+γ_κ1x_i1+…+γ_κ mx_im;                             ]

We let the means of the parameters vary for each treatment group by letting x_i1, …, x_im be dummy variables for treatment. The γ parameters are then coefficients that will represent treatment effect. 
In our application, with three treatment groups, m=2.
We complete the model by placing normal hyper-priors with means of zero and large variances on the hyperparameters γ, this is N(0,10^3). 
 We also choose vague priors for the variances:
σ^2,
σ^2_α,σ^2_β_1,σ^2_β_2,
σ^2_κInverse gamma(0.01,0.01).
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 
<cit.>.



 §.§ Model for Competing Risk Dropout

Competing risk is defined by 
<cit.> 
as the situation where one type of event precludes the occurrence or alters the probability of occurrence of another event 
<cit.>. 
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 d_ik denote the time-to-dropout for subject i due to risk k where k=1 for adverse event and k=2 for inefficacious treatment. Let c_ik denote the corresponding censoring indicators for subject i where c_ik=1 if subject i dropped out of the study due to risk k, 0 otherwise. Because dropout tends to occur in the earlier part of the study, we use a lognormal survival regression model for d_ik to allow for a humped-shaped hazard function,

    [                log(d_ik)∼N(θ_ik,ς^2_k); θ_ik=φ_k0+φ_k1x_i1+…+φ_kmx_im + ω_kν_i ]

Here φ_k0,…,φ_km are the cause-specific effects on dropout, x_i1, …, x_im are the dummy variables for treatment. The last term, ω_kν_i, is the component for the joint analysis where ν_i 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 
(<ref>), and
ω_k is the coefficient (or vector of coefficients). See next section for more
detail about ν_i.
  We use normal hyper-priors for φ_kj and ω_k with means of zero and large variances, say 1000. 
We also assume a vague prior for ς^2_k∼Inverse gamma(.01,.01).




§ APPLICATION TO THE TEMPO STUDY

We coded the joint model described in Section 2 in WinBUGS 
<cit.>. 
 The covariates in the model are x_i1 and x_i2 representing the two dummy variables required for three treatment groups. For example, x_i1=1 if subject i received methotrexate alone, 0 otherwise, and x_i2=1 if subject i received etanercept alone, 0 otherwise. Here, the combo-therapy is the reference group.

For ν_i, we consider different latent variables described in the literature. In the works of 
<cit.>, 

and 
<cit.>,  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. 
<cit.> 
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 <cit.>. 

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 ν_i=0, 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 ψ_i^*, ν_i=ψ_i^*. Model 3 considers dropout is associated with the disease activity at baseline, ν_i=α_i. In Model 4, dropout is associated with the initial change in disease activity prior to the change-point, ν_i=β_1i. In Model 5, dropout is associated with the change in disease activity following the change-point, ν_i=β_2i. Model 6 shares all parameters from the trajectory function where ν_i is now a vector with ν_𝐢=(α_i, β_1i, β_2i). Convergence and acceptable autocorrelation was met by all six models. Model 6 is preferred with the smallest DIC (Table <ref>).





Visually, Model 6 fits well to the subject-level data (Figure <ref>). Figure <ref> 
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 < 2.6). 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).




 










 
 
  
 
 










§ DISCUSSION

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.



















 


chicago
























 §.§ Appendix

WinBUGS code for Model 6: Joint model that shares parameters from the longitudinal trajectory function




