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 *S*_{i} 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 *Y*_{ij}, *T*_{ij}, *S*_{i}^{(c)}, and *R*_{ij} = *T*_{ij} − *S*_{i}^{(c)} be the *i*th patient's BDI score, trial time (months), starting time (months) of antidepressant use, and antidepressant duration time (months), respectively, at the *j*th visit. For all 1 ≤ *i* ≤ *n*, *T*_{i1} = 0, *T*_{ini} is the trial time at the last visit, and *λ*_{ij} = 1_{[Si < Tij]} is the indicator of antidepressants use at *T*_{ij}. The observed (*S*_{i}^{(c)},*δ*_{i}^{(c)}) is (*S*_{i}^{(c)} = *S*_{i},*δ*_{i}^{(c)} = 0) if the *i*th patient used antidepressants within the treatment period, (*S*_{i}^{(c)} = *T*_{ini},*δ*_{i}^{(c)} = 1) if the patient did not use antidepressant within the treatment period, and (*S*_{i} = 0,*δ*_{i}^{(c)} = 2) if the patient used antidepressants before baseline. When the linear models *μ*_{0}(*T*_{ij};*a*_{i0},*a*_{i1}) = *a*_{i0} + *a*_{i1}*T*_{ij} and *μ*_{1}(*R*_{ij};*b*_{i0},*b*_{i1}) = *b*_{i0} + *b*_{i1}*R*_{ij} are used, (*a*_{i0},*a*_{i1}) represents the intercept and slope of the *i*th subject's BDI trajectory before antidepressant use, and (*b*_{i0},*b*_{i1}) 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 *Y*_{ij} and *T*_{ij} 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

where (

*α*_{0},

*α*_{1},

*β*_{00},

*β*_{10})

^{T} is the unknown mean vector and Σ is the unstructured covariance matrix for the multivariate normal distribution

(·,·). 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} +

*β*_{10}*R*_{ij} indicates ignorable antidepressant effect on the mean BDI scores.

The left panel of 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

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

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.**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

*S*_{i}, (

*a*_{i0},

*a*_{i1}) and (

*b*_{i0},

*b*_{i1}) 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:

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 *σ*_{s}^{2}, 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 (*a*_{i0},*a*_{i1})^{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 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 , which gives

and

. Contrary to the interpretation from (4.1), we obtain from this submodel that

since 0 ≤

*R*_{ij} ≤ 6 and

, 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).