Agfdhyk

I read here, here, here, and elsewhere that linear probability model (LPM) might be used to get risk differences when the outcome variable is binomial.

LPM has some advantages such as ease of interpretation by simplifying the estimation of risk differences, which in certain fields might be preferable than odds ratio that is usually provided by logistic regression.

My concerns are however that "[u]sing the LPM one has to live with the following three drawbacks:

1. The effect ΔP(y=1∣X=x0+Δx) is always constant




2. The error term is by definition heteroscedastic




3. OLS does not bound the predicted probability in the unit interval"

Therefore, I would appreciate any idea on a better regression model for binomial data to get robust risk difference in R while avoiding these drawbacks of LPM.

The first "drawback" you mention is the definition of the risk difference, so there is no avoiding this.

There is at least one way to obtain the risk difference using the logistic regression model. It is the average marginal effects approach. The formula depends on whether the predictor of interest is binary or continuous. I will focus on the case of the continuous predictor.

Imagine the following logistic regression model:

$$lnbigg[frac{hatpi}{1-hatpi}bigg] = hat{y}^* = hatgamma_c times x_c + Zhatbeta$$

where $Z$ is an $n$ cases by $k$ predictors matrix including the constant, $hatbeta$ are $k$ regression weights for the $k$ predictors, $x_c$ is the continuous predictor whose effect is of interest and $hatgamma_c$ is its estimated coefficient on the log-odds scale.

Then the average marginal effect is:

$$mathrm{RD}_c = hatgamma_c times frac{1}{n}Bigg(sumfrac{e^{hat{y}^*}}{big(1 + e^{hat{y}^*}big)^2}Bigg)\$$

This is the average PDF scaled by the weight of $x_c$. It turns out that this effect is very well approximated by the regression weight from OLS applied to the problem regardless of drawbacks 2 and 3. This is the simplest justification in practice for the application of OLS to estimating the linear probability model.

For drawback 2, as mentioned in one of your citations, we can manage it using heteroskedasticity-consistent standard errors.

Now, Horrace and Oaxaca (2003) have done some very interesting work on consistent estimators for the linear probability model. To explain their work, it is useful to lay out the conditions under which the linear probability model is the true data generating process for a binary response variable. We begin with:

begin{align}
begin{split}
P(y = 1 mid X)
{}& = P(Xbeta + epsilon > t mid X) quad text{using a latent variable formulation for } y \
{}& = P(epsilon > t-Xbeta mid X)
end{split}
end{align}

where $y in {0, 1}$, $t$ is some threshold above which the latent variable is observed as 1, $X$ is matrix of $n$ cases by $k$ predictors, and $beta$ their weights. If we assume $epsilonsimmathcal{U}(-0.5, 0.5)$ and $t=0.5$, then:

begin{align}
begin{split}
P(y = 1 mid X)
{}& = P(epsilon > 0.5-Xbeta mid X) \
{}& = P(epsilon < Xbeta -0.5 mid X) quad text{since $mathcal{U}(-0.5, 0.5)$ is symmetric about 0} \
{}&=begin{cases}
0, & mathrm{if} Xbeta -0.5 < -0.5\
frac{(Xbeta -0.5)-(-0.5)}{0.5-(-0.5)}, & mathrm{if} Xbeta -0.5 in [-0.5, 0.5)\
1, & mathrm{if} Xbeta -0.5 geq 0.5
end{cases} quad text{CDF of $mathcal{U}(-0.5,0.5)$}\
{}&=begin{cases}
0, & mathrm{if} Xbeta < 0\
Xbeta, & mathrm{if} Xbeta in [0, 1)\
1, & mathrm{if} Xbeta geq 1
end{cases}
end{split}
end{align}

So the relationship between $Xbeta$ and $P(y = 1mid X)$ is only linear when $Xbeta in [0, 1]$, otherwise it is not. Horrace and Oaxaca suggested that we may use $Xhatbeta$ as a proxy for $Xbeta$ and in empirical situations, if we assume a linear probability model, we should consider it inadequate if there are any predicted values outside the unit interval.

As a solution, they recommended the following steps:

1. Estimate the model using OLS


2. Check for any fitted values outside the unit interval. If there are none, stop, you have your model.


3. Drop all cases with fitted values outside the unit interval and return to step 1

Using a simple simulation (and in my own more extensive simulations), they found this approach to recover adequately $beta$ when the linear probability model is true. They termed the approach sequential least squares (SLS). SLS is similar in spirit to doing MLE and censoring the mean of the normal distribution at 0 and 1 within each iteration of estimation, see Wacholder (1986).

Now how about if the logistic regression model is true? I will demonstrate in a simulated data example what happens using R:

# An implementation of SLS

s.ols <- function(fit.ols) {

  dat.ols <- model.frame(fit.ols)

  n.org <- nrow(dat.ols)

  fitted <- fit.ols$fitted.values

  form <- formula(fit.ols)

  while (any(fitted > 1 | fitted < 0)) {

    dat.ols <- dat.ols[!(fitted > 1 | fitted < 0), ]

    m.ols <- lm(form, dat.ols)

    fitted <- m.ols$fitted.values

  }

  m.ols <- lm(form, dat.ols)

  # Bound predicted values at 0 and 1 using complete data

  m.ols$fitted.values <- punif(as.numeric(model.matrix(fit.ols) %*% coef(m.ols)))

  m.ols

}



set.seed(12345)

n <- 20000

dat <- data.frame(x = rnorm(n))

# With an intercept of 2, this will be a high probability outcome

dat$y <- ((2 + 2 * dat$x + + rlogis(n)) > 0) + 0

coef(fit.logit <- glm(y ~ x, binomial, dat))

# (Intercept)           x 

#    2.042820    2.021912 

coef(fit.ols <- lm(y ~ x, dat))

# (Intercept)           x 

#   0.7797852   0.2237350 

coef(fit.sls <- s.ols(fit.ols))

# (Intercept)           x 

#   0.8989707   0.3932077 

We see that the RD from OLS is .22 and that from SLS is .39. We can also compute the average marginal effect from the logistic regression equation:

coef(fit.logit)["x"] * mean(dlogis(predict(fit.logit)))

#        x 

# 0.224426 

We can see that the OLS estimate is very close to this value.

How about we plot the different effects to better understand what they try to capture:

library(ggplot2)



dat.res <- data.frame(

  x = dat$x, logit = fitted(fit.logit),

  ols = fitted(fit.ols), sls = fitted(fit.sls))

dat.res <- tidyr::gather(dat.res, model, fitted, logit:sls)

ggplot(dat.res, aes(x, fitted, col = model)) +

  geom_line() + theme_bw()

From here, we see that the OLS results looks nothing like the logistic curve. OLS captures the average change in probability of y across the range of x (the average marginal effect). While SLS results in the linear approximation to the logistic curve in the region it is changing on the probability scale.

In this scenario, I think the SLS estimate better reflects the reality of the situation.

As with OLS, heteroskedasticity is implied by SLS, so Horrace and Oaxaca recommend heteroskedasticity-consistent standard errors.

• Horrace, W. C., & Oaxaca, R. L. (2003, January 1). New Wine in Old Bottles: A Sequential Estimation Technique for the Lpm. Retrieved from https://papers.ssrn.com/sol3/papers.cfm?abstract_id=383102
  


• Wacholder, S. (1986). Binomial regression in GLIM: estimating risk ratios and risk differences. American Journal of Epidemiology, 123(1), 174–184. https://doi.org/10.1093/oxfordjournals.aje.a114212
  

1. Every model has this problem. For example, logistic regression implies the constant log odd ratio.




2. For binomial distribution, the variance is $p(1-p)$ for one trial. So the different predict value of p implies the different variance. But in the model fitting process this problem is resolve by WLS (weighted least square).




3. For $hat p = Xhat beta$, it is possible for some $X$, the $hat p$ can go lower than 0, or higher than 1, especially when the model is used to predict the probability using the $X$s that is not in the dataset used to build the model.