The Variance of the Slope in a Regression Model

In my "applied linear models" exam, there was a tricky question (it was a multiple choice, so no details were asked). I was simply asking if the following statement was valid, or not

Consider a linear regression with one single covariate, y=β0+β1x1+ε and the least-square estimates. The variance of the slope is Var[β1] Do we decrease this variance if we add one variable, and consider y=β0+β1x1+β2x2+ε ?

For the exam, the expected answer was simply "no." In a nutshell, there are two cases where we should expect different changes:

• if x1 and x2 are highly correlated, then we should expect the variance to increase.
• if x1 and x2 are not correlated, then we should expect the variance to decrease.

We briefly observed (and discussed) those points on examples during the lecture... but I wanted to go a bit further since I couldn't find any analytical results. Let's generate a model, y=β0+β1x1+β2x2+ε, and then compare the variance, Var[β1], on the two fitted modes, depending on the correlation between x1 and x2.

library(mnormt)
n=200
s=function(r=0){
S=matrix(c(1,r,r,1),2,2)
X=rmnorm(n,c(0,0),S)
B=data.frame(y=-2+X[,1]+X[,2]+rnorm(n)/2,
x1=X[,1],
x2=X[,2])
reg12=lm(y~x1+x2,data=B)
reg1=lm(y~x1,data=B)
k=summary(reg12)$coefficients[2,2]/summary(reg1)$coefficients[2,2]
k}

Let's generate 500 samples for each value of the correlation, from -0.9 to _0.9. 

M=NULL
for(r in ((-9):9)/10) M=cbind(M,Vectorize(s)(rep(r,500)))

And let's plot the ratio of the two variances: 

plot(0:1,0:1,xlim=c(-1,1),ylim=c(0,2),col="white")
for(i in 1:19) points(rep((((-9):9)/10)[i],500),M[,i],col="light blue")
VM=apply(M,2,mean)
lines((((-9):9)/10),VM,col="red",lwd=2)
abline(h=1,lty=2)

If the ratio exceeds 1, the variance increases when adding a covariate.

Indeed, here, when the two variables are independent, the variance is divided by two. But when covariates are highly correlated, the variance is multiplied by two.

Now, what if, actually, x2 is not a real explanatory variable: the true model we generate is y=β0+β1x1+ε. In that case,

s=function(r=0){
S=matrix(c(1,r,r,1),2,2)
X=rmnorm(n,c(0,0),S)
B=data.frame(y=-2+X[,1]+rnorm(n)/2,
x1=X[,1],
x2=X[,2])
reg12=lm(y~x1+x2,data=B)
reg1=lm(y~x1,data=B)
k=summary(reg12)$coefficients[2,2]/summary(reg1)$coefficients[2,2]
k}

We get our samples just like we did earlier: 

M=NULL
for(r in ((-9):9)/10) M=cbind(M,Vectorize(s)(rep(r,500)))

And we plot those ratios:

plot(0:1,0:1,xlim=c(-1,1),ylim=c(0,2),col="white")
for(i in 1:19) points(rep((((-9):9)/10)[i],500),M[,i],col="light blue")
VM=apply(M,2,mean)
lines((((-9):9)/10),VM,col="red",lwd=2)
abline(h=1,lty=2)

In this case, we add a useless variable x2. Whatever the correlation with x1, it will always, on average, increase the variance of β1.