CODEWITHSUNDEEP
rlme4random-effectsmixed-models
Tomas
Tomas

Reputation: 59565

How to get coefficients and their confidence intervals in mixed effects models?

In lm and glm models, I use functions coef and confint to achieve the goal:

m = lm(resp ~ 0 + var1 + var1:var2) # var1 categorical, var2 continuous
coef(m)
confint(m)

Now I added random effect to the model - used mixed effects models using lmer function from lme4 package. But then, functions coef and confint do not work any more for me!

> mix1 = lmer(resp ~ 0 + var1 + var1:var2 + (1|var3)) 
                                      # var1, var3 categorical, var2 continuous
> coef(mix1)
Error in coef(mix1) : unable to align random and fixed effects
> confint(mix1)
Error: $ operator not defined for this S4 class

I tried to google and use docs but with no result. Please point me in the right direction.

EDIT: I was also thinking whether this question fits more to https://stats.stackexchange.com/ but I consider it more technical than statistical, so I concluded it fits best here (SO)... what do you think?

Upvotes: 45

Views: 100540

Answers (7)

Ben Bolker
Ben Bolker

Reputation: 226771

I'm going to add a bit here. If m is a fitted (g)lmer model (most of these work for lme too):

  • fixef(m) is the canonical way to extract coefficients from mixed models (this convention began with nlme and has carried over to lme4)
  • you can get the full coefficient table with coef(summary(m)); if you have loaded lmerTest before fitting the model, or convert the model after fitting (and then loading lmerTest) via coef(summary(as(m,"lmerModLmerTest"))), then the coefficient table will include p-values. (The coefficient table is a matrix; you can extract the columns via e.g. ctab[,"Estimate"], ctab[,"Pr(>|t|)"], or convert the matrix to a data frame and use $-indexing.)
  • As stated above you can get likelihood profile confidence intervals via confint(m); these may be computationally intensive. If you use confint(m, method="Wald") you'll get the standard +/- 1.96SE confidence intervals. (lme uses intervals(m) instead of confint().)

If you prefer to use broom.mixed:

  • tidy(m,effects="fixed") gives you a table with estimates, standard errors, etc.
  • tidy(as(m,"merModLmerTest"), effects="fixed") (or fitting with lmerTest in the first place) includes p-values
  • adding conf.int=TRUE gives (Wald) CIs
  • adding conf.method="profile" (along with conf.int=TRUE) gives likelihood profile CIs

You can also get confidence intervals by parametric bootstrap (method="boot"), which is considerably slower but more accurate in some circumstances.

Upvotes: 15

Tom Wenseleers
Tom Wenseleers

Reputation: 8019

There are two packages, lmerTest and emmeans, that can calculate 95% confidence limits for lmer and glmer output. Maybe you can look into those? And coefplot2, I think can do it too (though as Ben points out below, in a not so sophisticated way, from the standard errors on the Wald statistics, as opposed to Kenward-Roger and/or Satterthwaite df approximations used in lmerTest and emmeans)...

Upvotes: 15

user319487
user319487

Reputation:

I'd suggest tab_model() function from sjPlot package as alternative. Clean and readable output ready for markdown. Reference here and examples here.

For those more visually inclined plot_model() from the same package might come handy too.

Alternative solution is via parameters package using model_parameters() function.

Upvotes: 1

denis
denis

Reputation: 5673

To find the coefficient, you can simply use the summary function of lme4

m = lm(resp ~ 0 + var1 + var1:var2) # var1 categorical, var2 continuous
m_summary <- summary(m)

to have all coefficients : 

m_summary$coefficient

If you want the confidence interval, multiply the standart error by 1.96:

CI <- m_summary$coefficient[,"Std. Error"]*1.96
print(CI)

Upvotes: 1

jciloa
jciloa

Reputation: 1137

Not sure when it was added, but now confint() is implemented in lme4. For example the following example works:

library(lme4)
m = lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
confint(m)

Upvotes: 21

julieth
julieth

Reputation: 442

Assuming a normal approximation for the fixed effects (which confint would also have done), we can obtain 95% confidence intervals by

estimate + 1.96*standard error.

The following does not apply to the variance components/random effects. 

library("lme4")
mylm <- lmer(Reaction ~ Days + (Days|Subject),  data =sleepstudy)

# standard error of coefficient

days_se <- sqrt(diag(vcov(mylm)))[2]

# estimated coefficient

days_coef <- fixef(mylm)[2]

upperCI <-  days_coef + 1.96*days_se
lowerCI <-  days_coef  - 1.96*days_se

Upvotes: 9

Dieter Menne
Dieter Menne

Reputation: 10215

I suggest that you use good old lme (in package nlme). It has confint, and if you need confint of contrasts, there is a series of choices (estimable in gmodels, contrast in contrasts, glht in multcomp).

Why p-values and confint are absent in lmer: see http://finzi.psych.upenn.edu/R/Rhelp02a/archive/76742.html .

Upvotes: 8

Related Questions