---
title: "Comparison with Stata's 'margins' command"
date: "`r Sys.Date()`"
output:
  html_document:
    fig_caption: false
    toc: true
    toc_float:
      collapsed: false
      smooth_scroll: false
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Comparison with Stata's 'margins' command}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

**margins** is intended as a port of (some of) the features of Stata's `margins` command. This vignette compares output from Stata's `margins` command for linear models against the output of **margins**.

```{r}
library("margins")
options(width = 100)
```

## OLS marginal effects

### Stata

```
. quietly reg mpg cyl hp wt
. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OLS
Expression   : Linear prediction, predict()
dy/dx w.r.t. : cyl hp wt

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |  -.9416166   .5509165    -1.71   0.098    -2.070118    .1868846
          hp |  -.0180381   .0118763    -1.52   0.140    -.0423655    .0062893
          wt |  -3.166973    .740576    -4.28   0.000    -4.683975   -1.649972
------------------------------------------------------------------------------
```

### R

```{r}
library("margins")
x <- lm(mpg ~ cyl + hp + wt, data = mtcars)
summary(margins(x))
```

---


## OLS marginal effects with interaction

### Stata

```
. quietly reg mpg cyl c.hp##c.wt
. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OLS
Expression   : Linear prediction, predict()
dy/dx w.r.t. : cyl hp wt

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |  -.3652391   .5086204    -0.72   0.479    -1.408842    .6783638
          hp |  -.0252715   .0105097    -2.40   0.023    -.0468357   -.0037073
          wt |  -3.837584   .6730996    -5.70   0.000     -5.21867   -2.456498
------------------------------------------------------------------------------
```

### R

```{r}
x <- lm(mpg ~ cyl + hp * wt, data = mtcars)
summary(margins(x))
```

---

## OLS marginal effects with factor term

### Stata

```
. quietly reg mpg i.cyl hp wt
. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OLS
Expression   : Linear prediction, predict()
dy/dx w.r.t. : 6.cyl 8.cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |
          6  |  -3.359024    1.40167    -2.40   0.024    -6.235014   -.4830353
          8  |  -3.185884   2.170476    -1.47   0.154    -7.639332    1.267564
             |
          hp |  -.0231198   .0119522    -1.93   0.064    -.0476437    .0014041
          wt |  -3.181404   .7196011    -4.42   0.000    -4.657904   -1.704905
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.
```

### R

```{r}
x <- lm(mpg ~ factor(cyl) + hp + wt, data = mtcars)
summary(margins(x))
```


---

## OLS marginal effects with squared term

### Stata

```
. quietly reg mpg cyl c.hp##c.hp wt
. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OLS
Expression   : Linear prediction, predict()
dy/dx w.r.t. : cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |  -.3696041   .6163571    -0.60   0.554    -1.634264    .8950561
          hp |  -.0429018   .0178353    -2.41   0.023    -.0794969   -.0063066
          wt |  -2.873553   .7301251    -3.94   0.001    -4.371646    -1.37546
------------------------------------------------------------------------------

```

### R

```{r}
x <- lm(mpg ~ cyl + hp + I(hp^2) + wt, data = mtcars)
summary(margins(x))
```

---

## OLS marginal effects with squared term (but no first-order term)

### Stata

```
. gen hp2 = hp^2
. quietly reg mpg cyl hp2 wt
. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OLS
Expression   : Linear prediction, predict()
dy/dx w.r.t. : cyl hp2 wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |   -1.21919   .5030753    -2.42   0.022    -2.249693   -.1886869
         hp2 |   -.000028   .0000276    -1.01   0.320    -.0000846    .0000286
          wt |  -3.218637   .7570747    -4.25   0.000    -4.769435    -1.66784
------------------------------------------------------------------------------
```

### R

```{r}
x <- lm(mpg ~ cyl + I(hp^2) + wt, data = mtcars)
summary(margins(x))
```

---

## Logit effects on log-odds and probability scales

### Stata

```
. quietly logit am cyl hp wt
. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Pr(am), predict()
dy/dx w.r.t. : cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |   .0214527   .0469746     0.46   0.648    -.0706157    .1135212
          hp |   .0014339   .0006182     2.32   0.020     .0002224    .0026455
          wt |  -.4025475   .1154098    -3.49   0.000    -.6287466   -.1763484
------------------------------------------------------------------------------

. quietly logit am cyl hp wt
. margins, dydx(*) predict(xb)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Linear prediction (log odds), predict(xb)
dy/dx w.r.t. : cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |   .4875978   1.071621     0.46   0.649    -1.612741    2.587936
          hp |   .0325917   .0188611     1.73   0.084    -.0043753    .0695587
          wt |   -9.14947   4.153326    -2.20   0.028    -17.28984   -1.009101
------------------------------------------------------------------------------
```

### R

```{r}
x <- glm(am ~ cyl + hp + wt, data = mtcars, family = binomial)
# AME
summary(margins(x, type = "response"))
# AME and MEM equivalent on "link" scale
summary(margins(x, type = "link"))
```


---

## Logit effects with factor variable on log-odds and probability scales


### Stata

```
. quietly logit am i.cyl hp wt
. margins, dydx(*) predict(xb)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Linear prediction (log odds), predict(xb)
dy/dx w.r.t. : 6.cyl 8.cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |
          6  |   2.765754   3.156829     0.88   0.381    -3.421517    8.953025
          8  |  -8.388958   13.16745    -0.64   0.524     -34.1967    17.41878
             |
          hp |    .103209   .0960655     1.07   0.283    -.0850759    .2914939
          wt |  -10.67598   5.441998    -1.96   0.050     -21.3421   -.0098575
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.

. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Pr(am), predict()
dy/dx w.r.t. : 6.cyl 8.cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |
          6  |   .1197978   .1062873     1.13   0.260    -.0885214    .3281171
          8  |  -.3478575   .2067542    -1.68   0.092    -.7530883    .0573732
             |
          hp |   .0033268   .0029852     1.11   0.265    -.0025241    .0091777
          wt |  -.3441297   .1188604    -2.90   0.004    -.5770919   -.1111675
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.
```

### R

```{r}
x <- glm(am ~ factor(cyl) + hp + wt, data = mtcars, family = binomial)
# Log-odds
summary(margins(x, type = "link"))
# Probability with continuous factors
summary(margins(x, type = "response"))
```


---

## Logit with interaction on probability and Log-Odds scales

### Stata

```
. quietly logit am cyl c.hp##c.wt
. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Pr(am), predict()
dy/dx w.r.t. : cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |   .0215633   .0492676     0.44   0.662    -.0749994    .1181261
          hp |   .0026673   .0023004     1.16   0.246    -.0018414     .007176
          wt |  -.5157922   .2685806    -1.92   0.055    -1.042201    .0106162
------------------------------------------------------------------------------

. margins, dydx(*) predict(xb)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Linear prediction (log odds), predict(xb)
dy/dx w.r.t. : cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |   .5156396   1.169458     0.44   0.659    -1.776456    2.807735
          hp |   .0515116    .035699     1.44   0.149    -.0184571    .1214804
          wt |  -12.24264   7.678428    -1.59   0.111    -27.29208    2.806807
------------------------------------------------------------------------------
```

### R

```{r}
x <- glm(am ~ cyl + hp * wt, data = mtcars, family = binomial)
# AME
summary(margins(x, type = "response"))
# AME and MEM equivalent on "link" scale
summary(margins(x, type = "link"))
```


## Probit effects on latent and probability scales

### Stata

```
. quietly probit am cyl c.hp##c.wt
. margins, dydx(*) predict(xb)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Linear prediction, predict(xb)
dy/dx w.r.t. : cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |   .2974758   .6629205     0.45   0.654    -1.001825    1.596776
          hp |   .0277713   .0193121     1.44   0.150    -.0100797    .0656223
          wt |  -6.626949   4.096208    -1.62   0.106    -14.65537    1.401471
------------------------------------------------------------------------------

. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Pr(am), predict()
dy/dx w.r.t. : cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |    .022611   .0498253     0.45   0.650    -.0750447    .1202667
          hp |   .0025769   .0022607     1.14   0.254     -.001854    .0070077
          wt |   -.508829   .2625404    -1.94   0.053    -1.023399    .0057408
------------------------------------------------------------------------------
```

### R

```{r}
x <- glm(am ~ cyl + hp * wt, data = mtcars, family = binomial(link="probit"))
# AME (log-odds)
summary(margins(x, type = "link"))
# AME (probability)
summary(margins(x, type = "response"))
```


## Poisson effects on latent and probability scales

### Stata

```
. quietly poisson carb cyl c.hp##c.wt
. margins, dydx(*) predict(xb)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Linear prediction, predict(xb)
dy/dx w.r.t. : cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |  -.0993854   .1478936    -0.67   0.502    -.3892516    .1904808
          hp |   .0066519   .0024217     2.75   0.006     .0019054    .0113984
          wt |   .1225051   .2035185     0.60   0.547    -.2763837     .521394
------------------------------------------------------------------------------

. margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE    : OIM
Expression   : Predicted number of events, predict()
dy/dx w.r.t. : cyl hp wt
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |  -.2795214   .4169931    -0.67   0.503    -1.096813      .53777
          hp |   .0175935   .0067179     2.62   0.009     .0044267    .0307604
          wt |   .2075447   .4859868     0.43   0.669    -.7449719    1.160061
------------------------------------------------------------------------------
```


### R

```{r}
x <- glm(carb ~ cyl + hp * wt, data = mtcars, family = poisson)
# AME (linear/link)
summary(margins(x, type = "link"))
# AME (probability)
summary(margins(x, type = "response"))
```

