I have a logistic regression and one dummy and one continuous variable and their interaction term. What is the right way to obtain marginal effects, i.e. should one use log odds estimates or exponentiated into odds ratios?


migrated from stackoverflow.com Apr 4 at 12:33

This question came from our site for professional and enthusiast programmers.

  • $\begingroup$ This has been covered in several topics. Bottom line: form the contrasts you want, compute them, and get simultaneous confidence intervals for them. Example: difference in log odds between group B and group A at age=30. Also form double-difference (interaction) contrasts: the difference above minus the same difference at age=40. $\endgroup$ – Frank Harrell Apr 4 at 13:11
  • $\begingroup$ See stats.stackexchange.com/questions/246873/… $\endgroup$ – kjetil b halvorsen Apr 5 at 17:57
  • $\begingroup$ ok, I directly added some working examples to my answer. $\endgroup$ – Daniel Apr 11 at 14:04

One typical way is to compute predicted probabilities to investigate marginal effects. You can do this with eg the ggeffects package, see examples here, where you also find examples for interactions.

You find a concrete example for logistic regression with interaction between continuous and categorical predictors here.

Here is a code-example, marginal effects computed with different packages. The emmeans-package returns marginal effects on the link-scale by default. However, this is probably less intuitive to understand, and in this example I backtransformed the marginal effects.

To avoid redundance, I only show one plot. You'll see that all plots produced by this code-example are essentially identical.


# create dummy data

data <- data.frame(
  outcome = rbinom(100, 1, 0.5),
  var1 = rbinom(100, 1, 0.1),
  var2 = rnorm(100, 10, 7)

# fit example model
m <- glm(
  outcome ~ var1 * var2, 
  data = data, 
  family = binomial(link = "logit")

# with ggeffects-package, using "predict()
ggpredict(m, c("var2", "var1")) %>% plot()

# with ggeffects-package, using "effect()
ggeffect(m, c("var2", "var1")) %>% plot()

# with effects-package
eff <- as.data.frame(Effect(c("var1", "var2"), m, xlevels = list(var1 = c(0, 1))))

ggplot(eff, aes(x = var2, y = fit, colour = as.factor(var1))) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = as.factor(var1)), alpha = .1) +

# with emmeans
eff <- as.data.frame(emmeans(
  m, c("var1", "var2"), 
  at = list(var1 = c(0, 1), var2 = seq(-8, 30, 2))

# we get estimated marginal means on link-scale, 
# so get link-inverse function to back-transform to probabilities
linv <- insight::link_inverse(m)
eff$emmean <- linv(eff$emmean)
eff$asymp.LCL <- linv(eff$asymp.LCL)
eff$asymp.UCL <- linv(eff$asymp.UCL)

ggplot(eff, aes(x = var2, y = emmean, colour = as.factor(var1))) + 
  geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL, fill = as.factor(var1)), alpha = .1) +

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.