[R-sig-ME] Repeated measures ANOVA using lmer, nlme, and geepack
array chip
arrayprofile at yahoo.com
Wed Aug 10 20:15:15 CEST 2011
I am struggling to compare between R and SAS as well (see my other post), so I tried a number of things and was able to get close between the 2.
But with your questions below, first of all, anova() in R report Type I SS while in SAS output Type III SS (and Type II SS in R is hugely controversial especially with interaction term in the model. Prof Ripley has an article about it which you can search to find out). So the 2 anova is not comparable. However, the F value and P value for the last term in the model (interaction) should be comparable if model is specified equivalently. Another confounding factor, as far as I can understand, is how to determine the denominator degree of free in F test is controversial in mixed models, that's why lmer() doesn't even report a P value for fixed effect.
With my dataset (see my other post), I found the following R code is equivalent to the following SAS code for getting the same fixed effect coefficients, but not exactly the same p values for these coefficients because of the denominator DF issue.
summary(lme(y ~ age + time * trt, random=~1|pid/biopsy.site, data = test,correlation=corAR1(form= ~ time | pid/biopsy.site)))
proc mixed;
class pid trt time biopsysite;
model y = age trt time trt*time /s;
random biopsysite / type=cs subject=pid;
repeated / type=ar(1) subject=biopsysite(pid);
run;
With your example, can you try
summary(lme(ADASCOG ~ TREAT * MONTH, random = ~ 1 | PAT, correlation=corSymm(form = ~ 1 | PAT), data=unialz, na.action=na.omit))
vs
PROC MIXED DATA = UNIALZ;
CLASS TREAT MONTH PAT;
MODEL ADASCOG = TREAT MONTH TREAT*MONTH / s;
REPEATED / TYPE=UN SUBJECT=PAT RCORR;
RUN;
see if they gave the same coefficients? can similarly try corAR1() as well
John
----- Original Message -----
From: Paul Miller <pjmiller_57 at yahoo.com>
To: r-sig-mixed-models at r-project.org
Cc:
Sent: Wednesday, August 10, 2011 7:23 AM
Subject: [R-sig-ME] Repeated measures ANOVA using lmer, nlme, and geepack
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
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list