16

I have made a model that looks at a number of variables and the effect that has on pregnancy outcome. The outcome is a grouped binary. A mob of animals will have 34 pregnant and 3 empty, the next will have 20 pregnant and 4 empty and so on.

I have modelled this data using the glmer function where y is the pregnancy outcome (pregnant or empty).

mclus5 <- glmer(y~adg + breed + bw_start + year + (1|farm),
                data=dat, family=binomial)

I get all the usual output with coefficients etc. but for interpretation I would like to transform this into odds ratios and confidence intervals for each of the coefficients.

In past logistic regression models I have used the following code

round(exp(cbind(OR=coef(mclus5),confint(mclus5))),3)

This would very nicely provide what I want, but it does not seem to work with the model I have run.

Does anyone know a way that I can get this output for my model through R?

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
EmmaC
  • 169
  • 1
  • 1
  • 3
  • Type this: `?fixef` and `class(mclus5)` and `showMethods("fixef", includeDefs=TRUE)`. `glmer` is an S4 function. – IRTFM Oct 17 '14 at 04:56

3 Answers3

22

The only real difference is that you have to use fixef() rather than coef() to extract the fixed-effect coefficients (coef() gives you the estimated coefficients for each group).

I'll illustrate with a built-in example from the lme4 package.

library("lme4")
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               data = cbpp, family = binomial)

Fixed-effect coefficients and confidence intervals, log-odds scale:

cc <- confint(gm1,parm="beta_")  ## slow (~ 11 seconds)
ctab <- cbind(est=fixef(gm1),cc)

(If you want faster-but-less-accurate Wald confidence intervals you can use confint(gm1,parm="beta_",method="Wald") instead; this will be equivalent to @Gorka's answer but marginally more convenient.)

Exponentiate to get odds ratios:

rtab <- exp(ctab)
print(rtab,digits=3)
##               est 2.5 % 97.5 %
## (Intercept) 0.247 0.149  0.388
## period2     0.371 0.199  0.665
## period3     0.324 0.165  0.600
## period4     0.206 0.082  0.449

A marginally simpler/more general solution:

library(broom.mixed)
tidy(gm1,conf.int=TRUE,exponentiate=TRUE,effects="fixed")

for Wald intervals, or add conf.method="profile" for profile confidence intervals.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
6

I believe there is another, much faster way (if you are OK with a less accurate result).

From: http://www.ats.ucla.edu/stat/r/dae/melogit.htm

First we get the confidence intervals for the Estimates

se <- sqrt(diag(vcov(mclus5)))
# table of estimates with 95% CI
tab <- cbind(Est = fixef(mclus5), LL = fixef(mclus5) - 1.96 * se, UL = fixef(mclus5) + 1.96 * se)

Then the odds ratios with 95% CI

print(exp(tab), digits=3)
Gorka
  • 3,555
  • 1
  • 31
  • 37
3

Other option I believe is to just use package emmeans :

library(emmeans)
data.frame(confint(pairs(emmeans(fit, ~ factor_name,type="response"))))
Tom Wenseleers
  • 7,535
  • 7
  • 63
  • 103
  • 1
    No need to wrap the call in `data.frame`, though. In fact, doing so suppresses messages that may be useful. – Russ Lenth Jun 15 '16 at 13:34
  • Ha yes that's a good point - I usually make a data frame out of it to save it to a csv file though. But probably best to do that in two steps then... – Tom Wenseleers Jun 16 '16 at 00:32
  • 1
    note that this probably wants to switch to using `emmeans` now. – Ben Bolker Oct 23 '18 at 20:55
  • 1
    Yes of course one can also use the newer packae emmeans now, good point! – Tom Wenseleers Oct 25 '18 at 05:48
  • @BenBolker, the results from `emmeans` are different from the ones you provided in the [previous answer](https://stackoverflow.com/a/26417108/1060349), including the ones given by `broom.mixed`. – toto_tico Mar 03 '23 at 23:03
  • This answer is providing confidence intervals on pairwise contrasts. `emmeans()` is designed for providing estimates and CIs for expected marginal means, contrasts, etc. - different from the parameters themselves. If you want an explanation of the differences (which might be better for [CrossValidated](https://stats.stackexchange.com), please post a new question with a reproducible example and a clearly articulated question ... – Ben Bolker Mar 04 '23 at 00:09