Musings on Logistic Regression, CBPS, and Overlap Weights

Logistic regression via M-estimation

First, we specify the estimating equations. For logistic regression, we create the function s_bar_lr(), which takes in a vector of coefficient estimates, the treatment, and the covariates, and returns the value of the estimating equations. The first step is to compute \(p_i\) as \(\text{expit}(\mathbf{X}_i\beta)\) (in R, expit is called plogis()). Then, we supply that to the estimating equation and compute the means.

s_bar_lr <- function(b, A, X) {
  p <- plogis(drop(X %*% b))
  
  colMeans((A - p) * X)
}

When we find the values of b that make this vector equal to 0, we have found the maximum likelihood estimates of logistic regression coefficients. Let’s take a look at this in action. We’ll use rootSolve::multiroot() to find the roots. This takes in a function, some starting values, and other variables you need to supply to the function, and finds the values of the desired parameters that make the function return a vector of 0s. We’ll use the lalonde dataset as an example.

library("cobalt")

##  cobalt (Version 4.5.4.9001, Build Date: 2024-03-13)

data("lalonde")

X <- model.matrix(~age + educ + married + nodegree, data = lalonde)
A <- lalonde$treat

out <- rootSolve::multiroot(s_bar_lr,
                            start = rep(0, ncol(X)),
                            A = A, X = X)

setNames(out$root, colnames(X))

## (Intercept)         age        educ     married    nodegree 
## -2.54468907  0.01024968  0.12644332 -1.52238592  0.98034779

And we’ll do the same using logistic regression as implemented in glm() just to show we get the same estimates:

fit <- glm(treat ~ age + educ + married + nodegree, data = lalonde,
           family = binomial("logit"))
coef(fit)

## (Intercept)         age        educ     married    nodegree 
## -2.54468907  0.01024968  0.12644332 -1.52238592  0.98034779

all.equal(coef(fit), out$root, check.attributes = FALSE)

## [1] TRUE

And we can see that when we supply those values to the estimating equations, we get a vector of 0s:

s_bar_lr(out$root, A, X) |> round(7)

## (Intercept)         age        educ     married    nodegree 
##           0           0           0           0           0

Now let’s compute ATO weights from the propensity scores estimated from this model and compute the weighted mean differences (we’ll use col_w_smd() from cobalt, which computes weighted standardized mean differences column-wise for a matrix of covariates; we’ll omit the intercept here as well):

p <- plogis(drop(X %*% out$root))

w_ATO <- A * (1 - p) + (1 - A) * p

col_w_smd(X[,-1], A, w_ATO) |> round(7)

##      age     educ  married nodegree 
##        0        0        0        0

The ATO weights exactly balance the covariate means, as promised. If we compute the ATE weights, though, we find balance is good but not perfect.

w_ATE <- A / p + (1 - A) / (1 - p)

col_w_smd(X[,-1], A, w_ATE)

##          age         educ      married     nodegree 
## -0.066562991  0.059888054 -0.028829997 -0.004447331

CBPS via M-estimation

Now let’s use M-estimation to fit the CBPS model for the propensity scores. For this, we’ll use a different function that corresponds to the weighted difference in means for ATE weights. I’ll call this one s_bar_cb().

s_bar_cb <- function(b, A, X) {
  p <- plogis(drop(X %*% b))
  
  colMeans((A/p - (1-A)/(1-p)) * X)
}

We’ll again estimate the coefficients using M-estimation.

out <- rootSolve::multiroot(s_bar_cb,
                            start = rep(0, ncol(X)),
                            A = A, X = X)

setNames(out$root, colnames(X))

##  (Intercept)          age         educ      married     nodegree 
## -2.907727232  0.004142096  0.168217307 -1.535394932  1.110766225

The coefficients aren’t too different from the usual logistic regression coefficients. Let’s compare them to the coefficients estimated from CBPS::CBPS():

fit_cb <- CBPS::CBPS(treat ~ age + educ + married + nodegree,
                     data = lalonde, ATT = 0, method = "exact")

all.equal(out$root, drop(coef(fit_cb)), check.attributes = FALSE,
          tolerance = 1e-4)

## [1] TRUE

Looks like we got them right³! Let’s compute balance using the estimated ATE weights:

p <- plogis(drop(X %*% out$root))

w_ATE_cb <- A / p + (1 - A) / (1 - p)

col_w_smd(X[,-1], A, w_ATE_cb) |> round(7)

##      age     educ  married nodegree 
##        0        0        0        0

It turns out it’s pretty simple to make CBPS weights for the ATT or ATC, too; just apply the ATT or ATC weights formula to the estimating equations and then compute the weights from the resulting coefficients. I’ll demonstrate this without much commentary for the ATT below, noting that ATT weights are \(w_i^{ATT}=A + (1- A)\frac{p_i}{1-p_i}=p_i w_i^{ATE}\).

s_bar_cb_att <- function(b, A, X) {
  p <- plogis(drop(X %*% b))
  
  colMeans((A - (1 - A) * p / (1 - p)) * X)
}

out <- rootSolve::multiroot(s_bar_cb_att,
                            start = rep(0, ncol(X)),
                            A = A, X = X)

p <- plogis(drop(X %*% out$root))

w_ATT_cb <- A + (1 - A) * p / (1 - p)

col_w_smd(X[,-1], A, w_ATT_cb) |> round(7)

##      age     educ  married nodegree 
##        0        0        0        0

Overlap weights from canonical link GLM

I told you earlier that all generalized linear models with a canonical link have the same estimating equations as logistic regression, and I also told you that the estimating equations for logistic regression imply exact balance on the covariate means when using ATO weights. It turns out that this means if you use any generalized linear model with a canonical link, you also get exact mean balance when using ATO weights computed from its predicted values. I’ll demonstrate that below with Poisson and linear regression. This time I’ll skip the M-estimation and just use glm() to compute the propensity scores.

# Poisson PS
p <- glm(treat ~ age + educ + married + nodegree, data = lalonde,
         family = poisson("log")) |> fitted()

w_ATO_pois <- A * (1 - p) + (1 - A) * p

col_w_smd(X[,-1], A, w_ATO_pois) |> round(7)

##      age     educ  married nodegree 
##        0        0        0        0

# Linear PS
p <- glm(treat ~ age + educ + married + nodegree, data = lalonde,
         family = gaussian("identity")) |> fitted()

w_ATO_lin <- A * (1 - p) + (1 - A) * p

col_w_smd(X[,-1], A, w_ATO_lin) |> round(7)

##      age     educ  married nodegree 
##        0        0        0        0

Okay, but before I end, I want to show you something cool and maybe unexpected about the linear regression ATO weights. This was inspired by a finding in Hazlett and Shinkre (2024). Consider estimating a treatment effect using a linear regression of the outcome on the treatment and covariates but with no interaction between them. It turns out the coefficient on treatment is exactly equal to the weighted difference in means computed using the linear regression propensity score ATO weights. See below:

fit <- lm(re78 ~ treat + age + educ + married + nodegree,
          data = lalonde)
coef(fit)["treat"]

##    treat 
## 207.3495

# Computing weighted difference in means; not standardized
col_w_smd(lalonde$re78, A, w_ATO_lin, std = FALSE)

## [1] 207.3495

That’s pretty crazy! That means we can interpret this linear regression estimate as an overlap-weighted difference in means, where the propensity scores used in the overlap weights are fitted using linear regression (i.e., the linear probability model). Of course, the linear probability model is a pretty bad propensity score model since it can yield propensity scores greater than 1 and less than 0, which yield negative and possibly extreme ATO weights, unlike the logistic regression-based ATO weights, which are always positive and never extreme.