Multilevel Analysis – Fitting Slopes-as-Outcomes Linear Mixed Model Using lme4

lme4-nlmemultilevel-analysisrandom-effects-model

In education research there is a type of multilevel model referred to as a slopes-as-outcomes model by Raudenbush and Bryk (2002). Raudenbush and Bryk's HLM software uses its own notation in describing multilevel models, but it is often easy to convert into lme4 notation. The slopes-as-outcomes model they describe however, does not seem to have any easy analog. Model 4 at this link shows a worked example using the HLM software. This example fits the following model:

$$
MATHACH_{ij} = β_{0j} + β_{1j} SES + r_{ij}
$$

$$
β_{0j} = γ_{00} + γ_{01}(SCHTYPE) + γ_{02}(MEANSES) + u_{0j}
$$

$$
β_{1j} = γ_{10} + γ_{11}(SCHTYPE) + γ_{12}(MEANSES) + u_{1j}
$$

This model is saying that for each school j both an intercept term $B_{0j}$ and a slope for socioeconomic status (SES) $B_{1j}$ is estimated. However, both $B_{0j}$ and $B_{1j}$ are predicted by the variables MEANSES and school type (SCHTYPE).

In lme4 I would attempt to fit this model by using the following command:

library(lme4)
mod1 <- lmer(MATHACH ~ SES + MEANSES + SCHTYPE + (1 + SES | SCH)

In R I know this will return coefficients for SES, MEANSES, and SCHTYPE as fixed effects, and an intercept and SES slope for each school as part of the random effects. I also know that the slope and intercept for the schools are estimated conditional on the effects of SCHTYPE and MEANSES.

In HLM however, the reported outputs are different. As the example I linked to shows, HLM reports 2 coefficients for MEANSES and SCHTYPE — one each for the effect on the $j$ intercepts and $j$ SES slopes. So instead of an output where the call to fixef(mod1) looks like this:

fixef(fm1)
(Intercept)        SCHTYPE    MEANSES
251.40510         10.46729    9.7734

It would look like this:

fixef(fm1)
-- Intercepts --
(Intercept)        SCHTYPE    MEANSES
 251.410           10.49      9.74
-- Slope SES --
(Intercept)        SCHTYPE    MEANSES
  232.40            4.49      8.34

I am not sure how HLM derives these estimates. I am wondering if there is a way within lme4 to partial out the effect of a predictor like SCHTYPE between the random intercepts and the random slopes. Or if HLM is simply doing some post-estimation regression of SCHTYPE and MEANSES on the BLUPS for the slopes and intercepts and reporting those coefficients alongside the general results for the model.

I'd like to replicate this type of analysis but in R to take advantage of R's superior handling of missing data and binomial responses.

This question has previously been raised, and the model form suggested was this: 

lmer(mathach ~ sector + meanses + ses + sector:ses + meanses:ses + 
(1+ses|school),data=dat)

However, I'm not sure that using this model one could interpret the sector:ses interaction term as the prediction of the effect of sector on the school-level ses slopes. For the analysis I would like to do I'm interested in assessing the magnitude and direction of the impact of school-level indicators on the school intercept and slope estimates.

Best Answer

This may be a late answer but I will try to answer it however. In pages 117-130 of Raudenbush and Bryk (2002), as you rightly mentioned, the slope as outcome model has been represented as

$$ MATHACH_{ij} = β_{0j} + β_{1j} SES + r_{ij} $$

$$ β_{0j} = γ_{00} + γ_{01}(SCHTYPE) + γ_{02}(MEANSES) + u_{0j} $$

$$ β_{1j} = γ_{10} + γ_{11}(SCHTYPE) + γ_{12}(MEANSES) + u_{1j} $$

This translates to a a 2 level HLM model that can be written in lmer as 

lmer(mathach ~ sector + meanses + ses + sector:ses + meanses:ses + ...

The authors explain the significant coefficients in Page 128 which are nothing but the interaction effects coefficients.

I came across a webpage that describes Bayesian meta-regression where they model the slopes of the hierarchical model, which is what you had in mind perhaps.

Related Question