Model Sum of Squares from ANOVA table Background: I am studying a course on statistical experiments using the textbook by Douglas Montgomery on the analysis of experiments. This is an introductory course and so I am relatively new to the study of these types of models.
Question: Consider the model $Y_{i,j,k} = \mu + R_i + C_j + T_k + \epsilon_{i,j,k}$ where $R_i$ is the impact of from row $i$, $C_j$ is the impact from column $j$ and $T_k$ is the impact from treatment $k$ where there are exactly $n$ rows, $n$ columns, and $n$ treatments. We know that the row sum of squares $S^2_{\text{row}} = \sum _ {i=1}^n \frac{Y_{i \cdot \cdot}^2}{n} - \frac{Y_{\cdot \cdot \cdot}^2}{n^2}$ and similarly for $S^2_{\text{column}}$ and $S^2_{\text{treatment}}$ we have the same form but instead we sum over $j$ and $k$ respectively from $1$ to $n$ (we take the formulae for the sums of squares directly from the corresponding ANOVA table).
I have read in the textbook (without justification) the formula for the expected value $ \color{red}{\mathbb{E}(S^2_{\text{row}}\space)}$ and (again) the same applies for the columns and treatments by replacing $R_i^2$ with $C_j^2$ (for the columns) and $T_k^2$ (for the treatments).
$\text{How do we find the formula for the expectation?}$
Note: we are making the standard modelling assumptions for these types of models where the error terms are independent and normally distributed with mean $0$ and variance $\sigma^2$.
 A: If you are reading the textbook Design and Analysis of Experiments (8th edition) by Douglas Montgomery, although the model for the Latin squares design is formally written as (note that the notations in the cited textbook were kept unchanged in equation $(1)$. It is obvious that $\alpha, \beta, \tau, p$ in $(1)$ correspond to $R, C, T, n$ in your question):
\begin{align}
y_{ijk} = \mu + \alpha_i + \tau_j + \beta_k + \epsilon_{ijk}
\begin{cases}
i = 1, 2, \ldots, p \\
j = 1, 2, \ldots, p \\
k = 1, 2, \ldots, p
\end{cases}, \tag{1}
\end{align}
which seems no different from the main-effects only model for the three-factor factorial design (admittedly, this representation could have been made less ambiguous$^\dagger$), there is an important note below this equation:

Because there is only one observation in each cell, only two of the three subscripts $i, j, k$ are needed to denote a particular observation. This is a consequence of each treatment appearing exactly once in each row and column.

Therefore, instead of interpreting $y_{i..}$ as $\sum_{j = 1}^p\sum_{k = 1}^py_{ijk}$, it should be interpreted as row sum of the $i$-th row, meaning it is the sum of $p$ summands rather than $p^2$.  Keeping this in mind, as well as the constraints $\alpha_. = \beta_. = \tau_. = 0$, model $(1)$ entails:
\begin{align}
& y_{i..} = p\mu + p\alpha_i + \tau_. + \beta_. + \epsilon_{i..} 
= p\mu + p\alpha_i + \epsilon_{i..}, \\
& y_{...} = \sum_{i = 1}^p y_{i..} = p^2\mu + p\alpha_. + \epsilon_{...}
= p^2\mu + \epsilon_{...}.  \tag{2}
\end{align}
Since $\epsilon_{i..} \sim N(0, p\sigma^2), \epsilon_{...} \sim N(0, p^2\sigma^2)$, it follows that $E[\epsilon_{i..}^2] = p\sigma^2$, $E[\epsilon_{...}^2] = p^2\sigma^2$, hence $(2)$ implies that (all the cross-product terms between the random error and the fixed effect have expected value $0$ because $E[\epsilon_{i..}] = E[\epsilon_{...}] = 0$):
\begin{align}
& E(y_{i..}^2) = p^2\mu^2 + p^2\alpha_i^2 + 2p^2\mu\alpha_i + p\sigma^2, \\
& E(y_{...}^2) = p^4\mu^2 + p^2\sigma^2. \tag{3}
\end{align}
Substitute $(3)$ back to $E[SS_{\text{rows}}] = \frac{1}{p}\sum_{i = 1}^pE[y_{i..}^2] - \frac{1}{p^2}E[y_{...}^2]$ and use $\alpha_. = 0$, it follows that
\begin{align}
 & E[SS_{\text{rows}}] = \frac{1}{p}\sum_{i = 1}^p(p^2\mu^2 + p^2\alpha_i^2 + 2p^2\mu\alpha_i + p\sigma^2) - \frac{1}{p^2}(p^4\mu^2 + p^2\sigma^2) \\
=& p\sum_{i = 1}^p\alpha_i^2 + p\sigma^2 - \sigma^2 \\
=& (p - 1)\sigma^2 + p\sum_{i = 1}^p\alpha_i^2, 
\end{align}
as desired.

$^\dagger$ For example, a better representation can be:
\begin{align}
y_{ik} = \mu + \alpha_i + \tau_{L(i, k)} + \beta_k + \epsilon_{ik}
\begin{cases}
i = 1, 2, \ldots, p \\
k = 1, 2, \ldots, p
\end{cases}, 
\end{align}
where $L$ is the Latin square (which is an order $k$ matrix) selected for the experiment. See also this link.
A: The model you described is LSD.
Here $$ S^2_R = n\sum_i\left(\bar{y}_{i\cdot\cdot}- \bar{y}_{\cdot\cdot\cdot} \right)^2.$$
Then
\begin{align}\mathbb E\left[S^2_R\right] &= \mathbb E\left[n\sum_i\left(R_i + \bar{\varepsilon}_{i\cdot\cdot} - \bar{\varepsilon}_{\cdot\cdot\cdot}\right)^2\right]\\ &= n\left[\sum_iR_i^2+\mathbb E\left\{\sum_i(\bar{\varepsilon}_{i\cdot\cdot} - \bar{\varepsilon}_{\cdot\cdot\cdot})^2\right\} + 2\underbrace{\sum_i\alpha_i\mathbb E(\bar{\varepsilon}_{i\cdot\cdot} - \bar{\varepsilon}_{\cdot\cdot\cdot})}_{=0}\right]\\ &= n\left[\sum_iR_i^2+\mathbb E\left\{\sum_i\bar{\varepsilon}^2_{i\cdot\cdot} - n\cdot \bar{\varepsilon}^2_{\cdot\cdot\cdot}\right\} \right] \\ &= n\left[\sum_iR_i^2 + \sum_i\frac{\sigma^2} n - \frac{\sigma^2}n\right] \\ &= n\sum_i R_i^2 + (n-1)\sigma^2.\end{align}
