Convert mixed model with repeated measures from SAS to R

Question

I have been trying to convert a repeated measures model from SAS to R, since a collaborator will do the analysis but does not have SAS. We are dealing with 4 groups, 8 to 10 animals per group, and then 5 time points for each animal. The mock data file is available here https://drive.google.com/file/d/0B-WfycVUQyhaVGU2MUpuQkg4Mk0/edit?usp=sharing as a Rdata file and here https://drive.google.com/file/d/0B-WfycVUQyhaR0JtZ0V4VjRkTk0/edit?usp=sharing as an excel file:

The original SAS code (1) is :

proc mixed data=essai.data_test method=reml;
    class group time mice;
    model param = group time group*time / ddfm=kr;
    repeated time / type=un subject=mice group=group;
run;

Which gives :

    Type 3 Tests des effets fixes
               DDL     DDL     Valeur
Effet         Num.    Res.          F    Pr > F
group            3    15.8       1.58    0.2344
time             4    25.2      10.11    <.0001
group*time      12    13.6       1.66    0.1852

I know that R does not handle degrees of freedom in the same way as SAS does, so I am first trying to obtain results similar to (2) :

proc mixed data=essai.data_test method=reml;
    class group time mice;
    model param = group time group*time;
    repeated time / type=un subject=mice group=group;
run;

I have found some hints here Converting Repeated Measures mixed model formula from SAS to R and when specifying a compound symmetry correlation matrix this works perfectly. However, I am not able to obtain the same thing for a general correlation matrix.

With (2) in SAS, I obtain the following results :

   Type 3 Tests des effets fixes
              DDL     DDL     Valeur
Effet         Num.    Res.          F    Pr > F
group            3      32       1.71    0.1852
time             4     128      11.21    <.0001
group*time      12     128       2.73    0.0026

Using the following R code :

options(contrasts=c('contr.sum','contr.poly'))
mod <- lme(param~group*time, random=list(mice=pdDiag(form=~group-1)),
            correlation = corSymm(form=~1|mice),
            weights = varIdent(form=~1|group),
            na.action = na.exclude, data = data, method = "REML")
anova(mod,type="marginal")

I obtain:

            numDF denDF   F-value p-value
(Intercept)     1   128 1373.8471  <.0001
group           3    32    1.5571  0.2189
time            4   128   10.0628  <.0001
group:time     12   128    1.6416  0.0880

The degrees of freedom are similar, but not the tests on fixed effects and I don’t know where this comes from. Would anyone have any idea of what I am doing wrong here?

Answer 1

Your R code differs from the SAS code in multiple ways. Some of them are fixable, but I was not able to fix all the aspects to reproduce the SAS analysis.

1. The R code fits a mixed effects model with a random mice effect, while the SAS code fits a generalized linear model that allows correlation between the residuals, but there are no random effects (because there is no RANDOM statement). In R you would have to use the gls function from the same nlme package.

2. In the R code all observations within the same group have the same variance, while in the SAS code you have an unstructured covariance matrix, that is each time-point within each group has its own variance. You can achieve the same effect by using weights=varIdent(form=~1|group*time).

3. In the R code the correlation matrix is the same for every mouse regardless of group. In the SAS code each group has its own correlation matrix. This is the part that I don't know how to reproduce in R.

I have to note that the R model seems to be more meaningful - SAS estimates way too many variances and correlations (which, by the way, you can see meaningfully arranged using the R and RCORR options to the repeated statement).