# Why is the intercept of linear regression biased?

Out of curiosity, I conducted the following simulation (code below). Why is it that when the variance of the error term is large coefficient associated with the intercept is biased? Can you recommend some reference that discusses this? Or better yet is there a formal proof of that?

rm(list=ls())
set.seed(12345)
m  <- 10000
x1 <- runif(m,0,100)       # random numbers from uniform distribution
x2 <- 1:10000
u  <- rnorm(m,0,100)       # random numbers from standard normal distribution
y  <- 5*x1 + 2*x2 + 10 + u # generating y series
data  <- cbind(x1, x2, y)
beta1 <- c()
beta2 <- c()
beta3 <- c()
R2    <- c()
n       <- 1000 # number of loops
ksubset <- 100  # length of subset
for (i in 1:n){
datam <- data.frame(data[sample(nrow(data),ksubset), ])
ols   <- summary(lm(y~x1+x2, data=datam))
beta1 <- append(beta1, ols$coefficients[1]) beta2 <- append(beta2, ols$coefficients[2])
beta3 <- append(beta3, ols$coefficients[3]) R2 <- append(R2, ols$r.squared)
}

results <- c(mean(beta1), mean(beta2), mean(beta3))
results
## [1] 7.290909 5.027431 2.000345

The code takes thousand random samples of size 100 from the population and calculates a regression in each. Then I take the average of each of the estimated coefficients, which in theory should be equal to the original model. It works great for the slopes, but not for the intercept.

Update

I have noticed that with the population set at 10,000, the bias persists. If I increase the size of the population to 1,000,000, the bias disappears. The size of the target population was insufficient. Anyway, @Maarten's answer is a step in the right direction.

Update2

Entire population results:

> summary(lm(y~x1+x2))
Call:
lm(formula = y ~ x1 + x2)

Residuals:
Min      1Q  Median      3Q     Max
-367.24  -66.37   -0.33   66.64  385.26

Coefficients:
Estimate Std. Error  t value Pr(>|t|)
(Intercept) 7.0388165  2.6626801    2.644  0.00822 **
x1          5.0296554  0.0348480  144.331  < 2e-16 ***
x2          2.0002959  0.0003465 5772.764  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 100 on 9997 degrees of freedom
Multiple R-squared:  0.9997,  Adjusted R-squared:  0.9997
F-statistic: 1.667e+07 on 2 and 9997 DF,  p-value: < 2.2e-16

I updated the question, including simulation results (result). I also present the model results with the entire population, which in effect is biased due to the huge variance. The simulation, as expected, reproduce those results correctly. I repeat, my mistake is in the size of the target population.

• Could you add an (non-code) explanation about what the experiment does, and show the results? – Juho Kokkala May 11 '16 at 6:23
• There are some obvious problems in the code. For instance, x2 is not a vector of length 100, which it ought to be. As is well known and mathematically demonstrated, the coefficients are not biased. We have to conclude that you are asking us to debug this code. – whuber May 11 '16 at 16:38
• When I independently code the same (intended) simulation I cannot reproduce what you describe. Because your question relies entirely on the code, it evidently is purely a coding question. If you could express your question using a more general language, such as English and mathematics, we might be able to identify a statistical issue we could help you with. – whuber May 11 '16 at 18:07
• @MaartenBuis: sampling with replacement provides independent samples from your original sample, which are then dependent samples from the population of interest. What I think is confusing is that I was referring to the mistake of computing a t-statistic from mean(beta1 - 10) / (sd(beta1) / sqrt(n) ) (n is number of loops from code). You can see this is immediately flawed: fixing the original sample size (m), you can let n->infinity and your t-statistic is unbounded. – Cliff AB May 12 '16 at 16:27
• @CliffAB Absolutely agree. It's the answer I was expecting, I had left that detail aside. – Héctor Garrido May 12 '16 at 17:13

The role of the constant is to force the mean of the errors to be 0. In your program the mean of the error term is 0 in your population, but the variance is very high, so it will deviate a lot from 0 in many of your samples. This means that in each sample the constant can deviate a lot from the population value. So much so that a 1000 replications is not enough to find that on average you will get the right constant. Try and increase the number of replications to say 10,000 or 100,000 and see what happens.

Update

As promised I tried to replicate Héctor's results, but I gues I did something wrong the first time I tried, as I can't reproduce it anymore. Maybe I set the residual standard deviation at 1000 instead of 100, maybe I tried 100 replications instead of a 1000, or some silly typo like that. Anyhow here is my simulation with a 1000 replications in Stata, and a constant that is close enough to 10 given the spread.

. set seed 123456

. clear all

. program define sim
1.         drop _all
2.         set obs 100
3.         gen x1 = runiform() < .5
4.         gen x2 = rnormal()
5.         gen y  = 5*x1 + 2*x2 + 10 + rnormal(0,100)
6.         reg y x1 x2
7. end

. simulate b=_b[_cons], reps(1000) nodots : sim

command:  sim
b:  _b[_cons]

. sum

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
b |      1,000     9.94305    14.99877  -46.69164   60.36443
• $1000$ replications is more than enough replications for these parameters; a correct simulation will produce a nearly Normally distributed set of estimates of the intercept that is squarely centered on the specified value of $10$. – whuber May 11 '16 at 16:40
• @whuber I did a quick simulation myself and could replicate Hector's "results". The key was the huge residual variance. I will update the answer tomorrow with that simulation. – Maarten Buis May 11 '16 at 18:46
• Thank you, Maarten. Since you could replicate the results, I will vote to reopen the thread--but in the meantime I hope the OP will find a way to communicate the problem more clearly. – whuber May 11 '16 at 18:47
• @Cliff And that's exactly what an unbiased estimator of $10$ is supposed to do. – whuber May 11 '16 at 21:28
• @whuber: exactly. I was just offering a very easy way of informally showing this given the current code. – Cliff AB May 11 '16 at 21:33