[R-sig-ME] Repeated measures ANOVA using lmer, nlme, and geepack
Paul Miller
pjmiller_57 at yahoo.com
Wed Aug 10 16:23:30 CEST 2011
Hello All,
I sent the message below about two weeks ago and got no reply. I'm told that many people on the list were away at JSM in Miami around that time. My hope is that some people who otherwise would have replied were away and that my message simply fell through the cracks. So I'm going to try posting it again.
Hope no one will mind.
Thanks,
Paul
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