Understanding selection bias

Suppose we have data {(xi,yi):i∈{1,2,…,n}}$\left\{\right(x_i,y_i):i\in \{1,2,\dots ,n\left\}\right\}$ generated by the process yi=βxi+ui,$$y_i=\beta x_i+u_i,$$ where the ui$u_i$ are random errors with zero means, equal variances, and zero correlations with the xi$x_i$. This data generating process (DGP) satisfies the Gauss-Markov assumptions, so we can obtain an unbiased estimate ˆβ$\hat{\beta }$ of the coefficient β$\beta$ using ordinary least squares (OLS).

Now suppose we restrict our data to observations with xi≥0$x_i\ge 0$ or yi≥0$y_i\ge 0$. How will these restrictions change ˆβ$\hat{\beta }$?

To investigate, let’s create some toy data:

library(dplyr)

n <- 100
set.seed(0)
df <- tibble(x = rnorm(n), u = rnorm(n), y = x + u)

Here xi$x_i$ and ui$u_i$ are standard normal random variables, and yi=xi+ui$y_i=x_i+u_i$ for each observation i∈{1,2,…,100}$i\in \{1,2,\dots ,100\}$. Thus β=1$\beta =1$. The OLS estimate of β$\beta$ is ˆβ=Cov(x,y)Var(x),$$\hat{\beta }=\frac{\mathrm{Cov}(x,y)}{\mathrm{Var}\left(x\right)},$$ where x=(x1,x2,…,x100)$x=(x_1,x_2,\dots ,x_{100})$ and y=(y1,y2,…,y100)$y=(y_1,y_2,\dots ,y_{100})$ are data vectors, Cov$\mathrm{Cov}$ is the covariance operator, and Var$\mathrm{Var}$ is the variance operator. For these data, we have

cov(df$x, df$y) / var(df$x)

## [1] 1.138795

as our estimate with no selection.

Next, let’s introduce our selection criteria:

df <- df %>%
  tidyr::crossing(criterion = c('x >= 0', 'y >= 0')) %>%
  rowwise() %>%  # eval is annoying to vectorise
  mutate(selected = eval(parse(text = criterion))) %>%
  ungroup()

df

## # A tibble: 200 x 5
##        x       u      y criterion selected
##    <dbl>   <dbl>  <dbl> <chr>     <lgl>   
##  1 -2.22 -0.0125 -2.24  x >= 0    FALSE   
##  2 -2.22 -0.0125 -2.24  y >= 0    FALSE   
##  3 -1.56 -1.12   -2.68  x >= 0    FALSE   
##  4 -1.56 -1.12   -2.68  y >= 0    FALSE   
##  5 -1.54  0.577  -0.963 x >= 0    FALSE   
##  6 -1.54  0.577  -0.963 y >= 0    FALSE   
##  7 -1.44 -1.39   -2.83  x >= 0    FALSE   
##  8 -1.44 -1.39   -2.83  y >= 0    FALSE   
##  9 -1.43 -0.543  -1.97  x >= 0    FALSE   
## 10 -1.43 -0.543  -1.97  y >= 0    FALSE   
## # … with 190 more rows

Now df contains two copies of each observation—one for each selection criterion—and an indicator for whether the observation is selected by each criterion. We can use df to estimate OLS coefficients and their standard errors among observations with xi≥0$x_i\ge 0$ and yi≥0$y_i\ge 0$:

df %>%
  filter(selected) %>%
  group_by(criterion) %>%
  summarise(n = n(),
            estimate = cov(x, y) / var(x),
            std.error = sd(y - estimate * x) / sqrt(n))

## # A tibble: 2 x 4
##   criterion     n estimate std.error
##   <chr>     <int>    <dbl>     <dbl>
## 1 x >= 0       48    1.02      0.136
## 2 y >= 0       47    0.356     0.110

The OLS estimate among observations with xi≥0$x_i\ge 0$ approximates the true value β=1$\beta =1$ well. However, the estimate among observations with yi≥0$y_i\ge 0$ is much smaller than one. We can confirm this visually:

What’s going on? Why do we get biased OLS estimates of β$\beta$ among observations with yi≥0$y_i\ge 0$ but not among observations with xi≥0$x_i\ge 0$?

The key is to think about the errors ui$u_i$ in each case. Since the xi$x_i$ and ui$u_i$ are independent, selecting observations with xi≥0$x_i\ge 0$ leaves the distributions of the ui$u_i$ unchanged—they still have zero means, equal variances, and zero correlations with the xi$x_i$. Thus, the Gauss-Markov assumptions still hold and we still obtain unbiased OLS estimates of β$\beta$.

In contrast, the xi$x_i$ and ui$u_i$ are negatively correlated among observations with yi≥0$y_i\ge 0$. To see why, notice that if yi=xi+ui$y_i=x_i+u_i$ then yi≥0$y_i\ge 0$ if and only if xi≥−ui$x_i\ge -u_i$. So if xi$x_i$ is low then ui$u_i$ must be high (and vice versa) for the observation to be selected. Thus, among selected observations, we have ui=ρxi+εi,$$u_i=\rho x_i+{\epsilon }_i,$$ where ρ<0$\rho <0$ indexes (and, in this case, equals) the correlation between the xi$x_i$ and ui$u_i$, and where the residuals εi${\epsilon }_i$ are uncorrelated with the xi$x_i$. Our DGP then becomes yi=(β+ρ)xi+εi.$$y_i=(\beta +\rho )x_i+{\epsilon }_i.$$ The εi${\epsilon }_i$ have equal variances (equal to 1+ρ2$1+{\rho }^2$ in this case) and, again, are uncorrelated with the xi$x_i$. Therefore, the OLS estimate ˆρ=Cov(u,x)Var(x)$$\hat{\rho }=\frac{\mathrm{Cov}(u,x)}{\mathrm{Var}\left(x\right)}$$ of ρ$\rho$ is unbiased¹, and for our toy data equals ˆρ≈−0.644$\hat{\rho }\approx -0.644$ among observations with yi≥0$y_i\ge 0$. Subtracting ˆρ$\hat{\rho }$ from ˆβ$\hat{\beta }$ then gives ˆβ−ˆρ≈0.356−(−0.644)=1, recovering the true value β=1.

The table below reports 95% confidence intervals for ˆβ, ˆρ, and (ˆβ−ˆρ), estimated by simulating the DGP yi=xi+ui described above 100 times. The table confirms that the OLS estimate ˆβ of β=1 is unbiased among observations with xi≥0 but biased negatively among observations with yi≥0.

The estimate ˆβ always differs from β by ˆρ, which is significantly non-zero among observations with yi≥0. However, this pattern is not useful empirically because we generally don’t observe the ui and so can’t estimate ˆρ to back out the true value of β=ˆβ−ˆρ. Instead, we may use the Heckman correction to adjust for the bias introduced through non-random selection.

In empirical settings, selecting observations with xi≥0 may lead to biased estimates when (i) there is heterogeneity in the relationship between yi and xi across observations i, and (ii) OLS is used to estimate an average treatment effect.² In particular, if the xi are correlated with the observation-specific treatment effects then restricting to observations with xi≥0 changes the distribution, and hence the mean, of those effects non-randomly.