Multi-level data vs multi-level models

Question

I am confused about the applicability of linear multi-levels models (also called hierarchical linear models or mixed-effect models) on data with multiple levels. At first, I thought they are made for each other but after trying to use a multi-level model on my own data I got very confused.

By multi-level data I mean that the data is about different levels in an hierarchy. For example, in the popularity data used by Joop Hox there is data about the pupil (popularity, sex) and the teacher (experience and popularity). In my case I want to asses the relation between the student-level and teacher-level predictors on the popularity of the student.

What I do know now about multi-level modeling is that they assess the effects of predictors on the lowest level possibly grouped by higher levels of the hierarchy (at least that is what I assume). For example, you can use multi-level models to check if schools have a different intercept or if the effect of sex is different for schools.

The problem I am trying to solve is how to use data higher up in the hierarchy to predict values lower in the hierarchy. In my real-case example I want to predict the house-prices using data about the districts they are in. Since the number of districts is much lower than the number of houses I want to incorporate this information in the estimation of the variance of the coefficients.

Below is a parallel example with pupil-teacher data used by Joop Hox. In this example, I want to asses the effect of sex on the pupil level and the effect of teacher experience on the teacher level. The model should take into account that there are much less observations (100) of teacher experience than of student sex (2000). Note that there is only one teacher per school; for clarity I rename school to teacher.

Following Joop Hox's notation that would be the following in math:

$$\text{popularity}_{ij} = \beta_{0j} + \beta_{1j} \text{gender}_{ij} + e_{ij}$$

$$\beta_{0j} = \gamma_{00} + \gamma_{01} \text{experience}_{j} + u_{0j}$$

Ordinary linear model

library(foreign)
library(lme4)
library(arm)
library(tidyverse)
library(broom)
df <- read.dta("http://www.ats.ucla.edu/stat/stata/examples/mlm_ma_hox/popular.dta") %>%
  rename(teacher = school)
fit_lm <- lm(popular ~ sex + texp, df)
display(fit_lm)

This would underestimate the standard error of the effect of the teacher experience since it assumes there are 2000 observations of teacher experience instead of only 100 (number of teachers).

Stacked linear models

fit_lm_1 <- lm(popular ~ sex + factor(teacher) + 0, df)
teacher_effect <- tidy(fit_lm_1)[3:101,] %>%
  mutate(teacher = as.numeric(substr(term, 15,18))) %>%
  dplyr::select(teacher, estimate) %>%
  bind_rows(data.frame(teacher = 1, estimate = 0))
df_teacher <- df %>%
  group_by(teacher) %>%
  summarise(texp = first(texp)) %>%
  left_join(teacher_effect)
fit_lm_2 <- lm(estimate ~ texp, df_teacher)
display(fit_lm_2)

Takes into account that there is less data about teachers than about pupils. However, does not allow for dependencies between teacher-experience and gender (e.g. think of a private boys/girls school).

Multi-level model with teacher-experience per school

fit_lmer1 <- lmer(popular ~ sex + (texp + 0 | teacher), df)
display(fit_lmer1)
ranef(fit_lmer1)

This would give a separate slope for teacher experience for each teacher. Since there is only a single teacher per teacher (duh!) this is a strange thing to do.

Multi-level model with fixed teacher-experience and intercept per school

fit_lmer2 <- lmer(popular ~ sex + texp + (1 | teacher), df)
display(fit_lmer2)
fixef(fit_lmer2)

This results in a single estimate for teacher experience. According to this blog post lmer understands that texp and teacher are related and that there are only 100 observations of teacher experience.

Fine, you would say? Not totally, I do not want to model the teacher intercept but only the teacher experience, because what could the teacher experience possibly add to an intercept per teacher. Furthermore, I doubt that lmer really understands that texp and school belong together.

In short

I have data on two separate levels and I want to draw valid conclusions about there relationships with a dependent variable. Can I use multi-level modeling for this? How to do this in R with lme4 (or possibly nlme)

Answer 1

First of all, that there are only 100 different teachers for 2000 students will not lead to an underestimation of teacher experience on the popularity of the student (measured on the student level). There are not 100 observations of teacher experience but 2000! To explain this with an analogy: if we want to study the effect of gender on an outcome variable, say test score on an exam, we might have 2000 test scores from 2000 students, but only 2 different genders. This will not lead to underestimated standard errors for the effect of gender, because there are 2000 observations of the relation between gender and test score, even though there are only two different genders.

So I'd say that your ordinary regression model is fine from that perspective:

fit_lm <- lm(popular ~ sex + texp, df)

However, students are nested within schools and it might be the case that there are some differences among schools that might be captured by modeling schools as a random intercept:

fit_lmer <- lmer(popular ~ sex + (1 | school), df)

This model does not take teacher experience into account, and you say that you don't understand what teacher experience could add to a school-based intercept. To understand this, we might think about what might influence the outcome variable (popularity, whatever that is?) on a school level. Obviously teacher experience might do so, but also perhaps the environment in the school, the social class of the students, the personality of the teacher, the skill of the teacher (might be different from teacher experience) etc. All of those variables might be bundled together as a school-level random intercept.

Then why should we add teacher experience? We have information about teacher experience, so it might be used in addition to the random intercept of school, and it might provide a better fit. What if some very experienced teachers work in a very bad school (which might be associated with lower popularity)? Or a very inexperienced teacher works in a very good school? By using both teacher experience and school itself as a random variable, we take the rest of the random variance associated with teacher/school into account, while also taking teacher experience into account. This might provide a better fit for the data, and consequently more accurate estimates:

fit_lmer2 <- lmer(popular ~ sex + texp + (1 | school), df)

However, it is also possible that the influence of teacher experience on popularity varies across schools. In some schools, the environment might be so supportive, the education of the teachers so high etc, that the influence of the individual teacher's experience is lower. Conversely, some schools might be so dysfunctional that the individual teacher's experience might have a much higher influence. I don't work in the education system, but from my experience within the health care system, this is not a too far-fetched idea. Some clinics work very well, and all patients get more or less the same care regardless of clinician's experience, while other clinics are more dysfunctional and the level of care depends to a high degree on the individual clinician's expertise.

So this idea might be modeled as a random slope in addition to the random intercept of school:

fit_lmer3 <- lmer(popular ~ sex + texp + (1 + texp | school), df)

I would try model 2 and 3, and then compare them:

anova(fit_lmer2, fit_lmer3)

In summary: yes, you can use multilevel modeling to compare data on separate levels. There is no problem with using district-level data (perhaps income, education and other things) along with using the districts themselves as random variables when predicting house prices.

Answer 2

In the end I decided that the stacked linear model was the best option for my application. With more than a two-level hierarchy this could be difficult, but for me it is suitable.

In the first run the linear model distinguishes between the pupil and teacher effect. The pupil effect is directly explained by the pupil predictors (sex and intercept) and the teacher effect is a dummy coding of the teacher.

fit_lm_11 <- lm(popular ~ sex + factor(teacher) + 0, df)

Then I summarise the data on the teacher level, including the intercepts (coefficients for the dummy variables:

teacher_effect <- tidy(fit_lm_11)[3:101,] %>%
  mutate(teacher = as.numeric(substr(term, 16,19))) %>%
  dplyr::select(teacher, estimate) %>%
  bind_rows(data.frame(teacher = 1, estimate = 0))
df_teacher <- df %>%
  group_by(teacher) %>%
  summarise(texp = first(texp)) %>%
  left_join(teacher_effect)

And finally learn the influence of the teacher experience on his or hers intercept.

fit_lm_12 <- lm(estimate ~ texp, df_teacher)

This explains 43% of the variation between teachers:

> display(fit_lm_12)
lm(formula = estimate ~ texp, data = df_teacher)
            coef.est coef.se
(Intercept) -3.57     0.17  
texp         0.09     0.01  
---
n = 100, k = 2
residual sd = 0.71, R-Squared = 0.43

And the standard error of the teacher experience estimate is much more realistic (0.01086) than the one of the ordinary linear model (0.003299). About three times higher.