Skip to contents

Marginal effects

We can summarize the results of the comparisons() or slopes() functions using the modelsummary package.

library(modelsummary)
library(marginaleffects)

mod <- glm(am ~ wt + drat, family = binomial, data = mtcars)
mfx <- slopes(mod)

modelsummary(mfx)
 (1)
wt −0.217
(0.080)
drat 0.278
(0.168)
Num.Obs. 32
AIC 22.0
BIC 26.4
Log.Lik. −8.011
F 3.430
RMSE 0.28

The same results can be visualized with modelplot():

Contrasts

When using the comparisons() function (or the slopes() function with categorical variables), the output will include two columns to uniquely identify the quantities of interest: term and contrast.

dat <- mtcars
dat$gear <- as.factor(dat$gear)
mod <- glm(vs ~ gear + mpg, data = dat, family = binomial)

cmp <- comparisons(mod)
get_estimates(cmp)
#> # A tibble: 3 × 8
#>   term  contrast          estimate std.error statistic    p.value conf.low conf.high
#>   <chr> <chr>                <dbl>     <dbl>     <dbl>      <dbl>    <dbl>     <dbl>
#> 1 gear  mean(4) - mean(3)   0.0372    0.137      0.272 0.785       -0.230     0.305 
#> 2 gear  mean(5) - mean(3)  -0.340     0.0988    -3.44  0.000588    -0.533    -0.146 
#> 3 mpg   mean(+1)            0.0608    0.0128     4.74  0.00000218   0.0356    0.0860

We can use the shape argument of the modelsummary function to structure the table properly:

modelsummary(cmp, shape = term + contrast ~ model)
 (1)
gear mean(4) - mean(3) 0.037
mean(4) - mean(3) (0.137)
mean(5) - mean(3) −0.340
mean(5) - mean(3) (0.099)
mpg mean(+1) 0.061
mean(+1) (0.013)
Num.Obs. 32
AIC 26.2
BIC 32.1
Log.Lik. −9.101
F 2.389
RMSE 0.31

Cross-contrasts can be a bit trickier, since there are multiple simultaneous groups. Consider this example:

mod <- lm(mpg ~ factor(cyl) + factor(gear), data = mtcars)
cmp <- comparisons(
  mod,
  variables = c("gear", "cyl"),
  cross = TRUE)
get_estimates(cmp)
#> # A tibble: 4 × 9
#>   term  contrast_gear     contrast_cyl      estimate std.error statistic p.value conf.low conf.high
#>   <chr> <chr>             <chr>                <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
#> 1 cross mean(4) - mean(3) mean(6) - mean(4)    -5.33      2.77     -1.93 0.0542     -10.8  0.0953  
#> 2 cross mean(4) - mean(3) mean(8) - mean(4)    -9.22      3.62     -2.55 0.0108     -16.3 -2.13    
#> 3 cross mean(5) - mean(3) mean(6) - mean(4)    -5.16      2.63     -1.96 0.0500     -10.3  0.000166
#> 4 cross mean(5) - mean(3) mean(8) - mean(4)    -9.04      3.19     -2.84 0.00453    -15.3 -2.80

As we can see above, there are two relevant grouping columns: contrast_gear and contrast_cyl. We can simply plug those names in the shape argument:

modelsummary(
  cmp,
  shape = contrast_gear + contrast_cyl ~ model)
gear cyl  (1)
mean(4) - mean(3) mean(6) - mean(4) −5.332
mean(4) - mean(3) mean(6) - mean(4) (2.769)
mean(4) - mean(3) mean(8) - mean(4) −9.218
mean(4) - mean(3) mean(8) - mean(4) (3.618)
mean(5) - mean(3) mean(6) - mean(4) −5.156
mean(5) - mean(3) mean(6) - mean(4) (2.631)
mean(5) - mean(3) mean(8) - mean(4) −9.042
mean(5) - mean(3) mean(8) - mean(4) (3.185)
Num.Obs. 32
R2 0.740
R2 Adj. 0.701
AIC 173.7
BIC 182.5
Log.Lik. −80.838
F 19.190
RMSE 3.03

Marginal means

library("marginaleffects")
library("modelsummary")

dat <- mtcars
dat$cyl <- as.factor(dat$cyl)
dat$am <- as.logical(dat$am)
mod <- lm(mpg ~ hp + cyl + am, data = dat)
mm <- marginal_means(mod)

modelsummary(mm,
             title = "Estimated Marginal Means",
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             group = term + value ~ model)
Estimated Marginal Means
 (1)
cyl 6 18.960 (1.073)***
4 22.885 (1.357)***
8 19.351 (1.377)***
am TRUE 22.478 (0.834)***
FALSE 18.320 (0.785)***
Num.Obs. 32
R2 0.825
R2 Adj. 0.799
AIC 161.0
BIC 169.8
Log.Lik. −74.502
F 31.794
RMSE 2.48