Mixed Effects Model (3 level model?)

Question

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.

Answer 1

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.

And your "problem after" is:

• 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!