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

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?

• Are you trying to recreate logistic regression? Dec 23, 2015 at 23:50
• I am trying to estimate a linear model. I think that way I can get more interpretability from the model (the coefficients would give out the marginal effect of each independent variable on the probability). That's why I avoided using logistic regression models (and non-linear models in general). Dec 24, 2015 at 1:30
• Don't think of logistic regression as a "non-linear model in general". The predictors still enter the model linearly. The response variable is simply a transformation of the probabilities in which you are interested, in a way that avoids the weirdness that you observed, and which is quite interpretable once you think about it a while.
– EdM
Dec 24, 2015 at 2:52
• With a bounded response it is not possible to have constant marginal effects. If you use a model that imposes constant marginal effects, you can easily get predictions that are not valid. Dec 24, 2015 at 18:55
• A good place to start learning about binary response models would be Agresti's Introduction to Categorical Data Analysis which discusses logistic regression and its interpretations at great length, and in a very accessible way.
– Sycorax
Jan 4, 2016 at 12:40

## 1 Answer

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).

• that does make sense, I was expecting more robustness from the linear probability model. I know that some kinds of bias can be studied on the multi-variable linear regression model; is there an expression that shows the bias introduced in the linear probability model when the underlying model gives probabilities outside of the unit interval? Or could you point to any papers or other reading on that? Thanks huge lots. Dec 24, 2015 at 1:25
• I like the treatment of the topic in the book “Analysis of micro data“ by Winkelmann & Boes, see in particular Chapters 2 and 4. Dec 24, 2015 at 18:57