Write Several Normal Distributions into Matrix Variate Normal

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?

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?

• The reason is that I want to do Bayesian statistics : for $X$ and for $Y$.
– 温泽海
Commented Aug 3 at 16:53
• @温泽海 But isn't it doable using multivariate normal form? Commented Aug 3 at 16:54
• It is much messier. But I guess I have to do that.
– 温泽海
Commented Aug 3 at 22:45