15-marginal-effect.Rmd

# Marginal Effects

In cases without polynomials or interactions, it can be easy to interpret the marginal effect.

For example,

$$
Y = \beta_1 X_1 + \beta_2 X_2
$$

where $\beta$ are the marginal effects.

Numerical derivation is easier than analytical derivation.

-   We need to choose values for all the variables to calculate the marginal effect of $X$

Analytical derivation

$$
f'(x) \equiv \lim_{h \to 0} \frac{f(x+h) - f(x)}{h}
$$

E.g., $f(x) = X^2$

$$
\begin{aligned}
f'(x) &= \lim_{h \to 0} \frac{(x+h)^2 - x^2}{h} \\
&= \frac{x^2 + 2xh + h^2 - x^2}{h} \\
&= \frac{2xh + h^2}{h} \\
&= 2x + h \\
&= 2x
\end{aligned}
$$

For numerically approach, we "just" need to find a small $h$ to plug in our function. However, you also need a large enough $h$ to have numerically accurate computation [@gould2010maximum, chapter 1]

Numerically approach

One-sided derivative

$$
\begin{aligned}
f'(x) &= \lim_{h \to 0} \frac{(x+h)^2 - x^2}{h}  \\
& \approx \frac{f(x+h) -f(x)}{h}
\end{aligned}
$$

Alternatively, two-sided derivative

$$
f'_2(x) \approx \frac{f(x+h) - f(x- h)}{2h}
$$

Marginal effects for

-   discrete variables (also known as incremental effects) are the change in $E[Y|X]$ for a one unit change in $X$

-   continuous variables are the change in $E[Y|X]$ for very small changes in $X$ (not unit changes), because it's a derivative, which is a limit when $h \to 0$

|                  | Analytical derivation | Numerical derivation                                   |
|-----------------|-----------------|---------------------------------------|
| Marginal Effects | Rules of expectations | Approximate analytical solution                        |
| Standard Errors  | Rules of variances    | Delta method using the asymptotic errors (vcov matrix) |

## Delta Method

-   approximate the mean and variance of a function of random variables using a first-order Taylor approximation
-   A semi-parametric method
-   Alternative approaches:
    -   Analytically derive a probability function for the margin

    -   Simulation/Bootstrapping
-   Resources:
    -   Advanced: [modmarg](https://cran.r-project.org/web/packages/modmarg/vignettes/delta-method.html)

    -   Intermediate: [UCLA stat](https://stats.oarc.ucla.edu/r/faq/how-can-i-estimate-the-standard-error-of-transformed-regression-parameters-in-r-using-the-delta-method/)

    -   Simple: [Another one](https://www.alexstephenson.me/post/2022-04-02-standard-errors-and-the-delta-method/)

Let $G(\beta)$ be a function of the parameters $\beta$, then

$$
var(G(\beta)) \approx \nabla G(\beta) cov (\beta) \nabla G(\beta)'
$$

where

-   $\nabla G(\beta)$ = the gradient of the partial derivatives of $G(\beta)$ (also known as the Jacobian)

## Average Marginal Effect Algorithm

For one-sided derivative $\frac{\partial p(\mathbf{X},\beta)}{\partial X}$ in the probability scale

1.  Estimate your model
2.  For each observation $i$
    1.  Calculate $\hat{Y}_{i0}$ which is the prediction in the probability scale using observed values

    2.  Increase $X$ (variable of interest) by a "small" amount $h$ ($X_{new} = X + h$)

        -   When $X$ is continuous, $h = (|\bar{X}| + 0.001) \times 0.001$ where $\bar{X}$ is the mean value of $X$

        -   When $X$ is discrete, $h = 1$

    3.  Calculate $\hat{Y}_{i1}$ which is the prediction in the probability scale using new $X$ and other variables' observed vales.

    4.  Calculate the difference between the two predictions as fraction of $h$: $\frac{\bar{Y}_{i1} - \bar{Y}_{i0}}{h}$
3.  Average numerical derivative is $E[\frac{\bar{Y}_{i1} - \bar{Y}_{i0}}{h}] \approx \frac{\partial p (Y|\mathbf{X}, \beta)}{\partial X}$

Two-sided derivatives

1.  Estimate your model
2.  For each observation $i$
    1.  Calculate $\hat{Y}_{i0}$ which is the prediction in the probability scale using observed values

    2.  Increase $X$ (variable of interest) by a "small" amount $h$ ($X_{1} = X + h$) and decrease $X$ (variable of interest) by a "small" amount $h$ ($X_{2} = X - h$)

        -   When $X$ is continuous, $h = (|\bar{X}| + 0.001) \times 0.001$ where $\bar{X}$ is the mean value of $X$

        -   When $X$ is discrete, $h = 1$

    3.  Calculate $\hat{Y}_{i1}$ which is the prediction in the probability scale using new $X_1$ and other variables' observed vales.

    4.  Calculate $\hat{Y}_{i2}$ which is the prediction in the probability scale using new $X_2$ and other variables' observed vales.

    5.  Calculate the difference between the two predictions as fraction of $h$: $\frac{\bar{Y}_{i1} - \bar{Y}_{i2}}{2h}$
3.  Average numerical derivative is $E[\frac{\bar{Y}_{i1} - \bar{Y}_{i2}}{2h}] \approx \frac{\partial p (Y|\mathbf{X}, \beta)}{\partial X}$

```{r}
library(margins)
library(tidyverse)

data("mtcars")
mod <- lm(mpg ~ cyl * disp * hp, data = mtcars)
margins::margins(mod) %>% summary()

# function for variable
get_mae <- function(mod, var = "disp") {
    data = mod$model
    
    pred <- predict(mod)
    
    if (class(mod$model[[{
        {
            var
        }
    }]]) == "numeric") {
        h = (abs(mean(mod$model[[var]])) + 0.01) * 0.01
    } else {
        h = 1
    }
    
    data[[{
        {
            var
        }
    }]] <- data[[{
{
var
}
}]] - h

    pred_new <- predict(mod, newdata = data)

    mean(pred - pred_new) / h
}

get_mae(mod, var = "disp")
```

## Packages

### 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](https://stackoverflow.com/questions/72463092/estimate-marginal-effect-in-triple-interaction) 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

```{r}
library(marginaleffects)
library(tidyverse)
data(mtcars)

mod <- lm(mpg ~ hp * wt * am, data = mtcars)
predictions(mod) %>% head()
# for specific reference values
predictions(mod, newdata = datagrid(am = 0, wt = c(2, 4)))
plot_cap(mod, condition = c("hp","wt"))

```

```{r}
# Average Margianl Effects
mfx <- marginaleffects(mod, variables = c("hp", "wt"))
summary(mfx)

# Group-Average Marginal Effects
marginaleffects::marginaleffects(mod, by = "hp", variables = "am")

# Marginal effects at representative values
marginaleffects::marginaleffects(mod, 
                                 newdata = datagrid(am = 0, 
                                                    wt = c(2, 4)))

# Marginal Effects at the Mean
marginaleffects::marginaleffects(mod, newdata = "mean")
```

```{r}
# counterfactual
comparisons(mod, variables = list(am = 0:1)) %>% summary()
```

### 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).

```{r}
library(margins)

# examples by the package's authors
mod <- lm(mpg ~ cyl * hp + wt, data = mtcars)
summary(mod)
```

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

```{r}
summary(margins(mod))

# equivalently 
margins_summary(mod)

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**)

```{r}
margins(mod, at = list(hp = 150))

margins(mod, at = list(hp = 150)) %>% summary()
```

### mfx

Works well with [Generalized Linear Models]/`glm` package

For technical details, see the package [vignette](https://cran.rstudio.com/web/packages/mfx/vignettes/mfxarticle.pdf)

| Model             | Dependent Variable | Syntax       |
|-------------------|--------------------|--------------|
| Probit            | Binary             | `probitmfx`  |
| Logit             | Binary             | `logitmfx`   |
| Poisson           | Count              | `poissonmfx` |
| Negative Binomial | Count              | `negbinmfx`  |
| Beta              | Rate               | `betamfx`    |

```{r}
library(mfx)
data("mtcars")
poissonmfx(formula = vs ~ mpg * cyl * disp, data = mtcars)
```

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.