A Beginner‘s Guide to Mixed Effects Regression for Hierarchical Modeling

In many research settings, data has a hierarchical or nested structure. For example, students are grouped within classrooms, patients are nested within hospitals, and repeated measures are grouped within individuals. In these cases, observations within a group tend to be more similar to each other than to observations from different groups. This violates the independence assumption of standard statistical methods like linear regression.

Hierarchical linear modeling, also known as mixed effects modeling, extends standard regression to account for the grouping structure in the data. It allows modeling the effects of both individual-level and group-level variables on the outcome. Mixed effects regression has become an essential tool in many fields including social sciences, biology, and medical research for analyzing hierarchical data.

In this blog post, we‘ll dive into the basics of mixed effects models, work through an example in R, discuss some important considerations and limitations, and point to resources for learning more. Whether you‘re a practicing data scientist or an aspiring statistician, understanding mixed models will allow you to extract insights from a wide range of data structures. Let‘s get started!

Why Use Mixed Effects Models?

Before we get into the details of mixed effects regression, let‘s consider why we need it in the first place. Suppose we want to study the relationship between student test scores and family income. A standard regression model would estimate a single intercept and slope, assuming the relationship between income and test scores is the same across all schools:

TestScore = β0 + β1*Income + ε

However, this ignores the fact that students are nested within schools. Students in the same school are likely to have more similar test scores than students from different schools, due to shared teachers, resources, student demographics, etc. The standard errors of the regression coefficients will be underestimated because the model assumes all observations are independent.

A mixed effects model accounts for school-level clustering by estimating a separate intercept for each school. This models between-school variability in mean test scores. We can also allow the effect of income to vary across schools by estimating a random slope. The mixed model looks like:

TestScoreij = β0 + u0j + (β1 + u1j)*Incomeij + εij

• β0 is the fixed intercept, the overall mean test score across schools
• u0j is the random intercept, the deviation of school j‘s mean from the overall mean
• β1 is the fixed slope, the average effect of income on test scores
• u1j is the random slope, the deviation of school j‘s slope from the average
• εij is the residual error for student i in school j

The random intercepts uj and slopes u1j are assumed to come from a normal distribution with mean 0 and variance estimated from the data. This models the correlations among observations within each group.

By partitioning the variability into within-school and between-school components, mixed models provide more accurate estimates of the fixed effects and their standard errors. They can also estimate the variance and covariance of the school-level random effects, which is often of substantive interest. For example, we could test if there is significant variability in the income slope between schools.

Fitting Mixed Models in R

To demonstrate mixed effects models in action, we‘ll analyze a dataset of math test scores from students in 100 schools. The dataset has 3 variables:

• MathScore: The student‘s score on a standardized math test
• Income: The family income of the student (in thousands of dollars)
• SchoolID: A unique identifier for the school the student attends

Let‘s start by fitting a standard linear regression model predicting math scores from income, ignoring the school groupings:

lm_model <- lm(MathScore ~ Income, data=mathdata)
summary(lm_model)

Call:
lm(formula = MathScore ~ Income, data = mathdata)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.842 -11.290   0.736  12.428  39.649 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 51.94219    0.53548   97.01   <2e-16 ***
Income       3.88175    0.06943   55.91   <2e-16 ***
---
Signif. codes:  0 ‘***‘ 0.001 ‘**‘ 0.01 ‘*‘ 0.05 ‘.‘ 0.1 ‘ ‘ 1

Residual standard error: 16.04 on 4998 degrees of freedom
Multiple R-squared:  0.3845,    Adjusted R-squared:  0.3844 
F-statistic:  3126 on 1 and 4998 DF,  p-value: < 2.2e-16

The model estimates that for every $1000 increase in family income, math scores increase by 3.88 points on average. However, the residual standard error is quite large at 16 points. Let‘s see if accounting for schools can improve the model fit.

We‘ll use the lme4 package to fit a mixed effects model with random intercepts for schools:

library(lme4)
lmer_model <- lmer(MathScore ~ Income + (1|SchoolID), data=mathdata)
summary(lmer_model)

Linear mixed model fit by REML [‘lmerMod‘]
Formula: MathScore ~ Income + (1 | SchoolID)
   Data: mathdata

REML criterion at convergence: 38415.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4693 -0.6749  0.0288  0.6757  3.5665 

Random effects:
 Groups   Name        Variance Std.Dev.
 SchoolID (Intercept) 120.09   10.958  
 Residual             145.71   12.071  
Number of obs: 5000, groups:  SchoolID, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)  52.5148     1.1235   46.74
Income        3.8214     0.0579   66.04

Correlation of Fixed Effects:
       (Intr)
Income -0.266

The mixed model estimates a fixed income slope of 3.82, similar to the standard regression. But it also estimates the standard deviation of the school-level intercepts to be 10.96. This means 95% of schools are expected to have mean math scores within ±21.5 points of the overall mean, reflecting substantial between-school variability.

The residual standard deviation is reduced to 12.07, a 25% reduction compared to the standard model. This is because the mixed model separates the residual variance into within-school and between-school components. The intraclass correlation (ICC), the proportion of total variance due to schools, is 120.09 / (120.09 + 145.71) = 0.45.

We can also allow the income slope to vary across schools by adding a random slope term:

lmer_model2 <- lmer(MathScore ~ Income + (1 + Income|SchoolID), data=mathdata)

This model estimates a standard deviation of 0.86 for the school-level slopes. The correlation between the random intercepts and slopes is -0.61, suggesting schools with higher average scores tend to have flatter income gradients.

Comparing Models

To check if the random effects improve model fit, we can compare the models using likelihood ratio tests:

anova(lm_model, lmer_model)
anova(lmer_model, lmer_model2)

Adding random intercepts significantly improves fit compared to the standard model (p < 0.001). Adding random slopes further improves fit (p < 0.001), indicating there is significant between-school variation in the effect of income on achievement.

We can also compare models using information criteria like AIC and BIC, which balance model fit and complexity. Lower values indicate preferred models. The random intercepts and slopes model has the lowest AIC and BIC, confirming it provides the best fit to the data.

Assumptions and Diagnostics

Mixed models, like standard regression models, make certain assumptions about the data. These include:

1. Linearity: The relationship between the predictors and the outcome is linear.
2. Normality: The residuals are normally distributed with mean 0 and constant variance.
3. Independence: The observations are independent, conditional on the random effects.
4. Homoscedasticity: The residual variance is constant across levels of the predictors.

We can check these assumptions by examining residual plots and normal quantile plots of the random effects. The lme4 package provides several diagnostic functions for this purpose.

It‘s also important to consider the structure of the random effects. We should have at least 5-10 observations per random effect level to estimate the variance components reliably. If there are too few levels or observations per level, the model may not converge or may give unstable estimates.

Overly complex random effects structures can also lead to overfitting. It‘s recommended to start with a simple model and add complexity gradually, testing for improvement in fit at each step. Plotting the random effects can help identify outliers or unusual patterns that may drive results.

Conclusion

In this post, we‘ve covered the basics of mixed effects regression for analyzing hierarchical data. We saw how mixed models extend standard regression by estimating both fixed and random effects, providing a flexible framework to model group-level variability.

Mixed models are a powerful tool, but they also come with some caveats. They require careful specification of the random effects structure and attention to model assumptions. The estimation methods are more complex and computationally intensive than standard regression.

Despite these limitations, mixed effects models have become increasingly popular in many fields due to their ability to handle complex data structures. They can be extended to non-normal outcomes, crossed random effects, multiple levels of nesting, and other settings. We didn‘t have space to cover all of these topics here, but I encourage you to explore them further.

I hope this post has given you a foundation to start using mixed models in your own work. Remember, all models are wrong, but some are useful. The goal is not to find the "true" model, but to find a model that is a reasonable approximation to reality and helps answer your scientific questions. Happy modeling!

Further Reading

• Gelman, A., & Hill, J. (2006). Data Analysis Using Regression and Multilevel/Hierarchical Models (1st ed.). Cambridge University Press.
• Pinheiro, J., & Bates, D. (2006). Mixed-Effects Models in S and S-PLUS. Springer Science & Business Media.
• Galecki, A., & Burzykowski, T. (2013). Linear Mixed-Effects Models Using R: A Step-by-Step Approach. Springer Science & Business Media.
• Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1–48.