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

Doran, Harold HDoran at air.org
Tue Feb 17 18:12:04 CET 2009 

It is an attribute. You can grab it like this:

> example(lmer)

> qq <- VarCorr(fm1)
> str(qq)
List of 1
 $ Subject: num [1:2, 1:2] 612.1 9.6 9.6 35.1
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "(Intercept)" "Days"
  .. ..$ : chr [1:2] "(Intercept)" "Days"
  ..- attr(*, "stddev")= Named num [1:2] 24.74 5.92
  .. ..- attr(*, "names")= chr [1:2] "(Intercept)" "Days"
  ..- attr(*, "correlation")= num [1:2, 1:2] 1 0.0655 0.0655 1
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "(Intercept)" "Days"
  .. .. ..$ : chr [1:2] "(Intercept)" "Days"
 - attr(*, "sc")= Named num 25.6
  ..- attr(*, "names")= chr "sigmaREML"
> attr(qq, 'sc')
sigmaREML 
 25.59182  

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