How to Calculate the Expectation and Variance of this One-Way Random Effects Model? So, this is the model:
$Y_{ij} = \mu + T_{i} + \epsilon_{ij}$ for $i = 1, ..., t$ and $j = 1, ..., t$.
I understand that

*

*$T_{1}, T_{2}, \dots, T_{t} \stackrel{iid}\sim N(0, \sigma^{2}_{T})$;

*$\epsilon_{11}, \epsilon_{22}, \dots, \epsilon_{tn} \stackrel{iid}\sim N(0, \sigma^{2}_{T})$; and

*$T_{1}, T_{2}, \dots, T_{t}$ are independent of $E_{11}, E_{22}, \dots, E_{tn}$
My task is to find the expectation of $\bar Y_{++}$ aka, the mean across all values of $i$ and $j$.
Now, just from looking at it's definition: $\bar Y_{++} = \frac{1}{nt}\sum_{i = 1}^{t}\sum_{j = 1}^{n}Y_{ij}$
I can immediately recognize that this is the sample mean across all values of $i$ and $j$ and so recalling the fact that the sample mean is an unbiased estimator, mentally I'm aware of the fact that the expectation of $\bar Y_{++}$ should be just $\mu$, but undergoing the calculation here, I have:
$E[\bar Y_{++}] = E[\frac{1}{nt}\sum_{i = 1}^{t}\sum_{j = 1}^{n}Y_{ij}] = \frac{1}{nt}E[\sum_{i = 1}^{t}(\sum_{j = 1}^{n} Y_{ij})] = \frac{1}{nt}E[\sum_{i = 1}^{t}(n\bar y_{i+})] = \frac{1}{nt}E[nt\bar y_{ij}] = \frac{nt}{nt}E[\bar y_{ij}] = E[\bar y_{ij}] = \mu$
But now I have no idea how to compute the variance from this. I know the formula for variance (for some random arbitrary random variable $X$) is $E[X^2] - (E[X])^2$ but inputting the expression for $\bar Y_{++}$ for X just gets stupid messy. Am I just doing it wrong?
I was also thinking about calculating the mean and variance this way:
$E[\bar Y_{++}] = E[\mu + T_{i} + \epsilon_{ij}] = E[\mu] + E[T_{i}] + E[\epsilon_{ij}] = \mu + 0 + 0 = \mu$
$Var(\bar Y_{++}) = Var(\mu + T_{i} + \epsilon_{ij}) = Var(T_{i} + \epsilon_{ij}) = Var(T_{i}) + Var(\epsilon_{ij})$ (b/c independent) $ = \frac{\sigma^{2}_{T}}{t} + \frac{\sigma^{2}}{nt}$
Both computations for the expectation yield the same answer, and the variance computation seems correct, so I'm inclined to go with the latter method of calculating the two statistics, but is there a way to calculate the variance using that variance definition provided above with how $\bar Y_{++}$ is defined? Or is that just the wrong way to approach this problem and I should instead just go with the latter way of computing the two statistics like I did above?
Thanks!
 A: For the mean, you can apply the linearity of expectation a bit more aggressively than you do: since this is just a finite sum we know we can exchange expectation and summation so
$$
\text E\left(\sum_{ij} Y_{ij}\right) = \sum_{ij} \text E Y_{ij} = nt\mu
$$
and then $\text E \bar Y_{++}$ follows.
For the variance, the key step is remembering that the random intercepts induce a correlation within each group:
$$
\text{Var}(\bar Y_{++}) = (nt)^{-2} \text{Var}\left(\sum_{ij} Y_{ij}\right) \\
= (nt)^{-2} \sum_{ijkl} \text{Cov}(Y_{ij}, Y_{kl}).
$$
We can compute this covariance on a case-by-case basis:
$$
\text{Cov}(Y_{ij}, Y_{kl}) = \begin{cases} 0 & i \neq k \\ \sigma^2_T &  i=j,k\neq l \\ \sigma^2_T + \sigma^2_\epsilon & i=j, k=l \end{cases}
$$
which reflects how observations in different blocks are independent, observations within a block have a covariance of $\sigma^2_T$, and each individual observation has a variance of $\sigma^2_T + \sigma^2_\epsilon$.
Now we just need to count how many terms of $\sigma^2_T + \sigma^2_\epsilon$ we get and how many terms of just $\sigma^2_T$, add them up, and divide by $(nt)^2$ to get the answer.
Your derivation is incorrect because $\text{Var}(T_i) + \text{Var}(\epsilon_{ij}) = \sigma^2_T + \sigma^2_\epsilon$, not what you wrote. You need to account for both the individual variances but also the covariances.

$\newcommand{\one}{\mathbf 1}$Here's an alternative derivation with matrices and vectors which I personally prefer even though it's a bit more notation to set up.
We can rewrite the model as
$$
Y = \mu\one + ZT + \epsilon
$$
where $\one$ is a vector of all $1$s, $T = (T_1,\dots,T_t)^T$, $\epsilon$ collects all the errors, and $Z$ is a binary matrix that indicates which element of $T$ corresponds to each element of $Y$. For example, the first row of $Z$ will be $(1, 0, \dots, 0)$ so $(ZT)_1 = T_1$ which just picks out the correct element of $T$ for $Y_{11}$.
By the independence of $T$ and $\epsilon$ we have
$$
Y \sim \mathcal N(\mu\one, \sigma^2_T ZZ^T + \sigma^2_\epsilon I)
$$
so
$$
\bar Y_{++} = (nt)^{-1}Y^T\one \sim \mathcal N(\mu, \sigma^2_T(nt)^{-2}\one^TZZ^T\one + (nt)^{-2}\sigma^2_\epsilon \one^T\one ) \\
\stackrel{\text d}= \mathcal N(\mu, \sigma^2_T(nt)^{-2}\one^TZZ^T\one + (nt)^{-1}\sigma^2_\epsilon)
$$
since $\bar Y_{++}$ is just a linear transformation of $Y$ so it's still Gaussian.
In this case since every group is the same size we have $\one^TZ = (n, \dots, n)^T$ since each column of $Z$ has exactly $n$ 1s for each observation in each group. Then $\one^TZZ^T\one = n^2t$ so
$$
(nt)^{-2}\one^TZZ^T\one = t^{-1}
$$
and the rest can be worked out.
