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

Question

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?

Quinby

To get you started until someone like @BenBolker shows up (an expert): ?lmer lists methods fixef and ranef in addition to coef. Since your error says it's having trouble combining the two, the issue is likely that your model specification is somehow "unusual". – Cuttie

Thanks @joran. My model spec is maybe unusual in omitting the intercept - I want to do this, because otherwise the coefficients are nonsense. var1 is categorical and I want "group specific intercepts" for each its category. If I allow the intercept (remove 0 + from formula), coef runs but doesn't give what I expect. fixef works great, thanks! However the confint doesn't work at all. – Quinby

I would extract the data you need directly from the S4 object -- see this post's answers: #8527181 – Holsworth

Thanks @baha-kev, but are you sure the confidence intervals are in this object? I don't think so... – Quinby

Not directly, but you would just need to multiply the standard error by +/- 1.96 and add it to the coefficient estimates for a 95% confidence interval, for example. The S4 object has the standard errors and the coefficient estimates. – Holsworth

I am fixing the bug(let)? in coef in the r-forge versions of lme4 (lme4.0, the currently stable branch which corresponds to CRAN-lme4), and lme4, the development branch). confint is a bigger can of worms, as has been discussed, although the development branch of lme4 can calculate profile confidence intervals ... – Considerate

Answer 1

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.

Considerate

Hi Ben, thanks! I'm a bit confused, what does the standalone dot . mean? If it's just a model var name, why not use e.g. m? :-) – Quinby

I could use m. Sometimes I use . as a placeholder. – Considerate

Answer 2

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)

Africa

Answer 3

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)...

Cremator

coefplot2 does it very naively, by computing 1.96 times the Wald standard errors -- it doesn't address the very significant issues of finite-size corrections to the CIs – Considerate

Check also this post stats.stackexchange.com/questions/117641/… for a more detailed answer – Cremator

lmerTest is now nicely described in JoSS jstatsoft.org/article/view/v082i13 – Kyla

Note that many of these comments are now quite outdated. Using emmeans or lmerTest is the way to go, and there are plotting methods now. – Subject

Answer 4

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

Griceldagrid

Hi julieth, nice idea, however there is a difference between the real confidence intervals (computed by confint) and these .... Maybe the t-distribution would give the same result as confint (not sure about this though), but in this case I don't know the df which should be used. – Quinby

In other words, this is the reason why I prefer to use functions like confint etc. to do all this for me... (especially if I'm not sure about the normal distribution of coefficients). – Quinby

The t-distribution is asymptotically normal and the degrees of freedom for the error term in many multi-level designs is so high that the error distribution is normal at that point. Therefore, if you have a design with lots of degrees of freedom this is a perfectly reasonable confidence interval estimate. – Coulter

Answer 5

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 .

Rimple

Thanks Dieter, I will try the older package. The absence of p-value - and possibility to tell significance right away - also alarmed me! Doesn't make any sense to me, if I will be able to get confidence interval then I will simply look whether it contains zero - and have the significance anyway! Regards, – Quinby

I forget to mention that confint(glht... from package multcomp give asymptotic confidence intervals for lmer. Douglas Bates caveats still apply, but his bold move to leave the p-value out of lmer/gaussian certainly has stirred the soup. – Rimple

Dieter, what do you mean with "confint(glht" ? There's no confint function in multcomp package... – Quinby

Dieter, I tried the old package lme, nice, it has p-values. But my main concern is to get the confidence interval of fixed effect coefficients. How do I do that? confint returns some big matrix, glht seems too complicated.. – Quinby

There is confint.glht. See examples and cran.r-project.org/web/packages/multcomp/vignettes/… – Rimple

using intervals(mix1) will you give you asymptotic confidence intervals as in @julieth's answer below; intervals(mix1)$fixed extracts the fixed-effect intervals. These are based on the normal approximation, not the t distribution or anything more exotic ... – Considerate

You are right for the default table; I was thinking of non-default contrast. – Rimple

Answer 6

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)

Flaxseed

here the 1.96 factor is for 95% confidence interval of course – Flaxseed

Answer 7

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.

Kyla

or broom.mixed::tidy() – Considerate

@Ben Bolker's should be the top rated answer IMO :) – Ronironica