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}


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


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

  • 1
    $\begingroup$ 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"? $\endgroup$
    – JLee
    Apr 16, 2020 at 15:46
  • 1
    $\begingroup$ $\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) $$ $\endgroup$ Apr 16, 2020 at 19:46
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 16, 2020 at 20:00
  • $\begingroup$ 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? $\endgroup$
    – JLee
    Apr 19, 2020 at 12:34
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 19, 2020 at 13:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.