On 8/3/06, Frank Johannes <

[hidden email]> 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,

______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.