A Reminder about Taylor Series

The Taylor expansion is a useful tool because it allows us to restate a differentiable function, G(x), in terms of (an infinite sum of) the derivatives of G(x). To be more precise, an infinitely differentiable G(x) evaluated at a can be written as

$$G(x) = G(a) + \frac{G'(a)}{1!}(x - a) + \frac{G''(a)}{2!}(x-a)^2 + \frac{G'''(a)}{3!}(x-a)^3 + \dots$$

If we cut off the expansion after some number of terms (two is common), we can get a useful approximation of G(x).

Taylor Series and the Delta Method

In the case of predicted margins (levels), where the regression model has link function link², the column vector of predicted levels P_m at covariates X₁ is

P_m(X₁β) = link⁻¹(X₁β)

The predicted effects P_e of that same regression is a function of Xβ such that

$$P_e(X_1 \beta) = \frac{d(\text{link}^{-1}(X_1 \beta))}{d(X \beta)}$$

Depending on whether the effect is over a continuous or categorical variable, this may be an actual derivative (the instantaneous rate of change) or the subtraction of P(Xβ) calculated at one value of X from another (the first difference).

Using the Taylor expansion, we can approximate P, an arbitrary function of the random variable Xβ around the point X₁β, as³

$$P(X \beta) = P(X_1 \beta) + \frac{d(P(X_1 \beta))}{d(X\beta)}(X\beta - X_1 \beta)$$

Now we can substitute for the different margins we’ll care about. For predicted levels, the Taylor expansion is

$$P_m(X \beta) = \text{link}^{-1}(X_1 \beta) + \frac{d(\text{link}^{-1}(X_1 \beta))}{d(X \beta)}(X\beta - X_1 \beta)$$

For the predicted effects of categorical variables, we are trying to estimate the effects at P_e(X₁β − X₂β), which gives us

$$\begin{aligned} P_e(X \beta) &= \text{link}^{-1} (X_1 \beta - X_2 \beta) + \frac{d(\text{link}^{-1}(X_1 \beta - X_2 \beta))}{d(X \beta)}(X\beta - (X_1 \beta - X_2 \beta)) \\ % &= \text{link}^{-1} (X_1 \beta) - \text{link}^{-1} (X_2 \beta) + \frac{d(\text{link}^{-1})}{d(X \beta)}(X_1) - \frac{d(\text{link}^{-1})}{d(X \beta)}(X_2) \end{aligned} $$

For the predicted effects of continuous variables, the marginal effect is a derivative, so

$$\begin{aligned} P_e(X_1 \beta) &= \frac{d(\text{link}^{-1}(X_1 \beta))}{d(X \beta)} + \frac{d^2(\text{link}^{-1}(X_1 \beta))}{d(X \beta)}\left(\frac{d(\text{link}^{-1}(X \beta))}{d(X \beta)} - \frac{d(\text{link}^{-1}(X_1 \beta))}{d(X \beta)}\right) \end{aligned} $$

Point estimates

To get an estimate of the average levels of the outcome for different values of treatment, we can set treatment to 0 or 1 for the entire dataset, generate predicted outcomes for the dataset, and average the results (in other words, what would the outcome have been, if everybody was in the control / treatment group?). The linear predictors of our model (in this case, setting τ = t), are given by: ŷ_i = α + β₁ ⋅ t + β₂ ⋅ age_i + β₃ ⋅ (t ⋅ age_i) also known as ŷ_i = X_iβ. We can then transform the linear predictors to the scale of the outcome variable (predicted probabilities) with the inverse link function: $$f(z) = \frac{1}{1 + \exp(-z)}$$ Finally, we can find the average predicted level, by taking the average of these predicted probabilities. Combining these, we have: $$ P(X\beta) = \frac{1}{n} \sum_{i = 1}^n \frac{1}{1 + \exp(-X_i\beta)} $$ Now let’s do this in R: set treatment to 0 or 1, generate linear predictors, transform them to predicted outcomes, and take the mean.

# Extract the n x k matrix of data
x <- model.matrix(lg)

# Extract the coefficients from the model (a column vector with k entries)
beta <- matrix(lg$coefficients)

# CONTROL:

# Set treatment and treatment:age to 0 for all observations
x[, "treatment"] <- 0
x[, 'treatment:age'] <- x[, 'treatment'] * x[, 'age']

# Get linear predictors
pred_ctl <- x %*% beta

# Apply the inverse link function to get predicted probabilities
pp_ctl <- 1 / (1 + exp(-pred_ctl))

# Get the average predicted probability
mean_pp_ctl <- mean(pp_ctl)

# TREATMENT:

# Set treatment to 1 and treatment:age to age for all observations
x[, "treatment"] <- 1
x[, 'treatment:age'] <- x[, 'treatment'] * x[, 'age']

# Get linear predictors
pred_treat <- x %*% beta

# Apply the inverse link function to get predicted probabilities
pp_treat <- 1 / (1 + exp(-pred_treat))

# Get the average predicted probability
mean_pp_treat <- mean(pp_treat)

# RESULTS:

mean_pp_ctl

## [1] 0.1126685

mean_pp_treat

## [1] 0.208375

Variance

So far this is pretty straightforward, but we want to know the standard error of these estimates.

You might be tempted to just use the R command, predict(model, newdata, se.fit = TRUE), fix treatment to 0 or 1, and find the average of the resulting standard errors. Don’t do that! You’ll average the standard errors of the predicted outcome of each observation in your dataset - which is not the same thing as finding the standard error for a specific estimate, which is what we’re looking for.

As explained above, we can approximate the variance of the predicted margins: using J to represent the jacobian (the derivative of our transformation with respect to each β parameter), then

Var(P(Xβ)) = JVJ^T

Going back to our function from before, let’s start with the coefficient on treatment, β₁. First we just apply the quotient rule (and the chain rule in the last step):

Coming up for air for a second: what’s the last term in that final expression?

X_iβ = α + β₁τ_i + β₂age_i + β₃(τ_i ⋅ age_i), so the derivative of that with respect to β₁ is just τ_i. Therefore,

which is the first term of our jacobian.⁹

What about the other terms? Well, because the X_iβ term is additively separable, each term is the same, except for having a different term at the end. For example, the derivative with respect to β₂ is going to be

$$\frac{1}{n} \sum_{i = 1}^n \frac{\exp(-X_i\beta))}{(1 + \exp(-X_i\beta))^2} \cdot \text{age}_i$$.

So the full jacobian is $$ J = \left[ \begin{array}{l} \\ \frac{1}{n} \sum_{i = 1}^n \frac{\exp(-X_i\beta)}{(1 + \exp(-X_i\beta))^2} \cdot 1\\ \frac{1}{n} \sum_{i = 1}^n \frac{\exp(-X_i\beta)}{(1 + \exp(-X_i\beta))^2} \cdot \tau_i\\ \frac{1}{n} \sum_{i = 1}^n \frac{\exp(-X_i\beta)}{(1 + \exp(-X_i\beta))^2} \cdot \text{age}_i\\ \frac{1}{n} \sum_{i = 1}^n \frac{\exp(-X_i\beta)}{(1 + \exp(-X_i\beta))^2} \cdot \tau_i \cdot \text{age}_i \end{array} \right] $$

This can be rewritten as

$$ J = \frac{1}{n} \left[ \begin{array} \\ \frac{\exp(-X_1\beta)}{(1 + \exp(-X_1\beta))^2} & \frac{\exp(-X_2\beta)}{(1 + \exp(-X_2\beta))^2} & \cdots & \frac{\exp(-X_n\beta)}{(1 + \exp(-X_n\beta))^2} \end{array} \right] \left[ \begin{array}{cccc}\\ 1 & \tau_1 & \text{age}_1 & \tau_1 \cdot \text{age}_1 \\ 1 & \tau_2 & \text{age}_2 & \tau_2 \cdot \text{age}_2 \\ \vdots & \vdots & \vdots & \vdots\\ 1 & \tau_n & \text{age}_n & \tau_n \cdot \text{age}_n \\ \end{array} \right] $$

The first term applies the derivative of the inverse link function to every linear predictor in the model, and the second term is the covariate matrix from our data! How convenient! That’s going to be true for all general linear models. So we can rewrite this as

$$J = \frac{\left[\begin{array}\\ \frac{\exp(-X_1\beta)}{(1 + \exp(-X_1\beta))^2} & \frac{\exp(-X_2\beta)}{(1 + \exp(-X_2\beta))^2} & \cdots & \frac{\exp(-X_n\beta)}{(1 + \exp(-X_n\beta))^2}\end{array}\right] X}{n}$$

Recalling that our variance is JVJ^T, and that we’ve already calculated X_iβ for all i, this is relatively simple to do in R:

# Get the data
x <- model.matrix(lg)

# CONTROL ERROR

# Apply the derivative of the inverse link function to the linear predictors
deriv <- as.vector(exp(-pred_ctl) / (1 + exp(-pred_ctl))^2)

# Set treatment to 0
x[, 'treatment'] <- 0
x[, 'treatment:age'] <- x[, 'treatment'] * x[, 'age']

# Complete the chain rule by matrix-multiplying the derivatives by the data,
# now we have the jacobian
j <- deriv %*% x / nrow(x)

# The variance of our estimate is the cross product of the jacobian and the model's
# variance-covariance matrix
variance <- j %*% vcov(lg) %*% t(j)

# The error is the square root of that
se_ctl <- sqrt(diag(variance))

# TREATMENT ERROR: same logic

deriv <- as.vector(exp(-pred_treat) / (1 + exp(-pred_treat))^2)

x <- model.matrix(lg)

x[, 'treatment'] <- 1
x[, 'treatment:age'] <- x[, 'treatment'] * x[, 'age']

j <- deriv %*% x / nrow(x)

variance <- j %*% vcov(lg) %*% t(j)
se_treat <- sqrt(diag(variance))

se_ctl

## [1] 0.009334252

se_treat

## [1] 0.009089667

Review

OK, we’ve got our estimates and variance:

result <- data.frame(
  Label = c("treatment = 0", "treatment = 1"),
  Margin = c(mean_pp_ctl, mean_pp_treat),
  Standard.Error = c(se_ctl, se_treat)
)

result

##           Label    Margin Standard.Error
## 1 treatment = 0 0.1126685    0.009334252
## 2 treatment = 1 0.2083750    0.009089667

What do we get from modmarg?

marg <- modmarg::marg(mod = lg, var_interest = 'treatment')
marg[[1]][, c("Label", "Margin", "Standard.Error")]

##           Label    Margin Standard.Error
## 1 treatment = 0 0.1126685    0.009334252
## 2 treatment = 1 0.2083750    0.009089667

Hooray!