Very interesting - Estimating linear probability models -- what's wrong with this? (R)

Question

I am trying to estimate a model for an event modelled by probability of happening which is a linear function of x (distributed normally) plus an error term, u.

Then I simulate whether the event really happened for each X comparing the probability of it happening against a uniformly distributed random variable.

So, I wrote a little function that simulates this model for a given b0, b1, X (mean and sd.) and error term (mean = 0 and sd.):

SAMPLE_SIZE = 10000

underlying <- function(b0, b1, mean_x, sd_x, sd_u) {
  xs <- rnorm(SAMPLE_SIZE, mean_x, sd_x)
  us <- rnorm(SAMPLE_SIZE, 0.0, sd_u)
  ys <- b0 + (b1 * xs) + us
  ws <- runif(SAMPLE_SIZE) < ys  
  list(ws = ws, ys = ys, xs = xs, us = us)
}

It neatly returns both the probability of the event taking place in the ys component plus a simulation on the ws component.

I then tested whether I can correctly estimate b0 and b1 using linear regressions. And I got a very weird result.

This is how I simulated the samples and did the regressions:

b1s <- seq(from = 0.0, to = 1.0, length.out = 100)

datasets <- lapply(b1s, FUN = function(x) underlying(0.5, x, 1.0, 0.2, 0.05)) 
regs     <- lapply(datasets,  FUN = function(x) lm(data = x, ws ~ xs))
b0s_hat = sapply(regs, function(x) x$coefficients[[1]])
    b1s_hat = sapply(regs, function(x) x$coefficients[[2]])

So, for different b1s (and b0 = 0.5) I can plot the estimated b0 and b1 against the real b1:

plot(b1s, b0s_hat)
plot(b1s, b1s_hat)

And what we get for b1s_hat looks sigmoid-ish like a cumulative distribution function, and b0s_hat looks like a bell curve (like the density function).

I thought I could recover the coefficients using the linear regression. What exactly is smelling weird here?

Answer 1

Well, you get the usual problem that the linear probability model does not always lead to responses that are valid probabilities. Or in other words: The ys your underlying() function produces can be outside of 0 or 1. And then you will introduce a bias in the estimated intercept and slope. This is what leads to the weird patterns in your plot.

As a simple illustration consider the case where b1 = 1.0:

R> set.seed(0)
R> d <- as.data.frame(underlying(0.5, 1.0, 1.0, 0.2, 0.05))
R> summary(d)
##      ws                ys              xs               us            
##  Mode :logical   Min.   :0.735   Min.   :0.2594   Min.   :-0.1907833  
##  FALSE:3         1st Qu.:1.362   1st Qu.:0.8663   1st Qu.:-0.0338383  
##  TRUE :9997      Median :1.503   Median :1.0020   Median : 0.0002358  
##  NA's :0         Mean   :1.502   Mean   :1.0021   Mean   :-0.0002308  
##                  3rd Qu.:1.640   3rd Qu.:1.1378   3rd Qu.: 0.0338984  
##                  Max.   :2.246   Max.   :1.7136   Max.   : 0.1669680  

Note that moste ys are above 1 and hence only 3 ws are FALSE. There is no hope to consistently estimate intercept and slope from this.

For this (and other) reason(s), it is preferred to use a proper binary response model (logit or probit GLM etc.) that assures that all probabilities are actually in (0, 1).