Repeated measures ANOVA using lmer, nlme, and geepack

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;