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

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


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