[R-sig-ME] Repeated measures ANOVA using lmer, nlme, and geepack
Paul Miller
pjmiller_57 at yahoo.com
Thu Jul 28 18:01:36 CEST 2011
Hello Everyone,
I'm A SAS user who is in the process of learning R. I've read about a dozen R-related books and have been testing my understanding of how to perform various analyses by replicating the results from one of my SAS books. Most of the time, I can replicate the results from SAS using R. And when I can't, it's ususally clear that there are differences in the way SAS and R handle the analysis.
I've replicated almost all of the analyses in the of 22 chapters in my book, but am currently stuck on a couple of things in the repeated measures ANOVA chapter. Specifically, I'm having trouble with repeated measures ANOVA using mixed models and GEE.
Below is some R code and output along with corresponding code and output from SAS. I tried using both lmer and nlme and found I was able to come closest to replicating the mixed model results using the latter. The results in R and SAS look very different and this makes me suspect that I'm using R incorrectly.
Can anyone help me understand what I'm doing wrong?
Thanks,
Paul
###############################################################
#### Example 8.3: Disease Progression in Alzheimer's Trial ####
###############################################################
#### Data in Wide Format ####
connection <- textConnection("
TREAT PAT ADAS02 ADAS04 ADAS06 ADAS08 ADAS10 ADAS12
L 1 22 30 . 33 28 30
L 5 34 35 46 37 31 35
L 8 40 41 41 46 52 48
L 12 24 . 21 28 30 27
L 13 29 26 29 26 . 36
L 15 31 36 41 46 52 57
L 19 22 27 28 24 27 28
L 21 43 49 42 48 48 46
L 24 18 28 29 . 25 28
L 28 25 24 27 18 21 22
L 31 37 35 35 38 42 .
L 34 24 27 28 24 27 25
L 37 45 50 58 59 60 58
L 40 33 32 35 30 31 35
L 44 34 37 43 44 39 38
L 47 25 27 29 28 31 .
L 51 30 . 36 32 34 38
L 54 23 . 33 28 32 32
L 57 35 37 39 38 41 43
L 59 44 48 48 45 50 52
L 63 28 30 32 31 35 32
L 67 24 22 23 24 27 30
L 68 . 49 51 48 55 54
L 73 26 28 30 27 30 33
L 76 30 32 35 35 36 38
L 78 40 42 44 43 45 46
H 2 31 36 35 31 31 31
H 6 24 27 28 21 27 26
H 7 31 31 39 37 41 .
H 14 45 48 46 52 48 42
H 17 24 28 26 23 24 29
H 20 21 32 39 36 33 30
H 22 32 34 45 42 37 32
H 25 18 22 26 26 27 24
H 27 51 47 . 43 43 43
H 33 20 22 29 24 29 30
H 38 41 34 37 29 35 33
H 42 24 35 39 32 24 .
H 45 23 . 33 36 33 30
H 50 25 28 25 28 28 30
H 52 31 34 . 33 34 35
H 56 27 31 26 33 33 34
H 60 37 43 39 42 43 36
H 62 41 42 51 45 46 51
H 66 35 33 34 35 36 41
H 69 30 31 27 34 33 36
H 72 54 60 55 58 . 65
H 75 35 37 39 41 39 44
H 79 18 21 19 19 20 27
H 80 40 35 33 39 38 41
P 3 31 36 37 41 39 44
P 4 20 26 32 35 25 29
P 9 33 33 29 33 39 41
P 10 35 39 40 38 40 38
P 11 26 24 31 42 50 .
P 16 44 48 44 37 36 47
P 18 25 31 21 27 41 32
P 23 28 34 26 26 36 35
P 26 27 . 28 35 40 .
P 29 20 30 30 27 33 29
P 30 49 . 43 48 44 53
P 32 26 29 31 30 35 38
P 35 30 33 41 . 41 44
P 36 31 34 44 44 50 56
P 39 42 46 36 43 48 48
P 41 31 30 31 . 41 38
P 43 27 22 36 45 54 60
P 46 24 37 41 31 36 44
P 48 33 31 38 41 31 .
P 49 27 30 36 36 32 33
P 53 35 34 45 44 38 40
P 55 39 40 38 44 43 44
P 58 32 34 40 45 36 38
P 61 45 50 . 54 50 53
P 64 21 23 31 34 27 27
P 65 26 30 37 37 30 32
P 70 53 50 55 57 . .
P 71 32 34 27 30 36 35
P 74 . 50 52 56 52 54
P 77 24 32 31 37 35 30
")
alzhmrs <- read.table(connection, header=TRUE, sep="", na.strings=".")
alzhmrs$TREAT <- factor(alzhmrs$TREAT, labels=c("High", "Low", "Placebo"))
head(alzhmrs)
#### Reshape Data to Long Format ####
library(reshape)
unialz <- melt(alzhmrs, id=c("TREAT","PAT"))
names(unialz) <- c("TREAT", "PAT", "MONTH", "ADASCOG")
unialz$TREAT <- relevel(unialz$TREAT,"Placebo")
unialz$MONTH <- relevel(unialz$MONTH,"ADAS12")
unialz <- sort_df(unialz,vars=c("PAT","MONTH"))
head(unialz)
#############################################
#### Mixed Model Approach (Unstructured) ####
#############################################
#### R code and output ####
library("nlme")
model.un <- lme(ADASCOG ~ TREAT * MONTH, random = ~ TREAT | PAT,
correlation=corSymm(form = ~ 1 | PAT),
data=unialz, na.action=na.omit)
summary(model.un)
anova(model.un)
> anova(model.un)
numDF denDF F-value p-value
(Intercept) 1 359 1560.5256 <.0001
TREAT 2 77 0.8724 0.4220
MONTH 5 359 24.6385 <.0001
TREAT:MONTH 10 359 1.9795 0.0346
#### SAS code and output ####
ODS SELECT Mixed.RCorr Mixed.Dimensions Mixed.FitStatistics Mixed.Tests3;
PROC MIXED DATA = UNIALZ;
CLASS TREAT MONTH PAT;
MODEL ADASCOG = TREAT MONTH TREAT*MONTH;
REPEATED / TYPE=UN SUBJECT=PAT(TREAT) RCORR;
TITLE4 'PROC MIXED Using Unstructured Covariance (UN)';
RUN;
PROC MIXED Using Unstructured Covariance (UN)
The Mixed Procedure
Type 3 Tests of Fixed Effects
Num Den
Effect DF DF F Value Pr > F
TREAT 2 77 0.94 0.3967
MONTH 5 77 21.55 <.0001
TREAT*MONTH 10 77 2.08 0.0357
model.cs <- update(model.un, correlation=corCompSymm(form = ~ 1 | PAT))
summary(model.cs)
anova(model.cs)
####################################
#### Mixed Model Approach (AR1) ####
####################################
#### R code and output ####
model.ar1 <- update(model.un, correlation=corAR1(form = ~ 1 | PAT))
summary(model.ar1)
anova(model.ar1)
> anova(model.ar1)
numDF denDF F-value p-value
(Intercept) 1 359 1597.7569 <.0001
TREAT 2 77 1.1192 0.3318
MONTH 5 359 30.9931 <.0001
TREAT:MONTH 10 359 2.3784 0.0097
#### SAS code and output ####
ODS SELECT Mixed.RCorr Mixed.Dimensions Mixed.FitStatistics Mixed.Tests3;
PROC MIXED DATA = UNIALZ;
CLASS TREAT MONTH PAT;
MODEL ADASCOG = TREAT MONTH TREAT*MONTH;
REPEATED / TYPE=AR(1) SUBJECT=PAT(TREAT);
TITLE4 'PROC MIXED Using 1st-Order Auto-Regressive Covariance (AR(1))';
RUN;
Type 3 Tests of Fixed Effects
Num Den
Effect DF DF F Value Pr > F
TREAT 2 77 0.94 0.3952
MONTH 5 359 11.21 <.0001
TREAT*MONTH 10 359 1.43 0.1671
######################
#### GEE Approach ####
######################
#### R code and output ####
library("geepack")
model <- geeglm(ADASCOG ~ TREAT * MONTH, id = PAT,
data=unialz, family = gaussian, corstr = "ar1")
summary(model)
anova(model)
> anova(model)
Analysis of 'Wald statistic' Table
Model: gaussian, link: identity
Response: ADASCOG
Terms added sequentially (first to last)
Df X2 P(>|Chi|)
TREAT 2 3.7 0.15595
MONTH 5 95.0 < 2e-16 ***
TREAT:MONTH 10 30.7 0.00065 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#### SAS code and output ####
ODS SELECT Genmod.Type3;
PROC GENMOD DATA = UNIALZ;
CLASS TREAT MONTH PAT;
MODEL ADASCOG = TREAT MONTH TREAT*MONTH / DIST=NORMAL TYPE3;
REPEATED SUBJECT=PAT / TYPE=AR(1);
TITLE4 'GEE Analysis Using PROC GENMOD';
TITLE5 'Autoregressive Correlation (AR(1)) Working Correlation';
RUN;
GEE Analysis Using PROC GENMOD
Autoregressive Correlation (AR(1)) Working Correlation
The GENMOD Procedure
Score Statistics For Type 3 GEE Analysis
Chi-
Source DF Square Pr > ChiSq
TREAT 2 2.02 0.3639
MONTH 5 46.14 <.0001
TREAT*MONTH 10 21.76 0.0164
More information about the R-sig-mixed-models
mailing list