1
$\begingroup$

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
$\endgroup$
0

1 Answer 1

4
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.