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*}
Best Answer
You were on the right track. But the likelihood is simpler than what it looks in your post.
If $\boldsymbol X_i=(X_{ij})_{1\le j\le n}$ for every $i$, then $\boldsymbol X_i$'s are i.i.d $N_n(\mu\mathbf1_n,\Sigma)$ with
$$\Sigma=\sigma^2 I_n + \tau^2\mathbf1_n\mathbf1_n^T$$
So joint pdf of the $X_{ij}$'s is joint pdf of the $\boldsymbol X_i$'s, given by
\begin{align} f_{\boldsymbol\theta}(\boldsymbol x_1,\ldots,\boldsymbol x_p)&\propto \frac1{(\det\Sigma)^{p/2}}\exp\left\{-\frac12 \sum_{i=1}^p (\boldsymbol x_i-\mu\mathbf1_n)^T\Sigma^{-1}(\boldsymbol x_i-\mu\mathbf1_n)\right\} \\&=\frac1{(\det\Sigma)^{p/2}}\exp\left\{-\frac12 Q(\boldsymbol x_1,\ldots,\boldsymbol x_p;\boldsymbol\theta)\right\} \end{align}
For sufficient statistics of $\boldsymbol\theta=(\mu,\tau^2,\sigma^2)$, we only need to focus on the sum $Q$ in the exponent.
We have
$$\Sigma^{-1}=\frac1{\sigma^2}I_n-\frac{\tau^2}{\sigma^2(\sigma^2+n\tau^2)}\mathbf1_n\mathbf1_n^T$$
So,
$$ Q=\frac1{\sigma^2}\sum_{i,j} (x_{ij}-\mu)^2 - \frac{\tau^2}{\sigma^2(\sigma^2+n\tau^2)}\sum_i\left\{\sum_j (x_{ij}-\mu)\right\}^2 \tag{1} $$
Keep in mind that
$$\overline x_{i\bullet}=\frac1n\sum_{j=1}^n x_{ij} \quad, \quad \overline x_{\bullet\bullet}=\frac1p \sum_{i=1}^p \overline x_{i\bullet}$$
Recall the definition of sum of squares due to errors $(\text{SSE})$ and sum of squares due to a factor $A$ (say) corresponding to $\alpha_i$ (call this $\text{SSA}$). The right hand side of $(1)$ can be written in terms of $\text{SSE}$, $\text{SSA}$ and an additional term involving $\mu$. Once you do that, you would notice that $f_{\boldsymbol\theta}$ is a member of a regular (full-rank) exponential family and a sufficient statistic would be easy to identify.
Note that
\begin{align} \sum_{i,j}(x_{ij}-\mu)^2&=\sum_{i,j}\{(x_{ij}-\overline x_{i\bullet})+(\overline x_{i\bullet}-\mu)\}^2 \\&=\sum_{i,j}(x_{ij}-\overline x_{i\bullet})^2 + n\sum_i (\overline x_{i\bullet}-\mu)^2 \end{align}
And also
\begin{align} \sum_i (\overline x_{i\bullet}-\mu)^2&=\sum_{i}\{(\overline x_{i\bullet}-\overline x_{\bullet\bullet})+(\overline x_{\bullet\bullet}-\mu)\}^2 \\&=\sum_{i}(\overline x_{i\bullet}-\overline x_{\bullet\bullet})^2 + p (\overline x_{\bullet\bullet}-\mu)^2 \end{align}
Therefore continuing from $(1)$,
\begin{align} Q&=\frac1{\sigma^2}\sum_{i,j} (x_{ij}-\mu)^2 - \frac{n^2\tau^2}{\sigma^2(\sigma^2+n\tau^2)}\sum_i (\overline x_{i\bullet}-\mu)^2 \\&=\frac1{\sigma^2}\sum_{i,j}(x_{ij}-\overline x_{i\bullet})^2 + \frac{n}{\sigma^2+n\tau^2} \sum_i (\overline x_{i\bullet}-\mu)^2 \\&=\frac1{\sigma^2}\sum_{i,j}(x_{ij}-\overline x_{i\bullet})^2 + \frac{n}{\sigma^2+n\tau^2}\sum_i(\overline x_{i\bullet}-\overline x_{\bullet\bullet})^2 + \frac{np}{\sigma^2+n\tau^2}(\overline x_{\bullet\bullet}-\mu)^2 \end{align}
Expressing the likelihood in this form also helps to find MLE of $\boldsymbol \theta$.