Why variance of dependent variable differs from variance of error term in linear regression?

Question

Numerous books and lecture slides start to discuss regression analysis as follows:
$$Y=X\beta+\epsilon;\text{ where, } Y\sim N(X\beta, \sigma^2 ) \text{ and } \epsilon\sim N(0, \sigma^2)$$ George Seber wrote in his book (Linear Regression Analysis, Ed 2) at page 42: $$Var[Y] = Var[Y-X\beta] = Var[\epsilon]$$ I was intended to verify it as follows:

set.seed(123)
n=10000
int=rep(1,n)
x1=rnorm(n,5,3)
x2=rnorm(n, 20, 10)
x3=rnorm(n, -10, 2)
x=cbind(int,x1,x2, x3)
beta=c(10, 2, 0.5, 1.5)
err=rnorm(n, 0, 7)
y=x%*%beta+err
mean(err) # very close to zero
var(err) # very close to 49
mean(y) # very close to 15
var(y) # very close to 118

Somehow I triggered to check whether R square have a role in equating $Var[Y]$ and $Var[\epsilon]$. Then I checked:

var(y)*(1-(summary(lm(y~x1+x2+x3))$r.squared)) # very close to 49

Yes it is. So $Var[\epsilon] = Var[Y]\times (1-R^2)$. I do not know what I am missing here! Why George Seber and many others never mention it? Certainly, I am wrong, not George Seber. But what is my mistake?
Any lights on my confusion will be appreciated.
Thanks.

Answer 1

G. Seber assumes that $X\beta$ is constant. You do not assume that in your verification. If $X$ is random, then $Var(Y)=Var(X\beta+\varepsilon)=Var(X\beta)+Var(\varepsilon)$, with additional assumption that $X$ is independent from $\varepsilon$. In theory if $X$ is fixed then $Var(X\beta)=0$. Since you take a sample $X\beta$ clearly varies, hence $var(X\beta)$ (I made a distinction to denote empirical variance by $var$) is not zero.

Furthermore it should be noted that G. Seber is talking about theoretical model, so the variances are theoretical. This means that verifying it with single sample is bound to be troublesome, i.e. you cannot expect empirical variances to coincide, if theoretical variances do, unless your explicitly prove that they do. G. Seber does not make a claim that empirical variances coincide.

Concerning your last question, recall that definition of $R^2$ is the following:

$$R^2=1-\frac{SS_{res}}{SS_{tot}},$$

where $SS_{tot}$ happens to be $(n-1)var(y)$. So $(1-R^2)var(y)$ is the residual sum of squares divided by $n-1$. Now residual sum of squares is $\sum(y_i-x_i\hat\beta)^2$, where $\hat\beta$ is the estimate of true $\beta$. Since your model is the correct one, the $\hat\beta$ is close to the true $\beta$, so $y_i-x_i\hat\beta\approx y_i-x_i\beta=\varepsilon_i$. Since $var(\epsilon)=\frac{1}{n-1}(\sum_{i=1}\varepsilon_i-\bar\varepsilon)^2$ and $\bar\epsilon=\frac{1}{n}\sum\varepsilon\approx 0$ in your case you get your relationship. Which as you can see is not so much related to original question.