[R-sig-ME] Repeated measures ANOVA using lmer, nlme, and geepack

[Gabor Grothendieck]

On Thu, Jul 28, 2011 at 12:01 PM, Paul Miller <pjmiller_57 at yahoo.com> wrote:
> 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;

REPEATED in SAS refers to a structured error term, not to mixed
effects, so use gls (in nlme) instead of lme.  Also Type III is not
invariant to coding so try contr.sum or contr.SAS (and this also
implies you may need to omit the relevel lines above).

-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com