PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of biostsLink to Publisher's site
 
Biostatistics. Oct 2011; 12(4): 737–749.
Published online Jan 24, 2011. doi:  10.1093/biostatistics/kxq084
PMCID: PMC3202304
A shared parameter model for the estimation of longitudinal concomitant intervention effects
Colin O. Wu,* Xin Tian, and Wenhua Jiang
Office of Biostatistics Research, National Heart, Lung and Blood Institute, Bethesda, MD 20892, USA ; wuc/at/nhlbi.nih.gov
*To whom correspondence should be addressed.
Received February 27, 2010; Revised December 17, 2010; Accepted December 18, 2010.
We investigate a change-point approach for modeling and estimating the regression effects caused by a concomitant intervention in a longitudinal study. Since a concomitant intervention is often introduced when a patient's health status exhibits undesirable trends, statistical models without properly incorporating the intervention and its starting time may lead to biased estimates of the intervention effects. We propose a shared parameter change-point model to evaluate the pre- and postintervention time trends of the response and develop a likelihood-based method for estimating the intervention effects and other parameters. Application and statistical properties of our method are demonstrated through a longitudinal clinical trial in depression and heart disease and a simulation study.
Keywords: Change-point model, Concomitant intervention, Likelihood, Longitudinal study, Shared parameter model
Change-point models are often used for evaluating the intervention effects in longitudinal clinical trials or epidemiological studies. By treating the starting time of an intervention as a change-point time, the longitudinal trends of the response variable in these models are described using different regression models before and after the intervention. Estimation and inferences can be obtained using well-known procedures, such as the maximum likelihood estimation (MLE), restricted maximum likelihood estimation (REMLE), and generalized estimation equations (e.g. Verbeke and Molenberghs, 2000; Diggle and others, 2002). This approach works well when the intervention is randomly assigned. When an intervention is not randomly assigned, treating the intervention as a time-dependent covariate in a change-point model may lead to model misspecification and biased estimation of the intervention effects (Wall and others, 2005).
Concomitant interventions are generally not randomly assigned in biomedical studies. The Enhancing Recovery in Coronary Heart Disease (ENRICHD) Study is a typical example which involves a concomitant intervention in addition to the randomly assigned treatment regimens. This is a randomized clinical trial for evaluating the efficacy of a 6-month psychosocial treatment (PT) versus the usual cardiovascular care (UC) on survival, hospitalization, myocardial infarction, and depression severity in 2481 patients with depression and/or low perceived social support after acute myocardial infarction (ENRICHD, 2001, ENRICHD, 2003). Depression severity was measured by Hamilton Rating Scale for Depression (HRSD) at the baseline and 6-month visits and Beck Depression Inventory (BDI) at regularly scheduled visits during the treatment and follow-up periods, where higher HRSD and BDI scores indicate worsened depression. BDI scores for patients in the PT arm were repeatedly measured at weekly visits during the treatment and 4 yearly follow-up visits, while BDI scores for patients in the UC arm were only measured at baseline, the 6-month visit and the yearly follow-up visits. Patients in both treatment arms were eligible for pharmacotherapy with antidepressants as a concomitant intervention if they had high baseline BDI scores or nondecreasing BDI trends after enrollment. In addition, antidepressants were prescribed at the requests of the patients or their primary care physicians. Taylor and others (2005) investigated the effects of antidepressants on cardiovascular morbidity and mortality among 1834 depressed ENRICHD patients and found that pharmacotherapy improved survival for this patient population. Their results, however, did not address the question of whether pharmacotherapy was beneficial for lowering patients' depression severity. In an analysis using a subsample of 91 ENRICHD patients in the PT arm who have received pharmacotherapy within the treatment period, Wu and others (2008) showed that the naïve mixed-effects models were misspecified when a concomitant intervention was present and proposed a varying-coefficient mixed-effects model to reduce the potential bias associated with the estimated pharmacotherapy effects. A main drawback of Wu and others (2008) is the potential loss of information because their regression model cannot be applied to patients who have already received a concomitant intervention at baseline or have not received any concomitant intervention during the study.
We propose in this paper a comprehensive regression method for evaluating the concomitant intervention effects which is capable to incorporate information from all the patients in a longitudinal study. Based on the framework of shared parameter models (Follmann and Wu, 1995), we describe the covariate effects on the response variable through a change-point mixed-effects model and incorporate the random coefficients and the intervention starting time (change-point time) through a series of joint distributions. Patients who have received a concomitant intervention at baseline or have not received any concomitant intervention during the study are treated as censored. A likelihood-based method is established for statistical estimation and inferences and its computation is implemented through a 2-stage iteration procedure. Our application to the ENRICHD pharmacotherapy data and simulation results suggest that the proposed method leads to adequate estimates when a concomitant intervention is present, while the naïve mixed-effects model is likely misspecified under such situations.
We describe our shared parameter model and its interpretations in Section 2, develop the likelihood-based estimation and inference procedures in Section 3, and present the application to the ENRICHD pharmacotherapy data and the simulation results in Sections 4 and 5, respectively. In Section 6, we discuss some potential extensions of our method. Further theoretical and technical results can be found in the supplementary material available at Biostatistics online.
2.1. Data structure
Let T0 and T1 be the beginning and ending times of a study, and n the total number of randomly selected subjects. The ith subject has ni visits and observations (Tij,Yij,Xi) at the jth visit, where Tij, the study time, is the time elapsed from the beginning of the study to the jth visit, Xi is a time-invariant covariate vector, and Yij is the real-valued outcome variable. For simplicity, we assume in this paper that the study involves only one concomitant intervention, and each subject can change from “without intervention” to “intervention” only once during the study. Extension to multiple concomitant interventions and multiple intervention changes involves more notation and computational complexity, hence, is not discussed in the current manuscript. When the context is clear, we may simply refer a “concomitant intervention” as an intervention. Define Si to be the intervention starting time or change-point time, and λij = 0 or 1, if TijSi or Tij > Si, respectively, to be the intervention indicator for the ith subject. Since not every subject changes from without intervention to intervention during the study, the ith subject's change-point time is observed if Ti1SiTini. If Si < Ti1 or Si > Tini, the subject's change-point time is left or right censored, respectively. The indicator variable for censoring δi(c) is defined by δi(c) = 0 if Ti1SiTini, 1 if Si > Tini, and 2 if Si < Ti1. The observed change-point times are {Si(c) = (Si(c),δi(c));i = 1,…,n}, where Si(c) = Si if δi(c) = 0, Tini if δi(c) = 1, and Ti1 if δi(c) = 2.
2.2. Model formulation and rationale
The change-point model assumes that Yij follows different trajectories before and after an intervention. The ith subject's trajectory before the intervention is μ0(Tij,Xi;ai), which is parameterized by the subject-specific parameter aiT = (ai1,…,aid0)T, d0 ≥ 1. The change of the trajectory after the intervention is μ1(Tij,Xi,Rij;bi), which is parameterized by the subject-specific parameter bi = (bi1,…,bid1)T, d1 ≥ 1, and may depend on the “intervention duration time” Rij = TijSi as well as Tij and Xi. Following the usual mixed-effects model framework (Davidian and Giltinan, 1995; Verbeke and Molenberghs, 2000; Diggle and others, 2002), a naïve mixed-effects model for evaluating the pre- and postintervention trajectories is
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx1_ht.jpg
(2.1)
where (aiT,biT)T has unknown mean (αT,βT)T and covariance matrix Σ, ϵij are mean-zero random errors with cov(ϵij1,ϵij2) = σij1j2, and ϵi1j1 and ϵi2j2 are independent if i1i2. The joint distribution of ai and bi is usually assumed to be multivariate Gaussian [mathematical script N]{(αT,βT)T,Σ}. A positive (or negative) value for μ1(Tij,Xi,Rij;bi) would suggest that the intervention tends to increase (or decrease) the mean of Yij given (Tij,Xi,Rij). Wu and others (2008) showed that (2.1) could be a misspecified model even if μ0(·;ai), μ1(·;bi), ϵij, and G(·) were correctly specified because the “self-selectiveness” of the intervention was ignored.
If the change-point time Si depends on the preintervention trend of Yij, a natural extension of (2.1) is
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx2_ht.jpg
(2.2)
where Fa(·) is the cumulative distribution function (CDF) of ai, Fs(·|ai) and Fb(·|ai) are the conditional CDFs of Si and bi, respectively, given ai, and bi and Si are independent given ai. In contrast to the varying-coefficient model of Wu and others (2008), where the conditional means of ai and bi given Si are used, (2.2) incorporates Si through Fs(·|ai), which allows Si to be left or right censored. We assume for simplicity in (2.2) that bi does not depend on Si, although further generalizations may allow the distribution of bi to depend on (Si,ai). In practice, further specifications of Fa(·), Fs(·|ai), and Fb(·|ai) are usually needed to balance the computational feasibility and flexibility of the model.
2.3. Linear and additive shared parameter models
A useful and mathematically tractable specification for (2.2) is to assume that Si and bi depend on ai only through their conditional means. If the conditional means of Si and bi are linear models of ai, a linear shared parameter model for (2.2) is
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx3_ht.jpg
(2.3)
where α = (α1,…,αd0)T, αd[set membership]R, β = (β1T,…,βd1T)T, βl = (βl0,…,βld0)T, βld[set membership]R, γ = (γ0,…,γd0)T, and {ϵi = (ϵi1,…,ϵini)T,ei(a),ei(b),ei(s)} are independent mean-zero random errors with covariance matrices {Vy,Va,Vb,σs2}, respectively. Here, the unknown parameters are the mean components θ = (αT,β1T,…,βd1T,γT)T and the covariance structures V = {Vy,Va,Vb,σs2}.
When the relationship between Si and ai are unknown, a nonparametric model is Si = μ(s)(ai) + ϵi(s), where μ(s)(ai) = E(Si|ai) is a smooth function of ai. Since unstructured estimation of μ(s)(ai) could be difficult when ai is a high-dimensional vector, a simple additive model is to replace the relationship between Si and ai in (2.3) with
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx4_ht.jpg
(2.4)
where μd(s)(aid) are smooth functions of aid. Further generalizations of (2.4) are theoretically possible but at the expense of computational complexity.
3.1. General formulation and approximate log-likelihood for splines
The parameters in (2.2) can be estimated by a likelihood approach if the distributions can be specified. Let Yi = (Yi1,…,Yini)T, Ti = (Ti1,…,Tini)T, Di = (Ti,Xi), and fy(·), fb(·), fs(·), and fa(·) be the densities of Yij, bi, Si, and ai. The joint density of (bi,Si,ai) is f(bi,Si,ai) = fb(bi|ai)fs(Si|ai)fa(ai). Integrating over ai and bi, the conditional density of (Yi,Si) given Di = (Ti,Xi) is
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx5_ht.jpg
The conditional density fs(Si(c),δi(c)|ai) of Si(c) = (Si(c),δi(c)) given ai is fs(Si|ai) if δi(c) = 0, 1 − Fs(Tini|ai) if δi(c) = 1, and Fs(Ti1|ai) if δi(c) = 2. The log-likelihood function for (Yi,Si(c)) conditioning on Di, i = 1,…,n, is
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx6_ht.jpg
(3.1)
where f(y,s)(·|·) is given above, f(y,1)(·|Di,ai) and f(y,2)(·|Di,ai,bi) are the densities of Yi given {Di,ai,δi(c) = 1} and {Di,ai,bi,δi(c) = 2}, respectively, f(y,1)(Yi|Di) = ∫f(y,1)(Yi|Di,ai){1 − Fs(Tini|ai)}fa(ai)dai and
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx7_ht.jpg
When f(y,s)(·|·) belongs to a parametric family {f(y,s)(·;[var phi]|Di);[var phi] = (θ,V)}, we can estimate [var phi] = (θ,V) by maximizing the log-likelihood (3.1).
Estimation for the additive model (2.4) can be achieved by an approximate likelihood method. Under some mild smoothness conditions on μd(s)(·) (Huang and others, 2004), we can approximate μd(s)(·) by the following B-spline expansion:
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx8_ht.jpg
(3.2)
where {Bp(d)(·);1 ≤ pPd} is a spline basis, B(d)(aid) = (B1(d)(aid),…,BPd(d)(aid))T, and γ(d) = (γ1(d),…,γPd(d))T is a set of real-valued coefficients. It follows from (3.2) that Si≈∑d = 0d0(γ(d))TB(d)(aid) + ϵi(s). If we substitute ∑d = 0d0μd(s)(aid) of (2.4) with ∑d = 0d0(γ(d))TB(d)(aid), the parameters are θ = (αT,β1T,…,βd1T,γT)T, γ = ((γ(0))T,…,(γ(d0))T)T, and V = {Vy,Va,Vb,σs2}. Let fs*(·;γ,σs|ai) be the conditional density of ∑d = 0d0{(γ(d))TB(d)(aid)} + ϵi(s) given ai. If the distributions of {ϵij,ei(a),ei(b),ei(s)} are parameterized by V, we can approximate fs(Si|ai) by fs*(Si;γ,σs|ai) and obtain the following approximate log-likelihood function for (Yi,Si(c)) given Di,
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx9_ht.jpg
(3.3)
where [var phi] = (θ,V), f(y,s)*(·|Di), f(y,k)*(·|Di), k = 1,2, are given in (3.1) with fs(Si|ai) replaced by fs*(Si;γ,σs|ai). If Ld*([var phi]) satisfies the regularity conditions for MLEs, the approximate MLE An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx15_ht.jpg satisfy An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx15_ht.jpg.
3.2. A 2-stage estimation procedure
Since the computation of maximum likelihood (ML) or approximate ML estimators requires a nonlinear maximization over the parameter spaces of θ and V, a global maximization over θ and V simultaneously could be computationally unfeasible in practice. We propose a 2-stage procedure which combines REMLE with the Newton–Raphson algorithm:
Stage 1. Assume that {ϵij,ai,bi,Si} of (2.3) or (2.4) are independent random variables with covariance matrices V = {Vy,Va,Vb,σs2}, that is, the naïve mixed-effects model (2.1) holds. Compute An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx10_ht.jpg of V using the REMLE procedure.
Stage 2. Substitute V with An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx10_ht.jpg and maximize An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx11_ht.jpg with respect to θ using the Newton–Raphson procedure. The maximizer An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx12_ht.jpg is the approximate ML estimator for θ.
By (3.1) and the expressions of f(y,s)(Yi,Si|Di) and f(y,k)(Yi|Di) for k = 1,2, the Newton–Raphson algorithm for maximizing An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx11_ht.jpg at Stage 2 involves multidimensional integrations over the functions of ai, bi, and Si with respect to the joint distributions of ai and bi. Monte Carlo simulations are used to compute all the necessary quantities, including the log-likelihoods, and their gradients and Hessian matrices. Implementation of the Newton–Raphson algorithm can be computationally costly, since, in each iteration, large Monte Carlo samples are required to compute the gradient and the Hessian matrix. The computational complexity further increases with the number of subjects. If a suitable initial estimator is available, computation of the algorithm can be significantly reduced by a “1-step” Newton–Raphson procedure (Bickel, 1975; Cai and others, 2000). The estimators computed from the REMLE can be used as a natural candidate for the initial estimator An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx15_ht.jpg. Initial estimators of γ can be computed by fitting the regression model An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx13_ht.jpg using those subjects with Si observed (i.e. δi(c) = 0), where An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx14_ht.jpg is the predicted value for ai. Further details of this algorithm are provided in Appendix C of the supplementary material available at Biostatistics online.
3.3. Inferences
Inferences for the MLE of a parametric model (3.1) may be obtained using their asymptotic distributions. Let An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx15_ht.jpg be the MLE of [var phi]. The asymptotic normality of the MLE's (Serfling, 1980, Chapter 4) implies that, when n is large, An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx15_ht.jpg has approximately the multivariate normal distribution [mathematical script N]([var phi], Var([var phi])). An approximate [100×(1 − α)]th confidence interval (CI) for [ell]([var phi]), a linear combination of [var phi], is An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx16_ht.jpg, where An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx17_ht.jpg is the variance estimator and Zα/2 is the [100×(1 − α/2)]th percentile of the standard normal distribution. For the additive model (2.4), where nonparametric components are present, asymptotic distributions of the approximate MLE in (3.2) have not yet been developed. As a practical alternative, we can generate bootstrap samples by resampling the subjects with replacement and compute the corresponding bootstrap estimators. The estimates obtained from the original sample are natural choices for the initial estimates for the bootstrap samples. A [100×(1 − α)]th CI based on percentiles is given by the corresponding lower and upper [100×(α/2)]th percentiles (Lα/2,Uα/2) of the bootstrap estimators. This bootstrap approach has been used in Wu and others (2008) among others. We investigate in Section 5 the empirical coverage probabilities (CP) of our CIs in a simulation study.
4.1. The ENRICHD pharmacotherapy data
We applied our method to the ENRICHD Study and investigated whether the use of antidepressants as a pharmacotherapy could lower the BDI scores for patients during the PT period. Wu and others (2008) investigated this question by applying their varying-coefficient model to the subsample of 91 depressed patients in the ENRICHD PT arm who had their change-point time Si observed during the treatment period and concluded that the use of antidepressants was beneficial for lowering the BDI scores for patients who were eligible for pharmacotherapy. Our analysis was based on all 557 depressed patients who had their status of antidepressant use (yes or no) and exact dates of antidepressant starting time recorded and attended 5 or more PT sessions during the treatment period. In order to have a meaningful interpretation for this study, 3 types of patients were excluded from our analysis: patients in the usual care arm because of the lack of BDI scores observed between the baseline and the end of the treatment period, patients whose pharmacotherapy status and starting dates of antidepressant use were not recorded, and patients who had poor adherence to the required weekly PT sessions, that is, attended less than 5 (approximately 20%) of the sessions. Because antidepressant use for each patient was individually monitored and recorded as accurate as possibly by study psychiatrists (Taylor and others, 2005, p. 794), it is reasonable to assume that the missing records on antidepressant starting dates were missing at random. Within our sample, 11 patients used antidepressants before baseline, 92 started antidepressant during the treatment period, and 454 did not use antidepressants before and during the treatment period. The number of visits for these patients ranges from 5 to 36 and has a median of 12.
Let Yij, Tij, Si(c), and Rij = TijSi(c) be the ith patient's BDI score, trial time (months), starting time (months) of antidepressant use, and antidepressant duration time (months), respectively, at the jth visit. For all 1 ≤ in, Ti1 = 0, Tini is the trial time at the last visit, and λij = 1[Si < Tij] is the indicator of antidepressants use at Tij. The observed (Si(c),δi(c)) is (Si(c) = Si,δi(c) = 0) if the ith patient used antidepressants within the treatment period, (Si(c) = Tini,δi(c) = 1) if the patient did not use antidepressant within the treatment period, and (Si = 0,δi(c) = 2) if the patient used antidepressants before baseline. When the linear models μ0(Tij;ai0,ai1) = ai0 + ai1Tij and μ1(Rij;bi0,bi1) = bi0 + bi1Rij are used, (ai0,ai1) represents the intercept and slope of the ith subject's BDI trajectory before antidepressant use, and (bi0,bi1) is the intercept and slope of the change of the subject's BDI trajectory after antidepressant use.
4.2. Fitting the naïve mixed-effects model
In order to have a clear and parsimonious interpretation of the overall trend of BDI over the trial time, we examined some preliminary analysis of the data and approximated the relationship between Yij and Tij by a simple linear model. Using the naïve mixed-effects model (2.1), we first estimated the unknown parameters α0, α1, β00, and β10 of
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx18_ht.jpg
(4.1)
where (α0,α1,β00,β10)T is the unknown mean vector and Σ is the unstructured covariance matrix for the multivariate normal distribution [mathematical script N](·,·). Intuitively, α0 and α1 are the mean intercept and slope for BDI scores before antidepressant use, and β00 and β10 are the mean intercept and slope for the “correction term” after antidepressant use. Zero value of β00 + β10Rij indicates ignorable antidepressant effect on the mean BDI scores.
The left panel of Table 1 shows the estimates of α0, α1, β00, and β10, and their standard errors (SEs) and 95% CIs computed using REMLE with unstructured covariance matrix Σ. The negative estimate An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx19_ht.jpg suggests that the mean BDI score for these patients tends to decrease over the trial time since the start of the PT sessions, while the positive estimate An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx20_ht.jpg and its 95% CI seem to suggest that the use of antidepressants increase the patients' mean BDI scores. As discussed in Section 2.2, the self-selectiveness of the antidepressant use as a concomitant intervention is not considered in (4.1), so that the positive estimate of β00 may not reflect the real effect of pharmacotherapy on depression severity.
Table 1.
Table 1.
Parameter estimates, SEs and 95% CIs computed for the ENRICHD pharmacotherapy data based on the naïve linear mixed-effects model (4.1), the linear shared parameter model (4.2) and the submodel of (4.2) with β01, β02, β (more ...)
4.3. Fitting the shared parameter model
To account for the possible link between antidepressant starting time and the BDI trend before the use of antidepressants, we examined a series of preliminary plots among Si, (ai0,ai1) and (bi0,bi1) using patients in the PT arm who had started using antidepressants during the treatment period. These preliminary results (see Appendix D.1 of the supplementary material available at Biostatistics online) suggested that the following linear shared parameter model could give a parsimonious approximation to the variables considered:
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx21_ht.jpg
(4.2)
where ϵi(a) and ϵi(b) = (ϵi0(b),ϵi1(b))T are mean-zero bivariate normal random vectors with unstructured covariance matrices Σ(a) and Σ(b), respectively, ϵi(s) is a mean-zero normal random variable with variance σs2, and ϵi(a), ϵi(b), and ϵi(s) are independent. Similar to (4.1), (α0,α1) represents the mean intercept and slope for the BDI trend over the trial time before antidepressant use. Given (ai0,ai1)T, (β00,β01,β02,β10,β11,β12)T specifies the conditional mean structure for the BDI correction term after antidepressant use.
The middle panel of Table 1 shows the estimates of the coefficients in (4.2) computed using the 2-stage ML procedure with 10 Newton–Raphson iterations, and their SE and 95% percentile CIs computed from the bootstrap procedure. Since the 95% CIs for β01, β02, β10, and β11 in (4.2) include zero, we repeated the estimation procedure for the submodel of (4.2) with β01, β02, β10, and β11 set to zero and summarized the estimates, their SEs, and 95% percentile CIs from this submodel in the right panel of Table 1, which gives An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx22_ht.jpg and An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx23_ht.jpg. Contrary to the interpretation from (4.1), we obtain from this submodel that An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx24_ht.jpg since 0 ≤ Rij ≤ 6 and An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx25_ht.jpg, which suggests that pharmacotherapy has on average a beneficial effect for lowering a patient's depression severity. The 95% CIs for γ1 obtained from (4.2) and its submodel suggests that (4.1) is likely a misspecified model for this data set. These results are consistent with the analysis of Wu and others (2008).
5.1. Data generation
The simulation study has 2 objectives: (i) to demonstrate the biases of the naïve mixed-effects models; and (ii) to investigate the practical properties of our ML and approximate MLE procedures. In order to resemble a moderate sized longitudinal study with data structures similar to the ENRICHD pharmacotherapy data, each of our simulated samples contains n = 200 or 300 independent subjects who have unequal number of repeated measurements. For i = 1,…,n, the ith subject in the sample has a maximum of 30 “scheduled visits” at time points {Ti,1,…,Ti,30} = {0,0.2 + e1,…,5.8 + e29}, where {el} are independent “distortions” generated from the uniform U( − 0.1,0.1) distribution, but each scheduled visit has 40% probability to be skipped. The remaining observed visiting times are {Tij;j = 1,…,ni} with possibly unequal ni for different subjects. For each 1 ≤ in, the random parameter ai = (ai0,ai1)T is generated by adding (α0,α1)T = (20,10)T to a randomly generated ϵi(a) from the mean-zero bivariate normal distribution with covariance matrix diag(2.5,1), and bi = (bi0,bi1)T is generated by adding the conditional mean (β00 + β01ai0,β12ai1)T with (β00,β01,β12)T = ( − 4,0.01,0.02)T to a randomly generated ϵi(b) from the standard bivariate normal distribution. For each (aiT,biT)T, we generate Si from 2 different scenarios: (a) Si~N(γ1ai0 + γ2ai1,1) with γ = (γ1,γ2)T = (1, − 1.7)T; (b) Si~N(μ0(s)(ai0),0.04), where μ0(s)(ai0) = − 0.75(ai0 − 18)2 + 6 if ai0 ≤ 20, and 0.75(ai0 − 22)2 if ai0 > 20. The purpose of (b) is to evaluate the B-spline estimate when the parametric form of μ0(s)(ai0) is unknown. For each {Tij,Si,ai,bi}, we generate Yij from the subject-specific change-point model
An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx30_ht.jpg
(5.1)
where λij = 1[Tij > Si], Rij = max(TijSi,0), and ϵij~N(0,1). When Si < Ti1 or Si > Tini, we use censored Si and obtain {(Yij,Tij,λij,Si(c),δi(c));j = 1,…,ni}.
The simulation was separately repeated 200 times under scenarios (a) and (b) with n = 200 and 300. Within each simulation, we used μ0(Tij;ai) = ai0 + ai1Tij and μ1(Tij,Rij;bi) = bi0 + bi1Rij, and estimated the parameters α = (α0,α1)T and β = (β00,β01,β12)T using the naïve mixed-effects model (2.1) and the shared parameter model (2.2). For estimation under the naïve mixed-effects model (2.1), we ignored the relationship between ai and Si and computed the restricted ML estimators An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx26_ht.jpg.
5.2. Estimation for linear shared parameter models
For the Gaussian linear shared parameter model (2.3) with Si generated from (a), we computed the MLE An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx27_ht.jpg under (3.1) using the 2-stage estimation procedure with diagonal estimated covariance matrices and 1-step Newton–Raphson iteration. The REML estimators equation n1 and equation n2 were used as the initial estimators of α and β for the Newton–Raphson iteration. For the initial estimator equation n3 of γ, we used the subjects with Si observed (i.e. δi(c) = 0), computed the predicted values An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx28_ht.jpg and calculated equation n4 by the least squares estimate from the linear model An external file that holds a picture, illustration, etc.
Object name is biostskxq084fx29_ht.jpg.
Table 2 summarizes the averages of the estimates and their SEs, the mean squared errors (MSEs), and the empirical CP for the 95% bootstrap CIs based on the simulated samples with Si generated from (a). In each simulated sample, the SE and CP are computed using B = 200 bootstrap samples. The estimates of α and β obtained under the naïve mixed-effects model (2.1) have significantly larger biases and MSEs than the estimates obtained under the linear shared parameter model (2.3). The 95% bootstrap CIs constructed under (2.1) have low empirical CPs, while the CPs for (2.3) are close to the nominal level of 95%. These results suggest that, because of the potential model misspecification and estimation bias, the naïve mixed-effects models may lead to erroneous conclusions when a concomitant intervention is present, while the estimation bias can be dramatically reduced by a shared parameter model.
Table 2.
Table 2.
True parameter values, averages of the parameter estimates and their SEs, the MSEs, and the empirical CPs of the 95% CIs computed over the simulated samples generated from the Gaussian linear shared parameter model (2.3)
5.3. Estimation for additive shared parameter model
For samples with Si generated from (b), we used the additive shared parameter model (2.4) with Si = μ0(s)(ai0) + ϵi(s) and Gaussian errors. We approximated μ0(s)(ai0) using the quadratic B-splines with one interior knot at ai0 = 20 and estimated (αT,βT)T using the approximate likelihood (3.2) and the 2-stage estimation procedure with diagonal estimated covariance matrices and 1-step Newton–Raphson iteration. The estimated curve μ0(s)(ai0) was computed by the B-splines with the estimated spline coefficients. We computed the SEs and the 95% bootstrap CIs for (αT,βT)T using the “resampling-subject” bootstrap procedure with B = 200 replications. Pointwise 95% CIs for μ0(s)(ai0) were also computed using the same bootstrap procedure.
Table 3 summarizes the averages of the estimates and their SEs, the MSEs and the empirical CPs for the 95% bootstrap CIs based on the 200 simulated samples with Si generated from (b). Similar to Table 2, the estimates obtained under the naïve mixed-effects model (2.1) have mostly larger MSEs than the estimates obtained under the additive shared parameter model (2.4), and the CIs obtained under (2.1) have lower empirical CPs than those obtained under (2.4). The left panel of Figure 1 shows the true curve μ0(s)(ai0), and the spline estimated curve for μ0(s)(ai0) together with its 95% bootstrap CIs at 600 equally spaced ai0 points within [17,23] computed from one randomly selected simulated sample. The right panel of Figure 1 shows the empirical CPs of these pointwise 95% bootstrap CIs computed over all the simulated samples with n = 200, and the average of these CPs over the 600 equally spaced ai0 points within [17,23]. The spline estimated curve appears to be reasonably close to the true curve μ0(s)(ai0), and the CPs of the pointwise 95% bootstrap CIs are generally higher than the nominal level of 95% in the interior of the range of ai0.
Table 3.
Table 3.
True parameter values, averages of the parameter estimates and their SEs, the MSEs, and the empirical CPs of the 95% CIs computed over the simulated samples generated from the additive shared parameter model (2.4)
Fig. 1.
Fig. 1.
Left Panel: The thick solid line shows the true curve for μ0(s)(ai0). The thin solid and dashed lines are the spline estimated curve and the 95% pointwise bootstrap CIs obtained from one simulated sample. Right Panel: The solid line gives the (more ...)
Our approach suggests that concomitant interventions should not be treated the same way as a usual time-dependent covariate, and a shared parameter model may reduce the estimation bias by incorporating the self-selectiveness of the intervention. A number of extensions are worthwhile for further investigations. First, a patient may have multiple concomitant interventions in practice, so that a multiple change-point extension of (2.2) would be needed in such situations. Second, a worthy future development is to incorporate missing data and multiple concomitant interventions into a single framework of joint models. Third, asymptotic distributions of our estimators may provide useful insights and should be developed. Finally, since it may not be always clear whether an intervention is a concomitant intervention, model diagnostic methods for evaluating the appropriateness of a shared parameter model would be a valuable tool to be developed.
FUNDING
Financial support of the ENRICHD study was provided by contracts N01-HC-55140, N01-HC-55141, N01-HC-55142, N01-HC-55143, N01-HC-55144, N01-HC-55145, N01-HC-55146, N01-HC-55147, and N01-HC-55148 from the National Heart, Lung and Blood Institute.
Supplementary Material
Supplementary Data
Acknowledgments
We would like to thank Drs Myron Waclawiw and Minjung Kwak, and the reviewers for their thoughtful comments. Conflict of Interest: None declared.
  • Bickel PJ. One-step Huber estimates in linear models. Journal of the American Statistical Association. 1975;34:584–653.
  • Cai J, Fan J, Li R. Efficient estimation and inferences for varying-coefficient models. Journal of the American Statistical Association. 2000;95:888–902.
  • Davidian M, Giltinan DM. Nonlinear Models for Repeated Measurement Data. London: Chapman Hall; 1995.
  • Diggle PJ, Heagerty P, Liang K-Y, Zeger SL. Analysis of Longitudinal Data. 2nd edition. Oxford: Oxford University Press; 2002.
  • Investigators ENRICHD. Enhancing recovery in coronary heart disease patients (ENRICHD): study intervention rationale and design. Psychosomatic Medicine. 2001;63:747–755. [PubMed]
  • Investigators ENRICHD. Enhancing recovery in coronary heart disease patients (ENRICHD): the effects of treating depression and low perceived social support on clinical events after myocardial infarction. Journal of the American Medical Association. 2003;289:3106–3116. [PubMed]
  • Follmann D, Wu MC. An approximate generalized linear model with random effects for informative missing data. Biometrics. 1995;51:151–168. [PubMed]
  • Huang JZ, Wu CO, Zhou L. Polynomial spline estimation and inference for varying coefficient models with longitudinal data. Statistica Sinica. 2004;14:763–788.
  • Serfling RJ. Approximation Theorems of Mathematical Statistics. New York: Wiley & Sons; 1980.
  • Taylor CB, Youngblood ME, Catellier D, Veith RC, Carney RM, Burg MM, Kaufmann P, Shuster J, Mellman T, Blumenthal J. and others. Effects of antidepressant medication on morbidity and mortality in depressed patients after myocardial infarction. Archives of General Psychiatry. 2005;62:792–298. [PubMed]
  • Verbeke G, Molenberghs G. Linear Mixed Models for Longitudinal Data. New York: Springer; 2000.
  • Wall MM, Dai Y, Eberly LE. GEE estimation of a misspecified time-varying covariate: an example with the effect of alcoholism treatment on medical utilization. Statistics in Medicine. 2005;24:925–939. [PubMed]
  • Wu CO, Tian X, Bang H. A varying-coefficient model for the evaluation of time-varying concomitant intervention effects in longitudinal studies. Statistics in Medicine. 2008;27:3042–3056. [PubMed]
Articles from Biostatistics (Oxford, England) are provided here courtesy of
Oxford University Press