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. Any help? Thanks in advance! --Jake Michaelson
random interactions in lme
6 messages · Douglas Bates, Jake Michaelson, Ignacio Colonna
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)
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
1 day later
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.
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
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