How do you fit a mixed effects model in R?

Use the lme4 package: model <- lmer(outcome ~ fixed_predictors + (1 | group), data = df). The term in parentheses is the random effect that accounts for clustering. Load lmerTest first so the summary reports p-values, and refit with glmer() and a family when the outcome is binary or a count.

What is the difference between a random intercept and a random slope?

A random intercept lets the baseline outcome level vary across groups, written (1 | group). A random slope also lets the effect of a predictor vary across groups, written (1 + predictor | group). Add random slopes when theory suggests the relationship differs by cluster and the model still converges.

Why does lme4 not give p-values?

The lme4 authors omit p-values because the degrees of freedom for mixed models are not exactly defined. Loading the lmerTest package adds p-values to the summary using Satterthwaite degrees of freedom, which is the standard approach for reporting fixed effects in applied work.

What does a singular fit warning mean in lme4?

A singular fit means the random-effects structure is too complex for the data to support, often because a variance is estimated at or near zero. Simplify the model, usually by removing a random slope or a correlation term, until the warning disappears and the estimates are stable.

How to fit a mixed effects model in R

Q: What is the difference between a random intercept and a random slope?

A random intercept lets the baseline outcome level vary across groups, written (1 | group). A random slope also lets the effect of a predictor vary across groups, written (1 + predictor | group). Add random slopes when theory suggests the relationship differs by cluster and the model still converges.

Q: Why does lme4 not give p-values?

The lme4 authors omit p-values because the degrees of freedom for mixed models are not exactly defined. Loading the lmerTest package adds p-values to the summary using Satterthwaite degrees of freedom, which is the standard approach for reporting fixed effects in applied work.

Q: What does a singular fit warning mean in lme4?

A singular fit means the random-effects structure is too complex for the data to support, often because a variance is estimated at or near zero. Simplify the model, usually by removing a random slope or a correlation term, until the warning disappears and the estimates are stable.

To fit a mixed effects model in R, you use the lme4 package with lmer(outcome ~ fixed + (1 | group), data = df), where the term in parentheses is the random effect that accounts for clustering, such as repeated measures within people or students within schools. The model that holds up depends on choosing random intercepts versus random slopes correctly, getting p-values the right way, and checking the residuals. This guide explains the formula syntax and the decisions that go with it.

When you need a mixed model

A linear mixed effects model, also called a multilevel or hierarchical model, is for data with structure that violates the independence assumption of ordinary regression: measurements nested within participants, pupils within classrooms, patients within hospitals. Ignoring that nesting, as a plain linear regression in R would, understates the standard errors. For a purely within-subjects design with balanced data, a repeated measures ANOVA in R is a simpler option, but a mixed model tolerates the missing time points and unequal cluster sizes that ANOVA cannot.

Fixed effects and random effects

Fixed effects are the predictors whose specific influence you want to estimate and report. Random effects capture variation across the levels of a grouping variable you are treating as a sample from a larger population, such as the particular schools in your study. The split is what lets the model generalize beyond the exact clusters you observed.

Fitting a random-intercept model

Load lme4 and fit the model. The (1 | school) term gives each school its own intercept:

library(lme4)
library(lmerTest)   # adds p-values to the summary

model <- lmer(score ~ hours + ses + (1 | school), data = df)
summary(model)

Base lme4 deliberately omits p-values; loading lmerTest adds them using Satterthwaite degrees of freedom, which is the standard way to report fixed effects in a dissertation.

Random intercepts versus random slopes

A random intercept lets the baseline level vary by group. A random slope also lets the effect of a predictor vary by group, which you should include when there is reason to think the relationship itself differs across clusters:

# random intercept and random slope for hours, by school
model2 <- lmer(score ~ hours + ses + (1 + hours | school), data = df)

Add random slopes when theory supports them and the model still converges. Overloading the random structure on small data is the usual cause of the singular fit warning, which signals the model is too complex for the data to support.

Testing whether terms matter

Compare nested models with a likelihood-ratio test. To test fixed effects, fit both models with maximum likelihood (REML = FALSE) before comparing them:

m_full <- lmer(score ~ hours + ses + (1 | school), data = df, REML = FALSE)
m_null <- lmer(score ~ ses + (1 | school), data = df, REML = FALSE)
anova(m_null, m_full)   # does 'hours' improve the model?

For the final reported model, refit with the default REML = TRUE, which gives less biased variance estimates. The reasoning behind comparing models rather than chasing significance is the same one in understanding p-values.

How much variance the clustering explains

Quantify the clustering with the intraclass correlation, the share of variance sitting between groups, and the overall fit with a marginal and conditional R-squared:

library(performance)
icc(model)        # variance between groups vs total
r2(model)         # marginal (fixed) and conditional (fixed + random)
check_model(model)  # residual and assumption diagnostics

A near-zero ICC suggests the clustering adds little and a simpler model might suffice; a substantial ICC confirms the mixed model was the right call.

Generalized mixed models for non-normal outcomes

When the outcome is binary or a count rather than continuous, switch from lmer() to glmer() with a family, the multilevel counterpart of logistic regression in R:

model_bin <- glmer(passed ~ hours + (1 | school),
                   data = df, family = binomial)

Reporting a mixed effects model

Report each fixed effect as its estimate, standard error, and p-value with a confidence interval, state the random-effects structure in words (random intercepts for school, say), and give the ICC or variance components. Note the estimation method and the package. The formatting conventions are in reporting statistics in APA style, and the prose belongs in your paper's results section.