Why the variance of the ANOVA residuals increases as the sample size grows?

Question

Consider a standard oneway ANOVA model given by $$ Y_{ij}=\theta_i+\varepsilon_{ij} $$ for $i=1,\ldots,k$ and $j=1,\ldots,n_i$, where $\varepsilon_{ij}$ are iid $N(0,\sigma^2)$ random variables. The residuals are defined as $$ e_{ij}=Y_{ij}-\frac1{n_i}\sum_{k=1}^{n_i}Y_{ik} $$ for $i=1,\ldots,k$ and $j=1,\ldots,n_i$. The variance of the residuals is given by $$ \operatorname{Var}e_{ij} =\sigma^2\Bigl(1-\frac1{n_i}\Bigr). $$ This means that the variance of the residuals increases if we increase the sample size $n_i$. Since the residuals $e_{ij}$ essentially estimate the errors $\varepsilon_{ij}$, I would expect the variance of the residuals to be larger than $\sigma^2$ and get closer and closer to $\sigma^2$ as the sample size $n_i$ increases. So the expression of the $\operatorname{Var}e_{ij}$ seems counterintuitive to me.

Why the variance of the residuals increases as the sample size $n_i$ grows?

Any help is much appreciated!

Answer 1

$\newcommand{\one}{\mathbf 1}\newcommand{\e}{\varepsilon}$We don't really need multiple groups here. It's enough to just consider $y_i = \theta + \e_i$ and then we can apply this to each group separately. In this case we have $\hat\theta = \bar y$ and $e = y - \bar y\one = (I - n^{-1}\one\one^T)y$ so $\text{Var}(e) = \sigma^2 (I-H)$ where I've set $H = n^{-1}\one\one^T$. This does indeed mean $\text{Var}(e_i) = \sigma^2 (1 - 1/n)$.

I think the key idea is that the variation of the $y_i$ around the sample mean $\bar y$ is less than around the true mean of $\theta$, since the sample mean is the minimizer of $\sum_{i=1}^n (y_i - c)^2$ over $c \in \mathbb R$. This is why $\text{Var}(e_i) < \text{Var}(\e_i)$. But as $n\to\infty$ we'll have $\bar y \to_{\text{as}} \theta$ so the variation of $y_i$ around $\bar y$ approaches the variation of $y_i$ around $\theta$, which is precisely $\e$. That's why $\text{Var}(e_i)\nearrow\text{Var}(\e_i)$.

Another aspect is that the residuals are not independent as $\text{Cov}(e_i,e_j) = -\sigma^2/n$ for $i\neq j$. This is because they are constrained to a subspace (the space orthogonal to $\one$, in this case) so knowing what $e_1$ is, say, reduces the possible variability in $e_j$ for $j\neq 1$. This effect also decreases with $n$ and the residuals are increasingly independent (since they're multivariate Gaussian here; more generally, they'll be increasingly uncorrelated).