[R] random interactions in lme
Ignacio Colonna
iacolonn at uiuc.edu
Tue Apr 26 18:36:39 CEST 2005
The code below gives almost identical results for a split-block analysis in
lme and SAS proc mixed, in terms of variance components and F statistics. It
just extends the example in Pinheiro & Bates (p.162) to a split block
design.
I am including below the SAS code and the data in case you want to try it.
The only difference between both is in the df for the F denominator, which I
wasn't able to compute correctly in lme, but this may be my ignorance on how
to correctly specify the model. It is not a big issue though, as the F
values are identical, so you can compute the p-values if you know how to
obtain the correct DenDF.
# a split block design
spbl.an1<-lme(yield~rowspace*ordered(tpop),random=list(rep=pdBlocked(list(pd
Ident(~1),
pdIdent(~rowspace-1),pdIdent(~ordered(tpop)-1)))),data=spblock)
* SAS code
proc mixed data=splitblock method=reml;
class rep rowspace tpop;
model yield=rowspace tpop rowspace*tpop;
random rep rep*rowspace rep*tpop;
run;
# data
rowspace tpop rep plot yield
9 60 1 133 19
9 120 1 101 19.5
9 180 1 117 22
9 240 1 132 19.4
9 300 1 116 23.9
18 60 1 134 15.8
18 120 1 102 26.2
18 180 1 118 21.9
18 240 1 131 20
18 300 1 115 23.3
9 60 2 216 20.6
9 120 2 233 22
9 180 2 201 23.4
9 240 2 217 28.2
9 300 2 232 25.9
18 60 2 215 19.7
18 120 2 234 30.3
18 180 2 202 22.4
18 240 2 218 27.9
18 300 2 231 28.5
9 60 3 309 20.8
9 120 3 308 21.6
9 180 3 324 24.6
9 240 3 340 25.3
9 300 3 325 35.3
18 60 3 310 17.2
18 120 3 307 23.6
18 180 3 323 24.9
18 240 3 339 30.7
18 300 3 326 33
9 60 4 435 15.6
9 120 4 403 20.4
9 180 4 430 24.4
9 240 4 414 21
9 300 4 419 23.2
18 60 4 436 17.7
18 120 4 404 23.6
18 180 4 429 21.7
18 240 4 413 24.4
18 300 4 420 26.2
Ignacio
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Douglas Bates
Sent: Monday, April 25, 2005 6:40 PM
To: Jacob Michaelson
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] random interactions in lme
Jacob Michaelson wrote:
>
> On Apr 24, 2005, at 8:52 AM, Douglas Bates wrote:
>
>> Jacob Michaelson wrote:
>>
>>> Hi All,
>>> I'm taking an Experimental Design course this semester, and have
>>> spent many long hours trying to coax the professor's SAS examples
>>> into something that will work in R (I'd prefer that the things I
>>> learn not be tied to a license). It's been a long semester in that
>>> regard.
>>> One thing that has really frustrated me is that lme has an extremely
>>> counterintuitive way for specifying random terms. I can usually
>>> figure out how to express a single random term, but if there are
>>> multiple terms or random interactions, the documentation available
>>> just doesn't hold up.
>>> Here's an example: a split block (strip plot) design evaluated in SAS
>>> with PROC MIXED (an excerpt of the model and random statements):
>>> model DryMatter = Compacting|Variety / outp = residuals ddfm =
>>> satterthwaite;
>>> random Rep Rep*Compacting Rep*Variety;
>>> Now the fixed part of that model is easy enough in lme:
>>> "DryMatter~Compacting*Variety"
>>> But I can't find anything that adequately explains how to simply add
>>> the random terms to the model, ie "rep + rep:compacting +
>>> rep:variety"; anything to do with random terms in lme seems to go off
>>> about grouping factors, which just isn't intuitive for me.
>>
>>
>> The grouping factor is rep because the random effects are associated
>> with the levels of rep.
>>
>> I don't always understand the SAS notation so you may need to help me
>> out here. Do you expect to get a single variance component estimate
>> for Rep*Compacting and a single variance component for Rep*Variety?
>> If so, you would specify the model in lmer by first creating factors
>> for the interaction of Rep and Compacting and the interaction of Rep
>> and Variety.
>>
>> dat$RepC <- with(dat, Rep:Compacting)[drop=TRUE]
>> dat$RepV <- with(dat, Rep:Variety)[drop=TRUE]
>> fm <- lmer(DryMatter ~ Compacting*Variety+(1|Rep)+(1|RepC)+(1|RepV), dat)
>>
>>
>>
>
> Thanks for the prompt reply. I tried what you suggested, here's what I
> got:
>
> > turf.lme=lmer(dry_matter~compacting*variety+(1|rep)+(1|rc)+(1|rv),
> turf.data)
> Error in lmer(dry_matter ~ compacting * variety + (1 | rep) + (1 | rc) +
:
> entry 3 in matrix[9,2] has row 3 and column 2
>
> Just to see what the problem was, I deleted the third random term, and
> it didn't complain:
>
> > turf.lme=lmer(dry_matter~compacting*variety+(1|rep)+(1|rv), turf.data)
> > anova(turf.lme)
> Analysis of Variance Table
> Df Sum Sq Mean Sq Denom F value Pr(>F)
> compacting 5 10.925 2.185 36.000 18.166 5.68e-09 ***
> variety 2 2.518 1.259 36.000 10.468 0.0002610 ***
> compacting:variety 10 6.042 0.604 36.000 5.023 0.0001461 ***
>
> Now obviously this isn't a valid result since I need that third term;
> but interestingly, it didn't matter which term I deleted, so long as
> there were only two random terms. Any ideas as to what the error
> message is referring to?
>
> Thanks for the help,
>
> Jake Michaelson
Unfortunately, yes I do know what the error message is referring to - a
condition that should not happen. This is what Bill Venables would call
an "infelicity" in the code and others with less tact than Bill might
call a bug.
______________________________________________
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
More information about the R-help
mailing list