# Mixed Effects Model (3 level model?)

Consider the following problem. The dataset that I am considering has $$n=1800$$ units (high-end copying machines). Label the units $$i = 1,\dots,n$$. Unit $$i$$ has $$n_i$$ recordings. It is of interest to model the use-rate for these copying machines. All machines are in the same building.

The following linear mixed effects model is used:

\begin{equation} \begin{aligned} X_i(t_{ij}) &= m_i(t)+ \varepsilon_{ij} \\ &= \eta + z_i(t_{ij})w_i + \varepsilon_{ij}, \end{aligned} \end{equation}

where $$\eta$$ is the mean, $$z_i(t_{ij}) = [1, \log(t_{ij})]$$, $$w_i = (w_{0i}, w_{1i})^\top \sim N(0,\Sigma_w)$$, $$\varepsilon_{ij} \sim N(0, \sigma^2)$$, and

$$\begin{equation} \Sigma_w = \begin{pmatrix} \sigma^2_1& \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma^2_2 \end{pmatrix}. \end{equation}$$

I can write this model in matrix form. More specifically, I have the model (I write this out for a reason)

$$\begin{equation} X = 1\eta + Zw + \varepsilon, \end{equation}$$

where

$$\begin{equation} X = \begin{pmatrix} X_1\\ \vdots \\ X_n \end{pmatrix} \in \mathbb{R}^N, \varepsilon = \begin{pmatrix} \varepsilon_1\\ \vdots \\ \varepsilon_n \end{pmatrix} \in \mathbb{R}^N, 1 = \begin{pmatrix} 1_{n_1}\\ \vdots \\ 1_{n_n} \end{pmatrix} \in \mathbb{R}^{N \times p}, w = \begin{pmatrix} w_1\\ \vdots \\ w_n \end{pmatrix} \in \mathbb{R}^{2n}, \end{equation}$$

where $$N = \sum_{i=1}^n n_i$$. In addition,

$$\begin{equation} Z = \begin{pmatrix} Z_1 & 0_{n_1 \times 2} & \dots & 0_{n_1 \times 2} \\ 0_{n_2 \times 2} & Z_2 & \dots & 0_{n_2 \times 2} \\ \vdots & & \ddots & \vdots \\ 0_{n_n \times 2} & \dots & & Z_n \end{pmatrix} \in \mathbb{R}^{N \times 2n}, 0_{n_i \times 2} = \begin{pmatrix} 0 & 0 \\ \vdots& \vdots \\ 0 & 0 \end{pmatrix} \in \mathbb{R}^{2n_i}. \end{equation}$$

Furthermore, we have that

$$\begin{equation} \begin{bmatrix} w\\ \varepsilon \end{bmatrix} \sim N \begin{bmatrix} \begin{pmatrix} 0\\ 0 \end{pmatrix},&\sigma^2 \begin{pmatrix} G(\gamma) & 0 \\ 0 & R(\rho) \end{pmatrix} \end{bmatrix}, \end{equation}$$

where $$\gamma$$ and $$\rho$$ are $$r \times 1$$ and $$s \times 1$$ vectors of unknown variance parameters corresponding to $$w$$ and $$\varepsilon$$, respectively. Mathematically,

$$\begin{equation} G = \frac{1}{\sigma^2} \begin{pmatrix} \Sigma_w & \dots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \dots & \Sigma_w \end{pmatrix} \in \mathbb{R}^{2n \times 2n}, R = \begin{pmatrix} I_{n_1} & \dots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \dots & I_{n_n} \end{pmatrix} \in \mathbb{R}^{N \times N}, \end{equation}$$

where $$w_i \sim N(0, \Sigma_w)$$, and $$\varepsilon_i \sim N(0, \sigma^2I_{n_i})$$. Here $$\gamma = (\sigma_1, \sigma_2, \rho)^\top$$ and $$\rho = \sigma^2$$.

Imagine I now obtain a dataset for a new building with $$n$$ units. But now, unit $$i$$ is in the same room as unit $$i+1$$ for $$i = 1,3,5,\dots, n-1$$. How would I model the additional dependence between units in the same room? At first I thought to use the exact same model as above but changing $$G$$ to

$$\begin{equation} G = \frac{1}{\sigma^2} \begin{pmatrix} \Sigma_w & \Sigma_{1,2} & \dots &0& 0 \\ \Sigma_{1,2}& \Sigma_w & \dots &0& 0 \\ \vdots & \vdots& \ddots & \vdots& \vdots \\ 0 & 0& \dots& \Sigma_w & \Sigma_{1799,1800} \\ 0 & 0& \dots & \Sigma_{1799,1800}& \Sigma_w \end{pmatrix} \in \mathbb{R}^{2n \times 2n}, \end{equation}$$

where $$\Sigma_{i, i+1}$$ is the covariance matrix which models the dependence between units $$i$$ and $$i+1$$ for $$i = 1,3, \dots, 1799$$.

Is this a possible way to model the problem? I guess it would not be possible to use nlm in R to do it but it would be possible using an analytic solution.

What else could be done? I think a three level hierarchical model (instead of two level model) could also work, but I am not sure how to formulate a three level model.

Any advice on past modelling experiences and how to write down the three level model would be appreciated.

If I understand correctly, here's your "problem before":

• There are $$n$$ rooms in the building, rooms are indexed by, say, $$j=1:n$$, each room has one unit, and each unit has $$n_i$$ recordings, recordings are indexd by $$i=1:n_j$$. You want to group the recordings by unit and build a hierarchical model on it.

• There are $$n$$ rooms in the building, rooms are indexed by $$j=1:n$$, each room has $$m_j$$ units(in your case $$m_j\le 2$$, but it doesn't matter). Units in a room are indexed by $$k=1:m_j$$, the $$k$$th unit in room $$j$$ has $$n_{jk}$$ recordings, recordings are indexd by $$i=1:n_{jk}$$. You want to group the recordings by unit and room, and build a hierarchical model on it.

Your "problem before" is a basic Bayesian linear regression model with linear Gaussian observations and Gaussian priors. To simplify the representation, let's assume the recordings are already centered, and use $$\Sigma$$ to replace $$\Sigma_w$$. The CPDs will be: $$x_{ji} \sim N(z_{ji}w_j,\sigma^2), i=1:n_j, j=1:n \\ w_j \sim N(0,\Sigma), j=1:n$$ In your equations there's no prior distribution specified for $$\sigma^2$$ and $$\Sigma$$, so it's not really a "hierarchical model". In order to make it hierarchical, let's extend it with two prior distributions, say using inverse Gamma for $$\sigma^2$$ and inverse Wishart for $$\Sigma$$, the extended CPDs will be: $$x_{ji} \sim N(z_{ji}w_j,\sigma^2), i=1:n_j, j=1:n \\ w_{j} \sim N(0,\Sigma), j=1:n \\ \Sigma \sim InvWishart(v,S) \\ \sigma^2 \sim InvGamma(a,b)$$

With the hierarchical representation in mind, the CPDs for the "problem after" will be: $$x_{jki} \sim N(z_{jki}w_{jk},\sigma^2),i=1:n_{jk}, k=1:m_j j=1:n\\ w_{jk} \sim N(0,\Sigma_j), k=1:m_j,j=1:n \\ \Sigma_j \sim InvWishart(v,S) ,j=1:n\\ \sigma^2 \sim InvGamma(a,b)$$

The Markov blankets for the hidden random variables ($$w,\Sigma,\sigma^2$$) can be easily derived from the CPDs, Gibbs sampling will be a natural solution.

Here's an [article] shows how to run Gibbs sampling on a hierarchical Bayesian linear regression model in R. It's slightly different from your model but it provides all the materials you need to build one.

Good luck!

• I see how $w_{jk} \sim N(0, \Sigma_j)$ models correlation between units within the same room, however, there will be correlation between rooms. How does $w_{jk} \sim N(0, \Sigma_j)$ model this? Also, would the G matrix I proposed work for the "problem after"? – JLee Apr 16 at 15:46
• $\Sigma_j \sim InvWishart(v,S)$ captures the correlation of the variances of $w_jk,k=1:m_j$ in room $j$, if you also want to capture the correlation of the means, denote $\mu_j$, of $w_jk,k=1:m_j$ in room $j$, you should add a Gaussian prior to it. Further more, you can merge the Gaussian and Inverse-Wishart prior into a Normal-Inverse-Wishart prior for $\mu_j,\Sigma_j$: $$x_{jki} \sim N(z_{jki}w_{jk},\sigma^2),i=1:n_{jk}, k=1:m_j j=1:n\\ w_{jk} \sim N(\mu_j,\Sigma_j), k=1:m_j,j=1:n \\ \mu_j,\Sigma_j \sim NIW(m,k,v,S) ,j=1:n\\ \sigma^2 \sim InvGamma(a,b)$$ – Haotian Chen Apr 16 at 19:46
• Sorry I misread the question. You are right, to capture the correlation between rooms we need another layer, for example: $$x_{jki} \sim N(z_{jki}w_{jk},\sigma^2),i=1:n_{jk}, k=1:m_j j=1:n\\ w_{jk} \sim N(0,\Sigma_j), k=1:m_j,j=1:n \\ \Sigma_j \sim InvWishart(v_1,S_1) ,j=1:n\\ S_1 \sim InvWishart(v_2,S_2) \\ \sigma^2 \sim InvGamma(a,b)$$ Where $v_1$, $v_2$, $S_2$, $a$, and $b$ are hyperparameters. And yes I do think your $G$ can capture the correlation. But using Bayesian methods can be more robust, making the model less sensitive to outliers. – Haotian Chen Apr 16 at 20:00
• I think your answer (and the last comment) both suffice if I understand correctly. My understanding being that $\Sigma_j$ models the covariance between units in room $j$. In addition, if I wanted to model the recordings in a new room (say room $l$), I would need to sample a covariance matrix from the "population" of covariance matrices, which has a $InvWishart(v_1, S_1)$ distribution. Is this the correct interpretation of the model? – JLee Apr 19 at 12:34
• Yes that's correct. Keep in mind that $S_1$ is also a hidden random variable, so you also need to sample $S_1$ from $InvWishart(v_2,S_2)$ and the same for all other hidden random variables. – Haotian Chen Apr 19 at 13:50