Note : this question is a repost, as my previous question had to be deleted for legal reasons.
While comparing PROC MIXED from SAS with the function
lme from the
nlme package in R, I stumbled upon some rather confusing differences. More specifically, the degrees of freedom in the different tests differ between
PROC MIXED and
lme, and I wondered why.
Start from the following dataset (R code given below) :
- ind : factor indicating the individual where the measurement is taken
- fac : organ where measurement is taken
- trt : factor indicating the treatment
- y : some continuous response variable
The idea is to build the following simple models :
y ~ trt + (ind) :
ind as a random factor
y ~ trt + (fac(ind)) :
fac nested in
ind as a random factor
Note that the last model should cause singularities, as there's only 1 value of
y for every combination of
In SAS, I build the following model :
PROC MIXED data=Data; CLASS ind fac trt; MODEL y = trt /s; RANDOM ind /s; run;
According to tutorials, the same model in R using
nlme should be :
> require(nlme) > options(contrasts=c(factor="contr.SAS",ordered="contr.poly")) > m2<-lme(y~trt,random=~1|ind,data=Data)
Both models give the same estimates for the coefficients and their SE, but when carrying out an F test for the effect of
trt, they use a different amount of degrees of freedom :
SAS : Type 3 Tests of Fixed Effects Effect Num DF Den DF F Value Pr > F trt 1 8 0.89 0.3724 R : > anova(m2) numDF denDF F-value p-value (Intercept) 1 8 70.96836 <.0001 trt 1 6 0.89272 0.3812
Question1: What is the difference between both tests? Both are fitted using REML, and use the same contrasts.
NOTE: I tried different values for the DDFM= option (including BETWITHIN, which theoretically should give the same results as lme)
In SAS :
PROC MIXED data=Data; CLASS ind fac trt; MODEL y = trt /s; RANDOM fac(ind) /s; run;
The equivalent model in R should be :
In this case, there are some very odd differences :
- R fits without complaining, whereas SAS notes that the final hessian is not positive definite (which doesn't surprise me a bit, see above)
- The SE on the coefficients differ (is smaller in SAS)
- Again, the F test used a different amount of DF (in fact, in SAS that amount = 0)
SAS output :
Effect trt Estimate Std Error DF t Value Pr > |t| Intercept 0.8863 0.1192 14 7.43 <.0001 trt Cont -0.1788 0.1686 0 -1.06 .
R Output :
> summary(m4) ... Fixed effects: y ~ trt Value Std.Error DF t-value p-value (Intercept) 0.88625 0.1337743 8 6.624963 0.0002 trtCont -0.17875 0.1891855 6 -0.944840 0.3812 ...
(Note that in this case, the F and T test are equivalent and use the same DF.)
Interestingly, when using
lme4 in R the model doesn't even fit :
> require(lme4) > m4r <- lmer(y~trt+(1|ind/fac),data=Data) Error in function (fr, FL, start, REML, verbose) : Number of levels of a grouping factor for the random effects must be less than the number of observations
Question 2: What is the difference between these models with nested factors? Are they specified correctly and if so, how comes the results are so different?
Simulated Data in R :
Data <- structure(list(y = c(1.05, 0.86, 1.02, 1.14, 0.68, 1.05, 0.22, 1.07, 0.46, 0.65, 0.41, 0.82, 0.6, 0.49, 0.68, 1.55), ind = structure(c(1L, 2L, 3L, 1L, 3L, 4L, 4L, 2L, 5L, 6L, 7L, 8L, 6L, 5L, 7L, 8L), .Label = c("1", "2", "3", "4", "5", "6", "7", "8"), class = "factor"), fac = structure(c(1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L), .Label = c("l", "r"), class = "factor"), trt = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Cont", "Treat"), class = "factor")), .Names = c("y", "ind", "fac", "trt" ), row.names = c(NA, -16L), class = "data.frame")
Simulated Data :
y ind fac trt 1.05 1 l Treat 0.86 2 l Treat 1.02 3 l Treat 1.14 1 r Treat 0.68 3 r Treat 1.05 4 l Treat 0.22 4 r Treat 1.07 2 r Treat 0.46 5 r Cont 0.65 6 l Cont 0.41 7 l Cont 0.82 8 l Cont 0.60 6 r Cont 0.49 5 l Cont 0.68 7 r Cont 1.55 8 r Cont