interaction model and variance in lme

example(lmer)

qq <- VarCorr(fm1)
str(qq)

attr(qq, 'sc')

-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org 
[mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf 
Of David Springate
Sent: Tuesday, February 17, 2009 12:05 PM
To: 'Douglas Bates'
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] interaction model and variance in lme

Thanks for the help!

The ML model that seems to be appropriate for my work is

model = lmer(trait~treatment + (1|family) +
(1|family:treatment),data=d,REML=F)

as I am most concerned with determining the interaction term 
of family and treatment.

However, I still have a problem.  I need to get to the 
variance components of the random effects to do some further 
downstream analysis.
The part of the summary output I need is:

Random effects:
 Groups           Name        Variance Std.Dev.
 family:treatment (Intercept)  5.42852 2.32992 
 family           (Intercept)  0.79754 0.89305 
 Residual                     17.16665 4.14327 
Number of obs: 172, groups: family:treatment, 49; family, 25

I have used VarCorr to get the variances for family and 
family:treatment:
a = VarCorr(model)
family.var = a$family[1,1]

but VarCorr doesn't have the residual variance term.  

How can I get to this variable?  Do I need to put an 
additional error term in the model?

Thanks in advance!


-----Original Message-----
From: dmbates at gmail.com [mailto:dmbates at gmail.com] On Behalf 
Of Douglas Bates
Sent: 12 February 2009 14:03
To: David Springate
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] interaction model and variance in lme

On Wed, Feb 11, 2009 at 4:19 PM, David Springate 
<David.Springate at postgrad.manchester.ac.uk> wrote:

Hi,

I am new to R and I have been trying to build a model that I can 
extract

ML

variance components for the interaction from, but am struggling to 
make sense of the formula.

I am looking at the interaction of family (random) and two 
environmental treatments (fixed) on a trait.  So far I have tried:

model=lme(trait~treatment,data=dataset,random=~1|family/treatment)

But this seems to give me just the treatment nested in the family 
rather than family*treatment values.

Some of the difficulty here may be in nomenclature and 
notation.  The phrase "family*treatment values" means 
different things in SAS and in R.  It may help if you could 
describe verbally what you want to obtain rather than symbolically.

For a model like this I would recommend using lmer from the 
lme4 package as it has several enhancements relative to lme 
from the nlme package.

As you have seen, the model above, which would be written

lmer(trait ~ treatment + (1|family/treatment), dataset)

or, equivalently,

lmer(trait ~ treatment + (1|family) + (1|family:treatment), dataset)

provides fixed effects for the treatment (Intercept and an 
effect of one level) plus the random effects for each family 
plus the random effects for each family:treatment 
combination.  (In R an interaction term is written 
family:treatment; the notation family*treatment indicates 
crossing of fixed effects so that family*treatment expands to 
family + treatment + family:treatment.)

Another model incorporating random effects for each level of 
treatment within family is

lmer(trait ~ treatment + (treatment|family), dataset)

for which I prefer and alternative expression as

lmer(trait ~ treatment + (0+treatment|family), dataset)

These models are equivalent but have different parameterizations.
Suppose that the treatment levels are called A and B.  Then 
the default model matrix for the treatment term provides the 
intercept and the indicator of level B.  The model matrix for 
0+treatment suppresses the intercept and provides an 
indicator for A and an indicator for B.

The notation (0+treatment|family) produces a pair of random 
effects for each level of family, one for treatment A and one 
for treatment B, and the variance-covariance matrix for these 
random effects is a general positive-definite 2x2 matrix.  
(In SAS-speak this is an "unconstrained" variance-covariance 
matrix but the mathematician in me will not accept the 
concept of an unconstrained matrix that is subject to the 
constraints of being symmetric and positive definite.)

I hope this helps but I encourage you to follow up on your 
question if this did not answer it.

I also tried:

Model = lme(trait~family*treatment,data=dataset,random=~1|family)

which gives me an interaction MS and F value in an anova, 

but seems to
treat

each interaction between each treatment and family as a 

separate fixed 

factor.  Also VarCorr() doesn't seem to give a variance for the

interaction

term.

What am I doing wrong?

I am sure this should be simple, but the docs seem pretty 

unclear to 

me on modeling interactions.

Help please!

David Springate