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 ?
Rand got a warning message: – smillig Apr 13 '12 at 12:57Warning message: In mer_finalize(ans) : singular convergence (7)which usually happens when there are variance components estimated to be 0. I imagine that myRwarning is related to the SAS warning. You might want to graph the data. – smillig Apr 13 '12 at 13:09lme4and/orR? I'm usingR version 2.14.2andlme4_0.999375-42on a Mac (OS X 10.6.8). – smillig Apr 13 '12 at 16:05