random interactions in lme

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.