[R] random interactions in lme

Jacob Michaelson jjmichael at comcast.net
Tue Apr 26 19:56:53 CEST 2005


Thanks, Ignacio --

That was another thing I'd been wondering about (the DenDF in SAS vs.  
lme).  Your example will give me something to chew on as I continue to  
try and reconcile proc mixed and lme.

Thanks for the guidance.

Jake

On Apr 26, 2005, at 10:36 AM, Ignacio Colonna wrote:

> 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