STATISTICS IN MEDICINE, VOL. I I. 1825-1839 (1992)

AN OVERVIEW OF METHODS FOR THE ANALYSIS OF LONGITUDINAL DATA SCOTT L. ZEGER A N D KUNG-YEE LIANG Department of Biosratistics. Johns Hopkins University, 615 North Wolfe Streer. Baltimore. Maryland 21205. U . S . A .

SUMMARY This paper reviews statistical methods for the analysis of discrete and continuous longitudinal data. The relative merits of longitudinal and cross-sectional studies are discussed. Three approaches, marginal, transition and random effects models, are presented with emphasis on the distinct interpretations of their coefficients in the discrete data case. We review generalized estimating equations for inferences about marginal models. The ideas are illustrated with analyses of a 2 x 2 crossover trial with binary responses and a randomized longitudinal study with a count outcome.

INTRODUCTION Longitudinal data comprise repeated observations over time on each of many individuals. Longitudinal data are in contrast to cross-sectional data where only a single response is available for each person. The statistical analysis of longitudinal data presents special opportunities and challenges because the repeated outcomes for one individual tend to be correlated with one another. In many studies, the scientific question can be formulated in statistical terms as a regression which can be estimated from either longitudinal or cross-sectional data. Longitudinal studies typically have two important advantages, however: increased power, and robustness to model selection. The following trivial example illustrates the first advantage, increased power. Consider two simple studies to compare treatments A and B for reducing blood pressure. In design I, patients are randomized to treatment A or B, treated for two weeks and then blood pressure is measured for each person. The treatment effect is estimated by the difference between the group means. In design 11, patients are randomized to either treatment, their pressures are measured at baseline, treatment is administered for two weeks and B P is then measured a second time. The treatment effect is estimated by the average change in B P for group A minus the average change for group B. Note that design I is cross-sectional; design I1 is longitudinal. Design I1 has the same number of individuals as design I but twice as many observations. Suppose that BP measurements from distinct individuals are independent with variance a’ and that the correlation between repeated observations on a single individual is p. Then the variance of the cross-sectional estimate of the treatment effect is 2 a 2 / n where n is the number of persons in each group; the variance of the longitudinal estimate is 2a’(l - p ) / n . Hence the gain in efficiency by measuring the B P prior to treatment is 1/(1 - p ) . For p = 0.5, a common value for blood pressures two weeks apart, the longitudinal study is twice as efficient as the cross-sectional study 0277-6715/92/141825-15$12.50 0 1992 by John Wiley & Sons, Ltd.

1826

S. L. ZEGER A N D K.-Y. LIANG

with the same number of individuals. With p > 0.5, the longitudinal study is more efficient than a cross-sectional study with twice as many people and the same number of measurements. Heuristically, the longitudinal study allows each person to be used as his own control so that the ever-present heterogeneity among people is reduced. Longitudinal data also enable us to avoid some problems of model misspecification, a serious concern in non-experimental settings where treatments cannot be assigned at random. Suppose the ‘truth’ can be described by the equation Yif = x::,fl

+ z;cr +

Ei,

(1)

where xi, and ziare time-dependent and time-independent explanatory variables with regression coefficient fl and a, respectively, and is an error. Now suppose xi, can be observed but zi is missing. Given only y i l as in a cross-sectional study, fl cannot be estimated consistently unless z is uncorrelated with x , but with repeated observations, fl can be estimated without bias since Yi,- Y i , = ( x i , - x i l - c i l . The essential idea is that longitudinal designs permit analyses which are insensitive to omitted covariates that do not change with time. Inferences can therefore be made robust to some forms of selection bias common in observational research. Of course, there is an analytic cost associated with longitudinal data, namely, the analysis must take account of the correlation on repeated observations to obtain valid inferences about regression coefficients. Ignoring correlation when it exists results in two problems: inefficient estimates of regression parameters, and, more important, inconsistent estimates of precision. Figure 1 illustrates both problems for a simple regression in which a Gaussian response is a linear function of time and the errors follow a first-order autoregressive model = p i t- + a, where ail are independent Gaussian variates with zero mean and common variance. The figure shows three curves plotted against the degree of correlation p: (A) the variance of the efficient estimator of the slope, which takes account of the correlation; (B) the true variance of the ordinary least squares estimate of slope, which ignores the correlation; and (C) the incorrect estimate of the variance of the ordinary least squares slope one obtains when correlation is ignored. The inefficiency of ordinary least squares is seen by comparing curves A and B. Note that the ordinary least squares slope is reasonably efficient even when there is high correlation. See Bloomfield and Watson’ for details on the inefficiency of least squares. The inaccuracy of variance estimates based upon the independence assumptions can be seen by comparing curves A and C. When the correlation is positive (negative), ignoring correlation leads to an under(over)-estimate of the variance. In short, ignoring correlation when it exists can lead to grossly incorrect inferences about regression coefficients. Before turning to methods of analysis, we introduce three longitudinal data sets which are typical of the problems that arise in biomedical research. The statistical methods are illustrated throughout using these examples.

)’a +

Example 1: Children’s morbidity and vitamin A Vitamin A maintains epithelial cells, the first line of defence against infection. An important question in international health is whether vitamin A deficiency is associated with increased respiratory and diarrhoea1 disease. Sommer et al.’ have addressed this question with longitudinal data on over 3000 pre-school Indonesian children. Each child was visited quarterly for 18 months to assess whether he had respiratory infection and was vitamin A deficient. Explanatory variables include weight, height and age which are time-dependent as well as gender.

1827

OVERVIEW O F METHODS FOR ANALYSIS O F LONGITUDINAL DATA

2 0 1

.

1

, 1

: :

.’

z0 aD

880

8 -

..*

0

.”

-.*Her

-.**.---- - - - - - - - - - -,,,,.- ..-...“’ - - -‘“L_

ru

,,.......’ ......

8-

,

A: Best possible

,.. ... ,,......

I

I

I

I

I

First Lag Correlation

+

Figure 1. Variance of slope estimates on regression Y , = &, Dl r + cir, ci, = pi,-+ a, where a , are independent Gaussian (0, 1) variates r = - 3, - 2, . . . , 2, 3. Curve A is the variance of the efficient estimate (MLE); Curve B is the true variance of the ordinary least squares (OLS) estimate; Curve C is the incorrect variance reported by OLS

Example 2 Crossover trial Table I presents data from a 2 x 2 crossover trial3 to test an active treatment against a placebo for relieving cardiovascular deficiency. Half of the participants were randomized to receive the placebo followed by the active drug; half received the active treatment followed by the placebo. The response measure is whether or not there was a normal electrocardiogram. The frequency of each of the four possible bivariate outcomes is given in the table. Example 3 Epileptic seizures The final example is a clinical trial to compare progabide with placebo for reducing seizure rates in epileptics. The data are presented by Thall and Vai14 and recently reanalysed by Breslow and Clayton.’ They are displayed in Figure 2. The study began with an eight week run-in period during which seizures were recorded. The 59 patients were then randomized to treatment or placebo and followed for four additional two week periods. The number of seizures in each interval is presented in the figure; the run-in interval responses are divided by four to normalize units to seizures per two weeks.

1828

S. L. ZEGER AND K. Y. LlANG

Table I. Summary of regression coefficient interpretations when prevalence is small Model

Parameter

Interpretation

epl

ratio of population prevalences ratio of population incidences ratio of odds (risk) for individual subject

Marginal Transition Random effects

,fly

eflY

T I I

T

T

I I

c B 0Q c I I

I

I

I

T

I I

I I

I

T

T

T

I

I I

1

1

I

Plac Interval:

III

I

Trt

B

I

I

I I I

I

I I

I

I

I

I I I

1

1

1

1

1

1

1

1

Plac

Trt

Plac

Trt

Plac

Trt

Plac

Trt

1

2

I

3

4

Figure 2. Box-plots of seizures and rates (numberi2 weeks) for treatment and control group for baseline and four follow-up intervals. Data on 59 patients are presented by Thall and Vai14

In each of the examples, the scientific question can be formulated as a regression problem with the regression coefficients being of primary interest. The correlation must be taken into account but is not itself a focus. Examples 1 and 2 have binary responses. In Example 3, the outcome is a count. This paper reviews three approaches for analysing longitudinal data like those in these examples: marginal, transition and random effects models. After comparing the distinct approaches using Example 1, we focus on inference for the marginal approach. The ideas are illustrated with analyses of Examples 2 and 3. First, we review regression models for independent observations.

OVERVIEW OF METHODS FOR ANALYSIS OF LONGITUDINAL DATA

1829

REVIEW O F GENERALIZED LINEAR MODELS In biomedical applications, it is common to confront both discrete and cqntinuous outcome measures. Regression models for independent, discrete and continuous responses have been unified under the class of generalized linear models (GLM),6which has facilitated data analysis by providing a common set of methods regardless of the type of response. We will describe below extensions of GLMs for correlated outcomes. First, the salient features of the GLM class are reviewed. Consider the cross-sectional situation with response Yi and p x 1 vector of explanatory variables x i , i = 1 , . . . , K. The objective is to describe the dependence of the mean response p i = E(Yi) on the covariates. A GLM is a member of the exponential family with likelihood function of the form f ( y i ) = exp{(yiei + a(ei))/d + c ( y i , 4)) (2) where the mean pi equals the first derivative of the function and variance of Yi is proportional to the second derivative of a(&). The mean pi is related to the explanatory variables by g ( p i ) = x$ where g is the ‘link’ function. The variance is related to the mean by var( Y i )= ui = u ( p i ) 4 where u is called the ‘variance’ function. The GLM family includes linear, logistic, log-linear and some parametric survival regression models as special cases. For example, in logistic regression, logit(pi) = xi8 and = pi. var( Yi) = pi(l - pi); in log-linear models for count data, log(pi) = xi8 and var( Yi) In each GLM the regression coefficients B are estimated by solving the same estimating

The solution can be obtained by iterative weighted least squares and is a consistent estimate as long as g(pi)= x;B whether or not the variance function is correctly specified. See McCullagh and Nelder6 for a complete discussion. Wedderburn’ first pointed out that the estimating equation above provides a consistent estimate of p with a variety of link and variance functions whether or not they correspond to a member of the exponential family of distributions. The name ‘quasi-likelihood estimate’ was coined for the solution of equation (3) in the more general case since its integral does not necessarily constitute a proper likelihood function. STATISTICAL MODELLING OF CORRELATED RESPONSES We have two objectives for statistical models of longitudinal data:

1. to adopt the conventional regression tools, which relate the response variables to the explanatory variables; and 2. to account for the within subject correlation. In most longitudinal studies, the regression objective is of primary interest. While it is essential to account for the within-subject correlation, the nature of this correlation is often of secondary interest. To introduce the three approaches to longitudinal data analysis, let Yi = ( Yi,,. . . , Y,,, . . . , Y,,J

be a ni x 1 vector of repeated responses for the ith subject, i = 1, . . . , K. At t , one observes the response Yi, as well as a p x 1 vector of explanatory variables, xi,, thought to be related to Yit.

1830

S. L. ZEGER AND K.-Y.LIANG

Marginal models In this approach, the regression and within-subject correlation are modelled separately. Specifically, we assume: 1. the marginal expectation of Y i r ,E ( Yi,) = pit,is related to x i , by g(pit) = x ; t B >

where g is a known link function such as the logit function for binary responses; 2. the marginal variance is a function of the marginal mean, that is, var( YiO = u(pir)4> where u is a known function and 4 is the over-dispersion parameter which accounts for the variation of Y i , not explained by v(pi,); 3. the covariance between Y,, and Y i t ,s < t = 1,. . . , ni is a function of the marginal means and additional parameters a, that is, COV(

Yis, Yit) = C(pisr Piti a),

where c is a known function. Marginal regression coefficients have the same interpretation as coefficients from a crosssectional analysis. This modelling approach is a natural analogue for correlated responses of GLM's or quasi-likelihood for independent data. In the Indonesian study (Example I), we are interested in whether the risk of respiratory infection is related to vitamin A deficiency as manifested by the presence of xerophthalmia (an eye condition). In a marginal model, we might assume 1. logit pi, = Po + PI xeroi, + B2agei,, where pit = E ( Y,,)= Pr( Y,, = I), and Y,, and xeroit are indicator variables for the presence of respiratory infection and xerophthalmia, respectively, and agei, is the age in months for child i at the tth visit. 2. var( Y i t )= p i r ( l - pi,), and 3. corr( Yi,, Yi,) = als-tl, s < r, that is, that. the correlation coefficient between two repeated observations from the same subject depends only on the time between the two visits. Alternatively, one might assume log OR ( Yi,,Y i , )= als-rl, where

OR( Yis, Yi,) =

Pr( Yi, = Y i , = l)Pr( Yi, = Yit = 0) Pr(Yi, = 1, Y,, = O)Pr( Y i , = 0, Y i , = 1)'

(4)

In either case, the covariance between Y,, and Y,, is fully specified by assumptions 2 and 3 above. Remark I. In this example, eBIhas the interpretation as approximately the ratio of prevalences of respiratory infection from two populations of children, those with and without xerophthalmia at a given visit. In other words, g

prevalence of respiratory infection among children with xerophthalmia prevalence of respiratory infection among children without xerophthalmia

'

The marginal model regression coefficients have what we call 'population-averaged' interpretations' because they contrast odds of disease in the populations with and without the risk factor. Remark 2. Note that the effect of age, p2, in Example 1 can be estimated cross-sectionally if there is heterogeneity in the age distribution of the children sampled at one time. However, this

OVERVIEW OF METHODS FOR ANALYSIS OF LONGITUDINAL DATA

1831

cross-sectional age effect may be subject to bias if so-called 'cohort effects' (referred to as zi in equation (1)) are present in the study population. This bias can be detected and corrected with longitudinal data by expanding the model to include age at baseline and change in age from baseline as covariates. The coefficient for age at baseline is sensitive to cohort effects; that for change in age is not. It is of interest to examine if, and in what way, these two age coefficients are different. Remark 3. Because p describes the effect of covariates on the marginal expectation of the Y's, it has the same interpretation regardless of the number of repeated observations, ni, which may vary among subjects. Remark 4 . The third assumption is needed to account for within-subject correlation. However, the magnitude of this correlation, indexed by a, does not alter the interpretations of p. Transition (Markov) models

This approach and the next are different from the marginal model in that they attempt to address both the regression objective and the within-subject correlation simultaneously. In other words, parameters for the dependence on x and for correlation are introduced on the same scale in a common equation. Transition models assume 1'. the conditional expectation of Y i l ,p:l = E ( YilI Y i l - l ,. . . , Y i l ) ,depends on xi1 and past

responses as follows: V

g(pFt) = &p*

+ 1 aTfi( Yiz-

1t

. . . YilL 9

j= 1

where {fi,j= 1, . . . , v ) are known functions. 2'. the conditional variance of Yi, given the past is a function of p f l ; that is, var( Yit I Yit- 1, . . . Yi1) = v(pL&)4, 9

where u is a known function. In the Indonesian study, we might assume

1'. logit Pr( Yi, = 1 I Yil- . . . , Yil-) =

+

xeroil

+ p: ageil + a* yi1

- 1,

2'. var( Yil1 Yil- ) = pfl( 1 - &). In this special case, v = 1, so that Yi,is assumed to depend upon the past responses only through the immediately preceding response. Remark 1 . The table below gives the transition probability of respiratory infection according to 1' in Example 1.

1832

S. L. ZEGER AND K.-Y. LIANG

Suppose the disease duration is similar to the sampling interval. Then because of the adjustment for Y i t - in l’, efl: is approximately the ratio of incidences from two groups, one with xerophthalmia and one without. In other words, efl: ‘v

Incidence of respiratory infection among children with xerophthalmia Incidence of respiratory infection among children without xerophthalmia‘

Remark 2. Examining the expression 1’ closely, we note that the interpretations of P* changes as v and the functional forms change. This is to be contrasted with Remark 4 about the parameters in marginal models. Remark 3. Transition models are easily fit with standard software by treatingfl, . . , ,f v along with x, as the set of regressors.

If,}

Random effects models In the random effects model, the correlation among responses for an individual is assumed to arise from natural heterogeneity in regression coefficients across people. Given one person’s coefficients, the responses are assumed to be independent observations on a generalized linear model. The collection of individual coefficients are assumed to be like a sample from an underlying distribution. Formally, we assume: 1”. the conditional distribution of Yi, given bi satisfies a GLM with g ( E ( Yitl hi)) = .YI,P**

+ Zithi,

where zit is a q x 1 vector of covariates whose coefficients vary across people; zit is in general a subset of xi, so we can adopt the convention Ebi = 0. 2“. the conditional independence of Y i l , . . . , Yi,, given b i , 3”. a specific form for the distribution of hi, that is, F ( - ; 0). A traditional choice is the q-dimensional Gaussian distribution with mean zero and covariance matrix indexed by 6. In the Indonesian study, we might assume a random intercept model as follows: xeroir + P:* agei,. where zit = 1. 1”. logit Pr( Yi, = 1 Ibi) = Po** + hi + 2”. Var( YitI b i ) = E ( Yi,I h i ){ 1 - E( Yi,I h i ) ) . 3”. hi N(0, 6).

-

Remark I. In Example 1, eP;’ has interpretation as the odds ratio of respiratory infection for a child who has xerophthalmia versus the same child if he does not. In other words, it describes the change in the risk of respiratory infection for a child whose xerophthalmia status changes. Note the odds ratio is assumed in the model to be the same across subjects. Remark 2. The value of P** depends on the assumed distribution F foi the hi’s. For discrete longitudinal data with small ni, it is difficult to check the validity of F . It is therefore good practice, but difficult, to examine how sensitive the results are to the specified distribution for the hi’s. Table I summarizes the interpretation of the regression coefficients P, P* and fl** using the variable xerophthalmia to illustrate. Given different interpretations for coefficients in each model, their magnitudes are expected to be different as well. Assuming a random effects model with random intercept, Neuhaus er aL9 show that for any well-defined distribution F of the bi’s with var(bi) > 0. IPI < IP**I

OVERVIEW OF METHODS FOR ANALYSIS OF LONGITUDINAL DATA

with equality holding if and only if B**

= 0.

In particular, if bi

-

1833

N ( 0 , O), Zeger et aL8 show that

logit Pr( Yi, = 1) 2: x:,p**/(l + 0*350)1/2. Thus, the logit of the marginal probability is roughly a linear function of xir with regression coefficients, p, given by

p 2: p**/(i

+

0.350)1'2.

It is clear that p and p** are approximately equal if either p** = 0 or 8 = var(bi) = 0. The relationship between /3 and p* can be established in special cases. Returning to Example 1 and ignoring the age variable for simplicity, esl is the odds ratios relating the frequency of respiratory infection to xerophthalmia. On the other hand, e,;, is the odds ratio for respiratory infection and xerophthalmia controlling for the respiratory infection status at the previous visit. If there exists a positive correlation between Yr- and Y, as we expect, and there is a positive effect of xerophthalmia on the risk of respiratory effect, it can be shown that"

IP?l 6 IB11. STATISTICAL INFERENCES FOR MARGINAL MODELS Note that only the first two moments of the joint distribution of Yi are specified by assumptions (1H3).Likelihood inference requires additional assumptions as suggested for binary data by ' ~ likelihood is often intractable because it involves many Prentice and Zhao" or Liang et ~ 1 . The nuisance parameters. For this reason, we suggest use of the Generalized Estimating Equations (GEE) method, a multivariate analogue of quasi-likelihood, which we now briefly outline. In the absence of a likelihood function p can be estimated by solving a multivariate analogue of the quasi-score function7

where pi = (pil, . . . , pini)'.The idea is to choose jso that p ( j ) is close to Yi on average and to optimally weight each residual Yi - pi by the inverse of cov( Yi). The term dpi/d/3 transforms the equation from the data to the parameter space. In the multivariate case, there is the additional complication that S , depends on a as well as p since cov( Y i )= cov( Yi; B, a). The correlation parameters a may be estimated by simultaneously solving S , = 0 and

2 where wi = (rilli2, rilri3, . . . , rini-lrinr,&, . . . , rini)' rij = y i j - pij and qi = E(wi; b, a).13 The choice of the weight function Hi depends on the type of responses. For binary response, the last ni components of wi can be ignored since y: = yir. In this case, we can use

var(rilri2) Hi=(

0 var(rilri3)

0

.

*

var(rini- 1 r i n r )

which ensures that S, depends on E and a only. For counted data, we suggest the use of the n: x nr identity matrix for H i when solving S, = 0, here n: = [(ni/2) nil.

+

1834

S. L. ZEGER AND K.-Y. LIANG

One may argue that the suggested Hi are not optimal.“ Our experience, however, has been that the choice of H i has very little impact on inference for b. Furthermore, to use the most suitable weight, namely cov(wi), one needs to make further assumptions on the third and the fourth joint moments of Y i . See Liang et al.” for details in the case of binary responses. The solution (j?,Oi) of S , = 0 and S , = 0 is asymptotically normal” with variance consistently estimated by

evaluated at

(b,G), where 0

Bi=(

cov( Yi) 0

b

The proposed method of estimation enjoys three properties. First, is nearly efficient relaiive to the maximum likelihood estimates of /r in many practical situations.”*” Second, fl is consistent, as the number of subjects K -,00, even if the covariance structure of Yi is incorrectly specified. This robustness property is important because with small n i , it is difficult to check the reliability of the assumption on cov(Yi). Third, our experience has been that the assumption about cov( Yi) has very little impact on inference as long as the robust variance estimate of in equation (7) is used. In other words, slightly difTerent assumptions about cov( Yi) lead to similar estimates of fl and of their standard errors. When regression coefficients are the scientific focus as in the examples here, these comments indicate that one should invest the lion’s share of time modelling the mean structure rather than the covariance. We have focused on inference for marginal models. Transition models can also be fit using GEE. The reader is referred to Korn and Whittemore16 and Ware et al.” and references therein. Random effects GLMs are more difficult to fit because evaluation of the likelihood or even the first two moments requires numerical methods in most problems. One approach to avoid the numerical integrals is to approximate the integrands with simple expansions whose integrals have closed forms. This was the approach proposed from a Bayesian perspective by Stiratelli et al.’* and recently reviewed by Breslow and Clayton.’ These approximate techniques give effective estimates of the fixed effects but are somewhat biased for estimating random effects and the random effects variance matrix,’. l 9 especially when the variance is large. An alternative strategy based upon optimal estimating equations was proposed by Waclawiw and Liang.*’ A Bayesian approach for fitting random effects GLMs using Gibbs sampling has been proposed by Zeger and Karim.” This method is used in the crossover example in the next section. See Zeger and Karim’’ for a complete discussion of the methodology.

EXAMPLES Crossover trial

The data shown in Table I1 are from Jones and Kenward3 who reported results of a crossover trial on cerebrovascular deficiency in which a placebo (A) and an active drug (B) were compared. Only data from centre 2 are used here for illustration. Thirty-four patients received treatment A

1835

OVERVIEW O F METHODS FOR ANALYSIS OF LONGITUDINAL DATA

Table 11. Data from a 2 x 2 crossover trial on cerebrovascular deficiency3 Outcomes for (first and second) period (0, 1) (190) (1,1)

Group

(0,O)

AB BA

6 9

0 4

6 2

22 18

Period Total

1

2

34 33

28 20

22 22

A placebo B treatment

followed by treatment B; another 33 patients were treated in the reverse order. The response variable is either 0 or 1 indicating whether the patient has an abnormal or normal electrocardiogram reading. Thus, for example, 22 patients in the AB group responded normally on both treatments whereas 4 patients in the BA group resgonded abnormally to the first treatment (B) but normally to the second treatment (A).Also shown in Table I1 are the marginal distributions of responses for both periods for the AB and BA groups. Focusing on the data at period one, note that 82 per cent (28/34)of the patients receiving placebo (A) were normal as opposed to 61 per cent (20/33)of the patients receiving active drug (B). This gives (13x 28)/(20x 6)= 3.0 as an estimate of the treatment odds ratio with approximate 95 per cent confidence interval (097,9.1). The estimate of the treatment effect is more than one, indicating that treatment A is favoured. However, the confidence interval is wide including the null value 1.0.It may be natural then to combine the odds ratio estimate from period one with that for period two ( 12 x 22)/(22x 11)= 1.09for the sake of reducing the sampling error. Two problems arise. First, combining across periods ignores the possible treatment-by-period interaction or ‘carry-over effect’. Secondly, the two responses for one individual are likely to be correlated. In fact, both groups show strong degrees of within-subject dependence as the odds ratio defined in equation (4) is estimated as (22x 6)/(6x 0.5)= 44 and (18x 19)/(4x 2)= 20.3 for group AB and group BA, respectively. We now show that the GEE method can be used to estimate the treatment efficacy while alleviating both concerns. Conceptually, cross over trials can be viewed as longitudinal studies. In this example, K = 67 and ni = 2 for each i. The covanates are treatment assignment ( X , ) and period (X2),defined as X , = 1 if treatment B,0 if treatment A and X2 = 1 if period two and 0 if period one. We first fit a logistic regression model logit E( Y = 1) = Po

+ PIXl + j2X2 +

83xlx2

with a common odds ratio y for the association between repeated outcomes for individuals. Results from Table 111 indicate there is little evidence in favour of a treatment-by-period interaction = 1.01 & 0.98)while a strong and positive within-subject association is evident (f = e3‘54= 34.5).In other words, for those with an abnormal reading at the first period, the odds of an abnormal reading at the next period is about 35 times higher than those whose reading was normal at the first period. When dropping the interaction term from the model, we observe a statistically significant = - 058 f 0.23)assumed to be the same for both periods. Thus the overall treatment effect odds of normal electrocardiogram readings for patients with cerebrovascular deficiency is 44 per cent lower for the drug group than for the placebo group regardless of treatment order. We note that by incorrectly assuming there is no within-subject dependence (Model 3, Table 111),the treatment effect appears to be non-significant when using the model-based standard error

(b3

(s2

1836

S. L. ZEGER AND K.-Y. LIANG

Table 111. Results (estimate*,robust standard error,t model based standard error$) of fitting three marginal logistic regressions to the 2 x 2 crossover data in Table I using generalized estimating equations (GEE) to account for correlation between repeated responses for each person. Also given are results of fitting a random intercept model Marginal model 3

Random intercept model

1.24 0.30 0.32

- 1.22

- 4.99

- 1.11

- 0.58

- 0.56

0.57 0.57

0.23 0.23

0.23 0.38

-

0.29 0.23 0.23

- 0.27

-

Variable

1

Intercept

- 1.540*

2 -

0.450t 0.450$

Treatment ( X ) Period ( X , )

- 0.85

0.58 0.58 XI*X,

log y or var ( b i ) l i 2

-

0.23 0.37

- 3.58

2-11

2.77 2.03 3.32

1.02 0.98 0.97 3.55 0.82 0-87

0.30 0.34

-

3.32 3.56 0.83 0.87

4.92 -

1.84

(b2

= - 0.56 k 0.38). This mistake can be avoided if one uses the robust standard error from equation (7) (0.233) instead. fit the random intercept model in which logit We have also E ( Yi, I b i ) = j : * + /?:* X + /?:* X 2 + b**XI X 2 + bi and bi is assumed to be Gaussian with mean 0 and variance u 2 . Note in the last column of Table 111 that the random effects coefficients are roughly three times as large as those from the marginal models. The marginal coefficients describe the relative frequencies of normal electrocardiogram in the population of patients. The random intercept coefficients describe the change in odds for an individual patient assuming the random intercept model is valid. The substantial difference in the magnitudes of coefficients is due to the extreme heterogeneity in the propensity to normal electrocardiograms. See Zhao and Prentice” for a related analysis of these data.

Epileptic seizures We analyse a subset of data from a clinical trial of 59 epileptic patients previously discussed by .~ suffering from partial seizures were randomThall and Vai14 and Breslow and C l a y t ~ nPersons ized to receive either progabide ( K , = 31) or a placebo ( K 2 = 28). The number of seizures in a eight-week baseline period prior to treatment and in four contiguous two-week periods during treatment were recorded for each person. The average number of seizures in the baseline period was 31.6 and 30.8 for the treatment and placebo groups, respectively. Our analysis addresses two questions: ( I ) is change from baseline to follow-up in the average number of seizures greater in the treatment group than in the control group; and (2) are the group differences in seizure level associated with the baseline seizure rate? Table IV gives the mean seizure rate (seizures per two-week interval) stratified by treatment group, pre- or post-treatment interval (baseline versus intervals 1 4 ) and quartiles of the baseline

1837

OVERVIEW OF METHODS FOR ANALYSIS OF LONGITUDINAL DATA

Table IV. Mean seizure rate (seizures per two weeks) by treatment group, pre- or post-treatment and quartiles of baseline seizure count. The cross-product ratio is (progabide interval lM/progabide baseline)/( placebo interval lM/placebo baseline) Baseline quartile Treatment

Interval

< 12

13-22

2341

> 41

Cornbined

1.96

3.97

6.68

2.42

4.44

8.32

15.5 1 (12.56)* 19.80

7.96 (5.71)* 7.90

3.53 2.45

3.33 4.55

6.50 6.94

18.11 15.61

8.60 7.70

0.56

1.23

0.86

1.09 (0.55)*

0.90t (0.74)*

Progabide

Placebo

{

BaZne

Cross-product ratio?

t

Summary statistics when patient 207 is deleted 0 9 0 = (7.96/7.90)/(860/7.70)

seizure counts. The 25 per cent, 50 per cent and 75 per cent quartiles of the baseline counts in eight weeks are 12,22 and 41, respectively. Overall, there is very little reduction of seizure counts in the two week period for the treated group (7.96/7.90 = 101 per cent) while there is a small increase in the placebo group (8.60/7.70 = 112 per cent). The treatment effect, can be measured by the cross-product ratio which is reported in the last row of Table IV. Note it appears to vary by the baseline seizure counts. For example, the cross-product ratio for people with relatively low baseline counts is 0.56, considerably lower than one. Patient 207 appears to be influential; the overall cross-product ratio drops from 0.90 to 0.74 when the data from this patient was deleted. In fact, the direction of the treatment effect reverses for those people whose baseline seizure count is greater than 41. Patient 207, who received progabide at age 22, has an unusually high seizure count (151) during the baseline interval, and remains high in the subsequent intervals (102,65,72, 63). We now appeal to more formal inference. To address question one, estimating the overall treatment effect, we use the following model: log E( Y i j )= log mij + Do + Dl Tx

+ p2 Time + b3 Tx* Time,j = 0,1,2,3,4; i = 1, . . . , 59.

(9)

Here mij = 8 i f j = 0 and mij = 2 i f j = 1,2,3,4 to account for different observation periods for the baseline and subsequent intervals. Suppressing the double subscripts on the explanatory variables for notational convenience, Tx = 1 for progabide, 0 for placebo and Time = 1 if interval 1-4,O if baseline. The variable Tx is included in the model to allow for differential baseline seizure counts between the treated and placebo groups. The coefficient B2 is the log of the ratio of seizure rates during the treatment phase (visits 1 4 ) relative to baseline among patients receiving placebo. The coefficient of interest p3 represents the difference in the seizure rate ratios (log scale) between the progabide and placebo groups. A negative coefficient corresponds to a greater reduction of seizure counts for the progabide group. Results from Table V shows that, overall, there is very little difference between the treated and placebo groups in the reduction of seizure counts prior to and following randomization - 0.104 f 0.214). As a side remark, we note that a different conclusion would be drawn if the strong over-dispersion = 19.44) were ignored.

(p3:

(6

1838

S. L. ZEGER A N D K.-Y. LIANG

Table V. Log-linear regression results for equation (9).The model was fit using GEE with and without patient 207, who had a very high seizure rate Variable

Patient 207 deleted

Complete data

Intercept

1.35 0.16

Treatment (Tx)

0.027 0.22

Time

Tx Time

1.35 0.16 -011

0.19

0.1 1

0.1 1

012

0.16

- 0.10

-

0.21 Over-dispersion parameter

19.4

0.30 0.17

10.4

Table V1 Estimated difference in the logarithm of the rate ratio (standard error) between the progabide and group for each quartile of baseline seizure Baseline seizure level (seizures in 8 weeks)

Model coefficients

All data

Patient 207 deleted

83

- 0.60

- 0.78

< 12

Estimates

(0.32)

(0.35)

13-22

83

+ 8s

018 (038)

0.07 (0.37)

2 3 41

83

+ 89

0 01 (0.30)

- 0.10

0.25 (0.44

0.14 (0.44

> 41

83

+ 810

(0.30)

To address the second issue regarding the association of the treatment effect with the baseline count, we expand the model as follows

+ PITx + P2Time + 8, Tx *Time + f14 log(age) + fls Time * Base , + b,Time * Base2 + 8, Time * BaseJ + BBTx*Time * Base, + P9Tx * Time * Base, + PloTx *Time * Base,. Here Base, = 1 if 12 < Yo < 22,O otherwise; Base2 = 1 if 22 < Yo < 41,O otherwise; Base, log E ( Yij)

= log mij

= 1 if 41 < Y o , 0 otherwise. The reference group comprises those with less than 13 in the first eight weeks. The estimated treatment effects for each quartile of baseline seizure count are presented in Table VI. There is mild evidence that the drug is effective for those with the lowest seizure rates at baseline but little or no evidence that there is an advantage over all. REFERENCES

1. Bloomfield, P. and Watson, G . S. ‘The inefficiency of least squares’, Biornepika, 62, 121-128 (1975). 2. Sommer, A., Katz, J. and Tarwotjo, I. ‘Increased risk of respiratory infection and diarrhea in children with pre-existing vitamin A deficiency’. American Journal of Clinical Nutrition, 40, 1090-1095 (1984).

OVERVIEW O F METHODS FOR ANALYSIS O F LONGITUDINAL DATA

1839

3. Jones, B. and Kenward, M. Design and Analysis of Cross Over Trials, Chapman and Hall, London, 1989. 4. Thall, P. F. and Vail, S. C. ‘Some covariance models for longitudinal count data with overdispersion’, Biornetrics, 46, 657-671 (1990). 5. Breslow, N. E. and Clayton, D. G. ‘Approximate inference in generalized linear mixed models’, Technical report No. 106, Department of Biostatistics, University of Washington, 1991, Seattle, WA 99164, U.S.A. 6. McCullagh, P. and Nelder, J. A. Generalized Linear Models, 2nd edn, Chapman and Hall, London, 1989.

7. Wedderburn, R. W. M. ‘Quasilikelihood function, generalized linear models and the Gaussian method’, Biornetrika, 61, 4 3 9 4 7 (1974). 8. Zeger, S. L., Liang, K.-Y. and Albert, P. S. ‘Models for longitudinal data: a generalized estimating equation approach’, Biornetrics, 44, 1049-1060 (1988). 9. Neuhaus, J. M., Kalbfleisch, J. D. and Hauck, W. W. ‘A comparison of cluster-specific and populationaveraged approaches for analyzing correlated binary data’, International Statistical Review, 59, 25 -36 (1991). 10. Bishop, Y. M., Fienberg, S. E. and Holland, P. W. Discrete Multivariate Analysis: Theory and Practice, Charles Griffin and Company, 1975. 11. Prentice, R. L. and Zhao, L. P. ‘Estimating equations for parameters in means and covariances of multivariate discrete and continuous responses’, Biornetrics, 47, 825-840 (1991). 12. Liang, K.-Y., Zeger, S. L. and Qaqish, B. ‘Multivariate regression analyses for categorical data’ (with discussion), Journal O Jthe Royal Statistical Society, Series B, 54, 3 4 0 (1992). 13. Prentice, R. L. ‘Correlated binary regression with covariate specific to each binary observation’, Biornetrics, 44, 1022-1048 (1988). 14. Godambe, V. P. ‘An optimum property of regular maximum likelihood estimation’, Annals of Mathematical Statistics, 31, 909-917 (1960). 15. Liang, K.-Y. and Zeger, S. L. ‘Longitudinal data analysis using generalized linear models’, Biornetrika, 73, 13-22 (1986). 16. Korn, E. L. and Whittemore, A. S. ‘Methods for analyzing panel studies of acute health effects of air pollution’, Biornetrics, 35, 7 15-802 (1979). 17. Ware, J. H., Lipsitz, S. and Speizer, F. E. ‘Issues in the analysis of repeated categorical outcomes’, Statistics in Medicine, 7 , 95-107 (1988). 18. Stiratelli, R., Laird, N. and Ware, J. H. ‘Random effects models for serial observations with binary responses’, Biometrics, 40, 96 1-97 1 ( 1984). 19. Karim, M. R. Generalized Linear Models with Random Effects, Unpublished Ph.D. Thesis, Biostatistics Department, Johns Hopkins University, 1991, 615 N. Wolfe Street, Baltimore, M D 21205, U.S.A. 20. Waclawiw, M. A. and Liang, K.-Y. ‘Predication of random effects in the generalized linear model’, Technical Report No. 704, Department of Biostatistics, Johns Hopkins University, (1991) 615 N. Wolfe Street, Baltimore, M D 21205, U.S.A. 21. Zeger, S. L. and Karim, M. R. ‘Generalized linear models with random effects; a Gibbs sampling approach’, Journal of the American Statistical Association, 86, 79-86 (1991). 22. Zhao, L. P. and Prentice, R. L. ‘Correlated binary regression using a quadratic exponential model’, Biornetrika, 77, 642-648 (1990).