GEE and Mixed Models for

longitudinal data
Kristin Sainani Ph.D.
http://www.stanford.edu/~kcobb
Stanford University
Department of Health Research and Policy

Limitations of
rANOVA/rMANOVA

They assume categorical predictors.

They do not handle time-dependent covariates
(predictors measured over time).
They assume everyone is measured at the same
time (time is categorical) and at equally spaced
time intervals.
You dont get parameter estimates (just p-values)
Missing data must be imputed.
They require restrictive assumptions about the
correlation structure.
2

Example with timedependent, continuous

6 patients with depression are given a drug that increases levels of a
happy
chemical in the brain. At baseline, all 6 patients have similar
predictor
levels of this happy chemical and scores >=14 on a depression scale.
Researchers measure depression score and brain-chemical levels at
three subsequent time points: at 2 months, 3 months, and 6 months
post-baseline.
Here are the data in broad form:
id

time1

time2

time3

20

18

15

22

24

14

time4

chem1

chem2

chem3

chem4

20

1000

1100

1200

1300

18

22

1000

1000

1005

950

10

24

10

1000

1999

800

1700

38

34

32

34

1000

1100

1150

1100

25

29

25

29

1000

1000

1050

1010

30

28

26

14

1000

1100

1109

1500

Turn the data to long

form
data long4;
set new4;
time=0; score=time1;
time=2; score=time2;
time=3; score=time3;
time=6; score=time4;
run;

chem=chem1;
chem=chem2;
chem=chem3;
chem=chem4;

output;
output;
output;
output;

Note that time is being treated as a

continuous variablehere measured in
months.
If patients were measured at different times,
this is easily incorporated too; e.g. time can
be 3.5 for subject As fourth measurement
4
and 9.12 for subject Bs fourth

Data in long
form:

id

time

score

chem

20

1000

18

1100

15

1200

20

1300

22

1000

24

1000

18

1005

22

950

14

1000

10

1999

24

800

10

1700

38

1000

34

1100

32

1150

34

1100

25

1000

29

1000

25

1050

29

1010

30

1000

28

1100

26

1109

Graphically, lets see whats going on:

First, by subject.

All 6 subjects at once:

Mean chemical levels compared with

mean depression scores:

How do you analyze these

data?
Using repeated-measures ANOVA?
The only way to force a rANOVA here is
data forcedanova;
set broad;
avgchem=(chem1+chem2+chem3+chem4)/4;
if avgchem<1100 then group="low";
if avgchem>1100 then group="high";
run;
proc glm data=forcedanova;
class group;
model time1-time4= group/ nouni;
repeated time /summary;
run; quit;

Gives no
significant
results!

14

How do you analyze these

data?
We need more complicated models!
Todays lecture:
Introduction to GEE for longitudinal
data.
Introduction to Mixed models for
longitudinal data.

15

But firstnave analysis

The data in long form could be naively thrown into an

ordinary least squares (OLS) linear regression
I.e., look for a linear correlation between chemical
levels and depression scores ignoring the correlation
between subjects. (the cheating way to get 4-times as
much data!)
Can also look for a linear correlation between
depression scores and time.
In SAS:

proc reg data=long;

model score=chem time;
run;

16

Graphically
Nave linear regression here looks for significant slopes (ignoring
correlation between individuals):
Y= 24.90889 - 0.557778*time.

Y=42.44831-0.01685*chem

N=24as if we have 24 independent

observations!

17

The model
The linear regression model:

Yi 0 chem (chemi ) time (timei ) Errori

18

Results
The fitted model:

Yi 42.46803 .01704(chemi ) .07466(timei )

Parameter

Standard

Variable

DF

Estimate

Error

t Value

Pr > |t|

Intercept

42.46803

6.06410

7.00

<.0001

chem

-0.01704

0.00550

-3.10

0.0054

time

0.07466

0.64946

0.11

0.9096

1-unit increase in chemical is

associated with a .0174 decrease in
depression score (1.7 points per 100
units chemical)
Each month is associated only with a .
07 increase in depression score, after
correcting for chemical changes.

19

Generalized Estimating
Equations (GEE)

GEE takes into account the

dependency of observations by
specifying a working correlation
structure.
Lets briefly look at the model
(well return to it in detail later)

20

The model
Score1

Score 2

0
1
Score3

Score
4

Chem1
Chem2
2 (time) CORR Error
Chem3

Chem4

Measures linear correlation between chemical levels and depression

scores across all 4 time periods. Vectors!
Measures linear correlation between time and depression scores.
CORR represents the correction for correlation between observations.
A significant beta 1 (chem effect) here would mean either that people who
have high levels of chemical also have low depression scores (betweensubjects effect), or that people whose chemical levels change
correspondingly have changes in depression score (within-subjects effect),
21

SAS code (long form of

data!!)
Generalized Linear models (using
MLE)

proc genmod data=long4;

class id;
model score=chem time;
repeated subject = id / type=exch corrw;
run; quit;
Time is continuous (do not
place on class statement)!

The type of correlation

structure

Here we are modeling as a

linear relationship with score.

for time-dependent
predictors
NOTE,

--Interaction term with time (e.g. chem*time)

is NOT necessary to get a within-subjects
effect.
--Would only be included if you thought
22there

Results
Analysis Of GEE Parameter Estimates
Empirical Standard Error Estimates
Standard
Parameter Estimate

Error

95% Confidence
Limits

Z Pr > |Z|

Intercept

38.2431

4.9704

28.5013

47.9848

7.69

<.0001

chem

-0.0129

0.0026

-0.0180

-0.0079

-5.00

<.0001

time

-0.0775

0.2829

-0.6320

0.4770

-0.27

0.7841

In nave analysis,
standard error for
time parameter was:
0.64946 Its cut
by more than half
here.

In Nave analysis,
the standard error
for the chemical
coefficient was
0.00550 also cut
in half here.

23

Effects on standard
errors
In general, ignoring the dependency of the
observations will overestimate the standard errors
of the the time-dependent predictors (such as
time and chemical), since we havent accounted for
between-subject variability.
However, standard errors of the time-independent
predictors (such as treatment group) will be
underestimated. The long form of the data makes
it seem like theres 4 times as much data then there
really is (the cheating way to halve a standard error)!

24

What do the parameters

mean?

Time has a clear interpretation: .0775 decrease in score per

one-month of time (very small, NS).
Its much harder to interpret the coefficients from timedependent predictors:
Between-subjects interpretation (different types of people): Having a 100-unit
higher chemical level is correlated (on average) with having a 1.29 point lower
depression score.
Within-subjects interpretation (change over time): A 100-unit increase in chemical
levels within a person corresponds to an average 1.29 point decrease in depression
levels.
**Look at the data: here all subjects start at the same chemical level, but have
different depression scores. Plus, theres a strong within-person link between
increasing chemical levels and decreasing depression scores within patients (so
likely largely a within-person effect).

25

How does GEE work?

First, a naive linear regression analysis is carried

out, assuming the observations within subjects
are independent.
Then, residuals are calculated from the naive
model (observed-predicted) and a working
correlation matrix is estimated from these
residuals.
Then the regression coefficients are refit,
correcting for the correlation. (Iterative process)
The within-subject correlation structure is treated
as a nuisance variable (i.e. as a covariate)
26

OLS regression variancecovariance matrix

t1
t1

2
y/t
0

t2

t2

t3

y/t

y/t

Correlation structure
(pairwise correlations
between time points) is
Independence.

t3
Variance of scores is homogenous
across time (MSE in ordinary least
squares regression).

27

GEE variance-covariance
matrix
t1
t1

2
y/t
a

t2

t2

t3

y/t

y/t

Correlation structure must be

specified.

t3
Variance of scores is homogenous
across time (residual variance).

28

Choice of the correlation

structure within GEE
In GEE, the correction for within subject correlations is
carried out by assuming a priori a correlation structure
for the repeated measurements (although GEE is fairly
robust against a wrong choice of correlation matrix
particularly with large sample size)
Choices:

Independent (nave analysis)

Exchangeable (compound symmetry, as in rANOVA)
Autoregressive
M-dependent
Unstructured (no specification, as in rMANOVA)

We are looking for the simplest structure (uses up the fewest

degrees of freedom) that fits data well!

29

Independence
t1
t1

t2

t2

t3

0 0
0 0

0 0

t3
30

Exchangeable
t1

t1

t2
t3

t2

t3

Also known as compound

symmetry or sphericity. Costs 1 df
to estimate p.

31

Autoregressive
t1
t1

t2

t3

t2

t3

t4

3
2

t4
Only 1 parameter estimated.
Decreasing correlation for
farther time periods.

32

M-dependent
t1
t1

t2
2

t2

t3 0 2

t3

2
1

t4

0
2

t4
Here, 2-dependent. Estimate 2 parameters (adjacent
time periods have 1 correlation coefficient; time
periods 2 units of time away have a different
correlation coefficient; others are uncorrelated)

33

Unstructured
t1
t1
t2

1
2

t3 3

t2

1 2
5
5

4 6

t3

t4

3
4
6

t4
Estimate all correlations
separately (here 6)
34

How GEE handles missing

data
Uses the all available pairs method, in
which all non-missing pairs of data are
used in the estimating the working
correlation parameters.
Because the long form of the data are
being used, you only lose the
observations that the subject is
missing, not all measurements.

35

Back to our example

What does the empirical correlation matrix look like
for our data?
Pearson Correlation Coefficients, N = 6
Prob > |r| under H0: Rho=0

Independent?

time1

time2

time3

time4

time1

1.00000

0.92569
0.0081

0.69728
0.1236

0.68635
0.1321

time2

0.92569
0.0081

1.00000

0.55971
0.2481

0.77991
0.0673

time3

0.69728
0.1236

0.55971
0.2481

1.00000

0.37870
0.4591

time4

0.68635
0.1321

0.77991
0.0673

0.37870
0.4591

1.00000

Exchangeable?
Autoregressive?
M-dependent?
Unstructured?

36

Back to our example

I previously chose an exchangeable
correlation matrix
proc genmod data=long4;
class id;
model score=chem time;
repeated subject = id / type=exch corrw;
run; quit;
This asks to see
the working
37
correlation

Working Correlation Matrix

Working Correlation Matrix

Row1
Row2
Row3
Row4

Col1

Col2

Col3

Col4

1.0000
0.7276
0.7276
0.7276

0.7276
1.0000
0.7276
0.7276

0.7276
0.7276
1.0000
0.7276

0.7276
0.7276
0.7276
1.0000

Standard
Parameter Estimate

Error

95% Confidence
Limits

Z Pr > |Z|

Intercept

38.2431

4.9704

28.5013

47.9848

7.69

<.0001

chem

-0.0129

0.0026

-0.0180

-0.0079

-5.00

<.0001

time

-0.0775

0.2829

-0.6320

0.4770

-0.27

0.7841
38

Compare to
autoregressive
proc genmod data=long4;
class

id;

model score=chem time;

repeated subject = id / type=ar corrw;
run; quit;

39

Working Correlation Matrix

Working Correlation Matrix

Row1
Row2
Row3
Row4

Col1

Col2

Col3

Col4

1.0000
0.7831
0.6133
0.4803

0.7831
1.0000
0.7831
0.6133

0.6133
0.7831
1.0000
0.7831

0.4803
0.6133
0.7831
1.0000

Analysis Of GEE Parameter Estimates

Empirical Standard Error Estimates
Standard
Parameter Estimate
Error
Intercept
chem
time

36.5981
-0.0122
0.1371

4.0421
0.0015
0.3691

95% Confidence
Limits
28.6757
-0.0152
-0.5864

44.5206
-0.0092
0.8605

Z Pr > |Z|
9.05
-7.98
0.37

<.0001
<.0001
0.7104

40

Example tworecall
From rANOVA:
Within subjects effects,
but no between
subjects effects.
Time is significant.
Group*time is not
significant.
Group is not
significant.
This
is an example with
a binary timeindependent predictor.
41

Empirical Correlation
Pearson Correlation Coefficients, N = 6
Prob > |r| under H0: Rho=0
time1

time2

time3

time4

time1

1.00000

-0.13176
0.8035

-0.01435
0.9785

-0.50848
0.3030

time2

-0.13176
0.8035

1.00000

-0.02819
0.9577

-0.17480
0.7405

time3

-0.01435
0.9785

-0.02819
0.9577

1.00000

0.69419
0.1260

time4

-0.50848
0.3030

-0.17480
0.7405

0.69419
0.1260

1.00000

Independent?
Exchangeable?
Autoregressive?
M-dependent?
Unstructured?

42

GEE analysis
proc genmod data=long;
class group id;
model score=

group time group*time;

repeated subject = id / type=un corrw ;

run; quit;
NOTE,

for time-independent predictors

--You must include an interaction term with time to

get a within-subjects effect (development over time).

43

Working Correlation Matrix

Working Correlation Matrix
Col1
Col2
Row1
Row2
Row3
Row4

1.0000
-0.0701
0.1916
-0.1817

-0.0701
1.0000
0.1778
-0.5931

Col3

Col4

0.1916
0.1778
1.0000
0.5931

-0.1817
-0.5931
0.5931
1.0000

Group A is on average 8 points higher;

theres an average 5 point drop per
Analysis
Of GEE Parameter Estimates
time period for group B, and
an
Empirical
Standard Error Estimates
average 4.3 point drop more
for group
A.

Parameter
Intercept
group
A
group
B
time
time*group A

Standard
Estimate
Error
42.1433
7.8957
0.0000
-4.9184
-4.3198

6.2281
6.6850
0.0000
2.0931
2.1693

95% Confidence
Limits
29.9365
-5.2065
0.0000
-9.0209
-8.5716

54.3501
20.9980
0.0000
-0.8160
-0.0680

Comparable to within
effects for time and
time*group from
rMANOVA and rANOVA

Z Pr > |Z|
6.77
1.18
.
-2.35
-1.99

<.0001
0.2376
.
0.0188
0.0464

GEE analysis
proc genmod data=long;
class group id;
model score=

group time group*time;

repeated subject = id / type=exch corrw ;

run; quit;

45

Working Correlation Matrix

Working Correlation Matrix
Col1

Col2

Col3

Col4

Row1
1.0000
Row2
-0.0529
Row3
-0.0529
P-values areRow4
similar to rANOVA
-0.0529

-0.0529
1.0000
-0.0529
-0.0529

-0.0529
-0.0529
1.0000
-0.0529

-0.0529
-0.0529
-0.0529
1.0000

(which of course assumed

exchangeable, or compound Analysis Of GEE Parameter Estimates
symmetry, for the correlation
Empirical Standard Error Estimates
structure!)

Parameter
Intercept
group
A
group
B
time
time*group A

Standard
Estimate
Error
40.8333
7.1667
0.0000
-5.1667
-3.5000

5.8516
6.1974
0.0000
1.9461
2.2885

95% Confidence
Limits
29.3645
-4.9800
0.0000
-8.9810
-7.9853

52.3022
19.3133
0.0000
-1.3523
0.9853

Z Pr > |Z|
6.98
1.16
.
-2.65
-1.53

<.0001
0.2475
.
0.0079
0.1262

Introduction to Mixed Models

Return to our chemical/score example.

Ignore chemical for the moment, just ask if

theres a significant change over time in

47

Introduction to Mixed Models

Return to our chemical/score example.

48

Introduction to Mixed Models

Linear regression line for each person

49

Introduction to Mixed Models

Mixed models = fixed and random effects. For example,

Yit 0i ( random ) time ( fixed ) it

Residual
variance
:

Treated as a random variable with

a probability distribution.

0i ~ N ( 0 population ,

time constant

2
0

~ N (0, 2y / t )

This variance is comparable to the

between-subjects variance from
rANOVA.
Two parameters to estimate instead
of 1

50

Introduction to Mixed Models

What is a random effect?
--Rather than assuming there is a single intercept for the population, assume
that there is a distribution of intercepts. Every persons intercept is a
random variable from a shared normal distribution.
--A random intercept for depression score means that there is some average
depression score in the population, but there is variability between subjects.

0i ~ N ( 0 population ,

2
0

Generally, this is a
nuisance
parameterwe
have to estimate it for
making statistical
inferences, but we
dont care so much
about the actual
value.

51

Compare to OLS regression:

Compare with ordinary least squares regression (no
random effects):

Yit 0( fixed ) 1t ( fixed ) it

it ~ N (0, )
0 constant
time constant
2
y/t

Unexplained variability in Y.
LEAST SQUARES ESTIMATION
FINDS THE BETAS THAT MINIMIZE
THIS VARIANCE (ERROR)

52

RECALL, SIMPLE LINEAR REGRESSION:

The standard error of Y given T is the average variability around the
regression line at any given value of T. It is assumed to be equal at
all values of T.

y/t

y/t
y/t
y/t

y/t

y/t

All fixed effects

Yit 0( fixed ) 1t ( fixed ) it

it ~ N (0,

2
y/t

59.482929

0 constant

3 parameters to
estimate.

24.90888889

time constant
-0.55777778

54

The REG Procedure

Where to
find these
things in
OLS in SAS:

Model: MODEL1
Dependent Variable: score
Analysis of Variance
Sum of

Mean

DF

Squares

Square

F Value

Pr > F

Model

35.00056

35.00056

0.59

0.4512

Error

22

1308.62444

59.48293

Corrected Total

23

1343.62500

Source

Root MSE

7.71252

R-Square

0.0260

Dependent Mean

23.37500

Adj R-Sq

-0.0182

Coeff Var

32.99473
Parameter Estimates
Parameter

Standard

Variable

DF

Estimate

Error

t Value

Pr > |t|

Intercept

24.90889

2.54500

9.79

<.0001

time

-0.55778

0.72714

-0.77

0.4512

Introduction to Mixed Models

Adding back the random intercept term:

Yit 0i ( random ) 1t ( fixed ) it

0i ~ N ( 0 population ,

2
0

)
56

Meaning of random
intercept

Mean
population
intercept

Variation
in
intercepts

57

Introduction to Mixed Models

Yit 0i ( random ) 1t ( fixed ) it

it ~ N (0,

2
y/t

Residual variance:18.9264

0i ~ N ( 0 population ,
Same:24.90888889

time constant
Same:-0.55777778

2
0

4 parameters to
estimate.

Variability in intercepts
between subjects: 44.6121

58

Where to
find these
things in
from MIXED
in SAS:

Interpretation is the same

as with GEE: -.5578
decrease in score per
month time.
Effect

Covariance Parameter Estimates

Cov Parm

Subject

Variance

id

Estimate
44.6121

Residual

18.9264

Fit Statistics
-2 Res Log Likelihood

146.7

AIC (smaller is better)

152.7

AICC (smaller is better)

154.1

BIC (smaller is better)

152.1

44.6121
69%
18.9264 44.6121
69% of variability in
depression scores is
explained by the
differences between
subjects

Solution for Fixed Effects

Standard
Estimate

Error

DF

t Value

Pr > |t|

Intercept

24.9089

3.0816

8.08

0.0005

time

-0.5578

0.4102

17

-1.36

0.1916

Time coefficient is the same but standard error is nearly halved

(from 0.72714)..

With random effect for time,

but fixed intercept
Allowing time-slopes to be random:

Yit 0( fixed ) i ,time ( random ) it

i ,time ~ N ( time , population , )
2
t

60

Meaning of random beta

for time

61

With random effect for time,

but fixed intercept

Yit 0( fixed ) i ,time ( random ) it

it ~ N (0, y2 / t )

Residual variance:40.4937

i ,time ~ N ( time , population , )

2
t

Same: 24.90888889

0 constant

Variability in time slopes

between subjects: 1.7052

Same:-0.55777778

62

With both random

With a random intercept and random time-slope:

Yit 0i ( random ) i ,time ( random ) it

0i ~ N ( 0 population ,

2
0

i ,time ~ N ( time, population , )

2
t

63

Meaning of random beta

for time and random
intercept

64

With both random

With a random intercept and random time-slope:

Yit 0i ( random ) i ,time ( random ) it

16.6311

0i ~ N ( 0 population ,
24.90888889

2
0

53.0068

i ,time ~ N ( time, population , )

2
t

0.55777778

0.4162

Additionally, we have to
estimate the covariance of
the random intercept and
random slope:
here -1.9943
(adding random time therefore
cost us 2 degrees of freedom)
65

Choosing the best model

Aikake Information Criterion (AIC) : a fit statistic
penalized by the number of parameters

AIC = - 2*log likelihood + 2*(#parameters)

Values closer to zero indicate better fit

and greater parsimony.

Choose the model with the smallest AIC.

66

AICs for the four models

MODEL
AIC
All fixed
162.2
Intercept random 150.7
Time slope fixed
Intercept fixed
161.4
Time effect
random
All random

152.7
67

In SASto get model with

random intercept
proc mixed data=long;
class id;
model score = time /s;
random int/subject=id;
run; quit;

68

Model with chem (timedependent variable!)

proc mixed data=long;
class id;
model score = time chem/s;
random int/subject=id;
run; quit;
Typically, we take care of the repeated
measures problem by adding a random
intercept, and we stop therethough you can
try random effects for predictors and time.
69

Residual and
AIC are
reduced even
further due to
strong
explanatory
power of
chemical.

Interpretation is the same

as with GEE: we cannot
separate betweensubjects and withinsubjects effects of
chemical.
Effect

Cov Parm

Subject

Intercept

id

Residual

Estimate
35.5720
10.2504

Fit Statistics
-2 Res Log Likelihood

143.7

AIC (smaller is better)

147.7

AICC (smaller is better)

148.4

BIC (smaller is better)

147.3

Solution for Fixed Effects

Standard
Estimate

Error

DF

t Value

Pr > |t|

38.1287

4.1727

9.14

0.0003

time

-0.08163

0.3234

16

-0.25

0.8039

chem

-0.01283

0.003125

16

-4.11

0.0008

Intercept

New Example: timeindependent binary

From GEE:
predictor
Strong effect of time.
No group difference
Non-significant
group*time trend.

71

SAS code
proc mixed data=long ;
class id group;
model score = time group
time*group/s corrb;
random int /subject=id ;
run; quit;

72

Results (random intercept)

Fit Statistics
-2 Res Log Likelihood

138.4

AIC (smaller is better)

142.4

AICC (smaller is better)

143.1

BIC (smaller is better)

142.0

Solution for Fixed Effects

Standard
Effect

group

Estimate

Error

DF

t Value

Pr > |t|

Intercept

40.8333

4.1934

9.74

0.0006

time

-5.1667

1.5250

16

-3.39

0.0038

1.21

0.2444

group

7.1667

5.9303

16

group

time*group

-3.5000

2.1567

16

-1.62

time*group

.
0.1242
.

73

Compare to GEE results

Analysis Of GEE Parameter Estimates
Empirical Standard Error Estimates
Parameter
Intercept
group
A
group
B
time
time*group A

Standard
Estimate
Error
40.8333
7.1667
0.0000
-5.1667
-3.5000

5.8516
6.1974
0.0000
1.9461
2.2885

Same coefficient estimates.

Nearly identical p-values.

95% Confidence
Limits
29.3645
-4.9800
0.0000
-8.9810
-7.9853

52.3022
19.3133
0.0000
-1.3523
0.9853

Z Pr > |Z|
6.98
1.16
.
-2.65
-1.53

<.0001
0.2475
.
0.0079
0.1262

Mixed model with a random intercept is

equivalent to GEE with exchangeable
correlation(slightly different std. errors in
SAS because PROC MIXED additionally
allows Residual variance to change over
time.

Power of these models

Since these methods are based on generalized linear models,
these methods can easily be extended to repeated measures with a
dependent variable that is binary, categorical, or counts
These methods are not just for repeated measures. They are
appropriate for any situation where dependencies arise in the
data. For example,
Studies across families (dependency within families)
Prevention trials where randomization is by school, practice, clinic, geographical area, etc.
(dependency within unit of randomization)
Matched case-control studies (dependency within matched pair)
In general, anywhere you have clusters of observations (statisticians say that observations
are nested within these clusters.)
For repeated measures, our cluster was the subject.
75
In the long form of the data, you have a variable that identifies which cluster the observation

References

Jos W. R. Twisk. Applied Longitudinal Data Analysis for Epidemiology: A

Practical Guide. Cambridge University Press, 2003.

76