From: Douglas Bates <bates_at_stat.wisc.edu>

Date: Fri 04 Aug 2006 - 04:54:33 EST

*> (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)

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

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))

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-07number 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.
*