$\newcommand{\Var}{\operatorname{Var}}$ Let $y_{ij}$ denote the observed response of the $j$th experimental unit in the $i$th treatment group, and the $e_{ij}$ are i.i.d. $N(0,\sigma^2)$ experimental errors.
Assume that the $y_{ij}$ are independently distributed as $N(\mu_i, \sigma^2)$ r.v.'s.
The model is written as:
$$ y_{ij} = \mu_i + e_{ij},\text{ for }1 \leq i \leq a, 1 \leq j \leq n_i,$$
Four main assumptions underlie the model:
a) the errors $e_{ij}$ are normally distributed
b) the errors $e_{ij}$ are homoscedastic
c) the errors $e_{ij}$ are independently distributed
d) $E(e_{ij}) = 0$ or equivalently $E(y_{ij})=\mu_i$
Prove:
$$ \Var(\hat{e}_{ij}) = \sigma^2 \left(\frac{n_i-1}{n_i}\right)$$
Attempt:
$$ \Var(\hat{e}_{ij}) = \sigma^2 \left(\frac{n_i-1}{n_i}\right)$$
$$ \begin{align} \Var(\hat{e}_{ij}) & = \Var(y_{ij}-\bar{y}_i) \\ & = \Var\left(\mu_i -e_{ij}-\frac{\sum_{j=1}^{n_i}\mu_i -e_{ij}}{n_i}\right) \\ & = \Var\left(e_{ij}-\frac{\sum_{j=1}^{n_i}-e_{ij}}{n_i}\right) \end{align} $$
I'm basically stuck here. I know that the $n_i$ would be squared on the denominator, but I can't figure out how to deal with the summation on the numerator.
$\newcommand{\var}{\operatorname{var}}\newcommand{\cov}{\operatorname{cov}}$ \begin{align} \var(y_{ij}-\bar{y}_i) & = \var(y_{ij}) -2\cov(y_{ij},\bar{y}_i) + \var(\bar y_i) \\[8pt] & = \sigma^2 - 2\frac{\sigma^2}{n_i} + \frac{\sigma^2}{n_i}. \tag{$*$} \end{align} The covariance in the second term is $$ \cov\left( y_{ij}, \frac{y_{i1}+\cdots+y_{ij}+\cdots+y_{i,n_i}}{n_i}\right) = \cov\left( y_{ij}, \frac{y_{ij}}{n_i} \right) = \frac{\cov(y_{ij}, y_{ij})}{n_i} = \frac{\sigma^2}{n_i}. $$ The third term in $(*)$ is the usual variance of an average of i.i.d. random variables.