Random effects one-way layout analysis I am analyzing the model $X_{ij} = \mu + \alpha_i + \epsilon_{ij}$ where $\alpha_i \sim N(0, \tau^2)$ (i.i.d), $\epsilon_{ij} \sim N(0, \sigma^2)$ (i.i.d), and they're each independent of each other. Also, $i = 1, \dots, p$ and $j = 1, \dots, n$.
I am trying to find minimal sufficient statistics for $\mu, \sigma^2$ and $\tau^2$, check whether they're complete, and find MLEs.
I don't want to show too much work because I'm trying to avoid making the post too long. However, I will show some of the work I did to determine the likelihood.
Note: Below, $\mathbf P$ is a $p \times p$ matrix full of ones, and $\mathbf I$ is the $p \times p$ identity matrix.
My questions are: Am I on the right track? If so, how do I proceed from here? I have tried forming the likelihood ratio and also begun trying to find MLEs but everything is so messy that I can't figure out what I'm doing.
My Work:
The likelihood of $\mathbf X$ will be the product of the multivariate normal distributions of the blocks that comprise $\mathbf X$, i.e.
\begin{align*}
L(\mu, \sigma^2, \tau^2; \mathbf x) &= \prod_{j = 1}^p f_{\mathbf X_j}(\mathbf x_j).
\end{align*}
We already determined that $\mathbf X_i \sim \mathbf N_n(\boldsymbol{\mu}_i, \boldsymbol{\Sigma})$, and $\boldsymbol{\Sigma} = \sigma^2\mathbf I + \tau^2 \mathbf P$, so we need to write out the density $f_{\mathbf X_j}(\mathbf x_j)$.
In general, if $\mathbf Y \sim \mathbf N_k\big(\boldsymbol{\zeta}, \boldsymbol{\Omega}\big)$, then its density is
\begin{align*}
f_{\mathbf Y}(\mathbf y) &= (2\pi)^{-\frac{k}{2}}\det(\boldsymbol{\Omega})^{-\frac{1}{2}} \exp\Big[-\frac{1}{2}(\mathbf y - \boldsymbol{\zeta})^\top \boldsymbol{\Omega}^{-1} (\mathbf y - \boldsymbol{\zeta})\Big].
\end{align*}
In our case, using the given value for the determinant, the $i$th block has density
\begin{align*}
f_{\mathbf X_i}(\mathbf x_i) &= (2\pi)^{-\frac{n}{2}} \big[(\sigma^2)^{p - 1}(\sigma^2 + p\tau^2) \big]^{-\frac{1}{2}} \exp\Big[ -\frac{1}{2} (\mathbf x_i - \boldsymbol{\mu}_i)^\top \Big(-\frac{\tau^2}{\sigma^2(\tau^2n + \sigma^2)}\mathbf P + \frac{1}{\sigma^2}\mathbf I\Big) (\mathbf x_i - \boldsymbol{\mu}_i) \Big].
\end{align*}
The likelihood is therefore
\begin{align*}
L(\mu, \sigma^2, \tau^2; \mathbf x) &= \prod_{i = 1}^p f_{\mathbf X_i}(\mathbf x_i)\\
&= \prod_{i = 1}^p (2\pi)^{-\frac{n}{2}} \big[(\sigma^2)^{p - 1}(\sigma^2 + p\tau^2) \big]^{-\frac{1}{2}} \exp\Big[ -\frac{1}{2} (\mathbf x_i - \boldsymbol{\mu}_i)^\top \Big(-\frac{\tau^2}{\sigma^2(\tau^2n + \sigma^2)}\mathbf P + \frac{1}{\sigma^2}\mathbf I\Big) (\mathbf x_i - \boldsymbol{\mu}_i) \Big].
\end{align*}
If we let $\gamma := -\frac{\tau^2}{\sigma^2(\tau^2n + \sigma^2)}$ for simplicity's sake then we have
\begin{align*}
(\mathbf x_i - \boldsymbol{\mu}_i)^\top \Big(\gamma\mathbf P + \frac{1}{\sigma^2}\mathbf I\Big) (\mathbf x_i - \boldsymbol{\mu}_i) &= \gamma \mathbf x_i^\top \mathbf P \mathbf x_i + \frac{1}{\sigma^2} \mathbf x_i^\top \mathbf x_i - \gamma \boldsymbol{\mu}_i^\top \mathbf P \mathbf x_i - \frac{1}{\sigma^2}\boldsymbol{\mu}_i^\top \mathbf x_i\\
&- \gamma \mathbf x_i^\top \mathbf P \boldsymbol{\mu}_i - \frac{1}{\sigma^2}\mathbf x_i^\top \boldsymbol{\mu}_i + \gamma \boldsymbol{\mu}_i^\top \mathbf P \boldsymbol{\mu}_i + \frac{1}{\sigma^2}\boldsymbol{\mu}_i^\top \boldsymbol{\mu}_i\\ \\
&= \gamma \mathbf x_i^\top \mathbf P \mathbf x_i + \frac{1}{\sigma^2} \mathbf x_i^\top \mathbf x_i - \gamma \mu \mathbf 1_i^\top \mathbf P \mathbf x_i - \frac{1}{\sigma^2}\mu \mathbf 1_i^\top \mathbf x_i\\
&- \gamma \mu \mathbf x_i^\top \mathbf P \mathbf 1_i - \frac{1}{\sigma^2} \mu \mathbf x_i^\top \mathbf 1_i + \gamma \mu^2 \mathbf 1_i^\top \mathbf P \mathbf 1_i + \frac{1}{\sigma^2}\mu^2 \mathbf 1_i^\top \mathbf 1_i
\end{align*}
where we've written $\mu \mathbf 1_i = \boldsymbol{\mu}_i$.
This enables us to rewrite the likelihood in the form
\begin{align*}
L(\mu, \sigma^2, \tau^2; \mathbf x) &= (2\pi)^{-\frac{np}{2}} \big[(\sigma^2)^{p - 1}(\sigma^2 + p\tau^2) \big]^{-\frac{p}{2}}\\
&\cdot \exp \Bigg[ -\frac{1}{2}\Bigg( \gamma \sum_{i = 1}^p \mathbf x_i^\top \mathbf P \mathbf x_i + \frac{1}{\sigma^2} \sum_{i = 1}^p \mathbf x_i^\top \mathbf x_i - \gamma \mu \sum_{i = 1}^p \mathbf 1_i^\top \mathbf P \mathbf x_i - \frac{1}{\sigma^2}\mu \sum_{i = 1}^p \mathbf 1_i^\top \mathbf x_i\\
&- \gamma \mu \sum_{i = 1}^p \mathbf x_i^\top \mathbf P \mathbf 1_i - \frac{1}{\sigma^2} \mu \sum_{i = 1}^p \mathbf x_i^\top \mathbf 1_i + \gamma \mu^2 \sum_{i = 1}^p \mathbf 1_i^\top \mathbf P \mathbf 1_i + \frac{1}{\sigma^2}\mu^2 \sum_{i = 1}^p \mathbf 1_i^\top \mathbf 1_i \Bigg) \Bigg]\\ \\
&= (2\pi)^{-\frac{np}{2}} \big[(\sigma^2)^{p - 1}(\sigma^2 + p\tau^2) \big]^{-\frac{p}{2}}\\
&\cdot \exp \Bigg[ -\frac{1}{2}\Bigg( \gamma \sum_{i = 1}^p \mathbf x_i^\top \mathbf P \mathbf x_i + \frac{1}{\sigma^2} \sum_{i = 1}^p \mathbf x_i^\top \mathbf x_i - 2\gamma \mu \sum_{i = 1}^p \mathbf 1_i^\top \mathbf P \mathbf x_i - 2\frac{1}{\sigma^2}\mu \sum_{i = 1}^p \mathbf 1_i^\top \mathbf x_i\\
&+ \gamma \mu^2 \sum_{i = 1}^p \mathbf 1_i^\top \mathbf P \mathbf 1_i + \frac{1}{\sigma^2}\mu^2 \sum_{i = 1}^p \mathbf 1_i^\top \mathbf 1_i \Bigg) \Bigg].
\end{align*}
Edit:
Following StubbornAtom's help, I've rewritten the likelihood using $Q$ as described. I've also applied the standard trick of subtracting and adding a certain quantity, I can rewrite the sums in the exponent as below. But I don't know if this is helpful or not.
\begin{align*}
\sum_{i = 1}^{p} \sum_{j = 1}^n (x_{ij} - \mu)^2 &= \sum_{i = 1}^{p} \sum_{j = 1}^n (x_{ij} - \bar x_{\bullet \bullet})^2 + pn(\bar x_{\bullet \bullet} - \mu)^2, \text{ and}\\
\sum_{i = 1}^p \Big[ \sum_{j = 1}^n (x_{ij} - \mu) \Big]^2 &= \sum_{i = 1}^p \Big[ \sum_{j = 1}^n (x_{ij} - \bar x_{i \bullet}) \Big]^2 + \sum_{i = 1}^p \big[ n(\bar x_{i \bullet} - \mu) \big]^2.
\end{align*}

You have gotten slightly confused on the correct dimensions to use in parts of your work.  Here is a complete solution for finding the MLEs.  The reparametrization recommended also allows for a quick way to obtain their sufficient statistics.
Let $\boldsymbol{j}_n$ denote the n-dimensional vector of ones, ${\boldsymbol{J}_n = \boldsymbol{j}_n\boldsymbol{j}_n^{\prime}}$ represent the $n \times n$ matrix of ones, $\boldsymbol{P}_n = \boldsymbol{J}_n/n$, and $\boldsymbol{Q}_n = \boldsymbol{I}_n-\boldsymbol{P}_n$.  Furthermore, let $\boldsymbol{\Sigma} = \sigma^2 \boldsymbol{I}_n + \tau^2 \boldsymbol{J}_n$.  Then your model may be expressed as
\begin{eqnarray*}
\boldsymbol{x}_i = \begin{pmatrix}
x_{i1} \\
\vdots \\
x_{in}
\end{pmatrix} \overset{\mbox{i.i.d}}{\sim} N_n \left(\mu \boldsymbol{j}_n, \boldsymbol{\Sigma} \right) \quad \mbox{for} \quad i = 1,\cdots p.
\end{eqnarray*}
Hence, we can also express the model as
\begin{eqnarray}
\boldsymbol{x} = \begin{pmatrix}
\boldsymbol{x}_1 \\
\vdots \\
\boldsymbol{x}_p
\end{pmatrix} \sim N_{np} \left(\mu \boldsymbol{j}_{np},\boldsymbol{I}_p \otimes \boldsymbol{\Sigma} \right).
\end{eqnarray}
For convenience, let $\boldsymbol{z} = \boldsymbol{x} - \mu \boldsymbol{j}_{np}$, and note that the first differential of $\boldsymbol{z}$ is $\mbox{d}\boldsymbol{z} = - \boldsymbol{j}_{np} \mbox{d} \mu$.
Next, we must note that $\boldsymbol{\Sigma}$ may be re-written as
\begin{eqnarray*}
\boldsymbol{\Sigma} =  \sigma_{z1} \boldsymbol{P}_n + \sigma_{z2} \boldsymbol{Q}_n,
\end{eqnarray*}
where $\sigma_{z1}=n \tau^2 +\sigma^2$ and $\sigma_{z2} = \sigma^2$.  And that the determinant and inverse of $\boldsymbol{\Sigma}$ can be written as
\begin{eqnarray*}
\left|\boldsymbol{\Sigma}\right| &=& \sigma_{z1}\sigma_{z2}^{n-1} \\
\boldsymbol{\Sigma}^{-1} &=& \frac{\boldsymbol{P}_n}{\sigma_{z1}} + \frac{\boldsymbol{Q}_n}{\sigma_{z2}}.
\end{eqnarray*}
Clearly, we can invoke the invariance property of the MLEs to obtain the MLE of $\tau^2$ by using the formula
\begin{eqnarray*}
\tau^2 = \frac{\sigma_{z1}-\sigma_{z2}}{n}.
\end{eqnarray*}
Some final results include $\left|\boldsymbol{I}_p \otimes \boldsymbol{\Sigma}\right|=\left|\boldsymbol{\Sigma}\right|^p$ and $\left(\boldsymbol{I}_p \otimes \boldsymbol{\Sigma}\right)^{-1} = \boldsymbol{I}_p \otimes \boldsymbol{\Sigma}^{-1}$.
Up to an additive constant, the log-likelihood of $\boldsymbol{x}$ is
\begin{eqnarray*}
l \left(\mu,\sigma_{z1},\sigma_{z2}|\boldsymbol{x}\right) &=& -\frac{p}{2} \left[\ln \left(\sigma_{z1}\right) + (n-1)\ln \left(\sigma_{z2}\right) \right] -\frac{1}{2} \boldsymbol{z}^{\prime} \left[\boldsymbol{I}_p \otimes \left(\frac{\boldsymbol{P}_n}{\sigma_{z1}} + \frac{\boldsymbol{Q}_n}{\sigma_{z2}}\right) \right] \boldsymbol{z}.
\end{eqnarray*}
The 1st differential of the log-likelihood is
\begin{eqnarray*}
\mbox{d}l \left(\mu,\sigma_{z1},\sigma_{z2}|\boldsymbol{x}\right) &=& -\frac{p}{2} \left[\frac{\mbox{d}\sigma_{z1}}{\sigma_{z1}} + (n-1)\frac{\mbox{d}\sigma_{z2}}{\sigma_{z2}}\right] + \boldsymbol{z}^{\prime} \left[\boldsymbol{I}_p \otimes \left(\frac{\boldsymbol{P}_n}{\sigma_{z1}} + \frac{\boldsymbol{Q}_n}{\sigma_{z2}}\right) \right] \boldsymbol{j}_{np} \mbox{d}\mu \\ &+& \frac{1}{2} \boldsymbol{z}^{\prime} \left[\boldsymbol{I}_p \otimes \left(\frac{\boldsymbol{P}_n}{\sigma_{z1}^2}\mbox{d}\sigma_{z1} + \frac{\boldsymbol{Q}_n}{\sigma_{z2}^2}\mbox{d}\sigma_{z2}\right) \right] \boldsymbol{z}.
\end{eqnarray*}
We can simplify the term involving $\mbox{d}\mu$ by noting that $\boldsymbol{Q}_n \boldsymbol{j}_n = \boldsymbol{0}_n$, hence
\begin{eqnarray*}
\left[\boldsymbol{I}_p \otimes \left(\frac{\boldsymbol{P}_n}{\sigma_{z1}} + \frac{\boldsymbol{Q}_n}{\sigma_{z2}}\right) \right] \boldsymbol{j}_{np} &=& \frac{\boldsymbol{j}_{np}}{\sigma_{z1}}.
\end{eqnarray*}
Therefore, the MLEs of $\left(\mu,\sigma_{z1},\sigma_{z2}\right)$ are
\begin{eqnarray*}
\widehat{\mu} &=& \frac{\boldsymbol{j}^{\prime}_{np}\boldsymbol{x}}{np} \\
\widehat{\sigma_{z1}} &=& \frac{\boldsymbol{x}^{\prime}\left(\boldsymbol{Q}_p \otimes \boldsymbol{P}_n\right)\boldsymbol{x}}{p} \\
\widehat{\sigma_{z2}} &=& \frac{\boldsymbol{x}^{\prime}\left(\boldsymbol{I}_p \otimes \boldsymbol{P}_n\right)\boldsymbol{x}}{(n-1)p},
\end{eqnarray*}
the formulas for $\widehat{\sigma_{z1}}$ and $\widehat{\sigma_{z2}}$ follow upon simplifying $\frac{\widehat{\boldsymbol{z}}^{\prime}\left(\boldsymbol{I}_p \otimes \boldsymbol{P}_n\right)\widehat{\boldsymbol{z}}}{p}$ and $\frac{\widehat{\boldsymbol{z}}^{\prime}\left(\boldsymbol{I}_p \otimes \boldsymbol{Q}_n\right)\widehat{\boldsymbol{z}}}{(n-1)p}$, respectively, where $\widehat{\boldsymbol{z}} = \boldsymbol{x} - \widehat{\mu} \boldsymbol{j}_{np} = \boldsymbol{Q}_{np}\boldsymbol{x}$.
