[R] fitting a model with the nlme package

[Douglas Bates]

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,