On 8/3/06, Frank Johannes <fjohannes at fastmail.fm>
wrote:> Dear all,
> I am analyzing some data that requires a mixed model. I have been
> reading Pinheiro and Bates' book,
> but cannot find the notation to fit the following model:
>
> Suppose I have the dataset below. Here I am fitting variable p as a
> fixed effect, variable h
> as a random effect and variable t as nested within h. I would like to
> include the variable j as well
> as an independent (non-nested) random effect. I can't find the notation
> in the book to tell R
> to do this. Does anybody know?
>
> library(nlme)
>
> h<-c(1,1,1,2,2,2,3,3,3)
> j<-c(2,3,8,3,4,3,9,5,4)
> t<-c(1,2,3,1,2,3,1,2,3)
> f<-c(1,2,1,2,1,2,1,2,1)
> p<-c(45,32,35,12,23,12,2,9,12)
> set<-data.frame(p,h,j,t)
>
> out<-lme(p~ -1 + f,data=set, random=~1|h/t)
It is not easy to fit mixed models with non-nested grouping factors
for the random effects using lme. I would recommend using lmer from
the lme4/Matrix package instead. Even with lmer you can't fit the
model you want to fit because you have 9 levels of the interaction h:t
and only 9 observations. The version of lmer in the
soon-to-be-released Matrix_0.995-13 even produces a warning about
this.
> library(lme4)
Loading required package: Matrix
Loading required package: lattice
Loading required package: lattice> set <- data.frame(h = factor(c(1,1,1,2,2,2,3,3,3)),
+ j = factor(c(2,3,8,3,4,3,9,5,4)),
+ t = factor(c(1,2,3,1,2,3,1,2,3)),
+ f = c(1,2,1,2,1,2,1,2,1),
+ p = c(45,32,35,12,23,12,2,9,12))> (out <- lmer(p ~ f - 1 + (1|h/t) + (1|j), set))
Linear mixed-effects model fit by REML
Formula: p ~ f - 1 + (1 | h/t) + (1 | j)
Data: set
AIC BIC logLik MLdeviance REMLdeviance
39.21916 40.00806 -15.60958 36.17505 31.21916
Random effects:
Groups Name Variance Std.Dev.
t:h (Intercept) 1.7053e-22 1.3059e-11
j (Intercept) 1.0067e+02 1.0033e+01
h (Intercept) 1.2655e+02 1.1249e+01
Residual 3.4106e-13 5.8400e-07
number of obs: 9, groups: t:h, 9; j, 6; h, 3
Fixed effects:
Estimate Std. Error t value
f 12.000 4.899 2.4495
Warning messages:
1: Estimated variance for group(s) t:h is zero in:
"LMEoptimize<-"(`*tmp*`, value = list(maxIter = 200, tolerance
1.49011611938477e-08,
2: nlminb returned message singular convergence (7)
in: "LMEoptimize<-"(`*tmp*`, value = list(maxIter = 200, tolerance
1.49011611938477e-08,