The **margins** and **prediction** packages
are a combined effort to port the functionality of Stata’s (closed
source) `margins`

command to (open source) R. These tools provide ways of obtaining common
quantities of interest from regression-type models.
**margins** provides “marginal effects” summaries of models
and **prediction** provides unit-specific and sample
average predictions from models. Marginal effects are partial
derivatives of the regression equation with respect to each variable in
the model for each unit in the data; average marginal effects are simply
the mean of these unit-specific partial derivatives over some sample. In
ordinary least squares regression with no interactions or higher-order
term, the estimated slope coefficients are marginal effects. In other
cases and for generalized linear models, the coefficients are not
marginal effects at least not on the scale of the response variable.
**margins** therefore provides ways of calculating the
marginal effects of variables to make these models more
interpretable.

The major functionality of Stata’s `margins`

command -
namely the estimation of marginal (or partial) effects - is provided
here through a single function, `margins()`

. This is an S3
generic method for calculating the marginal effects of covariates
included in model objects (like those of classes “lm” and “glm”). Users
interested in generating predicted (fitted) values, such as the
“predictive margins” generated by Stata’s `margins`

command,
should consider using `prediction()`

from the sibling
project, **prediction**.

Stata’s `margins`

command is very simple and intuitive to
use:

```
. import delimited mtcars.csv
. quietly reg mpg c.cyl##c.hp wt
. margins, dydx(*)
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cyl | .0381376 .5998897 0.06 0.950 -1.192735 1.26901
hp | -.0463187 .014516 -3.19 0.004 -.076103 -.0165343
wt | -3.119815 .661322 -4.72 0.000 -4.476736 -1.762894
------------------------------------------------------------------------------
. marginsplot
```

With **margins** in R, replicating Stata’s results is
incredibly simple using just the `margins()`

method to obtain
average marginal effects and its `summary()`

method to obtain
Stata-like output:

```
library("margins")
x <- lm(mpg ~ cyl * hp + wt, data = mtcars)
(m <- margins(x))
```

`## Average marginal effects`

`## lm(formula = mpg ~ cyl * hp + wt, data = mtcars)`

```
## cyl hp wt
## 0.03814 -0.04632 -3.12
```

`summary(m)`

```
## 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
```

With the exception of differences in rounding, the above results
match identically what Stata’s `margins`

command produces. A
slightly more concise expression relies on the syntactic sugar provided
by `margins_summary()`

:

`margins_summary(x)`

```
## 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
```

Using the `plot()`

method also yields an aesthetically
similar result to Stata’s `marginsplot`

:

`plot(m)`

`margins()`

**margins** is intended as a port of (some of) the
features of Stata’s `margins`

command, which includes
numerous options for calculating marginal effects at the mean values of
a dataset (i.e., the marginal effects at the mean), an average of the
marginal effects at each value of a dataset (i.e., the average marginal
effect), marginal effects at representative values, and any of those
operations on various subsets of a dataset. (The functionality of
Stata’s command to produce *predictive* margins is not ported, as
this is easily obtained from the prediction
package.) In particular, Stata provides the following options:

`at`

: calculate marginal effects at (potentially representative) specified values (i.e., replacing observed values with specified replacement values before calculating marginal effects)`atmeans`

: calculate marginal effects at the mean (MEMs) of a dataset rather than the default behavior of calculating average marginal effects (AMEs)`over`

: calculate marginal effects (including MEMs and/or AMEs at observed or specified values) on subsets of the original data (e.g., the marginal effect of a treatment separately for men and women)

The `at`

argument has been translated into
`margins()`

in a very similar manner. It can be used by
specifying a list of variable names and specified values for those
variables at which to calculate marginal effects, such as
`margins(x, at = list(hp=150))`

. When using `at`

,
`margins()`

constructs modified datasets - using
`build_datalist()`

- containing the specified values and
calculates marginal effects on each modified dataset,
`rbind`

-ing them back into a single “margins” object.

Stata’s `atmeans`

argument is not implemented in
`margins()`

for various reasons, including because it is
possible to achieve the effect manually through an operation like
`data$var <- mean(data$var, na.rm = TRUE)`

and passing the
modified data frame to `margins(x, data = data)`

.

At present, `margins()`

does not implement the
`over`

option. The reason for this is also simple: R already
makes data subsetting operations quite simple using simple
`[`

extraction. If, for example, one wanted to calculate
marginal effects on subsets of a data frame, those subsets can be passed
directly to `margins()`

via the `data`

argument
(as in a call to `prediction()`

).

The rest of this vignette shows how to use `at`

and
`data`

to obtain various kinds of marginal effects, and how
to use plotting functions to visualize those inferences.

We can start by loading the **margins** package:

`library("margins")`

We’ll use a simple example regression model based on the built-in
`mtcars`

dataset:

`x <- lm(mpg ~ cyl + hp * wt, data = mtcars)`

To obtain average marginal effects (AMEs), we simply call
`margins()`

on the model object created by
`lm()`

:

`margins(x)`

`## Average marginal effects`

`## lm(formula = mpg ~ cyl + hp * wt, data = mtcars)`

```
## cyl hp wt
## -0.3652 -0.02527 -3.838
```

The result is a data frame with special class `"margins"`

.
`"margins"`

objects are printed in a tidy summary format, by
default, as you can see above. The only difference between a
`"margins"`

object and a regular data frame are some
additional data frame-level attributes that dictate how the object is
printed.

The default method calculates marginal effects for all variables
included in the model (ala Stata’s `, dydx(*)`

option). To
limit calculation to only a subset of variables, use the
`variables`

argument:

`summary(margins(x, variables = "hp"))`

```
## factor AME SE z p lower upper
## hp -0.0253 0.0105 -2.4046 0.0162 -0.0459 -0.0047
```

In an ordinary least squares regression, there is really only one way
of examining marginal effects (that is, on the scale of the outcome
variable). In a generalized linear model (e.g., logit), however, it is
possible to examine true “marginal effects” (i.e., the marginal
contribution of each variable on the scale of the linear predictor) or
“partial effects” (i.e., the contribution of each variable on the
outcome scale, conditional on the other variables involved in the link
function transformation of the linear predictor). The latter are the
default in `margins()`

, which implicitly sets the argument
`margins(x, type = "response")`

and passes that through to
`prediction()`

methods. To obtain the former, simply set
`margins(x, type = "link")`

. There’s some debate about which
of these is preferred and even what to call the two different quantities
of interest. Regardless of all of that, here’s how you obtain
either:

`x <- glm(am ~ cyl + hp * wt, data = mtcars, family = binomial)`

`## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred`

`margins(x, type = "response") # the default`

`## Average marginal effects`

`## glm(formula = am ~ cyl + hp * wt, family = binomial, data = mtcars)`

```
## cyl hp wt
## 0.02156 0.002667 -0.5158
```

`margins(x, type = "link")`

```
## Average marginal effects
## glm(formula = am ~ cyl + hp * wt, family = binomial, data = mtcars)
```

```
## cyl hp wt
## 0.5156 0.05151 -12.24
```

Note that some other packages available for R, as well as Stata’s
`margins`

and `mfx`

packages enable calculation of
so-called “marginal effects at means” (i.e., the marginal effect for a
single observation that has covariate values equal to the means of the
sample as a whole). The substantive interpretation of these is fairly
ambiguous. While it was once common practice to estimate MEMs - rather
than AMEs or MERs - this is now considered somewhat inappropriate
because it focuses on cases that may not exist (e.g., the average of a
0/1 variable is not going to reflect a case that can actually exist in
reality) and we are often interested in the effect of a variable at
multiple possible values of covariates, rather than an arbitrarily
selected case that is deemed “typical” in this way. As such,
`margins()`

defaults to reporting AMEs, unless modified by
the `at`

argument to calculate average “marginal effects for
representative cases” (MERs). MEMs could be obtained by manually
specifying `at`

for every variable in a way that respects the
variables classes and inherent meaning of the data, but that
functionality is not demonstrated here.

`at`

ArgumentThe `at`

argument allows you to calculate marginal effects
at representative cases (sometimes “MERs”) or marginal effects at means
- or any other statistic - (sometimes “MEMs”), which are marginal
effects for particularly interesting (sets of) observations in a
dataset. This differs from marginal effects on subsets of the original
data (see the next section for a demonstration of that) in that it
operates on a modified set of the full dataset wherein particular
variables have been replaced by specified values. This is helpful
because it allows for calculation of marginal effects for
*counterfactual* datasets (e.g., what if all women were instead
men? what if all democracies were instead autocracies? what if all
foreign cars were instead domestic?).

As an example, if we wanted to know if the marginal effect of
horsepower (`hp`

) on fuel economy differed across different
types of automobile transmissions, we could simply use `at`

to obtain separate marginal effect estimates for our data as if every
car observation were a manual versus if every car observation were an
automatic. The output of `margins()`

is a simplified summary
of the estimated marginal effects across the requested variable
levels/combinations specified in `at`

:

```
x <- lm(mpg ~ cyl + wt + hp * am, data = mtcars)
margins(x, at = list(am = 0:1))
```

`## Average marginal effects at specified values`

`## lm(formula = mpg ~ cyl + wt + hp * am, data = mtcars)`

```
## at(am) cyl wt hp am
## 0 -0.9339 -2.812 -0.008945 1.034
## 1 -0.9339 -2.812 -0.026392 1.034
```

Because of the `hp * am`

interaction in the regression,
the marginal effect of horsepower differs between the two sets of
results. We can also specify more than one variable to `at`

,
creating a potentially long list of marginal effects results. For
example, we can produce marginal effects at both levels of
`am`

and the values from the five-number summary (minimum,
Q1, median, Q3, and maximum) of observed values of `hp`

. This
produces 2 * 5 = 10 sets of marginal effects estimates:

`margins(x, at = list(am = 0:1, hp = fivenum(mtcars$hp)))`

`## Average marginal effects at specified values`

`## lm(formula = mpg ~ cyl + wt + hp * am, data = mtcars)`

```
## at(am) at(hp) cyl wt hp am
## 0 52 -0.9339 -2.812 -0.008945 2.6864
## 1 52 -0.9339 -2.812 -0.026392 2.6864
## 0 96 -0.9339 -2.812 -0.008945 1.9188
## 1 96 -0.9339 -2.812 -0.026392 1.9188
## 0 123 -0.9339 -2.812 -0.008945 1.4477
## 1 123 -0.9339 -2.812 -0.026392 1.4477
## 0 180 -0.9339 -2.812 -0.008945 0.4533
## 1 180 -0.9339 -2.812 -0.026392 0.4533
## 0 335 -0.9339 -2.812 -0.008945 -2.2509
## 1 335 -0.9339 -2.812 -0.026392 -2.2509
```

Because this is a linear model, the marginal effects of
`cyl`

and `wt`

do not vary across levels of
`am`

or `hp`

. The minimum and Q1 value of
`hp`

are also the same, so the marginal effects of
`am`

are the same in the first two results. As you can see,
however, the marginal effect of `hp`

differs when
`am == 0`

versus `am == 1`

(first and second rows)
and the marginal effect of `am`

differs across levels of
`hp`

(e.g., between the first and third rows). As should be
clear, the `at`

argument is incredibly useful for getting a
better grasp of the marginal effects of different covariates.

This becomes especially apparent when a model includes power-terms (or any other alternative functional form of a covariate). Consider, for example, the simple model of fuel economy as a function of weight, with weight included as both a first- and second-order term:

```
x <- lm(mpg ~ wt + I(wt^2), data = mtcars)
summary(x)
```

```
##
## Call:
## lm(formula = mpg ~ wt + I(wt^2), data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.483 -1.998 -0.773 1.462 6.238
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.9308 4.2113 11.856 1.21e-12 ***
## wt -13.3803 2.5140 -5.322 1.04e-05 ***
## I(wt^2) 1.1711 0.3594 3.258 0.00286 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.651 on 29 degrees of freedom
## Multiple R-squared: 0.8191, Adjusted R-squared: 0.8066
## F-statistic: 65.64 on 2 and 29 DF, p-value: 1.715e-11
```

Looking only at the regression results table, it is actually quite
difficult to understand the effect of `wt`

on fuel economy
because it requires performing mental multiplication and addition on all
possible values of `wt`

. Using the `at`

option to
margins, you could quickly obtain a sense of the average marginal effect
of `wt`

at a range of plausible values:

`margins(x, at = list(wt = fivenum(mtcars$wt)))`

`## Average marginal effects at specified values`

`## lm(formula = mpg ~ wt + I(wt^2), data = mtcars)`

```
## at(wt) wt
## 1.513 -9.8366
## 2.542 -7.4254
## 3.325 -5.5926
## 3.650 -4.8314
## 5.424 -0.6764
```

The marginal effects in the first column of results reveal that the
average marginal effect of `wt`

is large and negative except
when `wt`

is very large, in which case it has an effect not
distinguishable from zero. We can easily plot these results using the
`cplot()`

function to see the effect visually in terms of
either predicted fuel economy or the marginal effect of
`wt`

:

`cplot(x, "wt", what = "prediction", main = "Predicted Fuel Economy, Given Weight")`

`cplot(x, "wt", what = "effect", main = "Average Marginal Effect of Weight")`

A really nice feature of Stata’s margins command is that it handles
factor variables gracefully. This functionality is difficult to emulate
in R, but the `margins()`

function does its best. Here we see
the marginal effects of a simple regression that includes a factor
variable:

```
x <- lm(mpg ~ factor(cyl) * hp + wt, data = mtcars)
margins(x)
```

`## Average marginal effects`

`## lm(formula = mpg ~ factor(cyl) * hp + wt, data = mtcars)`

```
## hp wt cyl6 cyl8
## -0.04475 -3.06 1.473 0.8909
```

`margins()`

recognizes the factor and displays the
marginal effect for each level of the factor separately. (Caveat: this
may not work with certain `at`

specifications, yet.)

Stata’s `margins`

command includes an `over()`

option, which allows you to very easily calculate marginal effects on
subsets of the data (e.g., separately for men and women). This is useful
in Stata because the program only allows one dataset in memory. Because
R does not impose this restriction and further makes subsetting
expressions very simple, that feature is not really useful and can be
achieved using standard subsetting notation in R:

```
x <- lm(mpg ~ factor(cyl) * am + hp + wt, data = mtcars)
# automatic vehicles
margins(x, data = mtcars[mtcars$am == 0, ])
```

`## Average marginal effects`

`## lm(formula = mpg ~ factor(cyl) * am + hp + wt, data = mtcars)`

```
## am hp wt cyl6 cyl8
## 1.706 -0.03115 -2.441 -1.715 -1.586
```

```
# manual vehicles
margins(x, data = mtcars[mtcars$am == 1, ])
```

```
## Average marginal effects
## lm(formula = mpg ~ factor(cyl) * am + hp + wt, data = mtcars)
```

```
## am hp wt cyl6 cyl8
## 2.167 -0.03115 -2.441 -4.218 -2.656
```

Because a `"margins"`

object is just a data frame, it is
also possible to obtain the same result by subsetting the
*output* of `margins()`

:

```
m <- margins(x)
split(m, m$am)
```

`## $`0``

`## Average marginal effects`

`## lm(formula = mpg ~ factor(cyl) * am + hp + wt, data = mtcars)`

```
## am hp wt cyl6 cyl8
## 1.706 -0.03115 -2.441 -1.715 -1.586
##
## $`1`
```

```
## Average marginal effects
## lm(formula = mpg ~ factor(cyl) * am + hp + wt, data = mtcars)
```

```
## am hp wt cyl6 cyl8
## 2.167 -0.03115 -2.441 -4.218 -2.656
```

Using `margins()`

to calculate marginal effects enables
several kinds of plotting. The built-in `plot()`

method for
objects of class `"margins"`

creates simple diagnostic plots
for examining the output of `margins()`

in visual rather than
tabular format. It is also possible to use the output of
`margins()`

to produce more typical marginal effects plots
that show the marginal effect of one variable across levels of another
variable. This section walks through the `plot()`

method and
then shows how to produce marginal effects plots using base
graphics.

`plot()`

method for “margins” objectsThe **margins** package implements a `plot()`

method for objects of class `"margins"`

(seen above). This
produces a plot similar (in spirit) to the output of Stata’s
`marginsplot`

. It is highly customizable, but is meant
primarily as a diagnostic tool to examine the results of
`margins()`

. It simply produces, by default, a plot of
marginal effects along with 95% confidence intervals for those effects.
The confidence level can be modified using the `levels`

argument, which is vectorized to allow multiple levels to be specified
simultaneously.

There are two common ways of visually representing the substantive results of a regression model: (1) fitted values plots, which display the fitted conditional mean outcome across levels of a covariate, and (2) marginal effects plots, which display the estimated marginal effect of a variable across levels of a covariate. This section discusses both approaches.

Fitted value plots can be created using `cplot()`

(to
provide *conditional* predicted value plots or
*conditional* effect plots) and both the `persp()`

method and `image()`

method for `"lm"`

objects,
which display the same type of relationships in three-dimensions (i.e.,
across two conditioning covariates).

For example, we can use `cplot()`

to quickly display the
predicted fuel economy of a vehicle from a model:

```
x <- lm(mpg ~ cyl + wt * am, data = mtcars)
cplot(x, "cyl")
```