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*}

• What is $\mathbf I$ and $\mathbf P$? Nov 27 '21 at 21:41
• @StubbornAtom $\mathbf N$ is a $p \times p$ matrix full of ones, and $\mathbf I$ is the $p \times p$ identity matrix. I've added this to the question. Thanks. Nov 27 '21 at 21:44
• I meant $\mathbf P$, not $\mathbf N$. Sorry. Nov 29 '21 at 0:35
• There are closed-form solutions for the MLE of all parameters of this model. This may be found most simply by reparametrizing the variance parameters (Hint: obtain the eigen-decomposition of the Compound Symmetric (Exchangeable) covariance structure) Nov 30 '21 at 4:46
• @user277126 Thanks. Do you have any further information about finding the eigendecomposition of the covariance matrix for a given block? I figure it can't be done on a computer or via the characteristic polynomial. Dec 1 '21 at 3:25

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

• Thanks for your answer. I understand (1) but I'm not familiar enough with ANOVA to understand what you've written in the last paragraph. Any chance you could explain further or point me to something that explains further? Thanks. Nov 29 '21 at 2:25
• Applying 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*} Nov 29 '21 at 3:30
• $\text{SSE}=\sum_{i,j}(x_{ij}-\overline x_{i\bullet})^2$ and $\text{SSA}=n\sum_i(\overline x_{i\bullet}-\overline x_{\bullet\bullet})^2$, by definition. Use your standard trick of subtracting and adding a certain quantity to simplify $(1)$ in terms of $\text{SSE}$ and $\text{SSA}$. This term $\sum_{i = 1}^{p} \sum_{j = 1}^n (x_{ij} - \bar x_{\bullet \bullet})^2$ that you wrote is the total sum of squares, which equals $\text{SSE}+\text{SSA}$. Nov 29 '21 at 8:38
• Thanks for your help. I think I follow everything you've done to this point. But how do I manipulate the $(\bar x_{\bullet \bullet} - \mu)^2$ term to show that this is an exponential family? Thanks again. Dec 2 '21 at 17:00
• @Novice Actually this form of the likelihood is helpful in finding MLEs and unbiased estimators based on complete sufficient statistics. My argument here is same as showing $T_1=(\overline x,\sum (x_i-\overline x)^2)$ is complete sufficient for $(\mu,\sigma^2)$ in a $N(\mu,\sigma^2)$ population: First we show $T_2=(\sum x_i,\sum x_i^2)$ is complete sufficient from the structure of exponential family. Since $T_1$ is a one-to-one function of $T_2$, we say $T_2$ is also complete sufficient. Dec 3 '21 at 7:44

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