# Fitting a multilevel multivariate model in R with glmer

Background

I have a large dataset that contains three binary outcomes for individuals belonging to groups. I am interested in jointly modeling these binary outcomes because I have reason to believe they are positively correlated with one another. Most of my data are at the individual level; however, I also have some group-level information.

Because of the structure of my data, I am treating this as a 3-level logistic regression. The first level defines the multivariate structure through dummy variables that indicate for each outcome. Therefore level-1 accounts for the within-individual measurements. Level-2 provides the between-individuals variances and level-3 gives the between-group variance.

Hypothetical Example:

Suppose I have data for students in classrooms. I want to examine whether certain student-level characteristics are important predictors for passing three different pre-tests (math, history, and gym). The pre-tests are constructed such that about half of the students should pass each exam (no floor or ceiling effect). Since some students are better than others, I expect whether or not they pass their history exam to be correlated to their probability of passing their math and gym pre-tests. I also expect that students in the same classroom will perform more similarly than students across classrooms.

Here is my attempt at writing out the model

I use $h$ to index level-1, $i$ to index level-2, and $j$ to index level-3. Recall that level-1 corresponds to within-individual, level-2 corresponds to between-individuals, and level-3 corresponds to between-groups. So, I have $h$ measures for the $i^\text{th}$ individual in the $j^\text{th}$ group.

Let \begin{align} \alpha_{1ij} &= 1 \text{ if outcome}_1 = 1 \text{ and 0 otherwise} \\ \alpha_{2ij} &= 1 \text{ if outcome}_2 = 1 \text{ and 0 otherwise} \\ \alpha_{3ij} &= 1 \text{ if outcome}_3 = 1 \text{ and 0 otherwise} \end{align}

\begin{align} \text{log}\left( \frac{\pi_{hij}}{1 - \pi_{hij}} \right) &= \alpha_{0hij} + \alpha_{1hij}Z_{ij} + \eta_h \\ Z_{ij} &= \beta_{0ij} + X_{ij}\beta_{ij} + U_{j} + \epsilon_i \\ U_{j} &= \gamma_{0j} + X_{j}\gamma_{j} + \rho_j \end{align}

Trying to specify the model in R

I have been using R 3.0.1, but am open to solutions using other standard software (e.g. SAS, Stata, WinBUGs, etc.). Since I am modeling a binary response, I am using the glmer function in the lme4 package.

My data set is in long format, with one row per outcome per individual per group.

One of my current problems is correctly specifying a 3-level model. Given 5 individual-level measures and one group-level measure, I have tried to specify the model as:

> glmer(pi ~ outcome1:(x1 + x2 + x3 + x4 + u5) +
outcome2:(x1 + x2 + x3 + x4 + u5) +
outcome3:(x1 + x2 + x3 + x4 + u5) +
(1 + x1 + x2 + x3 + x4 | individual) +
(1 + u5 | group), family=binomial, data=pretest)


but this often produces warnings and does not converge.

My questions

1. Does my approach make sense? (i.e. does it make sense to model a multilevel multivariate model by specifying the multivariate structure in the first level?)

2. I have difficult with the proper notation and appreciate suggestions for clarifying my notation.

3. Am I using the right tools for this problem? Should I be using other packages or software?

• Where did you get the idea to model the multivariate responses this way? What is pi? – ndoogan Feb 22 '14 at 1:48
• I got the idea from a chapter in a book – which was only briefly described so I am not certain my interpretation is correct. pi is a binary response indicator; I didn't want to use y since it is sort of confusing with outcome. – Ellis Valentiner Feb 24 '14 at 16:14
• I think I understand. I believe outcome{1,2,3} are actually supposed to be indicators of the outcome represented by pi, not the outcome measurement itself--pi captures that. Then the parameter estimates for them are the outcome specific intercepts (thus, exclude the standard intercept with +0). Then to model the correlation between the outcomes, you need allow the intercepts to vary randomly and correlate with a term like (outcome1+outcome2+outcome3+0|individual). This is a comment because I am not comfortable enough with the procedure to make it an "answer". – ndoogan Feb 25 '14 at 15:51