Omitted variable bias: 3 correlated variables and 1 omitted (R simulation)

Question

I get weird results when trying to analyze omitted variable bias in R.

If I try to analyze the bias for coefficient $\beta_j$ of variable $x_j$ in case of one omitted variable from set of two correlated variables, the limit bias should be given by:

$$ BIAS_{j,(lim)} = \beta_k \frac{COV \left[ x_j, x_k \right]}{VAR \left[ x_j \right]} $$

This is something I can see in R, considering I have the following setup:

$corr(x_2, x_3) \neq 0$ for $y_{(true)} = \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \epsilon$ and $y_{(est)} = \gamma_1 x_1 + \gamma_2 x_2 + u$.

Then:

$$\gamma_2 = \beta_2 + \beta_3 \frac{COV \left[ x_2, x_3 \right]}{VAR \left[ x_2 \right]}$$

BUT when I do the scenario such that:

$corr(x_2, x_3) \neq 0 \land corr(x_2, x_4) \neq 0 \land corr(x_3, x_4) \neq 0$, while omitting only $x_4$, then I do not get the results above for $x_2$ and $x_3$. Why?

The code is below:

set.seed(100)

n=1000

corr12 = 0.7
corr23 = 0.5
corr13 = 0.2

sd1 = 20
sd2 = 5
sd3 = 7

covmat = matrix(c(sd1^2, corr12*sd1*sd2, corr13*sd1*sd3,
                  corr12*sd1*sd2,  sd2^2, corr23*sd2*sd3,
                  corr13*sd1*sd3, corr23*sd2*sd3, sd3^2), ncol = 3, byrow = T)



ostat = mvrnorm(n, c(100,80, 120), Sigma = covmat)



z = rnorm(n,100, 30)

x1 = rnorm(n,100,20)

x2 = ostat[,1]

x3 = ostat[,2]

x4 = ostat[,3]

y = 20 + 2*x1 + 4*x2 + 2*x3 + 5*x4 +  rnorm(n,0,1)


cor(cbind(x1,x2,x3,x4,y))

olsT = lm(y~x1 + x2 + x3 + x4)
summary(olsT)


ols1 = lm(y ~x1 + x2 + x3)
summary(ols1)

5*(cov(x2,x4)/var(x2)) ## this does not provide the right bias
5*(cov(x3,x4)/var(x3)) ## this does not provide the right bias

Answer 1

Your calculations are incorrect because in a regression with more than two explanatory variables the omitted-variable bias (OVB) for x3 is not a function of cov(x3, x4) and var(x3) only; you cannot ignore x1 and x2.

Actually, in your particular simulation, you can ignore x1 as it's uncorrelated with x2 and x3. However, this is a special case; let's write down a general solution for the omitted-variable bias using matrix notation.

Let $X_1$ and $X_2$ be the matrices of included and omitted predictors, respectively; here $X_1 = [\text{Intercept}, x_1,x_2, x_3]$ and $X_2 = [x_4]$. Also let $\beta_1,\beta_2$ be the corresponding parameters in the full model Y ~ X1 + X2 and $\gamma$ be the regression coefficients in the reduced model Y ~ X1. We can show:

$$ \begin{aligned} \operatorname{E}(\widehat{\gamma} | X_1^\vphantom{'}) &= \beta_1 + \left(X_1'X_1^\vphantom{'}\right)^{-1} X_1'X_2^\vphantom{'} \beta_2^\vphantom{'} \end{aligned} $$

The bias is the second term. For a derivation, see the wikipedia article on Omitted-variable bias as well as this answer by @YashaswiMohanty which explains the matrix math nicely (+1).

In short, you need to substitute cov(x2,x4) / var(x2) and cov(x3,x4) / var(x3) with a single matrix equation solve(t(X1) %*% X1) %*% t(X1) %*% (X2) where X1 is the matrix of included variables and X2 is the matrix of excluded variables. The term t(X1) %*% (X2) generalizes the covariances while the term t(X1) %*% X1 generalizes the variances. solve() computes the matrix inverse. (NB: In practice we wouldn't use solve() to fit a regression.)

And here are the calculations in R:

#> Coefficients (full model):
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 20.404477   0.642784   31.74   <2e-16 ***
#> x0           1.998913   0.001612 1240.16   <2e-16 ***
#> x1           4.002192   0.002158 1854.66   <2e-16 ***
#> x2           1.977233   0.009510  207.91   <2e-16 ***
#> x3           5.010879   0.005451  919.21   <2e-16 ***

#> Coefficients (reduced model):
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 314.98029   16.23879   19.40   <2e-16 ***
#> x0            1.86917    0.04679   39.95   <2e-16 ***
#> x1            3.53918    0.06115   57.88   <2e-16 ***
#> x2            6.55068    0.23619   27.73   <2e-16 ***

The omitted-variable bias for x1 is about 3.5 - 4 = -0.5 and that for x2 is 6.5 - 2 = 4.5.

Next we calculate the OVB in terms of the design matrices $X_1$ and $X_2$.

XT <- model.matrix(olsT)

X1 <- XT[, 1:4, drop = FALSE] # matrix of included predictors + intercept
X2 <- XT[, 5, drop = FALSE] # matrix of excluded predictors

solve(t(X1) %*% X1) %*% t(X1) %*% (X2) * 5 # beta2 = 5
#>                      x3
#> (Intercept) 293.9362880
#> x0           -0.1294633
#> x1           -0.4620058
#> x2            4.5635190