[R-sig-ME] interaction model and variance in lme

[David Springate]

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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>