15.3 Packages

15.3.1 MarginalEffects

MarginalEffects package is a successor of margins and emtrends (faster, more efficient, more adaptable). Hence, this is advocated to be used.

  • A limitation is that there is no readily function to correct for multiple comparisons. Hence, one can use the p.adjust function to overcome this disadvantage.

Definitions from the package:

  • Marginal effects are slopes or derivatives (i.e., effect of changes in a variable on the outcome)

    • margins package defines marginal effects as “partial derivatives of the regression equation with respect to each variable in the model for each unit in the data.”
  • Marginal means are averages or integrals (i.e., marginalizing across rows of a prediction grid)

To customize your plot using plot_cme (which is a ggplot class), you can check this post by the author of the MarginalEffects package

Causal inference with the parametric g-formula

  • Because the plug-in g estimator is equivalent to the average contrast in the marginaleffects package.

To get predicted values

library(marginaleffects)
library(tidyverse)
data(mtcars)

mod <- lm(mpg ~ hp * wt * am, data = mtcars)
predictions(mod) %>% head()
#> 
#>  Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
#>      22.5      0.884 25.4   <0.001 471.7  20.8   24.2
#>      20.8      1.194 17.4   <0.001 223.3  18.5   23.1
#>      25.3      0.709 35.7   <0.001 922.7  23.9   26.7
#>      20.3      0.704 28.8   <0.001 601.5  18.9   21.6
#>      17.0      0.712 23.9   <0.001 416.2  15.6   18.4
#>      19.7      0.875 22.5   <0.001 368.8  17.9   21.4
#> 
#> Columns: rowid, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, mpg, hp, wt, am
# for specific reference values
predictions(mod, newdata = datagrid(am = 0, wt = c(2, 4)))
#> 
#>  am wt Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %  hp
#>   0  2     22.0       2.04 10.8   <0.001  87.4  18.0   26.0 147
#>   0  4     16.6       1.08 15.3   <0.001 173.8  14.5   18.7 147
#> 
#> Columns: rowid, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, mpg, hp, am, wt
plot_cap(mod, condition = c("hp","wt"))

# Average Margianl Effects
mfx <- marginaleffects(mod, variables = c("hp", "wt"))
summary(mfx)
#> 
#>  Term    Contrast Estimate Std. Error     z Pr(>|z|)   2.5 % 97.5 %
#>    hp mean(dY/dX)  -0.0381     0.0128 -2.98  0.00291 -0.0631 -0.013
#>    wt mean(dY/dX)  -3.9391     1.0858 -3.63  < 0.001 -6.0672 -1.811
#> 
#> Columns: term, contrast, estimate, std.error, statistic, p.value, conf.low, conf.high

# Group-Average Marginal Effects
marginaleffects::marginaleffects(mod, by = "hp", variables = "am")
#> 
#>  Term          Contrast  hp Estimate Std. Error      z Pr(>|z|)   S  2.5 %
#>    am mean(1) - mean(0)  52    3.976       5.20  0.764    0.445 1.2  -6.22
#>    am mean(1) - mean(0)  62   -2.774       2.51 -1.107    0.268 1.9  -7.68
#>    am mean(1) - mean(0)  65    2.999       4.13  0.725    0.468 1.1  -5.10
#>    am mean(1) - mean(0)  66    2.025       3.48  0.582    0.561 0.8  -4.80
#>    am mean(1) - mean(0)  91    1.858       2.76  0.674    0.500 1.0  -3.54
#>    am mean(1) - mean(0)  93    1.201       2.35  0.511    0.609 0.7  -3.40
#>    am mean(1) - mean(0)  95   -1.832       1.97 -0.931    0.352 1.5  -5.69
#>    am mean(1) - mean(0)  97    0.708       2.04  0.347    0.728 0.5  -3.28
#>    am mean(1) - mean(0) 105   -2.682       2.37 -1.132    0.258 2.0  -7.32
#>    am mean(1) - mean(0) 109   -0.237       1.59 -0.149    0.881 0.2  -3.35
#>    am mean(1) - mean(0) 110   -0.640       1.57 -0.407    0.684 0.5  -3.73
#>    am mean(1) - mean(0) 113    4.081       3.94  1.037    0.300 1.7  -3.63
#>    am mean(1) - mean(0) 123   -2.098       2.10 -0.998    0.318 1.7  -6.22
#>    am mean(1) - mean(0) 150   -1.429       1.90 -0.753    0.452 1.1  -5.15
#>    am mean(1) - mean(0) 175   -0.416       1.56 -0.266    0.790 0.3  -3.48
#>    am mean(1) - mean(0) 180   -1.381       2.47 -0.560    0.576 0.8  -6.22
#>    am mean(1) - mean(0) 205   -2.873       6.24 -0.460    0.645 0.6 -15.11
#>    am mean(1) - mean(0) 215   -2.534       6.95 -0.364    0.716 0.5 -16.16
#>    am mean(1) - mean(0) 230   -1.477       7.07 -0.209    0.835 0.3 -15.34
#>    am mean(1) - mean(0) 245    1.115       2.28  0.488    0.625 0.7  -3.36
#>    am mean(1) - mean(0) 264    2.106       2.29  0.920    0.358 1.5  -2.38
#>    am mean(1) - mean(0) 335    4.027       3.24  1.243    0.214 2.2  -2.32
#>  97.5 %
#>   14.18
#>    2.14
#>   11.10
#>    8.85
#>    7.26
#>    5.80
#>    2.02
#>    4.70
#>    1.96
#>    2.87
#>    2.45
#>   11.79
#>    2.02
#>    2.29
#>    2.64
#>    3.46
#>    9.36
#>   11.09
#>   12.39
#>    5.59
#>    6.59
#>   10.38
#> 
#> Columns: term, contrast, hp, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted

# Marginal effects at representative values
marginaleffects::marginaleffects(mod, 
                                 newdata = datagrid(am = 0, 
                                                    wt = c(2, 4)))
#> 
#>  Term Contrast am wt Estimate Std. Error      z Pr(>|z|)   S   2.5 %   97.5 %
#>    am    1 - 0  0  2   2.5465     2.7860  0.914   0.3607 1.5 -2.9139  8.00694
#>    am    1 - 0  0  4  -2.9661     3.0381 -0.976   0.3289 1.6 -8.9207  2.98852
#>    hp    dY/dX  0  2  -0.0598     0.0283 -2.115   0.0344 4.9 -0.1153 -0.00439
#>    hp    dY/dX  0  4  -0.0309     0.0187 -1.654   0.0981 3.3 -0.0676  0.00572
#>    wt    dY/dX  0  2  -2.6762     1.4194 -1.885   0.0594 4.1 -5.4582  0.10587
#>    wt    dY/dX  0  4  -2.6762     1.4199 -1.885   0.0595 4.1 -5.4591  0.10676
#> 
#> Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, am, wt, predicted_lo, predicted_hi, predicted, mpg, hp

# Marginal Effects at the Mean
marginaleffects::marginaleffects(mod, newdata = "mean")
#> 
#>  Term Contrast Estimate Std. Error      z Pr(>|z|)    S  2.5 %  97.5 %
#>    am    1 - 0  -0.8086    1.52383 -0.531  0.59568  0.7 -3.795  2.1781
#>    hp    dY/dX  -0.0323    0.00956 -3.375  < 0.001 10.4 -0.051 -0.0135
#>    wt    dY/dX  -3.7959    1.21310 -3.129  0.00175  9.2 -6.174 -1.4183
#> 
#> Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, mpg, hp, wt, am
# counterfactual
comparisons(mod, variables = list(am = 0:1)) %>% summary()
#> 
#>  Term          Contrast Estimate Std. Error      z Pr(>|z|) 2.5 % 97.5 %
#>    am mean(1) - mean(0)  -0.0481       1.85 -0.026    0.979 -3.68   3.58
#> 
#> Columns: term, contrast, estimate, std.error, statistic, p.value, conf.low, conf.high

15.3.2 margins

  • Marginal effects are partial derivative of the regression equation with respect to each variable in the model for each unit in the data
  • Average Partial Effects: the contribution of each variable the outcome scale, conditional on the other variables involved in the link function transformation of the linear predictor

  • Average Marginal Effects: the marginal contribution of each variable on the scale of the linear predictor.

  • Average marginal effects are the mean of these unit-specific partial derivatives over some sample

margins package gives the marginal effects of models (a replication of the margins command in Stata).

prediction package gives the unit-specific and sample average predictions of models (similar to the predictive margins in Stata).

library(margins)

# examples by the package's authors
mod <- lm(mpg ~ cyl * hp + wt, data = mtcars)
summary(mod)
#> 
#> Call:
#> lm(formula = mpg ~ cyl * hp + wt, data = mtcars)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.3440 -1.4144 -0.6166  1.2160  4.2815 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 52.017520   4.916935  10.579 4.18e-11 ***
#> cyl         -2.742125   0.800228  -3.427  0.00197 ** 
#> hp          -0.163594   0.052122  -3.139  0.00408 ** 
#> wt          -3.119815   0.661322  -4.718 6.51e-05 ***
#> cyl:hp       0.018954   0.006645   2.852  0.00823 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 2.242 on 27 degrees of freedom
#> Multiple R-squared:  0.8795, Adjusted R-squared:  0.8616 
#> F-statistic: 49.25 on 4 and 27 DF,  p-value: 5.065e-12

In cases where you have interaction or polynomial terms, the coefficient estimates cannot be interpreted as the marginal effects of X on Y. Hence, if you want to know the average marginal effects of each variable then

summary(margins(mod))
#>  factor     AME     SE       z      p   lower   upper
#>     cyl  0.0381 0.5999  0.0636 0.9493 -1.1376  1.2139
#>      hp -0.0463 0.0145 -3.1909 0.0014 -0.0748 -0.0179
#>      wt -3.1198 0.6613 -4.7176 0.0000 -4.4160 -1.8236

# equivalently 
margins_summary(mod)
#>  factor     AME     SE       z      p   lower   upper
#>     cyl  0.0381 0.5999  0.0636 0.9493 -1.1376  1.2139
#>      hp -0.0463 0.0145 -3.1909 0.0014 -0.0748 -0.0179
#>      wt -3.1198 0.6613 -4.7176 0.0000 -4.4160 -1.8236

plot(margins(mod))

Marginal effects at the mean (MEM):

  • Marginal effects at the mean values of the sample
  • For discrete variables, it’s called average discrete change (ADC)

Average Marginal Effect (AME)

  • An average of the marginal effects at each value of the sample

Marginal Effects at representative values (MER)

margins(mod, at = list(hp = 150))
#>  at(hp)    cyl       hp    wt
#>     150 0.1009 -0.04632 -3.12

margins(mod, at = list(hp = 150)) %>% summary()
#>  factor       hp     AME     SE       z      p   lower   upper
#>     cyl 150.0000  0.1009 0.6128  0.1647 0.8692 -1.1001  1.3019
#>      hp 150.0000 -0.0463 0.0145 -3.1909 0.0014 -0.0748 -0.0179
#>      wt 150.0000 -3.1198 0.6613 -4.7175 0.0000 -4.4160 -1.8236

15.3.3 mfx

Works well with Generalized Linear Models/glm package

For technical details, see the package vignette

Model Dependent Variable Syntax
Probit Binary probitmfx
Logit Binary logitmfx
Poisson Count poissonmfx
Negative Binomial Count negbinmfx
Beta Rate betamfx
library(mfx)
data("mtcars")
poissonmfx(formula = vs ~ mpg * cyl * disp, data = mtcars)
#> Call:
#> poissonmfx(formula = vs ~ mpg * cyl * disp, data = mtcars)
#> 
#> Marginal Effects:
#>                    dF/dx   Std. Err.       z  P>|z|
#> mpg           1.4722e-03  8.7531e-03  0.1682 0.8664
#> cyl           6.6420e-03  3.9263e-02  0.1692 0.8657
#> disp          1.5899e-04  9.4555e-04  0.1681 0.8665
#> mpg:cyl      -3.4698e-04  2.0564e-03 -0.1687 0.8660
#> mpg:disp     -7.6794e-06  4.5545e-05 -0.1686 0.8661
#> cyl:disp     -3.3837e-05  1.9919e-04 -0.1699 0.8651
#> mpg:cyl:disp  1.6812e-06  9.8919e-06  0.1700 0.8650

This package can only give the marginal effect for each variable in the glm model, but not the average marginal effect that we might look for.