# SAS fails to fit a mixed model

Here are some simulated data:

    library(mvtnorm)
I <- 3 # positions (fixed factor)
J <- 4 # tubes (random factor)
K <- 4 # repeats
n <- I*J*K
set.seed(123)
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)) )
tube <- rep(tube, each=K)
position <- rep(position, each=K)
Mu_ij <- rep(Mu_ij, each=K)
dat <- data.frame(tube, position, Mu_ij)
sigmaw <- 2
dat$y <- rnorm(n, dat$Mu_ij, sigmaw)
dat$tube <- factor(dat$tube)

> str(dat)
'data.frame':   48 obs. of  4 variables:
$tube : Factor w/ 4 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 2 2 2 2 3 3 ...
$Mu_ij : num 2.44 2.44 2.44 2.44 6.13 ...$ y       : num  3.24 2.66 1.33 6.01 7.12 ...
tube position    Mu_ij        y
1    1        A 2.439524 3.241067
2    1        A 2.439524 2.660890
3    1        A 2.439524 1.327842
4    1        A 2.439524 6.013351
5    1        B 6.129288 7.124989
6    1        B 6.129288 2.196053


I fit a mixed model with R, it works well:

> library(lme4)
> lmer(y ~ position + (0+position|tube), data=dat)
Linear mixed model fit by REML
Formula: y ~ position + (0 + position | tube)
Data: dat
AIC   BIC logLik deviance REMLdev
212.6 231.3  -96.3    194.8   192.6
Random effects:
Groups   Name      Variance Std.Dev. Corr
tube     positionA 0.30123  0.54885
positionB 0.68317  0.82654  -0.695
positionC 1.66666  1.29099  -0.408  0.940
Residual           3.14003  1.77201
Number of obs: 48, groups: tube, 4

Fixed effects:
Estimate Std. Error t value
(Intercept)   3.3533     0.5211   6.435
positionB     3.1098     0.8923   3.485
positionC     5.6138     1.0144   5.534

Correlation of Fixed Effects:
(Intr) postnB
positionB -0.753
positionC -0.651  0.744


But the same model does not work well with SAS:

PROC MIXED DATA=dat ;
CLASS POSITION TUBE ;
MODEL y = POSITION ;
RANDOM POSITION / subject=TUBE type=UN G GCORR ;
RUN; QUIT;


gives

 Estimated G matrix is not positive definite.
Estimated G Matrix

Row    Effect      position    tube        Col1        Col2        Col3

1    position    A           1        0.08895     -0.5823     -0.1545
2    position    B           1        -0.5823      0.1455      1.2431
3    position    C           1        -0.1545      1.2431      1.4835


Is it possible to remedy this failure ?

-
I'm not sure about the SAS part, but I ran this model in R and got a warning message: –  smillig Apr 13 '12 at 12:57
What does the warning say ? –  Stéphane Laurent Apr 13 '12 at 13:04
Oops, my message got cut off. My warning message says: Warning message: In mer_finalize(ans) : singular convergence (7) which usually happens when there are variance components estimated to be 0. I imagine that my R warning is related to the SAS warning. You might want to graph the data. –  smillig Apr 13 '12 at 13:09
That's strange: I don't have this warning. If you type head(dat), you get the same as me ? –  Stéphane Laurent Apr 13 '12 at 14:03
Yes, it's the same. Could it be the version of lme4 and/or R? I'm using R version 2.14.2 and lme4_0.999375-42 on a Mac (OS X 10.6.8). –  smillig Apr 13 '12 at 16:05
show 2 more comments