Tell me more ×
Cross Validated is a question and answer site for statisticians, data analysts, data miners and data visualization experts. It's 100% free, no registration required.

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 ...
> head(dat)
  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 ?

share|improve this question
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

Know someone who can answer? Share a link to this question via email, Google+, Twitter, or Facebook.

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.