Re: [R] fitting a model with the nlme package

From: Douglas Bates <bates_at_stat.wisc.edu>
Date: Fri 04 Aug 2006 - 04:54:33 EST

On 8/3/06, Frank Johannes <fjohannes@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,



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Fri Aug 04 05:01:54 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Fri 04 Aug 2006 - 06:18:23 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.