I have a dilemma with respect to the included (decomposition) between bias and variance in the calculation of the Mean square error (MSE) for the OLS estimator with the equation: MSE = bias ^ 2 + variance
I calculated with R software the bias, the variance and the MSE. As you will see I run the code many times (replications = 1000 times).
And it is clear that MSE is not equal to bias ^ 2 + variance. However, I do not think I was wrong in the 3 formulas (bias, variance and MSE) neither in my R code in general. Since equality is not respected, can I then deduce a strong/extra noise in the data ? If not, how to explain this inequality ?
You can copy and paste the code below.
###########################
# Data
PIB.hab<-c(12000,34000,25000,43000,12500,32400,76320,45890,76345,90565,76580,45670,23450,34560,65430,65435,56755,87655,90755,45675)
ISQ.2018<-c(564,587,489,421,478,499,521,510,532,476,421,467,539,521,478,532,449,487,465,500)
Dataset=data.frame(ISQ.2018,PIB.hab)
#plot
plot(ISQ.2018,PIB.hab)
plot(ISQ.2018,PIB.hab, main="Droite de régression linéaire", xlab="Score ISQ 2018", ylab="PIB/hab")
#OLS fit
fit1<-lm(PIB.hab~ISQ.2018)
lines(ISQ.2018, fitted(fit1), col="blue", lwd=2)
# Create a list to store the results
lst<-list()
# This statement does the repetitions (looping)
for(i in 1 :1000)
{
n=dim(Dataset)[1]
p=0.667
sam<-sample(1 :n,floor(p*n),replace=FALSE)
Training <-Dataset [sam,]
Testing <- Dataset [-sam,]
fit2<-lm(PIB.hab~ISQ.2018)
ypred<-predict(fit2,newdata=Testing)
y<-Dataset[-sam,]$PIB.hab
MSE <- mean((y-ypred)^2)
biais <- mean(ypred-y)
variance <-mean((ypred- mean(ypred))^2)
lst[[i]] <- c(MSE = MSE,
biais = biais,
variance = variance)
# lst[i]<-MSE
# lst[i]<-biais
# lst[i]<-variance
}
# convert to a matrix
x <- as.matrix(do.call(rbind, lst))
colMeans(x)
###########################