16.3 Packages
16.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.adjustfunction 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)
marginspackage 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
marginaleffectspackage.
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
#> Type: response
# Create a data frame with specific reference values
newdata <- datagrid(am = 0, wt = c(2, 4), model = mod)
# Plot predictions for 'hp' and 'wt' with specific reference values
marginaleffects::plot_predictions(mod, newdata = newdata, condition = c("hp", "wt"))
# Average Margianl Effects
mfx <- marginaleffects(mod, variables = c("hp", "wt"))
summary(mfx)
#> rowid term estimate std.error
#> Min. : 1.00 Length:64 Min. :-8.48333 Min. :0.01329
#> 1st Qu.: 8.75 Class :character 1st Qu.:-3.06560 1st Qu.:0.01665
#> Median :16.50 Mode :character Median :-0.07093 Median :0.47802
#> Mean :16.50 Mean :-1.98858 Mean :0.91167
#> 3rd Qu.:24.25 3rd Qu.:-0.03768 3rd Qu.:1.79509
#> Max. :32.00 Max. : 0.63503 Max. :4.88798
#> statistic p.value s.value conf.low
#> Min. :-4.6022 Min. :0.0000042 Min. : 0.1574 Min. :-12.09618
#> 1st Qu.:-2.9093 1st Qu.:0.0036223 1st Qu.: 3.3772 1st Qu.: -7.48128
#> Median :-2.0666 Median :0.0387763 Median : 4.6888 Median : -1.83083
#> Mean :-2.2166 Mean :0.1182531 Mean : 5.9584 Mean : -3.77543
#> 3rd Qu.:-1.6631 3rd Qu.:0.0965502 3rd Qu.: 8.1089 3rd Qu.: -0.06772
#> Max. : 0.1299 Max. :0.8966318 Max. :17.8678 Max. : -0.03297
#> conf.high predicted_lo predicted_hi predicted
#> Min. :-4.870481 Min. :12.05 Min. :12.05 Min. :12.05
#> 1st Qu.:-0.023276 1st Qu.:16.02 1st Qu.:16.02 1st Qu.:16.02
#> Median :-0.003027 Median :19.22 Median :19.22 Median :19.22
#> Mean :-0.201739 Mean :20.09 Mean :20.09 Mean :20.09
#> 3rd Qu.: 0.065697 3rd Qu.:22.71 3rd Qu.:22.70 3rd Qu.:22.70
#> Max. :10.215292 Max. :33.16 Max. :33.15 Max. :33.15
#> mpg hp wt am
#> Min. :10.40 Min. : 52.0 Min. :1.513 Min. :0.0000
#> 1st Qu.:15.43 1st Qu.: 96.5 1st Qu.:2.581 1st Qu.:0.0000
#> Median :19.20 Median :123.0 Median :3.325 Median :0.0000
#> Mean :20.09 Mean :146.7 Mean :3.217 Mean :0.4062
#> 3rd Qu.:22.80 3rd Qu.:180.0 3rd Qu.:3.610 3rd Qu.:1.0000
#> Max. :33.90 Max. :335.0 Max. :5.424 Max. :1.0000
# 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
#> Type: response
# 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.4 -0.0676 0.00571
#> wt dY/dX 0 2 -2.6762 1.4199 -1.885 0.0595 4.1 -5.4591 0.10676
#> wt dY/dX 0 4 -2.6762 1.4206 -1.884 0.0596 4.1 -5.4605 0.10816
#>
#> Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, am, wt, predicted_lo, predicted_hi, predicted, hp, mpg
#> Type: response
# 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.5238 -0.531 0.59568 0.7 -3.7952 2.1781
#> hp dY/dX -0.0422 0.0133 -3.181 0.00147 9.4 -0.0683 -0.0162
#> wt dY/dX -2.6762 1.4193 -1.886 0.05935 4.1 -5.4579 0.1055
#>
#> Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, hp, wt, am, mpg
#> Type: response# counterfactual
comparisons(mod, variables = list(am = 0:1)) %>% summary()
#> rowid term contrast estimate
#> Min. : 1.00 Length:32 Length:32 Min. :-2.87275
#> 1st Qu.: 8.75 Class :character Class :character 1st Qu.:-1.73488
#> Median :16.50 Mode :character Mode :character Median :-0.57928
#> Mean :16.50 Mean :-0.04811
#> 3rd Qu.:24.25 3rd Qu.: 1.30608
#> Max. :32.00 Max. : 4.08073
#> std.error statistic p.value s.value
#> Min. :1.568 Min. :-1.13223 Min. :0.2139 Min. :0.1823
#> 1st Qu.:1.975 1st Qu.:-0.64610 1st Qu.:0.3608 1st Qu.:0.5996
#> Median :2.320 Median :-0.27781 Median :0.5345 Median :0.9036
#> Mean :2.894 Mean :-0.08309 Mean :0.5321 Mean :1.0028
#> 3rd Qu.:3.176 3rd Qu.: 0.56094 3rd Qu.:0.6604 3rd Qu.:1.4709
#> Max. :7.073 Max. : 1.24288 Max. :0.8813 Max. :2.2249
#> conf.low conf.high predicted_lo predicted_hi
#> Min. :-16.163 Min. : 1.957 Min. :10.76 Min. : 9.553
#> 1st Qu.: -6.222 1st Qu.: 2.480 1st Qu.:15.90 1st Qu.:15.012
#> Median : -5.016 Median : 4.038 Median :19.03 Median :16.955
#> Mean : -5.720 Mean : 5.624 Mean :19.32 Mean :19.269
#> 3rd Qu.: -3.507 3rd Qu.: 8.082 3rd Qu.:22.62 3rd Qu.:22.881
#> Max. : -2.324 Max. :14.176 Max. :29.18 Max. :33.155
#> predicted mpg hp wt
#> Min. :12.05 Min. :10.40 Min. : 52.0 Min. :1.513
#> 1st Qu.:16.02 1st Qu.:15.43 1st Qu.: 96.5 1st Qu.:2.581
#> Median :19.22 Median :19.20 Median :123.0 Median :3.325
#> Mean :20.09 Mean :20.09 Mean :146.7 Mean :3.217
#> 3rd Qu.:22.70 3rd Qu.:22.80 3rd Qu.:180.0 3rd Qu.:3.610
#> Max. :33.15 Max. :33.90 Max. :335.0 Max. :5.424
#> am
#> Min. :0.0000
#> 1st Qu.:0.0000
#> Median :0.0000
#> Mean :0.4062
#> 3rd Qu.:1.0000
#> Max. :1.000016.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-12In 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.7175 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.7175 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.823616.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.8650This 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.