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

Douglas Bates bates at stat.wisc.edu
Thu Feb 12 15:03:04 CET 2009

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

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

More information about the R-sig-mixed-models mailing list