R help - Jan 2005 - Specification of factorial random-effects model

I want to specify two factors and their interaction as random effects using
the function lme().   This works okay when I specify these terms using the
function Error() within the function aov(), but I can't get the same model
fitted using lme().   The code below illustrates the problem.

 

a <- factor(rep(c(1:3), each = 27))

b <- factor(rep(rep(c(1:3), each = 9), times = 3))

c <- factor(rep(rep(c(1:3), each = 3), times = 9))

y <- c(74.59,75.63,76.7,63.48,63.17,65.99,64,66.35,64.5,

   46.57,44.16,47.96,35.09,36.14,35.16,36.4,34.72,34.58,

   41.82,47.35,45.74,33.33,36.8,33.38,34.13,34.39,34.48,

   89.73,85.24,90.86,82.5,79.44,81.65,77.74,77.02,81.62,

   59.32,62.29,60.7,55.42,55.5,51.17,50.54,53.54,51.85,

   64.5,63.6,65.19,55.07,50.26,53.73,54.57,47.8,48.8,91.56,

   94.49,92.17,82.14,83.16,81.31,83.58,78.63,77.08,60.53,

   60.79,58.57,51.28,52.9,51.54,49.15,48.97,51.61,59.44,

   60.07,60.07,51.94,52.2,50.2,49.45,50.75,49.56)

anovamodel <- aov(y ~ c + Error(a*b))

summary(anovamodel)

lmemodel <- lme(y ~ c, random = ~ 1|a*b)

anova(lmemodel)

 

N.B. You have to load the package 'nlme' before lme() will run at all.
When it does run, it gives the error

 

Error in getGroups.data.frame(dataMix, groups) : 

        Invalid formula for groups

 

I think there is something I haven't understood about the syntax of the
argument 'random' in lme() .   Can anyone explain it to me?

 

Nick Galwey


	[[alternative HTML version deleted]]

If you read the Help file for lme (!), you'll see that ~1|a*b is certainly
incorrect.

Briefly, the issue has been discussed before on this list: the current
version of lme() follows the original Laird/Ware formulation for **nested**
random effects. Specifying **crossed** random effects is possible but
difficult, and the fitting algorithm is not optimized for this. See p. 163
in Bates and Pinheiro for an example.

-- Bert Gunter
Genentech Non-Clinical Statistics
South San Francisco, CA
 
"The business of the statistician is to catalyze the scientific learning
process."  - George E. P. Box
 
 
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Nicholas Galwey
> Sent: Wednesday, January 26, 2005 1:45 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] Specification of factorial random-effects model
> 
> I want to specify two factors and their interaction as random 
> effects using
> the function lme().   This works okay when I specify these 
> terms using the
> function Error() within the function aov(), but I can't get 
> the same model
> fitted using lme().   The code below illustrates the problem.
> 
>  
> 
> a <- factor(rep(c(1:3), each = 27))
> 
> b <- factor(rep(rep(c(1:3), each = 9), times = 3))
> 
> c <- factor(rep(rep(c(1:3), each = 3), times = 9))
> 
> y <- c(74.59,75.63,76.7,63.48,63.17,65.99,64,66.35,64.5,
> 
>    46.57,44.16,47.96,35.09,36.14,35.16,36.4,34.72,34.58,
> 
>    41.82,47.35,45.74,33.33,36.8,33.38,34.13,34.39,34.48,
> 
>    89.73,85.24,90.86,82.5,79.44,81.65,77.74,77.02,81.62,
> 
>    59.32,62.29,60.7,55.42,55.5,51.17,50.54,53.54,51.85,
> 
>    64.5,63.6,65.19,55.07,50.26,53.73,54.57,47.8,48.8,91.56,
> 
>    94.49,92.17,82.14,83.16,81.31,83.58,78.63,77.08,60.53,
> 
>    60.79,58.57,51.28,52.9,51.54,49.15,48.97,51.61,59.44,
> 
>    60.07,60.07,51.94,52.2,50.2,49.45,50.75,49.56)
> 
> anovamodel <- aov(y ~ c + Error(a*b))
> 
> summary(anovamodel)
> 
> lmemodel <- lme(y ~ c, random = ~ 1|a*b)
> 
> anova(lmemodel)
> 
>  
> 
> N.B. You have to load the package 'nlme' before lme() will run at
all.
> When it does run, it gives the error
> 
>  
> 
> Error in getGroups.data.frame(dataMix, groups) : 
> 
>         Invalid formula for groups
> 
>  
> 
> I think there is something I haven't understood about the 
> syntax of the
> argument 'random' in lme() .   Can anyone explain it to me?
> 
>  
> 
> Nick Galwey
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at 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
>