$\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.