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. 
 A: 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!
