I use simple linear regression and I want to find the decomposition of MSE, that is as a sum of the bias, the variance and the variance of the error term. I have the following code:
n <- 100
Sigma <- 2
x <- rnorm(n, 0, 3)
y <- x + rnorm(n, 0, sqrt(Sigma))
sim_data <- data.frame(y,x)
model <- lm(y ~ x, data = sim_data)
mse1 <- mean(summary(model)$residuals^2)
fitted <- predict(model,data.frame(x))
mse2 <- mean((fitted-y)^2)
# decomposition of MSE into bias, variance and variance of error
mse3 <- mean((x-mean(fitted))^2) + var(fitted) + Sigma
mse1
and mse2
are equal, but mse3
is not the correct MSE. What am I doing wrong in the decomposition of MSE into bias, variance and variance of the error terms?