Write Several Normal Distributions into Matrix Variate Normal

Question

For $i = 1, \ldots, k$ and $j = 1,\ldots, n_i$, let $X_{i,j} \in \mathbb{R}^{p \times 1}$'s be independent identically distributed $N_p(\mu, \Sigma)$ random variables, and, independent to $X_{i,j}$'s, let $\epsilon_{i,j}$'s be independent identically distributed $N(0, \sigma_{\epsilon}^2)$ random variables. Write $N = \sum_{i=1}^kn_i$. For each $i = 1, \ldots, k$, let $u_i$ be independent identically distributed $N(0, \sigma_u^2)$ random variables which is independent to both $X_{i,j}$'s and $\epsilon_{i,j}$'s. Fix $a \in \mathbb{R}$ and $b \in \mathbb{R}^p$. For each $i,j$, put :

\begin{equation*} Y_{i,j} = a + b^T X_{i,j} + u_i + \epsilon_{i,j} \end{equation*}

That is, linear mixed model with random intercept on $k$ classes . Let $Z$ be the $(p+1) \times N$ matrix:

\begin{equation*} Z = \begin{bmatrix} X_{1,1}&X_{1,2}&\cdots& X_{1,n_1} & \cdots & X_{k, n_k}\\ Y_{1,1}&Y_{1,2}&\cdots& Y_{1,n_1} & \cdots & Y_{k, n_k} \end{bmatrix} \end{equation*}

The question is how to write the distribution of $Z$ as a matrix variate normal distribution. Write $X_{i,j}=[X_{i,j}^{(1)}, \ldots, X_{i,j}^{(p)}]^T$. We have :

\begin{equation*} \text{Vec}(Z^T) = \begin{bmatrix} X_{1,1}^{(1)}\\ X_{1,2}^{(1)}\\ \vdots \\ X_{1, n_1}^{(1)}\\ \vdots \\ X_{k, n_k}^{(1)}\\ \vdots \\ X_{1,1}^{(p)} \\ \vdots \\ X_{k, n_k}^{(p)}\\ Y_{1,1} \\ \vdots \\ Y_{k, n_k} \end{bmatrix} \end{equation*}

It seems that the $X$ part of $\text{Vec}(Z^T)$ has covariance matrix $\Sigma \otimes I_N$. However, the $Y$ part seems to need a matrix like $E$, where $E$ is block diagonal matrix with $1_{n_i}1_{n_i}^T$ as blocks. Here $1_n \in \mathbb{R}^{n \times 1}$ denotes the vector of entry $1$ everywhere. So that, it seems impossible to write that $Z$ follows $N_{p+1,N}$ distribution with some covariance matrix in format $A \otimes B$, where $A$ is $(p+1)\times (p+1)$ and $B$ is $N \times N$. I am probably missing something obvious.

If matrix variate normal does not work out for $Z$, is there any way to write a matrix containing all $X$'s and $Y$'s so that the result is a matrix normal distribution?

Answer 1

I doubt if it is possible for unbalanced random blocks. I seems to me that $\otimes$ isn't flexible to handle different block sizes. Maybe I am wrong.

In case of balanced design with $n_1=n_2=\cdots=n_k=n$, an obvious way is to put all data within a block in a column, i.e., $$\tilde{Z}=\begin{bmatrix} X_{11} & X_{21} & \dots & X_{k1} \\ Y_{11} & Y_{21} & \dots & Y_{k1} \\\hline X_{12} & X_{22} & \dots & X_{k2} \\ Y_{12} & Y_{22} & \dots & Y_{k2} \\\hline \vdots & \vdots & \vdots & \vdots \\\hline X_{1n} & X_{2n} & \dots & X_{kn} \\ Y_{1n} & Y_{2n} & \dots & Y_{kn} \\ \end{bmatrix}.$$ Then $\text{vec}(\tilde{Z})$ has covariance $I_k \otimes V_{\text{block}}$, where $$ V_{\text{block}} = I_n \otimes V_{x,y-u} + (1_n 1_n^T ) \otimes V_u, $$ with $$V_{x,y-u}=\begin{bmatrix} \Sigma & b^T\Sigma \\ \Sigma^Tb & b^T\Sigma b+ \sigma_{\epsilon}^2 \end{bmatrix}$$ and $$V_u=\begin{bmatrix} 0_{p\times p} & 0_{p\times 1} \\ 0_{1\times p} & \sigma_u^2 \\ \end{bmatrix}.$$

PS: What's the benefit to write it in matrix variate normal form, compared to multivariate normal form?