Week 3: Generalized Linear Models

Previously

• Likelihood: a function of the parameters \(\theta\) evaluated at the observed data

  \[\mathcal{L}(\theta| \mathbf{y}) = f(\mathbf{y}| \theta) \qquad \ell_n(\theta) = \sum_{i=1}^n \log f(y_i|\theta)\]

• Maximum likelihood estimator: \(\hat{\theta} = \argmax_{\theta \in \Theta} \ \ell_n(\theta)\)

  • Consistent: \(\hat{\theta} \overset{p}{\to} \theta_0\)
  • Asymptotically normal: \(\sqrt{n}(\hat{\theta} - \theta_0) \overset{d}{\to} \text{Normal}\big(0, \mathcal{I}(\theta_0)^{-1}\big)\)
  • Asymptotically efficient: attains the Cramér-Rao lower bound

This week

• More properties of the MLE
  • Invariance and the delta method
• How does numerical optimization actually work?
  • Newton-Raphson, Fisher scoring, gradient ascent
• Generalized linear models
  • Extending the classic linear form to functions of the CEF
  • \(g(\mathbb{E}[Y_i|X_i]) = X_i^{\prime}\beta\)
• Binary outcome models
  • Logit and probit; deriving and estimating the likelihood
  • Interpretation: log-odds, marginal effects, and the delta method
• An application: propensity scores and inverse probability weighting

More on the MLE

Invariance

• We often want to do inference on some function of \(\theta\)
  • (e.g.) we estimate \(\beta\) but want a predicted probability, a risk ratio, or a marginal effect.
• The MLE has a convenient property: invariance.
  • If \(\hat{\theta}\) is the MLE of \(\theta\), then \(g(\hat{\theta})\) is the MLE of \(g(\theta)\) – for any function \(g\).
• Combined with consistency, \(g(\hat{\theta}) \overset{p}{\to} g(\theta_0)\) for continuous \(g\) (continuous mapping theorem).
• But how do we get the standard errors (asymptotically)?
  • The delta method

Delta method

• We know the asymptotic variance of \(\hat{\theta}\)

  • It’s the inverse Fisher information, estimated by the inverse negative Hessian at the MLE.

• To get \(\text{Var}(g(\hat{\theta}))\), start with a first-order Taylor approximation around the true value \(\theta_0\)

  \[g(\hat{\theta}) \approx g(\theta_0) + \big[\nabla g(\theta_0)^{\prime}\big](\hat{\theta} - \theta_0)\]

• Take the variance

  \[\text{Var}\big(g(\hat{\theta})\big) \approx \text{Var}\bigg(g(\theta_0) + \big[\nabla g(\theta_0)^{\prime}\big]\hat{\theta} - \big[\nabla g(\theta_0)^{\prime}\big]\theta_0\bigg)\]

• \(g(\theta_0)\) and \(\nabla g(\theta_0)\) are constants, so the variance of the first and third terms is zero

  \[\text{Var}\big(g(\hat{\theta})\big) \approx \text{Var}\bigg(\big[\nabla g(\theta_0)^{\prime}\big]\hat{\theta}\bigg)\]

Delta method

• Pull out and “square” the constant. We get an expression in terms of the original variance-covariance matrix

  \[\text{Var}\big(g(\hat{\theta})\big) \approx \big[\nabla g(\theta_0)^{\prime}\big]\text{Var}(\hat{\theta})\big[\nabla g(\theta_0)\big]\]

• And so our asymptotic approximation is

  \[\sqrt{n}\big(g(\hat{\theta}) - g(\theta_0)\big) \overset{d}{\to} \text{Normal}\bigg(0, \ \big[\nabla g(\theta_0)^{\prime}\big]\mathcal{I}(\theta_0)^{-1}\big[\nabla g(\theta_0)\big]\bigg)\]

• We can use a plug-in estimator at the MLE: evaluate the gradient at \(\hat{\theta}\) and use our existing estimator \(\widehat{\text{Var}}(\hat{\theta})\).

How optimizers work

Finding the maximum

• How does maxLik() actually find a solution?

• Our estimator is defined as

  \[\hat{\theta} = \argmax_{\theta \in \Theta} \ \ell_n(\theta)\]

• For well-behaved problems the log-likelihood is concave, so finding the maximum amounts to finding the unique solution to \(S_n(\theta) = 0\).

• Closed form solutions are rare (the Normal linear model is a special case).

  • Generally, we’ll use an iterative algorithm to find this solution.

Newton-Raphson

• We start at some initial guess \(\hat{\theta}^{(0)}\)

• Let \(\hat{\theta}^{(t)}\) denote the “current” value and \(\hat{\theta}^{(t+1)}\) our update – we keep iterating until convergence.

• Our goal is to solve for a zero of \(S_n(\theta)\). Take a first-order Taylor approximation around our current guess

  \[S_n(\theta) \approx S_n(\hat{\theta}^{(t)}) + \nabla S_n(\hat{\theta}^{(t)}) \left(\theta - \hat{\theta}^{(t)}\right)\]

• What’s \(\nabla S_n(\hat{\theta}^{(t)})\)?

  • It’s the Jacobian of the gradient… or the matrix of second-order partial derivatives of the log-likelihood… or the Hessian \(H_n(\hat{\theta}^{(t)})\)!

Newton-Raphson

• Our next value of \(\hat{\theta}\) is the value that sets this approximation to zero

  \[0 = S_n(\hat{\theta}^{(t)}) + H_n(\hat{\theta}^{(t)}) \left(\hat{\theta}^{(t+1)} - \hat{\theta}^{(t)}\right)\]

• Multiply through by the inverse Hessian and rearrange

  \[\hat{\theta}^{(t+1)} = \hat{\theta}^{(t)} - H_n^{-1}(\hat{\theta}^{(t)}) S_n(\hat{\theta}^{(t)})\]

• Recall that the negative Hessian at \(\hat\theta\) is the observed Fisher information \(\mathcal{J}_n(\hat\theta) = -H_n(\hat\theta)\)

  • So the update is the score scaled by the inverse negative hessian.

• An alternative algorithm, Fisher scoring, substitutes the expected Fisher information \(\mathcal{I}_n\) for the observed one.

  • Identical for GLMs with the canonical link – we’ll see why later.

Gradient ascent

• Compare the Newton-Raphson update to the classic “gradient descent/gradient ascent” algorithm

  \[\hat{\theta}^{(t+1)} = \hat{\theta}^{(t)} + \gamma_t S_n(\hat{\theta}^{(t)})\]

• Gradient ascent chooses some sequence of \(\gamma_t\) (the “learning rate” or “step size”) at iteration \(t\)

  • The gradient tells us which direction to move in and \(\gamma_t\) tells us how far to move.

• Contrast with Newton-Raphson, where “how far to move” comes from computing and inverting the Hessian

  • More information per step, but each step is computationally expensive (inverting a matrix is computationally hard)

• Trade-off: Newton-Raphson converges quickly but each step is costly; gradient ascent steps are cheap but need more of them.

  • Also typically need to tune the learning rate \(\gamma_t\) - too large and we overshoot the true optimum

Illustration

• Let’s illustrate with a simple one-parameter likelihood.

  \[y_i \underset{\text{i.i.d}}{\sim} \text{Bernoulli}(\pi), \qquad \pi = \frac{1}{1 + \exp(-\theta)}\]

  • One parameter \(\theta\), the log-odds of success.

• The average log-likelihood, its score, and its Hessian all have closed forms

  \[\bar{\ell}_n(\theta) = \bar{y}\,\theta - \log\big(1 + e^{\theta}\big), \qquad \bar{S}_n(\theta) = \bar{y} - \pi, \qquad \bar{H}_n(\theta) = -\pi(1-\pi)\]

• The Newton-Raphson step can be written in closed form

  \[\hat{\theta}^{(t+1)} = \hat{\theta}^{(t)} + \frac{\bar{y} - \pi^{(t)}}{\pi^{(t)}\big(1 - \pi^{(t)}\big)}\]

• And the true MLE is also known in closed form: \(\hat{\theta} = \log\frac{\bar{y}}{1-\bar{y}}\)

Coding the two algorithms

• Simulate some data and write down the three components

set.seed(60637)
n <- 100
y <- rbinom(n, 1, 1/(1 + exp(-0.8)))   # true theta = 0.8
ybar <- mean(y)

loglik <- function(th) ybar*th - log(1 + exp(th))      # average log-likelihood
score  <- function(th) ybar - 1/(1 + exp(-th))         # average score
hess   <- function(th){ p <- 1/(1 + exp(-th)); -p*(1-p) }  # average Hessian

theta_hat <- log(ybar/(1 - ybar))                      # the analytic MLE
theta_hat

[1] 0.895

Coding the two algorithms

newton <- function(th, iter){                 # step = score / curvature
  path <- th
  for (t in 1:iter){ th <- th - score(th)/hess(th); path <- c(path, th) }
  path
}
grad_ascent <- function(th, gamma, iter){     # step = score * fixed learning rate
  path <- th
  for (t in 1:iter){ th <- th + gamma*score(th); path <- c(path, th) }
  path
}

Coding the two algorithms

• Start both at \(\hat{\theta}^{(0)} = -1\)

nr_path <- newton(-1, iter = 4)
ga_path <- grad_ascent(-1, gamma = 1, iter = 12)
round(nr_path, 6)

[1] -1.000  1.243  0.863  0.895  0.895

round(ga_path, 6)

 [1] -1.0000 -0.5589 -0.2127  0.0503  0.2477  0.3961  0.5083  0.5939  0.6597
[10]  0.7105  0.7500  0.7808  0.8049

• Newton-Raphson is at the MLE (0.895) by the fourth iteration.
• Gradient ascent has taken twelve steps and still isn’t there.

Newton-Raphson

Gradient ascent

Intro to GLMs

Motivation: Propensity Scores

• Researchers wanting to estimate causal effects from observational designs often use a weighting estimator to account for non-random treatment assignment.
  • Observe treatment \(D_i\), confounders \(X_i\), outcome \(Y_i\)
  • Need to estimate the propensity score \(\pi_i = \Pr(D_i = 1 | X_i)\) to construct “inverse propensity of treatment weights”
• Key problem: with many continuous covariates, it’s hard to estimate \(\Pr(D_i = 1 | X_i)\) non-parametrically
  • One solution: assume a parametric model for \(D_i\).
• Propensity scores are predictions from the model.

The LaLonde dataset

• LaLonde (1986) asked: can selection-on-observables designs recover the answer from a randomized experiment?
  • The National Supported Work Demonstration (NSW) randomly assigned disadvantaged workers to a job training program.
  • Outcome: real earnings in 1978 (re78).
  • Because assignment was random, the experiment gives a credible benchmark: roughly $1,794 for the subsample we’ll use (Dehejia and Wahba 1999).
• LaLonde created an observational dataset by replacing the experimental controls with a survey comparison group drawn from the Panel Study of Income Dynamics (PSID).
  • Now treatment is not randomly assigned – the comparison group is much older, better educated, and far richer at baseline.
  • If our adjustment methods work, they should recover something near $1,794 anyway.
• LaLonde (1986) found typical OLS regression worked poorly
  • Dehejia and Wahba (1999) tackle this with propensity score methods (specifically matching, but we’ll do weighting).

The LaLonde dataset

• The dataset is in a lot of R packages, such as causalsens

library(causalsens)
data(lalonde.psid)
lalonde <- lalonde.psid   # shorter name for the slides
glimpse(lalonde)

Rows: 2,675
Columns: 12
$ age       <dbl> 37, 22, 30, 27, 33, 22, 23, 32, 22, 33, 19, 21, 18, 27, 17, …
$ education <dbl> 11, 9, 12, 11, 8, 9, 12, 11, 16, 12, 9, 13, 8, 10, 7, 10, 13…
$ black     <dbl> 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ hispanic  <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ married   <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, …
$ nodegree  <dbl> 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, …
$ re74      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ re75      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ re78      <dbl> 9930, 3596, 24910, 7506, 290, 4056, 0, 8472, 2164, 12418, 81…
$ u74       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ u75       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ treat     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …

• \(n = 2{,}675\): 185 NSW treated units and 2,490 PSID comparison units.

The experimental benchmark

• The experimental controls are not in that dataset – LaLonde replaced them. But they ship alongside it.

data(lalonde.exp)   # 185 treated + 260 *randomized* controls
lalonde.exp %>% group_by(treat) %>%
  summarize(n = n(), age = mean(age), education = mean(education),
            married = mean(married), re74 = mean(re74))

# A tibble: 2 × 6
  treat     n   age education married  re74
  <int> <int> <dbl>     <dbl>   <dbl> <dbl>
1     0   260  25.1      10.1   0.154 2107.
2     1   185  25.8      10.3   0.189 2096.

• Randomization did its job – the two arms look alike on every pre-treatment covariate.

exp_att <- lm_robust(re78 ~ treat, data = lalonde.exp, se_type = "HC2")
coef(summary(exp_att))["treat",]

  Estimate Std. Error    t value   Pr(>|t|)   CI Lower   CI Upper         DF 
  1.79e+03   6.71e+02   2.67e+00   7.77e-03   4.76e+02   3.11e+03   4.43e+02 

• Training raised 1978 earnings by about $1,794. That is the number our adjustment methods have to recover.
  • The 185 treated units are the same people in both datasets; only the comparison group changes.
  • Note the benchmark is itself an estimate, with a 95% CI of roughly \([\$476,\ \$3{,}113]\).

The naive comparison

• What happens if we just compare treated and control earnings directly?

lalonde %>% group_by(treat) %>% 
  summarize(n = n(), mean_re78 = mean(re78), mean_re74 = mean(re74), 
            mean_age = mean(age), mean_married = mean(married))

# A tibble: 2 × 6
  treat     n mean_re78 mean_re74 mean_age mean_married
  <dbl> <int>     <dbl>     <dbl>    <dbl>        <dbl>
1     0  2490    21554.    19429.     34.9        0.866
2     1   185     6349.     2096.     25.8        0.189

• Complete ignorability is implausible – the PSID men are nine years older, nearly two years better educated, four times as likely to be married, and earned nine times as much in 1974.

The naive comparison

• Regressing 1978 earnings on treatment alone:

naive_fit <- lm_robust(re78 ~ treat, data=lalonde, se_type="HC2")
summary(naive_fit)

Call:
lm_robust(formula = re78 ~ treat, data = lalonde, se_type = "HC2")

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper   DF
(Intercept)    21554        312    69.1  0.00e+00    20943    22165 2673
treat         -15205        657   -23.1 3.93e-108   -16493   -13916 2673

Multiple R-squared:  0.0609 ,   Adjusted R-squared:  0.0606 
F-statistic:  535 on 1 and 2673 DF,  p-value: <2e-16

• The training program appears to lower earnings by about $15,205.

Generalized Linear Models

• Generalized linear models (GLMs) (Nelder and Wedderburn 1972) have three components:

  1. A parametric distribution for \(Y_i|X_i\) (“stochastic component”)

  2. A linear predictor: \(\eta_i = X_i^{\prime}\beta = \beta_0 + \beta_1X_{i1} + \beta_2X_{i2} + \dotsc + \beta_kX_{ik}\) (“systematic component”)

  3. A link function \(g()\) applied to the CEF \(\mathbb{E}[Y_i|X_i]\) that yields the linear predictor

     \[g\big(\mathbb{E}[Y_i | X_i]\big) = \eta_i\]

• Alternatively, we can write the CEF in terms of the “inverse link” \(g^{-1}()\) applied to the linear predictor

  $$\mathbb{E}[Y_i | X_i] = g^{-1}(\eta_i)$$

Exponential Family

• The types of probability distributions permitted for \(Y_i\) are quite general: the exponential family

  • This contains the Normal as well as many other common distributions – the Bernoulli, Poisson, exponential, gamma, and more.

• Exponential family distributions have density functions of the form

  \[f(Y|\theta) = h(Y) \exp\bigg\{b(\theta) \cdot T(Y) - A(\theta)\bigg\}\]

• \(h(Y)\), \(b(\theta)\), \(T(Y)\) and \(A(\theta)\) are known functions

  • \(T(Y)\) is the sufficient statistic, \(b(\theta)\) the canonical (natural) parameter, \(A(\theta)\) the log-normalizer.

• The key intuition: exponential family distributions factorize in a convenient way

  • The parameter and the data interact only through the product \(b(\theta) \cdot T(Y)\).
  • If \(b(\theta) = \theta\), the distribution is in “canonical” form.

Example: Bernoulli

• Consider the Bernoulli PMF

  \[f(Y_i|\pi_i) = \pi_i^{Y_i}(1-\pi_i)^{1-Y_i}\]

• Take the log, then the exponent

  \[f(Y_i|\pi_i) = \exp\bigg\{\log\bigg[\pi_i^{Y_i}(1-\pi_i)^{1-Y_i}\bigg]\bigg\}\]

• Properties of logs

  \[f(Y_i|\pi_i) = \exp\bigg\{Y_i\log(\pi_i) + (1-Y_i)\log(1-\pi_i)\bigg\}\]

Example: Bernoulli

• Rearranging and using properties of logs again

  \[f(Y_i|\pi_i) = \exp\bigg\{Y_i\log\bigg(\frac{\pi_i}{1-\pi_i}\bigg) + \log(1-\pi_i)\bigg\}\]

• So our exponential form is

  • \(h(Y) = 1\)
  • \(T(Y) = Y_i\)
  • \(b(\pi_i) = \log\bigg(\frac{\pi_i}{1-\pi_i}\bigg)\)
  • \(A(\pi_i) = -\log(1-\pi_i)\)
    • remember the exponential form carries a minus sign in front of \(A\)

• Critically, this is where we get a good link function

  • The canonical parameter is \(\log\bigg(\frac{\pi}{1-\pi}\bigg)\) – the log-odds.
  • The canonical link is the function that equates this parameter with the linear predictor: \(\log\bigg(\frac{\pi_i}{1-\pi_i}\bigg) = X_i^{\prime}\beta\)

Binary outcome models

Logistic regression

• The “logit” or “logistic” GLM models the log-odds of a binary outcome as a function of the linear predictor \(X_i^{\prime}\beta\)

  \[Y_i \underset{\text{i.i.d}}{\sim} \text{Bernoulli}(\pi_i)\]

  \[\mathbb{E}[Y_i | X_i] = \Pr(Y_i = 1 | X_i) = \pi_i\]

  \[\log\bigg(\frac{\pi_i}{1-\pi_i}\bigg) = X_i^{\prime}\beta\]

• Alternatively, this is written in terms of the “inverse link” function (the logistic function) relating \(\pi_i\) to \(g^{-1}(X_i^{\prime}\beta)\)

  \[\pi_i = \frac{\exp(X_i^{\prime}\beta)}{1 + \exp(X_i^{\prime}\beta)} = \frac{1}{1 + \exp(-X_i^{\prime}\beta)}\]

Inverse-link functions

• The logistic function maps inputs on \(\mathbb{R}\) to \((0, 1)\)

Inverse-link functions

• Another link is the “probit”, whose inverse link is the Normal CDF \(\Phi()\)

Inverse-link functions

• When do we use probit vs. logit?
  • Computational convenience – the logit is the canonical link, so it’s better behaved.
  • Probit has some good properties for Bayesian inference (we’ll come back to this).
• Can you tell the difference?

Deriving the logit likelihood

• Start from the Bernoulli log-likelihood

  \[\ell_n(\beta) = \sum_{i=1}^n Y_i\log(\pi_i) + (1-Y_i)\log(1 - \pi_i)\]

• Collecting the \(Y_i\) terms gives us the log-odds, which is just \(X_i^{\prime}\beta\)

  \[\ell_n(\beta) = \sum_{i=1}^n Y_i \log\bigg(\frac{\pi_i}{1-\pi_i}\bigg) + \log(1 - \pi_i)\]

• And from the inverse link,

  \[1 - \pi_i = 1 - \frac{1}{1 + \exp(-X_i^{\prime}\beta)} = \frac{\exp(-X_i^{\prime}\beta)}{1 + \exp(-X_i^{\prime}\beta)} = \frac{1}{1 + \exp(X_i^{\prime}\beta)}\]

• Substituting in

  \[\ell_n(\beta) = \sum_{i=1}^n \bigg[Y_i X_i^{\prime}\beta - \log\bigg(1 + \exp(X_i^{\prime}\beta)\bigg)\bigg]\]

Example: Estimating the propensity score

• We want to predict assignment to training, treat, from the pre-treatment covariates
  • age - Age in years
  • education - Years of schooling
  • black, hispanic - Race/ethnicity indicators
  • married - Married; binary indicator
  • nodegree - No high school degree; binary indicator
  • re74 - Real earnings in 1974 (pre-treatment)
  • re75 - Real earnings in 1975 (pre-treatment)

Example: Estimating the propensity score

• Make the design matrix \(\mathbf{X}\)

X_mat <- model.matrix(treat ~ age + education + black + hispanic + married + nodegree + re74 + re75, data=lalonde)
head(X_mat) # View the top of the matrix

  (Intercept) age education black hispanic married nodegree re74 re75
1           1  37        11     1        0       1        1    0    0
2           1  22         9     0        1       0        1    0    0
3           1  30        12     1        0       0        0    0    0
4           1  27        11     1        0       0        1    0    0
5           1  33         8     1        0       0        1    0    0
6           1  22         9     1        0       0        1    0    0

• How many observations?

n_obs <- nrow(X_mat)
n_obs

[1] 2675

Example: Estimating the propensity score

• Put the log-likelihood into code (this returns a vector of the log-likelihood for each observation)

• We derived the score analytically, so let’s supply it too – one row per observation

Example: Estimating the propensity score

• Now let’s optimize it to get the MLE

• What did we get?

est <- coef(logit_mle) # Our optimization routine
names(est) <- colnames(X_mat)
est

(Intercept)         age   education       black    hispanic     married 
   0.637053   -0.078996    0.039806    2.200219    2.008432   -1.627318 
   nodegree        re74        re75 
   0.774187   -0.000109   -0.000261 

Example: Estimating the propensity score

• Compare to the built-in R routine

logit_Rglm <- glm(treat ~ age + education + black + hispanic + married + nodegree + re74 + re75, 
                  data=lalonde, family=binomial(link="logit"))
coef(logit_Rglm)

(Intercept)         age   education       black    hispanic     married 
   0.637053   -0.078996    0.039806    2.200219    2.008432   -1.627318 
   nodegree        re74        re75 
   0.774187   -0.000109   -0.000261 

• Let’s obtain our (asymptotic) variance-covariance matrix – the inverse of the observed information

logit_vcov <- solve(-maxLik::hessian(logit_mle))

• Square root of the diagonal gives our SEs

logit_SEs <- sqrt(diag(logit_vcov))

Example: Estimating the propensity score

• Let’s get our test statistics and p-values, using the usual asymptotic normal approximation

results <- rbind(coef(logit_mle), logit_SEs, 
                     coef(logit_mle)/logit_SEs, 
                     2*pnorm(-abs(coef(logit_mle)/logit_SEs)))
colnames(results) <- colnames(X_mat)
rownames(results) <- c("Estimate", "Std. Error", "Test statistic", "p-value")
round(results, 4)

               (Intercept)     age education black hispanic married nodegree
Estimate             0.637 -0.0790    0.0398 2.200   2.0084  -1.627   0.7742
Std. Error           1.180  0.0148    0.0697 0.318   0.5592   0.263   0.3529
Test statistic       0.540 -5.3453    0.5709 6.915   3.5914  -6.195   2.1941
p-value              0.589  0.0000    0.5681 0.000   0.0003   0.000   0.0282
                  re74    re75
Estimate       -0.0001 -0.0003
Std. Error      0.0000  0.0000
Test statistic -3.7472 -6.4661
p-value         0.0002  0.0000

Interpreting logit coefficients

• How do we interpret the \(\beta\) parameters substantively?

  \[\log\bigg(\frac{\pi_i}{1-\pi_i}\bigg) = \beta_0 + \beta_1X_{i1} + \beta_2X_{i2} + \dotsc + \beta_k X_{ik}\]

• Take the partial derivative w.r.t. \(X_{ik}\)

  \[\frac{\partial}{\partial X_{ik}} \log\bigg(\frac{\pi_i}{1-\pi_i}\bigg) = \beta_k\]

• So \(\beta_k\) captures the change in the log-odds for a one-unit change in \(X_{k}\)

  • Descriptively, it’s the difference in log-odds between two observations that differ in \(X_{k}\) by one unit.
  • Exponentiating, \(\exp(\beta_k)\) is an odds ratio. This is not a risk ratio, and it is not a difference in probabilities.

Interpreting logit coefficients

• On the “log-odds” scale, the change due to an increase in \(X_{k}\) does not depend on the values of the other \(X\) variables (unless we explicitly specify an interaction).

  • But thinking on the log-odds scale is hard! We think in terms of probabilities.
  • This additivity does not hold when we take \(\frac{\partial}{\partial X_{i1}} \pi_i\)

• Working through the chain rule,

  \[\frac{\partial \pi_i}{\partial X_{ik}} = \beta_k \, \pi_i(1 - \pi_i)\]

  • The effect on the probability scale depends on \(\pi_i\), and therefore on every covariate.

• Logit models implicitly encode interactions with respect to the CEF \(\mathbb{E}[Y_i | X_i]\).

Interpreting logit coefficients

• Remember: A one-unit change in the linear predictor corresponds to different changes in probability depending on your baseline.

Inference on transformations

Transformed quantities

• We have \(\hat{\beta}\), but we actually want the propensity score \(\hat{\pi}_i = \Pr(D_i = 1 | X_i)\)

  • Just apply the inverse link to get the quantity we want

  \[\hat{\pi}_i = \frac{1}{1 + \exp(-X_i^{\prime}\hat{\beta})}\]

• By invariance, this is the MLE of \(\pi_i\), and it’s consistent for the true propensity score (under our modeling assumptions).

• But if we want a confidence interval for \(\hat{\pi}_i\), we use the delta method.

Example: a predicted propensity

• Let’s get the propensity of receiving training for a typical NSW participant – the median covariate profile among the treated

(Intercept)         age   education       black    hispanic     married 
          1          25          11           1           0           0 
   nodegree        re74        re75 
          1           0           0 

• Construct our prediction function

      [,1]
[1,] 0.888

Example: a predicted propensity

• Let’s do the gradient numerically (we could calculate it analytically though!)

• Applying the delta method

• Making our 95% asymptotic CI for \(\hat{\pi}(X)\)

[1] 0.838 0.939

marginaleffects

• The marginaleffects R package provides a straightforward implementation of the delta method for a lot of transformed quantities
  • (e.g.) the marginal change in the probability of \(D_i=1\) for a one-unit change in \(X\)
• avg_predictions() will calculate point estimates + CIs for the expected outcome just like we did above.
• avg_comparisons() will do the implied “marginal effects” for a unit change in \(X\) given a set baseline.
• We have to tell it which covariate profile to evaluate at. We’ll use the treated units – the NSW participants are the group we care about.

library(marginaleffects)
treated <- subset(lalonde, treat == 1)

grid_median <- datagrid(model = logit_Rglm, newdata = treated,
                        FUN_numeric = median, FUN_integer = median, FUN_binary = median)
grid_mean   <- datagrid(model = logit_Rglm, newdata = treated,
                        FUN_numeric = mean,   FUN_integer = mean,   FUN_binary = mean)

marginaleffects

• The median profile is the one we just built by hand; the mean profile is the other common convention

rbind(median = unlist(grid_median[, -1]), mean = unlist(grid_mean[, -1]))

        age black education hispanic married nodegree re74 re75
median 25.0 1.000      11.0   0.0000   0.000    1.000    0    0
mean   25.8 0.843      10.3   0.0595   0.189    0.708 2096 1532

• Note the mean profile is not a person anyone could be – it is 84% Black and 71% without a degree.
  • The median profile at least corresponds to a possible respondent.

marginaleffects: predicted propensity

avg_predictions(logit_Rglm, newdata = grid_median)   # matches our hand calculation

 Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
    0.888     0.0258 34.4   <0.001 859.9 0.838  0.939

Type: response

avg_predictions(logit_Rglm, newdata = grid_mean)

 Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
    0.645     0.0342 18.8   <0.001 260.6 0.578  0.712

Type: response

• Evaluating at the mean moves us from \(\hat{\pi} \approx 0.89\) down to \(\approx 0.64\) – the profile is different, so the point on the curve is different.

marginaleffects: marginal effect of married

avg_comparisons(logit_Rglm, variables = "married", newdata = grid_median)

 Estimate Std. Error     z Pr(>|z|)    S  2.5 % 97.5 %
   -0.278     0.0612 -4.55   <0.001 17.5 -0.398 -0.158

Term: married
Type: response
Comparison: 1 - 0

avg_comparisons(logit_Rglm, variables = "married", newdata = grid_mean)

 Estimate Std. Error     z Pr(>|z|)    S  2.5 % 97.5 %
   -0.385     0.0566 -6.81   <0.001 36.5 -0.496 -0.274

Term: married
Type: response
Comparison: 1 - 0

• Being married lowers the probability of NSW participation by 28 points at the median profile, 39 at the mean profile.
  • Same \(\hat{\beta}\), same data – only the baseline \(\pi_i\) differs. This is the \(\beta_k\pi_i(1-\pi_i)\) scaling at work.

Inverse propensity of treatment weighting

• Now we have \(\hat{\pi}_i\) for every unit.

• The control group is unrepresentative of the treated group.

  • So we reweight it. Upweight controls who look like treated units; downweight the ones who don’t.

• Under conditional ignorability (\(\{Y_i(1), Y_i(0)\} \perp D_i \mid X_i\)) and positivity (\(0 < \pi_i < 1\)), the IPTW estimator identifies a causal contrast (Rosenbaum and Rubin 1983).

• For the average treatment effect on the treated (ATT), the weights are

  \[w_i = D_i + (1 - D_i)\frac{\hat{\pi}_i}{1 - \hat{\pi}_i}\]

  • Treated units get weight 1 – we keep the treated sample as-is and re-weight controls to match.
  • Controls get the odds of treatment: a control who looked very likely to be treated counts for more.

Inverse propensity of treatment weighting

• The ATT estimator is then a difference of weighted means

  \[\widehat{\text{ATT}} = \frac{1}{n_1}\sum_{i: D_i = 1} Y_i \ - \ \frac{\sum_{i: D_i = 0} w_i Y_i}{\sum_{i: D_i = 0} w_i}\]

• Dividing by the sum of weights rather than \(n_0\) gives the Hájek (stabilized) estimator

• For the ATE instead, we would weight every unit by the inverse probability of the treatment it actually received

  \[w_i^{\text{ATE}} = \frac{D_i}{\hat{\pi}_i} + \frac{1 - D_i}{1 - \hat{\pi}_i}\]

• We’ll target the ATT, since that’s what the experimental benchmark of $1,794 estimates.

Checking overlap first

• Before weighting anything, look at the estimated propensity scores by treatment status

lalonde$pscore <- as.vector(pred_prob(coef(logit_mle), X_mat))

• A large mass of PSID controls has \(\hat\pi \approx 0\), not close to the treated group.

Checking overlap first

• Effectively, weighting discards most of the PSID sample.

lalonde$w_att <- with(lalonde, ifelse(treat == 1, 1, pscore/(1 - pscore)))
c(max_control_weight = max(lalonde$w_att[lalonde$treat==0]),
  sum_control_weights = sum(lalonde$w_att[lalonde$treat==0]))

 max_control_weight sum_control_weights 
               10.8               184.2 

Did the weights fix the imbalance?

• The point of the weights is balance. After weighting, the covariate distributions should look alike.

• The standard diagnostic is the standardized mean difference

  \[\text{SMD}_k = \frac{\bar{X}_{k,\text{treated}} - \bar{X}_{k,\text{control}}^{\,w}}{s_{k,\text{treated}}}\]

  • Scaling by the treated-group SD puts every covariate on a common, unit-free scale.
  • Common heuristic for observational studies (from biostats/epi): \(|\text{SMD}| < 0.1\) is “balanced.”

• The cobalt package does this.

Did the weights fix the imbalance?

library(cobalt)
bt <- bal.tab(treat ~ age + education + black + hispanic + married + nodegree + re74 + re75,
              data = lalonde, weights = lalonde$w_att,
              method = "weighting", estimand = "ATT",
              un = TRUE, binary = "std", thresholds = c(m = .1))
bt$Balance

             Type Diff.Un Diff.Adj        M.Threshold
age       Contin.  -1.263  0.09045     Balanced, <0.1
education Contin.  -0.881  0.00646     Balanced, <0.1
black      Binary   1.630 -0.02417     Balanced, <0.1
hispanic   Binary   0.114 -0.02169     Balanced, <0.1
married    Binary  -1.729  0.07538     Balanced, <0.1
nodegree   Binary   0.886 -0.04782     Balanced, <0.1
re74      Contin.  -3.547 -0.10173 Not Balanced, >0.1
re75      Contin.  -5.446 -0.06331     Balanced, <0.1

Did the weights fix the imbalance?

• And the effective sample sizes:

bt$Observations

           Control Treated
Unadjusted  2490.0     185
Adjusted      42.6     185

• The unweighted imbalances are enormous – re74 starts at \(-3.5\) standard deviations, married at \(-1.7\).
  • After weighting nearly everything is inside or roughly at the \(0.1\) band.
• Note the effective sample size: the 2,490 controls are worth about 43 after weighting.
  • This is what our histogram shows - weighting puts a large amount of weight on a small number of units.
  • Low effective sample size = higher variance estimator.
• Don’t forget that we by definition can’t check balance on things we didn’t control for.
  • Balance on the covariates we included says nothing about the confounders we omitted.

The Love plot

• Common to present this in a Love plot

love.plot(treat ~ age + education + black + hispanic + married + nodegree + re74 + re75,
          data = lalonde, weights = lalonde$w_att,
          method = "weighting", estimand = "ATT",
          binary = "std", thresholds = c(m = .1),
          var.names = c(age = "Age", education = "Years of education", married = "Married",
                        nodegree = "No HS degree", re74 = "Earnings 1974",
                        re75 = "Earnings 1975", black = "Black", hispanic = "Hispanic"),
          sample.names = c(Unadjusted = "Unweighted", Adjusted = "IPTW (ATT)"),
          colors = c("indianred", "dodgerblue"), shapes = c("circle", "triangle"),
          title = "Covariate balance before and after weighting")

• Red circles are the unweighted sample, blue triangles the weighted one. Dashed lines mark \(\pm 0.1\).

The IPTW estimate

• Now compute the weighted difference in means

att_hajek <- with(lalonde,
  mean(re78[treat == 1]) - weighted.mean(re78[treat == 0], w_att[treat == 0]))
att_hajek

[1] 1759

The IPTW estimate

• Equivalently, as a weighted least squares regression

summary(lm_robust(re78 ~ treat, data = lalonde, weights = w_att, se_type = "HC2"))

Call:
lm_robust(formula = re78 ~ treat, data = lalonde, weights = w_att, 
    se_type = "HC2")

Weighted, Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper   DF
(Intercept)     4590        724    6.34 2.75e-10   3169.8     6011 2673
treat           1759        927    1.90 5.79e-02    -58.9     3577 2673

Multiple R-squared:  0.0144 ,   Adjusted R-squared:  0.014 
F-statistic:  3.6 on 1 and 2673 DF,  p-value: 0.0579

The IPTW estimate

• Compare to the alternatives:

naive_est <- with(lalonde, mean(re78[treat==1]) - mean(re78[treat==0]))
reg_fit <- lm_robust(re78 ~ treat + age + education + black + hispanic + married + nodegree + re74 + re75,
                     data=lalonde, se_type="HC2")
reg_est <- coef(reg_fit)["treat"]
data.frame(estimator = c("Naive difference", "Regression adjustment", "IPTW (ATT)", "Experimental benchmark"),
           estimate  = c(naive_est, as.numeric(reg_est), att_hajek, coef(exp_att)["treat"]))

               estimator estimate
1       Naive difference   -15205
2  Regression adjustment      752
3             IPTW (ATT)     1759
4 Experimental benchmark     1794

Estimating the standard errors

• For the standard errors, we want to account for estimation uncertainty in the propensity scores.
• Two options:
  • Bootstrap the whole procedure – refit the logit inside each bootstrap iteration.
  • Stack the score functions of both stages and apply the sandwich estimator
• We’ll implement the bootstrap

att_fn <- function(dd){
  Xb <- model.matrix(treat ~ age + education + black + hispanic + married + nodegree + re74 + re75, data=dd)
  fit <- suppressWarnings(glm.fit(Xb, dd$treat, family = binomial()))
  ps <- as.vector(1/(1 + exp(-Xb %*% fit$coefficients)))
  w  <- ifelse(dd$treat == 1, 1, ps/(1 - ps))
  mean(dd$re78[dd$treat == 1]) - weighted.mean(dd$re78[dd$treat == 0], w[dd$treat == 0])
}
set.seed(53703)
boot_att <- replicate(2000, att_fn(lalonde[sample(nrow(lalonde), replace = TRUE),]))
c(bootstrap_SE = sd(boot_att), quantile(boot_att, c(.025, .975)))

bootstrap_SE         2.5%        97.5% 
         966         -170         3599 

Results

• Collect the four estimators, each with its own standard error

grab <- function(f) coef(summary(f))["treat", c("Estimate", "Std. Error", "CI Lower", "CI Upper")]
results_tab <- data.frame(
  estimator = c("Naive difference", "Regression adjustment", "IPTW (ATT)", "Experimental benchmark"),
  rbind(grab(naive_fit), grab(reg_fit),
        c(att_hajek, sd(boot_att), unname(quantile(boot_att, c(.025, .975)))),
        grab(exp_att)),
  se = c("HC2", "HC2", "bootstrap", "HC2"),
  row.names = NULL)
names(results_tab) <- c("estimator", "estimate", "std_error", "ci_lower", "ci_upper", "se_type")

Results

results_tab

               estimator estimate std_error ci_lower ci_upper   se_type
1       Naive difference   -15205       657   -16493   -13916       HC2
2  Regression adjustment      752       785     -788     2292       HC2
3             IPTW (ATT)     1759       966     -170     3599 bootstrap
4 Experimental benchmark     1794       671      476     3113       HC2

• Adjustment moved us an enormous distance: \(-\$15{,}205 \rightarrow +\$1{,}759\), against a benchmark of \(\$1{,}794\).
  • Regression adjustment alone only reaches about \(+\$752\)
• But the bootstrap interval is wide and includes zero.
  • With the comparison group weighted down to an effective 43 units, there simply isn’t much information left.
• IPTW estimators are extremely noisy when propensities are close to \(0\) or \(1\).
  • Doubly-robust approaches (e.g. AIPW) combine the weights with an outcome model, which helps reduce variance (Kang and Schafer 2007).

Next week

• Latent variables and discrete choice
  • The latent-variable representation of logit and probit
  • Ordinal logit and proportional odds; multinomial logit and IIA
• Random utility and conditional logit
  • Modelling the attributes of the choices
  • Analyzing “conjoint” experiments with model-based approaches.
• Duration models
  • The survival function, the hazard, and censoring
  • Parametric (exponential, Weibull), semi-parametric (Cox), and non-parametric (Kaplan-Meier) approaches
• Readings: Chapter 15 of Gelman, Hill, and Vehtari (2020); Clark et al. (2003); Stensrud and Hernán (2020)