I'm having trouble fitting a general (not generalized...) linear mixed-effect model using the lme4 package. Can anybody help?

Question

I need to fit a linear mixed model in the "Laird and Ware" framework. This type of model is usually specified by (as you may know):

$\mathbf{y}_i = X_i \beta + Z_i \mathbf{b}_i + \mathbf{\epsilon}_i $, where

• $\mathbf{y}_i$ is the response for group $i$
• $X_i$ is the design matrix for the fixed effects, with coefficients $\beta$
• $Z_i$ is the design matrix for the random effects, with coefficients $\mathbf{b}_i$
• $\mathbf{\epsilon}_i$ is normally distributed error, mean 0 and covariance structure $R_i$

This specific type of model is usually quite easy to fit in R using the lmer function in lme4.

If we were to stack the $\mathbf{y}_i$ vectors into a larger vector $\mathbf{y}$ that contains all profile information, we would be able to write our model as

$\mathbf{y} = \begin{bmatrix} X_1 \\ \vdots \\ X_n \end{bmatrix} \beta + \begin{bmatrix} Z_1 & 0 & \cdots & 0 \\ 0 & Z_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & Z_n \end{bmatrix} \begin{bmatrix} \mathbf{b}_1 \\ \mathbf{b}_2 \\ \vdots \\ \mathbf{b}_n \end{bmatrix} + \mathbf{\epsilon}$

Now, my problem is that I do not necessarily have this nice form for my model. I have something like the following:

$\mathbf{y}_i = X_i \beta + S_i \mathbf{u} + Z_i \mathbf{b}_i + \mathbf{\epsilon}_i $, where $\mathbf{u}$ is a random effect, but is not indexed by group.

This leads my model to be something like this:

$\mathbf{y} = \begin{bmatrix} X_1 \\ \vdots \\ X_n \end{bmatrix} \beta + \begin{bmatrix} S_1 & Z_1 & 0 & \cdots & 0 \\ S_2 & 0 & Z_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ S_n & 0 & 0 & \cdots & Z_n \end{bmatrix} \begin{bmatrix} \mathbf{u} \\ \mathbf{b}_1 \\ \mathbf{b}_2 \\ \vdots \\ \mathbf{b}_n \end{bmatrix} + \mathbf{\epsilon}$

This type of model is seen, for example, on page 6 of http://www.stat.vt.edu/research/Technical_Reports/TechReport10-2.pdf.

I want to know how I would produce this type of model in R, with or without the lme4 package. The main issue I am having is that to specify random effects in lmer, we are required to group by a certain variable. However, I need to specify a random effect without grouping by any variable. Can anyone help? Thanks

Answer 1

The model

$$ y_{ij} = \beta_0 + \beta_1 x_{ij} + b_j + \varepsilon_{ij} $$

means that

$$ b_j \sim \mathcal{N}(\mu, \gamma) \ \text{and} \ \varepsilon_{ij} \sim \mathcal{N}(0, \sigma) $$

so there is grouping effect for $b_j$. If $b$'s were randomly distributed they would be redundant with $\beta_0 + \varepsilon_{ij}$ terms, so it would not provide any additional information and in fact, would be unidentifiable. So this model definition does not seem to make much sens.