# Multiple correlated random non-nested intercepts in R

I am trying to estimate a longitudinal model in R in which there are several random intercepts that are correlated with each other, and the data are non-nested. For example, consider a simple longitudinal model in which student test scores are regressed on three random intercepts, one random school effect, one random teacher effect, and one random teacher-school match effect. The data are non-nested because teachers frequently switch schools in the longitudinal data. I am primarily interested in the random effects parameters themselves, and so I want to specify all three random effects and allow them to be correlated with each other. There are no random slopes, only random intercepts.

I tried:

model1 <- lmer(test_score ~ FEs + (1|schoolid) + (1|teacherid) + (1|matchid), data)


but this appears to force the random effects to be uncorrelated with each other. Is there a way, using lmer or any other R package to estimate a similar model in which the random intercepts are allowed to be correlated with each other (using non-nested data)?

More specifically, the above model imposes the covariance structure: $$Cov \left[ \begin{array} \\ \theta_{s} \\ \theta_{t} \\ \theta_{m} \\ \end{array} \right] = \left[ \begin{array} \\ \sigma^{2}_{\theta_{s}}I_{s} & 0 & 0 \\ 0& \sigma^{2}_{\theta_{t}}I_{t} & 0\\ 0&0&\sigma^{2}_{\theta_{m}}I_{m}\\ \end{array} \right]$$

What I would like to estimate instead is a random effects model with the following covariance structure:

$$Cov \left[ \begin{array} \\ \theta_{s} \\ \theta_{t} \\ \theta_{m} \\ \end{array} \right] = \left[ \begin{array} \\ \sigma^{2}_{\theta_{s}}I_{s} & \sigma^{2}_{\theta_{st}}C_{st} & \sigma^{2}_{\theta_{sm}}C_{sm} \\ \sigma^{2}_{\theta_{ts}}C_{ts} & \sigma^{2}_{\theta_{t}}I_{t} & \sigma^{2}_{\theta_{tm}}C_{tm}\\ \sigma^{2}_{\theta_{ms}}C_{ms} &\sigma^{2}_{\theta_{mt}}C_{mt} &\sigma^{2}_{\theta_{m}}I_{m}\\ \end{array} \right]$$ where $C_{st}$, for example, is a matrix of zeros and ones, with $C(i,j)=1$ where there is an observation of the pair $(school=i, teacher=j)$ in the data, and 0 if no such pair is observed.

• If your interest is in the intercepts themselves, you may want to use fixed effects. W/ random effects you can't get estimates of the paramters, only predictions. Are you wanting to test the variances? I'm not sure this makes sense. Sep 8, 2014 at 19:37
• The model is not identifiable with fixed effects unless the matchid is assumed to be uncorrelated with schoolid and teacherid, which is precisely the assumption I am trying to avoid. The matchid would absorb all of the variation from schoolid and teacherid in this case, leaving nothing to identify those parameters. Sep 8, 2014 at 19:45
• In response to question about variances: Yes, by parameters I am referring to variance parameters. I want to know what fraction of the variation in test scores is attributable to schools vs teachers vs school-teacher pairs (where there is reason to believe that certain teachers have comparative advantages at certain schools). For further motivation on these types on problems see, for example, link Sep 8, 2014 at 19:45
• Hmmm. Do you know if there's an existing LMM framework (MLwin, HLM, etc.) that handles this structure? At a quick glance, the linked document estimates independent variances ... Sep 9, 2014 at 23:05