[R-sig-ME] estimating contrasts and their standard errors from mixed effects models in lme4

Nicholas Lewin-Koh nikko at hailmail.net
Sun Nov 22 20:15:38 CET 2009


Hi Eva,
Be careful when estimating contrasts when you have interactions in the
model. When doing contrasts on
the fixed effects you need to account for the interactions when looking
at the main effects.
Since I invariably screw up the contrast matrices, I like to use both
the contrast package to get the contrast matrix,
and then multcomp to calculate the contrasts and adjust for multiple
comparisons. Here is 
an example

library(contrasts)
library(multcomp)
library(mlmRev)
library(lme4)
data(egsingle)
 
fm1 <- lmer(math~year*size+female+(1|childid)+(1|schoolid), egsingle) ##
The mixed model
fm2 <- lm(math~year*size+female, data=egsingle) ## Main effects to get
contrast matrix, there is no method for lme4 objects

cc<-contrast(fm2,a=list(year=c(.5,1.5,2.5),size=380, female="Male"),
b=list(year=c(.5,1.5,2.5),size=800, female="Male"))

summary(glht(fm1, linfct = cc$X)) ## Plug in the contrast matrix from
contrast, multcomp will use the variance covariance matrix

	 Simultaneous Tests for General Linear Hypotheses

Fit: lmer(formula = math ~ year * size + female + (1 | childid) + 
    (1 | schoolid), data = egsingle, verbose = 1)

Linear Hypotheses:
       Estimate Std. Error z value Pr(>|z|)  
1 == 0  0.12774    0.08018   1.593   0.1271  
2 == 0  0.15323    0.08065   1.900   0.0661 .
3 == 0  0.17871    0.08177   2.185   0.0341 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- single-step method)


Note,the degrees of freedom in these tests is still a contentious issue.
But I think this is what
you are looking for.

Nicholas


> Message: 1
> Date: Fri, 20 Nov 2009 08:05:24 -0500
> From: eva petkova <evap120 at gmail.com>
> Subject: Re: [R-sig-ME] estimating contrasts and their standard errors
> 	from	mixed effects models in lme4
> To: adik at ilovebacon.org
> Cc: r-sig-mixed-models at r-project.org
> Message-ID:
> 	<754b10c00911200505u378741c5y3d01c5a0c236fa23 at mail.gmail.com>
> Content-Type: text/plain
> 
> Thank you Adam,
> 
> i am actually interested in estimating various contrasts without having
> to
> re-parametrize the factors in order to estimate different contrasts.  in
> your example, the parametrization of the effect factor you are using will
> estimate the contrasts  hea vs.cri or ach vs.cri, but if you wanted to
> estimate hea vs. ach, you will have to reparametrize the effect factor
> and
> refit the model.  in your case it it just one more reparametrization and
> fitting of the model but with factors that have more levels and
> interactions
> between them it is quite an elaborate process.  so i wondered if there is
> a
> function that can use the fit from one model and take as an input a
> linear
> combination of the regression coefficients and estimate the standard
> error
> of teh linear combination -- something similar to the "estimate"
> statement
> in Proc Mixed in SAS does.
> 
> Thanks again
> 
> e
> 
> On Thu, Nov 19, 2009 at 9:18 PM, Adam D. I. Kramer
> <adik at ilovebacon.org>wrote:
> 
> > Hi Eva,
> >
> >        If your data frame contains the factor variables you are interested
> > in analyzing, use the contrasts() function to get and set which contrasts
> > you would like.
> >
> >        Then, just use the factor variable in your fixed-effects formula.
> > lmer will automatically use the contrasts you provided, and you will get a
> > effect line for each contrast (including estimate and standard error).
> >
> > So, for example, I just completed this analysis and set it to my advisor.
> > We're looking at the perception of personal risk of an adverse outcome
> > occurring. There were three adverse outcomes, which I had contrast-coded in
> > the following manner:
> >
> >  contrasts(h3a2$effect)
> >>
> >    heaVcri achVcri
> > cri       0       0
> > hea       1       0
> > ach       0       1
> >
> > ...so when I use the h3a2$effect variable using lmer:
> >
> >  summary( lmer(risk ~ effect*order + income + (1|subj), data=h3a2) )
> >>
> > ... I get these fixed effects:
> >
> > Fixed effects:
> >                    Estimate Std. Error t value
> > (Intercept)          0.92254    0.21431   4.305
> > effectheaVcri       -0.53837    0.24859  -2.166
> > effectachVcri       -0.40662    0.25000  -1.627
> > order               -0.31194    0.07823  -3.987
> > income              -0.08551    0.03618  -2.363
> > effectheaVcri:order  0.27001    0.11736   2.301
> > effectachVcri:order  0.19991    0.11764   1.699
> >
> > ...output lines 2 and 3 are the two contrasts (labeled with the name of the
> > variable and then the name of the contrast: effectheaVcri is really
> > "effect""heaVcri"), with their estimates, and standard errors.  Lines 6 and
> > 7 are the same contrasts' interactions with the "order" variable.
> >
> > Really, the order effect is what I'm interested in here--order predicted
> > less risk, meaning that people who rated a given disorder later on in the
> > experiment thought they were less at risk for the disorder (even after we
> > control for what the disorder is, so for me the contrasts are more or less
> > control variables).
> >
> > Hope this helps!
> >
> > --Adam
> >
> >
> > On Thu, 19 Nov 2009, eva petkova wrote:
> >
> >  This must have come up before, but i did not find it in the help archive.
> >>
> >> in a mixed effects model fitted with lmer, i have an interaction term
> >> between two factors, each with more than two levels and would like to
> >> estimate and test  various contrasts between different combinations of the
> >> factor levels. i need the point estimates and the standard errors for
> >> these
> >> many contrasts.  does anyone know if there is a function that calculates
> >> the
> >> standard errors of contrasts?
> >>
> >> thank you
> >>
> >>        [[alternative HTML version deleted]]
> >>
> >> _______________________________________________
> >> R-sig-mixed-models at r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >>
> >>
> 
> 
> -- 
> Eva,
> 
> 	[[alternative HTML version deleted]]
> 
> 
> 
> ------------------------------




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