Estimating Rates of Change in Randomized Clinical Trials Nan M. Laird, PhD, and Fong Wang, ScD Harvard University, School of Public Health, Department of Biostatistics, Boston, MA
ABSTRACT: This article deals with the extension of the pretest-posttest clinical trial to the
longitudinal data setting. We assume that a baseline (or pretest) measurement is taken on all individuals, who are then randomized, without regard to baseline values, to a treatment group. Repeated measurements are taken postrandomization at specified times. Our objective is to estimate the average rate of change (or slope) in the expe r imental groups and the differences in the slopes. Our focus is on the optimal use of the baseline measurements in the analysis. We contrast two different approaches:-a multivariate one that regards the entire vector of responses (including the baseline) as random outcomes and a univariate one that uses each individual's least squares slope as an outcome. Our multivariate approach is essentially a generalization of Stanek's Seemingly Unrelated Regression (SUR) estimator for the pretest-posttest design (Am Stat 42:178-183, 1988). The multivariate approach is natural to apply in this setting, and optimal if the assumed model is correct. However, the most efficient estimator requires assuming that the baseline mean parameters are the same for all experimental groups. Although this assumption is reasonable in the randomized setting, the resulting multivariate estimator uses postrandomization data as a covariate; if the assumed linear model is not correct, this can lead to distortions in the estimated treatment effect. We propose instead a reduced form multivariate estimator that may be somewhat less efficient, but protects against model misspecification.
INTRODUCTION M a n y clinical trials are d e s i g n e d to estimate the rate of c h a n g e in a variable o v e r time a n d c o m p a r e the rates of change for two or m o r e groups. App r o a c h e s to the p r o b l e m of m o d e l i n g a n d estimating the rate of c h a n g e (and the effect of t r e a t m e n t on change) m a y differ in the w a y the baseline or p r e t r e a t m e n t m e a s u r e m e n t s are handled. In simplest p r e t e s t - p o s t t e s t designs (only one m e a s u r e m e n t is m a d e after treatment), baseline m e a s u r e m e n t s are u s e d b o t h to adjust posttest o u t c o m e s for p r e t r e a t m e n t differences a n d to reduce the variance of the e s t i m a t e d effects. The r e c o m m e n d a t i o n is often m a d e to use analysis of variance (ANOVA) on the differences in n o n r a n d o m i z e d studies a n d analysis of covariance (ANCOVA) on the posttest r e s p o n s e (or on posttest m i n u s pretest difference), using the pretest as a covariate in r a n d o m i z e d trials. The A N C O V A a p p r o a c h
Address reprint requests to: Nan M. Laird, PhD, Harvard University, School of Public Health, Department of Biostatistics, 677 Huntington Avenue, Boston, MA 02115. Received December 12, 1989; revised August 10, 1990. Controlled Clinical Trials 11:405-419 (1990) © N.M. Laird and F. Wang, 1990 655 A v e n u e of the Americas, N e w York, N e w York 10010
405
406
N.M. Laird and F. Wang is generally more efficient, and makes the proper adjustment for baseline differences in randomized trials by correcting posttest outcomes for regression to a common pretest mean [1,2]. In the absence of randomization, the groups may have different underlying distributions for the baseline measurements, and thus posttest measurements will not, in general, regress to a common mean. Thus a difference analysis is usually preferred in nonrandomized trials [3]. Stanek [4] considers a multivariate, Seeming!y Unrelated Regression (SUR) approach to estimation in the pretest-posttest design, which leads to the ANCOVA estimators. The SUR approach does yield a different estimated variance, essentially because the baseline values are treated as random variables rather than fixed as in ANCOVA, and thus the variances of the estimates of change are not conditional on baseline values. We consider here the extension of the pretest-posttest design to randomized longitudinal clinical trials where repeated measurements are taken postrandomization, and the main interest is in estimating the mean rate of change (slope) and its difference among the treatment groups. We assume that randomization is made without regard to baseline measurements, so that the population mean of the baseline measurements does not depend on group. The general approach that we take can be considered an extension of Stanek's SUR approach. Specifically, we will be concerned with three different analysis strategies: (1) fit a multivariate growth curve model to each individual's data vector, using generalized least squares (GLS) or maximum likelihood (ML) under a normality assumption: (2) reduce the data of each individual to least squares estimates of slope and intercept and use multivariate regression to estimate the treatment effects (reduced form multivariate): and (3) reduce the data as in (2), but analyze only the slopes using univariate methods (i.e., OLS and ANCOVA). We first consider the case where measurements are obtained on all individuals at the same set of time points. This assumption has the virtue of permitting closed form solutions for different estimators, allowing us to make direct comparisons among them. Although this is a typical study design, it often happens that, due to attrition and/or missed examinations, the data are incomplete. We use the term unbalanced design to denote this case, where an individual may be observed at an arbitrary subset of the prespecified times. This type of unbalanced data can also arise when individuals are measured at different time points after baseline. We will also discuss extensions to this unbalanced data setting. For simplicity, we consider designs involving only two groups, but our results are generalizable to any number of groups.
ESTIMATION WITH COMPLETE A N D BALANCED DATA: TWO-STAGE MULTIVARIATE MODEL
We will use a two-stage growth curve model [5,6] to develop our multivariate estimators. Let k = 1,2 index the groups, i = 1. . . . . nk index the individuals in each group, and j = 0 . . . . T index the occasions of measurement. For the first stage, the observation for the ith subject in the kth group on the jth occasion is modeled as Y~j = ~1~ + ~2~tj + e~/,
(1)
407
Estimating Rates of Change in RCTs
w h e r e to = 0. The (~k~, 132a) and the e~/s are mutually i n d e p e n d e n t sets of i n d e p e n d e n t r a n d o m variables, with
E(~sti ) = %~, cov(~gki, 13s'a) = dgg', cov(13xk,', ~g'k',',) = 0
E(ek,i) = O,
ki ¢= k'i',
var(ea/) = o?,
cov(ea, ek'ig') = 0
kij ¢= k'i'j',
for g,g' = 1,2, k,k' = 1,2, j,j' = 0 . . . . . T, and i,i' = 1. . . . . n~. The multivariate normality assumptions on (131a, ~z~) (second stage of the model) imply a final model of the form Ykij = al~ + c~2ktj + rkij,
(2)
where E(raj) = 0 var(raj) = o-2 + dH + ~d22 + 2tA12 and cov (rki/'ra/,)
= dll + titi,dz2 + (tj+tj,)d12 ki ¢= k'i'.
coy (rai,rk,&,) = 0,
The resulting covariance structure allows the variance to vary over time and the covariance to d e p e n d u p o n the time b e t w e e n m e a s u r e m e n t s . Model 2 implies that each Ya = (Yao. . . . . yhvc) T is multivariate normal with m e a n Z4xa, w h e r e Z is a ( T + 1) x 2 design matrix, and a covariance matrix that has a r a n d o m effects structure. Alternatively, we m a y express the model in terms of each individuals's ordinary least squares estimate of slope and intercept: is bivariate normal with m e a n "k, and variance-covariance matrix oa(ZTZ)-~ + D, w h e r e
g,g' = 1,2.
D = {da~,},
We refer to this latter representation as the reduced form. Using the r e d u c e d form, we can write the multivariate model as b
=
2n×l
X
a
2nxpp×l
+
8, 2nxl
where bT
=
( b l l 1. . . . .
bll,l,b121 .... b12,2,b2n ..... b21.1,b221 ..... b2a.2),
E(8) = 0, var (8) = I~ (~) I, and 1~ = D + o a (ZtZ) -~,
w h e r e ~ (~) I
[
2x2
nxn
o-nI
0"12I1
0"211
o2~lJ
denotes a (2n) x (2n) matrix of the form
(3)
408
N.M. Laird and F. Wang The design matrix X and the parameter vector ot d e p e n d u p o n the model assumptions below. Model E: Because we set to = 0 in model 2, setting ~Xll = ~x~2implies the means of the baseline m e a s u r e m e n t s are the same in both groups. U n d e r this a s s u m p t i o n of equal baseline m e a n s in the two groups, the vector ot can thus be written O~T
((X.1C, O~21,O~22),
=
w h e r e a~c is a c o m m o n intercept of the two groups and a2k is the slope for the kth group. The design matrix X u n d e r model E is 1
X=
00l
-}
,o, [i'
0 i 0 0
Ol 1 01) ! nl 10I 0 11}
i
~
xl
nx2
-=
Lnxl
X2 [ nx2
J
n2
901j
Model U: If the intercepts of the two groups are not equal, the vector o¢ is O¢T =
(Otll,Ot12,~21,a22),
w h e r e (C~k,~2k) are the intercept and slope for the kth group. For this case the design matrix is
0]
x 2n × 4
X2
MULTIVARIATE ESTIMATION ASSUMING KNOWN COVARIANCE PARAMETERS We n o w turn to estimation of cx and the slope difference, A = ~ 2 1 - 0t.22. It is well k n o w n that for model U, the GLS and ML estimators for a are the same as the ordinary least squares (OLS) estimator. Thus, we have for model U,
{~ML
:
~IGLS
=
~I~OI,S =.
(XTX)-IXFb
=
nk
w h e r e big. = E %dnk. i=l
The estimated difference of the slopes b e t w e e n the two groups and its variance are thus /~U ~___ b21.
-
b22.
409
Estimating Rates of Change in RCTs and var(~U) = ( 1
+ 1 ) 0.22,
(4)
where cr22denotes the (2,2) element of Z, as defined in eq. 3. Assuming o.2 and D are known, the GLS and ML estimators of ol are again equivalent for model E, and given by a ~
= a <~s =
_~,.-
\b22.
(0.,j0.,,)(_~,,.
-
(0.,2/0.,,)(b,2.
}_,..)
,
(5)
b,..)i
where-b1.. = (ni-bli. + n2bi2.)ln and or,i is the (i,j)th element of ]~. The estimate of A under model E and its variance are thus XE
=
~,.
_ ~..
_ (o-,do',,)(~,,.
-
~,~.)
(6)
and var (~ff)= ( 1
+
1)
0"9.2 (1
--
p2),
(7)
where ~)2 =
o.22/0.1i0.22.
U N K N O W N VARIANCE COMPONENTS
Although we have derived b E assuming that o2 and D are known, eqs. 5 and 6 still give the expression for the ML (and GLS) estimators of c~ and for the case where ~ is unknown and ML (or moment) estimators of the 0.ijs are u 3ed; eq. 7 also gives the asymptotic variance in this case. Closed-form expressions for ML estimators for model E are not available, although unbiased estimators are. Neglecting the functional dependence of ~ on 02 and D, the following estimators are unbiased for ]~: 2
nk
~u = ~_, ~ ( b k , - bk.)(bk,- bk.)V/(n--2), k=l
(8)
i=l
where nk
bk. = ~, bk,ln~, i=1
and ~E = ]~u except for the (1,1) element where we have ( n - l ) @lE1 = (n--2) @u + ~ nk (bik. -- bl..)2,
(9)
Ic
where &El and 6l~ are the (1,1) elements of ~E and ~u, respectively. Since /liE is the GLS under model E, we can use the (1,1) and (1,2) elements from ]~E to obtain a GLS estimator of di when ~ is unknown. We remark that these same estimators of ~ under the two models will be obtained if we use the
410
N.M. Laird and F. Wang moment unbiased estimators of o~ and D instead of X [7,8], provided that the estimate of D is nonnegative definite. In this case the ML and GLS estimators should be very dose. Note, however, that simply estimating I~ rather than o.2 and D has the advantage that the estimated variance of AE will be valid even if the random effects correlation structure for the data does not hold.
UNIVARIATE ESTIMATION From eqs. 5 and 6 we see that bE has the same form as a simple univariate ANCOVA estimator using b2~ as the response and bl~a as the covariate under the model b2~ = ~2k + ~(bl~i - bl..) + e~,
(10)
where the ekjs are independently distributed with zero mean and common variance or22 ( 1 - p2). However, under the univariate ANCOVA model, we would calculate the variance of the estimate of A as var(~a'4c) = ( 1 + 1 )
c r 2 2 ( 1 - p2)K,
(11)
where
z _ a ~ /k
Since K > 1, the variance of ~ is bigger under the univariate ANCOVA model than it is under the multivariate model E because we condition on the LS intercept in the former model, but treat it as a random variable in the latter. This result generalizes Stanek's result for the SUR estimator in a pretestposttest design. However, under model E, K is distributed as I plus (F~,~-2)/(n - 2), hence for large n, K should be close to 1. For the univariate ANCOVA analysis, ~/equals 0~u2/0~, so that the estimated slope on b,~ differs only slightly for the ANCOVA and GLS analyses. It follows from eq. 9 that ~/is always larger in absolute value than 0{2/~ (except possibly when n is small and group intercepts are nearly equal): thus the ANCOVA estimator makes a bigger adjustment for differences in intercepts than does the GLS, although the difference should be slight in randomized studies.
MULTIVARIATE GLS VERSUS UNIVARIATE ANCOVA: SOME TRADEOFFS A multivariate analysis (either full or reduced form) without the assumption of equal intercepts will generally be much less efficient than either the univariate ANCOVA or the multivariate under model E because it fails to make any covariate adjustment for baseline measurements. The multivariate GLS and the univariate ANCOVA approaches should give similar results in randomized trials, although the ANCOVA standard errors may be slightly larger. However, one drawback with the multivariate approach is that it requires
Estimating Rates of Change in RCTs
411
using post'treatment data (LS intercept) as a covariate in evaluating treatment effects. If the true growth pattern is not linear, then the estimated intercepts could depend quite strongly on the posttreatment responses. In the univariate model, this problem is solved by using y0~ as the covariate as opposed to blk~. Likewise, for the reduced form of the multivariate model, we can replace blh. by Y0~ X for the reduced form model becomes X~ = D + O2Wy, where the (1,2) (2,1), and (2,2) elements of Wy are the same as the corresponding elements of (ZTZ)-~, and the (1,1) element is 1. Since ~r~l = dn + o2
T-~
+ -t2/Y(tj _ ~)2 , < %11 it follows that
%n %=
o'no'z2
so that using the baseline rather than the estimated intercept will be less effident (especially if T is large) if the random effects model is correct, but it may be more desirable because of robustness to nonlinearity. It is natural to ask whether this version of the reduced form analysis (with Y~k replacing bin) corresponds to the appropriate analysis under some corresponding full multivariate model. The answer is no, although we can find a full multivariate model that yields the same ML or GLS estimator of A, A~ = b2~. - b=. - (%~/%1,)(,~01. - y0~.), but has variance equal to the corresponding univariate ANCOVA model variance using y0~ as a covariate, i.e., the analogous version of eq. 11 substituting py for p and y0~ for b~. This is accomplished by fitting the full multivariate model to Ya under the two-stage model, where at stage-two we assume E ( ~ ) = % + ~/y0ki. This is a bit counterintuitive, since y0~ actually appears both as a response (as the first element of y~ and as the covariate (at stage 2). An alternative approach would be to specify this multivariate model for the T vector of posttreatment responses only, but this could be quite inefficient, especially if T is small, as we lose leverage by omitting the baseline in the calculation of slope. For this reason, the univariate analysis should also use the LS slope calculated with the baseline value included.
GENERAL MULTIVARIATE MODEL
We have derived the multivariate results assuming a variance-components structure for X. We now briefly consider a general model, in which we assume normality with E(y~) = Z(x~ and var(y~) = F.
412
N.M. Laird and F. Wang Under model U, ML estimation is no longer equivalent to OLS, but rather is equivalent to ANCOVA on the LS slopes, using a ( T - 1) linearly i n d e p e n d e n t subset of ( I - z ( z T z ) -1ZT)Yki) as covariates [9]. By analogy with the two-stage model, u n d e r model E we would also include (blki - - b l ) in the covariate set for the optimal full multivariate analysis. Thus, the conclusion that ANCOVA on the LS slopes using b ~ as a covariate yields a fully efficient estimator u n d e r the multivariate model is not correct under the assumption of an arbitrary covariance structure. However, for the same reason that we might be reluctant to use blkJ as a covariate to analyze posttreatment differences, using the random-effects structure m a y be preferable to the more general one, since the optimal analysis u n d e r the arbitrary structure model uses postrandomization data as covariates and thus will be sensitive to model misspecification.
EXTENSION TO UNEQUAL SLOPES The usual analysis of covariance model assumes that the regression coefficient on the covariate is the same in all treatment groups. It is sometimes suggested that this assumption be tested before using the model, although interpretation of treatment effects in an ANCOVA model with different regression coefficients for each treatment group is difficult. The multivariate model offers a natural interpretation of the unequal slope model. Suppose we assume that the treatment m a y also affect variances, so that instead of a single D, we n o w have Dk, k = 1,2. To be consistent with the assumption under model E that the distribution of response at baseline does not d e p e n d u p o n k, we shall also assume that the (1,1) element is the same for all Dk and cr does not d e p e n d u p o n k, so o~n is the same for all k. The GLS (or ML) estimate of a u n d e r model U is still OLS. Deriving the GLS estimate of ot under model E, we now find ~
= b~.
-
~k(blk.
-
bl..),
where "Y~ = ~2/o'11,
k = 1,2.
Hence for the two-group problem with nl = n2, /~E becomes a~
=
a2~
-
~,~
=
(b2~.
-
bz;.)
-
~.(bl,.
-
b~2.),
with ~. = ('Yl+ "y2)/2" Further, the conditional distribution of b2k~given b~k~is normal, with mean azk - (~2/Crll)(blki- alk) and variance (r~2(1- PS): hence the ANCOVA model that parallels the multivariate model E assumes b2k, = a2k + "yk(blki -- b l ) + eki, where the eki are N(0,~).
Estimating Rates of Change in RCTs
413
EXAMPLE
To illustrate what the magnitude of the differences in these approaches might be in a realistic setting, we consider three examples. The first data set, taken from Potthoff and Roy [10], is dental measurements of 11 girls and 16 boys taken at four different ages (8, 10, 12 and 14 years). We code the ages as 0, 2, 4, and 6, respectively, and estimate I~ (pooling across sexes) as ]~ = ~4.435
L
-0.287]
0.118j"
Suppose we plan a clinical trial, randomizing children into two groups at age 8 (blocking on sex) and following them biannually until age 12; eqs. 4 and 7 allow us to calculate the variance ratio of GLS model E versus GLS model U as (1 - pa) = 84%. This also gives the expected relative variance of A N C O V A versus OLS on the slopes, where we use the estimated intercept as a covariate. If we use the baseline value as a covariate, the estimated expected relative variance is only 93%. As a second example, we use a simulated data set to compare the actual GLS versus A N C O V A estimated variances (eq. 7 and 11 with appropriate estimates substituted for (~2~ and p) w h e n the multivariate growth model is correct. The data are simulated using a model for a study of antiproteolytic replacement therapy among individuals with PiZ e m p h y s e m a developed by Wu and Carroll [11]. Here we take oa = 0.024, D = /'0.1521
\
-0.0124"~ 0.0083,/
and assume 13 equally spaced observations over 3 years. With nl = n2 100, u n d e r this model, (1 - p2) = 0.853. We find ~E = (0.1589
-----
-0.0171~ 0.0093,]
and ~ u = (0.1583
-0.0171~. 0.0093]'
hence ~E = 0.445, ~u = 0.446 and K = 1.0143, giving the following estimated variance ratios: GLSE vs GLSu = 81.8%, GLSE vs A N C O V A = 97.9%. This example confirms that there is little difference in the A N C O V A and GLS estimators w h e n n is large and we have randomized treatment assignments.
414
N.M. Laird and F. Wang If we use yok~ as a covariate in A N C O V A , the estimated variance ratio for GLSE versus A N C O V A d r o p s to 89.0%. A third e x a m p l e [12] uses results f r o m a s t u d y of the effects of e x p o s u r e to o z o n e versus clean air on forced expiratory v o l u m e in 1 second (FEV1). FEV1 w a s m e a s u r e d at baseline a n d at the e n d of each of six 50-minute periods. There are ten individuals in each group. (In fact, the s a m e ten individuals are u s e d in b o t h conditions in a crossover design, but w e ignore that complication here for illustrative p u r p o s e s only.) The data are s h o w n in Table 1, a n d are plotted in Figures 1 a n d 2. We use these data to illustrate two c o m m o n problems: d e p a r t u r e from linearity a n d unequal variances, b o t h evident in Figures 1 a n d 2. We note that a m o r e detailed analysis w o u l d consider t r a n s f o r m a t i o n s to linearity or m o r e complex models, a n d should also a d d r e s s the issue of outliers. In this setting the difference in baseline values is small (Yor - 702. = - 2 3 . 4 ml) relative to the difference in the estimated intercepts (bn. - b12. = - 7 2 . 7 mL), suggesting s o m e nonlinearity in the data. Table 2 s h o w s the estimates a n d their estimated variances a s s u m i n g equal variances. Again, the m o d e l estimates from the GLS m o d e l E a n d A N C O V A a p p r o a c h e s are v e r y close, but the A N C O V A estimated variances are bigger t h a n both the GLS variances a n d the OLS variance (particularly w h e n the estimated intercept is u s e d as a covariate) because the correlations are small a n d baseline differences are not. The a d j u s t m e n t for initial differences is less using yogi as a covariate rather than b~ki because of the smaller difference in Table 1
O z o n e Data Time
Individual
6 7 8 9 10 Average 1 2 3 4 5 6 7 8 9 10
Average
Group
0
1
2
3
4
5
6
C C C C C C C C C C
3971a 4598 3891 4283 4712 4432 3720 4076 4338 4362 4238
4162 4449 3956 4430 4685 4508 3622 3930 4110 4458 4231
4:222 4428 3869 4314 4557 4544 3762 3956 4008 4208 4187
4121 4544 3846 4380 4520 4416 3802 4059 4226 4145 4206
4160 4563 3912 4646 4602 4517 3854 4222 4074 3965 4252
4043 4618 3913 4555 4732 4766 3895 4193 4169 4171 4306
4165 4497 3933 4625 4555 4907 3826 4398 4016 4202 4312
T T T T T T T T T T
4132 4458 3730 4365 4599 4611 3853 4355 4339 4175 4262
4124 4370 3645 4397 4599 4685 3636 4120 4155 3834 4157
4027 4090 3784 4338 4533 4685 3636 4039 4134 3743 4101
3678 4328 3829 4412 4715 4647 3540 3928 4092 3362 4053
3690 4283 3664 4483 4647 4647 2644 4042 3939 3561 3960
3417 4276 3573 4393 4707 4472 2547 3923 3578 3056 3794
3590 4322 3483 4516 4641 4342 2874 3935 3319 2189 3721
"Response is FEVI(ml). Source: Follinsbee et al. [12].
415
Estimating Rates of Change in RCTs
Ix
iool
m
IX
-i
,
,
I
'
'
'
I
'
'
'
I
'
' ~+
O i 0)
si~+
45 ~
,
,
,
I
'
'
'
I
'+ '
'
,
~
I
'
'
45 41
39
35 )! ~6
27
2+
,
,
,
0
I
,
,
2
L
I 4
i
~
,
I 8
,
21 I
~ 0
,
,
I
,
J
,
2
I 4
,
I 6
till
tim
Figure 1 Individual FEV1 Measurements: Control Group
Figure 2 Individual FEV1 Measurements: Exposure Group
group means. Notice that the estimated variances show that analyses using y ~ as a covariate are m o r e efficient than those using bt~. This is a reflection of two features: (1) the nonlinearity has induced a difference in intercept means larger than that in baseline means, and (2) a larger correlation between slope and baseline than between slope and intercept suggests the randomeffects covariance structure [which implies corr(bt~,b2~) > corr(y0~,b2~)] does not accurately reflect the empirical covariance. Table 3 shows the same results, now assuming D and hence ~ are different in the groups. Although estimates of o~ do differ, the OLS and GLS results are nearly unchanged for estimating ix. The estimated ANCOVA variances are somewhat smaller, since the correlation increases when estimated separately for each group. Table
2
The OLS, GLS, and ANCOVA Estimates" of the Mean Slope for Model U and Model E: Ozone Data
OLS (Model U) GLS (Model E) B1
Y0 ANCOVA B1 Y0
Control
Exposure
Difference
Estimated Variance
15.57
-88.83
104.40
1168.40
16.94 16.24
-90.21 -89.50
107.15 105.74
1144.82 1110.32
16.89 16.21
- 90.16 -89.47
107.05 105.68
1233.84 1180.70
"Assuminga commonvariancestructure in each group.
416
N.M. Laird and F. Wang Table 3
The OLS, GLS, and ANCOVA Estimates a for the Mean Slope for Model U and Model E" Ozone Data Control
Exposure
Difference
Estimated Variance
15.57
- 88.83
104.40
1167.23
13.88 15.11
- 93.60 - 90.83
107.49 105.91
1141.88 1102.61
14.21 15.21
- 93.19 - 90.63
107.41 105.84
1191.06 1076.86
OLS (Model U) GLS (Model E) B1 Y1 ANCOVA B1 Y1
aAssuming variance and covariance parameters depend upon group.
UNBALANCED DATA SETS
The general conclusion from the balanced data situation is that there is little advantage to using a multivariate approach rather than a simple ANCOVA estimator. The situation is more complicated w h e n we have unbalanced designs, as will h a p p e n w h e n we have missing data or attrition. The multivariate approaches generalize readily to this setting, although the computations are more complex. If we assume o-2 and D are k n o w n , then one can show that for model U, the GLS or ML estimate of slope for the kth group is a weighted combination of both the slopes (b2~) and the intercepts (blk3 for individuals in the kth group. For model E, it is a weighted combination of individual slopes and intercepts from all the groups. Although one can derive algebraic expressions for the estimators of a and its variance u n d e r both models, they are not easily interpreted or amenable to comparison. W h e n o-2 and D are u n k n o w n , they must both be estimated for a multivariate analysis, since we require them to determine the weight matrix Wk, = (ZTZki)1 03 + D associated with bki. Here Zkt gives the design matrix for the kith individual. Maximum likelihood computations are complex, but simple unbiased estimators exist in closed form [7]; thus estimators can easily be compared via simulation. The univariate analysis of covariance approach is n o w more complicated, since we no longer have homoscedastic errors for the ANCOVA model given eq. 10. If we assume that the underlying two-stage multivariate growth model is correct, the variance of ek, for model 10 is a function of both o'2 and D (as well as wk,), hence to implement ANCOVA properly would still require estimation of both o'2 and D, and thus computationally would be no simpler to implement than GLS. An ad hoc approach to univariate ANCOVA in this setting would be to use weighted ANCOVA u n d e r model 10, assuming ea has variance OZWk~ + d, where wa is the (2,2) element of Wk~. This is similar to an estimator proposed by Wu and Bailey [13] for a s o m e w h a t different setting. The estimated parameters can be obtained by a weighted least squares analysis of model 10, plugging in m e t h o d of m o m e n t estimates for 03 and d. (For details, see the Appendix.)
417
Estimating Rates of Change in RCTs
Table 4
Simulated Data Multivariate Model, GLS Estimators oa and D Known
var (~u) vat (bE) Ratio
oa and D Estimated
No Attrition
50% Dropouts
No Attrition
50% Dropouts
0.000221 0.000186 84%
0.000253 0.000220 87%
0.000221 0.000183 83%
0.000270 0.000242 90%
To study the effect of using model U versus model E in the multivariate setting for unbalanced data, we used the simulation model previously described, except we assume that approximately 50% of the participants drop out of the study in a random fashion (uniform dropouts over time). The results based on 600 simulated data sets are given in Table 4. The variance ratio is 90%, using variances estimated via simulation. Had there been no attrition, we would expect a variance ratio of 0.83. To see how much of the loss in efficiency is due to estimation of 0-2 and D as opposed to attrition, we also computed the variance of/~ under models E and U, assuming 0-2 and D are known. This gives a variance ratio of 87% (Table 4), compared to 84% when there is no attrition. The variance for the weighted ANCOVA (with 50% attrition and oa and d estimated) is 0.000269, so the GLSE versus ANCOVA variance ratio is 90%. Now the ANCOVA approach does not compare favorably with the GLSE approach because the approach is using incorrect weights for the simulated data. In clinical trials it often happens that data are missing for nonrandom reasons, e.g., treatment failure, adverse events, etc. In general, moment type estimates will give biased results in this case, although the nature and extent of bias clearly depends upon the reasons for the missing data. Maximum likelihood estimators are readily available for the models we propose, and will give protection against bias under some types of nonrandom missingness processes, namely those where the probability of missingness depends only on observed characteristics, e.g. previous responses, or covariates included in the model. See Laird [14] for a comprehensive discussion of missing data in longitudinal studies. DISCUSSION Our comparison of the multivariate and ANCOVA approaches to estimation of treatment effects in linear growth curve models leads to the conclusion that a simple ANCOVA on the least squares slopes, using the estimated intercept as a covariate, is nearly optimal, provided we have done a randomized trial and are confident in the assumed linear model. Using a multivariate model for analysis (under model E), will give a smaller estimate of the variance, since the intercepts are treated as random in computing the variance. Although it is optimal from the point of view of efficiency to use the estimated intercepts as a covariate, rather than just the baseline values, using the baseline values protects against possible nonlinearity and avoids using postrandomization data as a covariate. In the absence of randomization, these approaches are only appropriate if
418
N.M. Laird and F. Wang we are willing to assume the intercept means are equal. Although this m a y seem s o m e w h a t counterintuitive, the interpretation is related to Lord's paradox, namely that covariance adjustment is inappropriate if individuals in different groups regress to different pretreatment means. Hence, the reco m m e n d a t i o n is sometimes m a d e that one substitute a difference analysis for a covariance analysis in n o n r a n d o m i z e d pretest-posttest studies (see Bock [2], for example). This would correspond to using model U in the multivariate analysis. W h e n T= 1, /~u becomes the difference in the mean change (time 1 minus baseline) between treatment and control group. The advantage of a multivariate approach over the ANCOVA is that it is straightforward to extend w h e n we have unbalanced data. The appropriate ANCOVA model for this case appears no less complex to implement than GLS; a simpler ad h o c - w e i g h t e d ANCOVA model appears to lose efficiency relative to the multivariate approach. Comparing the multivariate approach and OLS, our limited simulation results suggest that the efficiency gains m a y be slightly diminished with unbalanced data sets. The advantage of the reduced form multivariate model, i.e., assuming a linear model for b~ with cov(bk,) = O'2(ZTkiZk,)-1 + D, rather than a full multivariate model for Yk~, is that we can easily adapt the approach w h e n the response vector is (Y0~,b2k~) rather than (blk.,b2~). We recommend this adaptation because it should protect better against nonlinearity, and avoids using post-randomization data as a covariate, even t h o u g h there is some loss of efficiency if the two-stage growth model is correct. We have limited our discussion to models which assume linearity in the response; obviously tj in models 1 and 2 can be replaced by any arbitrary function of time, e.g., ln(1 + t/) and the results still hold. Our models can also be extended in a straightforward way to allow polynomial models in tj (or a suitable transformation). The same general results hold concerning the equivalence and relative efficiency of the three approaches. Even t h o u g h the models no longer specify linearity, it is still the case that estimated intercepts are sensitive to model misspecification, so that using Y0,k as a covariate is to be preferred over b0ik. This research was sponsored by grant GM29745 from the National Institutes of Health, Department of Biostatistics, Harvard School of Public Health, 677 Huntington Avenue, Boston, Massachusetts 02115.
REFERENCES 1. Laird N" Further comparative analyses of pretest-posttest research designs. Am
Star 37: 329-330, 1983 2. Bock RD: Multivariate Statistical Methods in Behavioral Research. New York,
McGraw-Hill, 1975, pp 493-495 3. Lord FM: Lord's paradox. In Anderson SB, et al. (eds): Encyclopedia of Educational Evaluation. San Francisco, CA, Jossey Bass, 1973 4. Stanek EJ: Choosing a pretest posttest analysis. Am Star 42: 178-183, 1988 5. Rao CR: The theory of least squares when the parameters are stochastic and its applications to the analysis of growth curves. Biometrika 52: 447-458, 1965 6. Lange N, Laird NM: Random effects and growth curve modeling for balanced and complete longitudinal data. J Am Star Assoc 84:241-247, 1989
Estimating Rates of Change in RCTs
419
7. Reinsel GC: Mean squared error properties of empirical Bayes estimators in a multivariate random effects general linear model. J Am Star Assoc 80: 642-650, 1985 8. Vonesh ER, Carter RL: Efficient inference for random-coeffident growth curve models with unbalanced data. Biometrics 43: 617-628, 1987 9. Grizzle JE, Allen DM: Analysis of growth and dose response curves. Biometrics 25: 357-381, 1969 10. Potthoff RF, Roy SN: A generalized multivariate analysis of variance model useful especially for growth curve problems. Biometrika 51: 313-326, 1964. 11. Wu MC, Carroll RJ: Estimation and comparison of changes in the presence of informative right censoring in modeling the censoring process. Biometrics 44: 175188, 1988 12. Follinsbee LJ, McDonnell WF, Horstman DH: Pulmonary function and symptom responses after 6.6-hour exposure to 0.12 ppm ozone with moderate exercise. J Air Pollut Control Assoc 38: 28-35, 1988 13. Wu MC, Bailey RJ: Analysing changes in the presence of informative right censoring caused by death and withdrawal. Stat Med 7: 337-346, 1988 14. Laird NM: Missing data in longitudinal studies. Star Med 7: 305-315, 1988
APPENDIX
In this a p p e n d i x w e give the estimation formulas for the weighted ANCOVA model b2~ = o,2k + ~(bl~ - bl..) + e~, w h e r e we use weights u~ - (o2w~ + d ) - L
(A1)
We first estimate o 2 by pooling the residual sums of squares from each individual regression. To estimate d, we first p e r f o r m an u n w e i g h t e d A N C O V A to obtain provisional estimates ot2t a n d ~ and t h e n form MS = ~ ~ (b2~ - &2k - ~/(b~ - ~..))2 k
i
and define ~/ = max{[MS - 0-2 1~ X k
wkjn],
0}.
i
We t h e n replace o.2 and d by their estimates in eq. A1 to d e t e r m i n e ~k~, and minimize Z (b2~ k
a2, -
v(bl~, -
k
to d e t e r m i n e 0,2k and y as ~
k
i
k
i
and a~ = ~
-
~(~y'~ - ~ , ) .
H e r e ~/t. : ~ b/~ Oe/~, ilk..
~.))2~