# Equivalent mixed models yielding different results in SAS

Below I show two equivalent ways to write a mixed model in R and SAS. The two R models as well as the two SAS models yield the same estimates of the random and fixed effects and the same standard errors of the estimates of the fixed effects. But the two SAS models do not give the same confidence intervals of the fixed effects. My question is: which are the correct confidence intervals ?

Here are simulated data:

library(mvtnorm)
I <- 3
J <- 6
K <- 5
n <- I*J*K
set.seed(666)
tube <- rep(1:J, each=I)
position <- rep(LETTERS[1:I], times=J)
Mu_i <- 3*(1:I)
Mu_ij <- c(t(rmvnorm(J, mean=Mu_i, sigma=diag(I)+2)) )
tube <- rep(tube, each=K)
position <- rep(position, each=K)
Mu_ij <- rep(Mu_ij, each=K)
dat <- data.frame(tube, position)
sigmaw <- 2
dat$y <- rnorm(n, Mu_ij, sigmaw) dat$tube <- factor(dat$tube) > str(dat) 'data.frame': 90 obs. of 3 variables:$ tube    : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
$position: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 2 2 2 2 2 ...$ y       : num  2.76 5.5 2.54 1.56 6.46 ...
tube position        y
1    1        A 2.759443
2    1        A 5.496689
3    1        A 2.540150
4    1        A 1.558261
5    1        A 6.462050
6    1        B 4.239749


The corresponding nlme model is the following:

> # firstly set position C as the "intercept" for concordance with SAS
> dat$position <- relevel(dat$position, "C")
> library(nlme)
> # fit the model
> ( fit1 <- lme(y ~ position, data=dat, random= list(tube = pdCompSymm(~ 0+position ))) )
Linear mixed-effects model fit by REML
Data: dat
Log-restricted-likelihood: -199.0196
Fixed: y ~ position
(Intercept)   positionA   positionB
8.526544   -4.800535   -3.322507

Random effects:
Formula: ~0 + position | tube
Structure: Compound Symmetry
StdDev   Corr
positionC 1.892433
positionA 1.892433 0.082
positionB 1.892433 0.082 0.082
Residual  1.932750

Number of Observations: 90
Number of Groups: 6


This model is equivalent to a 2-way ANOVA with mixed effets (in the sense that the marginal models are the same), which is more easy to fit with lme4:

> library(lme4)
> lmer(y ~ position + (1|tube) + (1|tube:position), data=dat)
Linear mixed model fit by REML
Formula: y ~ position + (1 | tube) + (1 | tube:position)
Data: dat
AIC BIC logLik deviance REMLdev
410 425   -199    402.5     398
Random effects:
Groups        Name        Variance Std.Dev.
tube:position (Intercept) 3.28587  1.81270
tube          (Intercept) 0.29543  0.54354
Residual                  3.73552  1.93275
Number of obs: 90, groups: tube:position, 18; tube, 6

Fixed effects:
Estimate Std. Error t value
(Intercept)   8.5265     0.8493  10.039
positionA    -4.8005     1.1595  -4.140
positionB    -3.3225     1.1595  -2.866

Correlation of Fixed Effects:
(Intr) postnA
positionA -0.683
positionB -0.683  0.500


Below I check that the results are indeed the same:

> # same standard erros
> sqrt(diag(fit1\$varFix))
(Intercept)   positionA   positionB
0.8493533   1.1594505   1.1594505
> # the total variance in the second model is given in the first model:
> sqrt(1.81270^2+ 0.54354^2)
[1] 1.892437


Well. Now here are the two equivalent SAS models:

/* First model */
PROC MIXED DATA=dat ;
CLASS POSITION TUBE ;
MODEL y = POSITION ;
RANDOM POSITION / subject=TUBE type=CS ;
RUN; QUIT;

/* Second model */
PROC MIXED DATA=dat ;
CLASS POSITION TUBE ;
MODEL y = POSITION ;
RANDOM TUBE TUBE*POSITION ;
RUN; QUIT;


Results are identical to the R results. But SAS assigns different degrees of freedom for the fixed effects and consequently gives different confidence intervals, as shown below.

The first model gives degrees of freedom 5, 10, 10:

Effect  position    Estimate    StandardError   DF  Alpha   Lower   Upper
Intercept       8.5265  0.8494  5   0.05    6.3432  10.7099
position    A   -4.8005 1.1595  10  0.05    -7.384  -2.2171
position    B   -3.3225 1.1595  10  0.05    -5.9059 -0.7391
position    C   0   .   .   .   .   .


whereas the second model gives degrees of freedom 15, 15, 15:

Effect  position    Estimate    StandardError   DF  Alpha   Lower   Upper
Intercept       8.5265  0.8494  15  0.05    6.7162  10.3369
position    A   -4.8005 1.1595  15  0.05    -7.2718 -2.3292
position    B   -3.3225 1.1595  15  0.05    -5.7938 -0.8512
position    C   0   .   .   .   .   .

• There are different ways to compute degrees of freedom in a mixed model. You can change the method SAS uses with the ddfm option after the model statement. Maybe using a different method produces the same results as R?
– Glen
Apr 16 '12 at 19:55
• Yes but: 1) This is a balanced 2-way ANOVA with mixed effects, so I believe there is a natural choice of the degrees of freedom ; 2) The first method gives equal degrees of freedom (15, 15, 15), not the second one (5,10,10), that sounds strange, doesn't it ? And finally R does not provide such confidence intervals (nlme provides "narrow" confidence intervals in which I'm not interested) Apr 17 '12 at 4:58
• Are you sure about the "narrow" confidence intervals? I've not heard about that. Would you share a reference or link? Thanks. Jun 1 '12 at 4:27
• Maybe I'm doing a confusion. But I am almost sure about one fact: the nlme intervals are always narrower than the SAS intervals. Jun 1 '12 at 5:35
• @StéphaneLaurent: I haven't noticed that, but I don't think I've ever looked at the CIs between the two. I am pretty sure that the standard errors are the same between SAS and R, though, so would be surprised that the CIs would be different (except for differences due to in the calculation of degrees of freedom). Jun 3 '12 at 18:43

By default, SAS uses a very simple method to compute the degrees of freedom; it looks to see if the fixed effect appears in a term in the random effects; if so, it uses that random effect as the degrees of freedom.

In this case, the model with random effects of TUBE and TUBE*POSITION do the right thing as this method can tell that POSITION should use TUBE*POSITION. But the model with POSITION / subject=TUBE cannot. Instead, one should tell it to compute the degrees of freedom using another method, usually the Satterthwaite (satterth) for models with only random effects and the Kenward-Roger method (kr) for models with repeated effects structures.

PROC MIXED DATA=dat ;
CLASS POSITION TUBE ;
MODEL y = POSITION /s ddfm=satterth;
RANDOM POSITION / subject=TUBE type=CS;
RUN;

• Indeed, since the day I asked my question, I have found that the Kenward & Roger degrees of freedom are correct for this model (I have not tried Satterthwaite) May 31 '12 at 14:13
• I give this link where you also provide useful related information : stats.stackexchange.com/questions/29551/… May 31 '12 at 19:33
• Hi Stéphane, I was wondering if you could elaborate on why you say that the Kenward-Roger degrees of freedom is the correct one.
– jkd
May 31 '12 at 20:03
• @jkd: This is a balanced design, so the df should match the df in the traditional least squares approach, which the K-R method does. May 31 '12 at 21:03
• Yes. Exactly as Aaron say in the link above. Jun 1 '12 at 5:33