Multilevel Regression (HLM)

The Basics

When to Use It

Use it when you have one continuous outcome measured on people who come in groups (students in classrooms, clients in therapy groups, employees in teams, patients in clinics), with predictors measured on the people, on the groups, or both. It's the regression to use for clustered data (Laird & Ware, 1982; Raudenbush & Bryk, 2002): it estimates each predictor's coefficient while treating the groups as a sample, so group-level predictors are judged against the number of groups rather than the number of people. For example:

Are students less anxious about math in classrooms where teachers give more growth-mindset feedback?
Do clients in therapy groups with stronger group cohesion report fewer depression symptoms at the end of treatment, beyond their own starting scores?
Do nurses on wards with more supportive supervisors report less burnout on the Maslach Burnout Inventory, after accounting for their own workload?

The model on this page gives each classroom its own intercept but one shared grade slope. A random slope lets that slope differ by classroom too, and a classroom-level predictor of the slope (is the link between grades and anxiety weaker where feedback is high?) enters as a cross-level interaction, grade × feedback. Leaving out a slope variance that's really there makes that slope's standard error too small, so test it by comparing the model with and without it. Here it was effectively zero (χ²(2) = 0.04, p = .98, roughly, since a variance tested at zero sits on a boundary), so we kept the simpler model.

Vocabulary You Need

Definition Fixed Effect (

b

)

A coefficient that a multilevel (mixed-effects) model estimates as one value for everyone: the overall intercept, or the slope of a predictor averaged across groups. Fixed effects are read like ordinary regression coefficients (the difference in the predicted outcome per unit of a predictor, holding the others constant), and they're what research questions are usually about. A predictor measured on the groups themselves (a teacher's practice, a clinic's size) can only have a fixed effect, since it doesn't vary inside a group. HLM notation writes them γ (“gamma”); many papers and programs write b or β. On this site β means a standardized coefficient. They contrast with random effects, which describe how groups vary around these averages. Glossary →

Definition Random Effect

In a multilevel model, a quantity allowed to differ from group to group, such as each classroom's own intercept or each person's own slope, modeled as a random draw from a population of groups rather than estimated as a separate parameter per group. The model estimates the variance of that distribution (a variance component) instead of dozens of group-specific coefficients, which lets results generalize to groups like the ones sampled, borrows strength across groups, and lets predictors measured on the groups be tested alongside it. Grouping factors are usually treated as random when their levels are a sample (classrooms, therapists, participants) and there are enough of them to estimate a variance; a handful of levels chosen on purpose (three specific clinics) are usually better entered as fixed predictors. Glossary →

Definition Random Intercept (

u_{0j}

)

A random effect that lets each group have its own intercept: its own level of the outcome (on the model's scale, such as the log-odds in a logistic model), after the predictors are taken into account. Group j's intercept is the overall intercept plus u₀ⱼ (often written u_j), its deviation, and the model estimates how much those deviations vary (their variance, written τ₀₀ or τ²). A random-intercept model assumes the predictors' slopes are the same in every group, so the groups' regression lines are parallel, shifted up or down. It's the simplest multilevel model; if a predictor's slope varies between groups, that predictor usually also needs a random slope for its standard error to be right. Glossary →

Definition Random Slope (

u_{1j}

)

A random effect that lets a lower-level predictor's slope differ from group to group: in one classroom prior grades might predict anxiety strongly, in another hardly at all. Group j's slope is the fixed (average) slope plus u₁ⱼ, and the model estimates the variance of those deviations and usually their correlation with the random intercepts. A classroom-level predictor of why slopes differ enters the model as a cross-level interaction. In a growth model the groups are people and the predictor is time, so each person gets their own rate of change, and a person-level predictor of that rate, such as treatment, is the time × predictor interaction. Leaving out a slope variance that really exists can make the standard error of that predictor's fixed effect too small; including one the data can't support often produces a warning that the fit is singular (a variance estimated at zero or a correlation at ±1). Glossary →

Definition Variance Component

In a multilevel model, the outcome's variance is split into parts, one per source of variation: in a two-level model, the variance of the group intercepts (τ₀₀, how much groups differ) and the residual variance within groups (σ², how much individuals differ from their own group's average), plus any slope variances. Each part is a variance component. From the model with no predictors, the groups' component divided by the total is the intraclass correlation, the share of variation that lies between groups. Adding predictors usually shrinks the components, and how much each one shrinks is a common, if rough, measure of variance explained at that level. Variance components can't be negative; an estimate of exactly zero means the data show no variation at that level beyond what chance would produce. Glossary →

Definition Intraclass Correlation (Clustering) (

\text{ICC}

)

In clustered data (students in classrooms, employees in teams, weekly scores within people), the between-cluster variance divided by the total variance: τ₀₀ ÷ (τ₀₀ + σ²) from a model with no predictors. It is also the expected correlation between two members of the same cluster. Even a small ICC matters when clusters are large, because the design effect, 1 + (m − 1) × ICC for clusters of size m, grows with cluster size; ignoring it makes standard errors too small. It describes how the data are structured, not how reliable a measure is, so reliability labels such as “poor” or “excellent” don't apply to it. The formula is the same as the one-way reliability ICC, but the question is different. Glossary →

Definition Design Effect (

\text{deff}

)

The factor by which clustering inflates the sampling variance of an estimate compared with a simple random sample of the same size. For a mean, with an average cluster size m and an intraclass correlation ρ, it's about 1 + (m − 1)ρ. Dividing the sample size by it gives the effective sample size: how many independent observations the clustered data are worth. Its square root says how much too small a standard error is when the clustering is ignored. It grows with cluster size as well as with the intraclass correlation, so even a small ρ matters when clusters are large. The same inflation applies to predictors measured at the group level; predictors that vary within groups are affected much less. Glossary →

Definition Restricted Maximum Likelihood (

\text{REML}

)

A version of maximum likelihood for mixed-effects models that estimates the variance components from the part of the data left after the fixed effects are accounted for. Plain maximum likelihood treats the fixed effects as if they were known, so it underestimates variances, much as dividing by n instead of n − 1 does; REML corrects that, which matters most when there are few groups. It's the default in lme4, SPSS MIXED and most programs. The catch: REML likelihoods from models with different fixed effects can't be compared, so likelihood ratio tests or AICs for fixed effects need models refitted by ordinary maximum likelihood. Comparing models that differ only in their random effects is fine under REML. Glossary →

Assumptions

1. Independent Observations The people (or other units) in your study don't influence each other.

What It Means: Each person's data comes from a separate, unconnected source, so knowing one person's scores tells you nothing about anyone else's. In designs that measure people more than once, their own scores are meant to be related; it's different people who must be independent. Methods built for clustered data, such as multilevel models and GEE, also allow members of the same group to be related, as long as the analysis names the grouping; then it's the groups that must be independent.
Applies To: Observations
How to Check: By design, not by a test. Did each person contribute their own data, with no clusters such as roommates, couples, therapy groups or classrooms whose members influence each other? If there are clusters, does the analysis model every one of them?
How Much It Matters: The model handles the dependence it's told about: students in the same classroom. It still needs the classrooms to be independent of each other, and students' residuals to be independent once their classroom is known. If classrooms come from a handful of schools, the classrooms in a school are alike too, and the classroom-level standard errors come out too small, just as ignoring classrooms made them too small here.

If It's Violated

Correction Fit a three-level model (students in classrooms in schools), or enter the schools as fixed predictors if there are only a few.
Design Change Draw classrooms from different schools, one or two per school.

2. Normally Distributed Random Effects The groups' deviations from the average are roughly bell-shaped.

What It Means: A multilevel model treats each group's intercept (and any random slope) as a draw from a normal distribution with a mean of zero and a variance the model estimates. That's what lets it estimate a variance component from a sample of groups and pull unreliable group estimates toward the average. Groups whose deviations are skewed, bimodal or include a few extreme groups break it.
Applies To: Residuals
How to Check: Plot the model's predicted group effects (the empirical Bayes estimates, ranef() in R) in a histogram or Q-Q plot and look for strong skew, two humps or a group far from the rest; then look at what makes any outlying group different. With few groups, the plot says little, so reason about how the groups were sampled too.
How Much It Matters: Usually little for the fixed effects and their tests; a simulation study found mixed-model estimates robust to skewed and bimodal random effects, with slightly inflated variance estimates when groups really fall into two clusters (Schielzeth et al., 2020). What matters more is a single extreme classroom, which with only 40 can pull the classroom-level slope around. Our 40 estimated classroom effects look roughly symmetric, with no classroom far from the rest.

If It's Violated

Robust Estimation Report cluster-robust ("sandwich") standard errors for the fixed effects, which work from how much the classrooms' residuals actually vary, or bootstrap whole classrooms. Either is a check that doesn't lean on the shape.
Proceed with Caveat Refit without any classroom that stands far from the rest, and report both if the conclusion changes.

3. Normality The scores the test relies on are roughly bell-shaped, or your sample is big enough that it doesn't matter.

What It Means: The test's model assumes the scores it relies on are roughly normal: the scores within each group (or the scores themselves, for a one-sample test), the difference scores in a paired design, or the residuals (the errors left around a model's predictions) in a regression. For p-values and confidence intervals, what really matters is that the estimates (such as a mean, a difference between means or a model's coefficients) are roughly normal across repeated samples. A large enough sample usually provides that; how large depends on how skewed the scores are and how extreme their outliers.
Applies To: Residuals
How to Check: Look at a histogram or Q–Q plot for strong skew or extreme outliers. Don't rely on a normality test: with small samples it misses real problems, and with large ones it flags trivial ones.
How Much It Matters: Little, with 800 students. The assumption is about the students' residuals within classrooms, not about raw anxiety scores, and the coefficients' sampling distributions are close to normal with this many students. For feedback, a classroom-level predictor, the count that matters is the 40 classrooms, not the 800 students. Our scaled residuals run from −2.41 to 3.02, roughly symmetric.

If It's Violated

Transformation Log-transform a strongly right-skewed outcome such as reaction time. Coefficients are then on the log scale.
Robust Estimation Bootstrap whole classrooms for confidence intervals that don't assume a residual shape.

4. Equal Spread of Residuals Your scores scatter about equally around the predictions, whatever the prediction.

What It Means: The residuals have about the same spread at every predicted value, so the model misses by about as much for people predicted to score low as for people predicted to score high. The opposite, heteroscedasticity, often looks like a funnel. If the model's form is right, unequal spread doesn't bias the coefficients, though it makes them less precise than they could be; it does make their standard errors, and so the p-values and confidence intervals, unreliable.
Applies To: Residuals
How to Check: Plot the residuals against the predicted values (in the SPSS REGRESSION procedure, *ZRESID against *ZPRED; for MIXED, save them with /SAVE=PRED RESID and plot those). A band of even width is fine; a funnel or fan that widens or narrows is not. Don't use a test such as Breusch–Pagan as a gatekeeper: with small samples it misses real problems, and with large ones it flags trivial ones.
How Much It Matters: Moderate. The model assumes students vary as much around their classroom's average in every classroom and at every predicted level. If some classrooms are far more spread out than others, the standard errors can be off, most for the predictors those classrooms differ on. Plot the residuals against the predicted values and compare each classroom's residual spread.

If It's Violated

Correction Let the residual variance differ by a grouping variable (nlme's weights argument in R), or add a random slope if the spread fans out with a predictor.
Robust Estimation Cluster-robust standard errors don't assume equal spread.

5. Linear Relationship The relationship between your variables follows a straight line, not a curve.

What It Means: As one variable goes up, the other tends to change by about the same amount across its whole range. A curved relationship, one that levels off or a U shape, breaks this.
Applies To: Relationship
How to Check: Look at a scatterplot of each pair of variables the analysis relates: a roughly oval cloud tilted one way is fine; a clear bend, a ceiling or a U is not. In a regression with several predictors, plot the residuals against the predicted values instead: a flat, even band is fine; a curve means the model is missing a bend. A component-plus-residual plot for each predictor (R: car::crPlots) shows whether one of them bends once the others are held constant.
How Much It Matters: A lot, for what the slopes mean. Each fixed effect is a straight-line slope; if anxiety drops steeply at low grades and levels off at high ones, one slope misdescribes both ends. Check the residuals against each predictor, at both levels.

If It's Violated

Correction Add a squared term for the predictor, centered first.
Transformation Log-transform a predictor whose effect levels off.

Limitations

This Is Not Your Test If…

Your outcome is yes/no, like passing or failing Multilevel Logistic Regression →
You measured the same people repeatedly to model change Multilevel Growth-Curve Modeling →
You want average effects, not variation between groups Generalized Estimating Equations (GEE) →
Your outcome is a latent variable measured by several items Multilevel SEM →
You have only a handful of groups, like three clinics Multiple Linear Regression →
Your people aren't grouped at all Multiple Linear Regression →

What It Doesn't Tell You

Why. Classrooms weren't assigned to more or less feedback. Teachers who give growth-mindset feedback may differ in many other ways (warmth, experience, the schools they work in), and any of those could lower anxiety instead.
Whether it matters. A point of feedback goes with 0.82 fewer points of anxiety on a 9–45 scale, and 95% of the variance in anxiety lies beyond what the fixed effects predict. Significance doesn't make that worth a training program.
Who else. The result describes classrooms like these 40, and with only 40 the classroom-level estimates are rough. It says nothing about how feedback works for any one student, since feedback was measured on classrooms, not students.
That classroom mean grade doesn't matter. Its non-significant coefficient (p = .907) means these 40 classrooms show no clear link, with an interval from −0.15 to 0.16 points per grade point; it isn't evidence that a class's overall ability is irrelevant.

Common Mistakes

Don'tIgnore the classrooms and run an ordinary regression on all the students.

InsteadFit a multilevel model. Here, ignoring classrooms shrank feedback's standard error from 0.35 to 0.19 and its p from .025 to below .001.
Don'tTreat the number of students as the sample size for a classroom-level predictor.

InsteadCount classrooms. With 40 classrooms, feedback's test has about 37 degrees of freedom, however many students there are.
Don'tEnter a student-level predictor uncentered and read its slope as a within-classroom effect.

InsteadCenter it within groups and add the group means, or say that the slope blends within- and between-group relationships.
Don'tCompare models with different fixed effects using their REML likelihoods.

InsteadRefit both by maximum likelihood for a likelihood ratio test, or use the fixed effect's t test from the REML fit.
Don'tReport “the R²” of a multilevel model without saying which one.

InsteadName it (marginal R², or the reduction in each variance component) and say that different definitions give different answers.
Don'tAdd random slopes for every predictor because the software allows it.

InsteadAdd a random slope when theory says a slope should vary, and drop it if the fit is singular (a variance estimated at zero or a correlation at ±1) or the slope variance is near zero.
Don'tTreat a grouping variable as random when it has only a handful of levels.

InsteadWith very few groups, the variance between them is estimated badly; enter the groups as fixed predictors instead.
Don'tConclude clustering doesn't matter because the intraclass correlation is small.

InsteadWork out the design effect. With 30 people per group, an ICC of .05 already inflates a mean's variance 2.45 times.
Don'tRead a p-value as the chance that the null hypothesis is true.

InsteadRead it as: if the null were true, how often data at least this extreme would turn up.

Result Decoder

Significance Guidelines

Benchmarks Significance (p Value)

p < .001
Very strong evidence against the null hypothesis.
.001 ≤ p < .01
Strong evidence against the null hypothesis.
.01 ≤ p < .05
Statistically significant at α = .05. Moderate evidence against the null hypothesis.
.05 ≤ p < .10
Not significant. Weak evidence; sometimes called a “trend.” Treat it as inconclusive.
p ≥ .10
Not significant. Little evidence against the null, which is not the same as evidence of no effect.

These evidence labels are one common convention for describing p in graded language rather than as significant or not; the boundaries are rough guides, not rules, and other authors draw them differently. Muff et al. (2022) argue for this kind of language of evidence over a single cutoff. α is chosen before the analysis, and p never says how big or important an effect is.

More about statistical significance ↓

How It's Reported

Equation Form

Example

b = −0.82, 95% CI [−1.53, −0.11], β = −0.15, t(36.85) = −2.34, p = .025

Table 1. Fixed Effects Predicting Math Anxiety

Example

	$b$	$SE$	$df$	$t$	$p$	95% CI
Intercept	21.93	0.420	37.00	52.21	< .001	[21.08, 22.78]
Prior grade, relative to classmates	−0.11	0.023	759.17	−4.76	< .001	[−0.16, −0.06]
Classroom mean prior grade	0.01	0.077	36.93	0.12	.907	[−0.15, 0.16]
Growth-mindset feedback	−0.82	0.351	36.85	−2.34	.025	[−1.53, −0.11]

Note. N = 800 students in 40 classrooms. Random-intercept model fit by REML (R's lme4 2.0.6 with lmerTest 3.2.1); df by Satterthwaite's method. Prior grade is centered within each classroom; the classroom mean grade and feedback are centered at the average classroom (75.30 and 3.86). Feedback is the effect worked through in the Napkin Math.

Table 2. Variance Components

Example

	Empty Model	Full Model	Full Model $SD$	Proportion Reduced
Between classrooms ( $\tau_{00}$ )	5.78	5.19	2.28	.10
Students within classrooms ( $\sigma^2$ )	37.56	36.52	6.04	.03

Note. Empty model: no predictors; its classroom share, 5.78 ÷ (5.78 + 37.56) = .13, is the intraclass correlation. Proportion reduced: (empty − full) ÷ empty, a rough measure of how much variance at each level lines up with the predictors. Marginal R² (fixed effects only) = .05, conditional R² (fixed and random effects together) = .16 (Nakagawa & Schielzeth, 2013).

Full Text (APA Style)

Example

Students' math anxiety averaged 22.00 on the 9–45 scale (SD = 6.58), and 13% of its variance lay between classrooms (ICC = .13, design effect = 3.53 with 20 students per classroom), so we fit a two-level random-intercept model by REML (lme4 and lmerTest in R), with Satterthwaite degrees of freedom. Prior grade was centered within classrooms, and each classroom's mean grade was entered at the classroom level. Students in classrooms with more growth-mindset feedback were less anxious: each point on the 1–7 feedback rating went with 0.82 points lower average anxiety, 95% CI [−1.53, −0.11], β = −0.15, holding prior grades constant, t(36.85) = −2.34, p = .025. Within classrooms, students with higher grades than their classmates were less anxious (b = −0.11 per grade point, 95% CI [−0.16, −0.06], t(759.17) = −4.76, p < .001), while a classroom's mean grade showed no clear link with its average anxiety (b = 0.01, 95% CI [−0.15, 0.16], p = .907). The predictors reduced the classroom variance by 10% (from 5.78 to 5.19) and the student variance by 3%; together the fixed effects predicted 5% of the total variance (marginal R² = .05). Table 1 gives every fixed effect and Table 2 the variance components.

Give the ICC from the empty model, the number of groups and the group sizes, so readers can judge how much the clustering matters.
Say how the model was fit (REML or ML), which random effects it has, how the predictors were centered, and how the degrees of freedom were found (Satterthwaite or Kenward-Roger).
Report the fixed effects in a table (b, SE, df, t, p, 95% CI) and the variance components in another, as in Tables 1 and 2.
Name the "variance explained" measure you use; there are several, and they disagree.

Reading Each Piece

The equation form above, piece by piece, in this order: how big, how certain, then how surprising.

Marginal R² = .05	The fixed effects together predict 5% of the variance in anxiety, students and classrooms combined. Most of what separates anxious from calm students lies within classrooms, and a classroom-level predictor like feedback can only line up with differences between classrooms.
b = −0.82, 95% CI [−1.53, −0.11]	Classrooms one point apart in growth-mindset feedback differ by 0.82 points in average anxiety, holding prior grades constant; effects from 0.11 points, barely anything on a 9–45 scale, to 1.53 points are all compatible with the data at the 95% level.
β = −0.15	Classrooms one SD apart in feedback, with the same mean grade, differ by 0.15 of a student SD in average anxiety, about 1 point on the 9–45 scale. That's small for any one student, but it's 0.41 of the SD between classrooms' true averages; whether a point of anxiety per classroom is worth a feedback program is a question about schools, not statistics.
t(36.85) = −2.34	Feedback's coefficient is 2.34 standard errors below zero. The 36.85 is the Satterthwaite degrees of freedom, set mostly by the 40 classrooms.
p = .025	If feedback were unrelated to classroom anxiety (holding prior grades constant), a coefficient at least this far from zero, in either direction, would turn up about 2.5% of the time. That's below .05, so it's statistically significant, but with 40 classrooms, not by much. Benchmarks ↑

Calculations

The Running Example

Example

The math-anxiety study: 800 students in 40 classrooms, math anxiety (9–45) predicted by each student's prior grade and by how often their teacher gave growth-mindset feedback (1–7).

Null Hypothesis:: the population fixed effect of feedback is 0
Alternative (Two-Sided):: the population fixed effect of feedback is not 0

The fixed effect of feedback is how much average anxiety differs between classrooms one point apart in feedback, holding prior grades constant, among classrooms like these 40. Each fixed effect in Table 1 gets the same kind of test; the prior-grade effects are there to separate how students compare with their classmates from how classrooms compare with each other.

What Your Data Should Look Like

One row per person, in long format: a column for the group each person belongs to, a column for the outcome, and a column for each predictor. A group-level predictor, such as the teacher's feedback rating, repeats the same value on every row from that group. Groups don't need to be the same size, and a person missing the outcome drops out without taking their group with them.

Example

student	classroom	math_anxiety	prior_grade	growth_feedback
1	1	12	91	4.0
2	1	31	89	4.0
3	1	25	75	4.0
4	1	20	78	4.0
5	1	16	69	4.0
6	1	22	51	4.0
… 794 more rows

	M	SD	1	2
1. Math Anxiety	22.00	6.58
2. Prior Grade	75.33	10.73	−.14
3. Growth-Mindset Feedback	3.87	1.20	−.15	.03

N = 800.

Correlations are Pearson r.

Simulated data.

Download all 800 rows (data.csv)

Napkin Math

The whole test, one idea at a time. Under each heading, a line says why the step matters, and the green note reads the equation in plain English with the math-anxiety study's numbers. Collapse all to read just the headings: they tell the whole story on their own.

1. Find How Much of the Anxiety Differences Lie between Classrooms

Before adding any predictors, split students' anxiety into two parts, how much classroom averages differ and how much students differ within their classroom. The classrooms' share is the intraclass correlation.

$\rho = \dfrac{\tau_{00}}{\tau_{00} + \sigma^2}$ = $\dfrac{5.78}{5.78 + 37.56}$

= .13

$\tau_{00}$ = 5.78
$\sigma^2$ = 37.56

An empty model, with no predictors, estimates $\tau_{00}$ , how much the classrooms' true average anxiety varies, as 5.78 (the 40 observed averages vary a little more, because each also carries sampling error), and $\sigma^2$ , the variance of students around their own classroom's average, as 37.56. 5.78 ÷ 43.34 = .13: about 13% of the differences in anxiety are differences between classrooms.
2. Work Out How Much Less Information Clustered Students Carry

Students in the same classroom are partly repeating each other. The design effect says how much that inflates the uncertainty of an average, and dividing by it gives how many independent students the data are really worth (Raykov & DiStefano, 2022).

$\text{deff} = 1 + (\bar{m} - 1)\,\rho$ = $1 + (20 - 1) \times .133$

= 3.53

$N_{\text{eff}}$ = 226 students

$\bar{m}$ is our average classroom size, 800 ÷ 40 = 20 students, and $\rho$ is the .133 from step 1. 1 + 19 × .133 = 3.53, so the 800 students carry about as much information about average anxiety as 800 ÷ 3.53 = 226 independent students would.
3. Separate Each Student's Standing in Class from the Class's Standing

A prior grade mixes two things: how a student compares with classmates, and how strong the whole class is. Subtracting each classroom's mean grade gives the first; the classroom means themselves go in as a separate, classroom-level predictor (Enders & Tofighi, 2007).

$x_{ij} - \bar{x}_j$ = $91 - 73.09$

= 17.91 points

$\bar{x}_1$ = 73.09

Student 1 had a 91 last year. Their classroom's 23 students averaged $\bar{x}_1$ = 73.09, so their centered grade is 91 − 73.09 = 17.91: nearly 18 points above their classmates. Their classroom's mean, 73.09, is 2.22 below the average classroom's 75.30, and that goes in the model separately.
4. Fit One Model for Students and Classrooms Together

The model predicts each student's anxiety from predictors at both levels, plus a random intercept that lets each classroom sit a little above or below what its predictors say. Software estimates all of it at once by restricted maximum likelihood (REML; Patterson & Thompson, 1971), so we take its numbers.

$y_{ij} = b_0 + b_1 (x_{ij} - \bar{x}_j) + b_2 (\bar{x}_j - \bar{x}) + b_3 (w_j - \bar{w}) + u_{j} + e_{ij}$ = $21.93 - 0.11\,(x_{ij} - \bar{x}_j) + 0.01\,(\bar{x}_j - 75.30) - 0.82\,(w_j - 3.86) + u_j + e_{ij}$

$b_0$ = 21.93 points
$b_1$ = −0.11
$b_2$ = 0.01
$b_3$ = −0.82

$x_{ij}$ is a student's prior grade, $\bar{x}_j$ their classroom's mean grade and $w_j$ their teacher's feedback rating; $u_j$ is classroom $j$ 's random intercept, how far it sits above or below what its predictors say, and $e_{ij}$ is the student's own residual. $b_0$ = 21.93 is the predicted anxiety for a student at their class average, in an average classroom. $b_1$ = −0.11: each grade point above classmates goes with 0.11 points less anxiety. $b_2$ = 0.01: classrooms with higher mean grades weren't clearly more or less anxious; its interval, −0.15 to 0.16 per grade point, is too wide to rule out a difference either way. $b_3$ = −0.82: classrooms one point apart in feedback differ by 0.82 points in average anxiety, holding mean grade constant. $b_2 - b_1$ = 0.01 − (−0.11) = 0.12 is the contextual effect: how two students with the same grade differ when one is in a class averaging a point higher. Centering the grade at the grand mean instead, with the class mean still in the model, makes $b_2$ this contextual effect. HLM papers write the model as one equation per level: $y_{ij} = \beta_{0j} + \beta_{1j}(x_{ij} - \bar{x}_j) + e_{ij}$ for students, and $\beta_{0j} = \gamma_{00} + \gamma_{01}(\bar{x}_j - 75.30) + \gamma_{02}(w_j - 3.86) + u_j$ and $\beta_{1j} = \gamma_{10}$ (no random slope) for classrooms. Substituting the second into the first gives the equation above, with their γ's as our b's (their β's are coefficients, not standardized ones).
5. See How Much of Each Level's Variation Lines Up with the Predictors

With predictors in the model, the variance components shrink by the part the predictors predict. Comparing them with the empty model gives a rough "variance explained" at each level, where "explained" means predicted statistically, not caused; marginal R squared gives one for the whole model.

$\dfrac{\tau_{00}^{\text{empty}} - \tau_{00}^{\text{full}}}{\tau_{00}^{\text{empty}}}$ = $\dfrac{5.78 - 5.19}{5.78}$

= .10

$\tau_{00}^{\text{full}}$ = 5.19
$\sigma^2_{\text{full}}$ = 36.52
$\text{Reduction}_{\text{students}}$ = .03
$R^2_m$ = .05

The classroom variance $\tau_{00}$ fell from 5.78 in the empty model to 5.19, so (5.78 − 5.19) ÷ 5.78 = 10% of the differences between classrooms line up with feedback and the classroom mean grade. The student variance fell from 37.56 to 36.52, just 3%. Counting both levels, the fixed effects' predictions vary by 2.01, so marginal $R^2$ = 2.01 ÷ (2.01 + 5.19 + 36.52) = .05.
6. Count How Much Information There Is about Classrooms

Feedback varies only from classroom to classroom, so its test rests on 40 classrooms, not 800 students. Satterthwaite's approximation turns that into degrees of freedom by asking how uncertain feedback's standard error is.

$\nu = \dfrac{2\,\left(SE_{b}^{2}\right)^{2}}{\widehat{\text{Var}}\!\left(SE_{b}^{2}\right)}$

= 36.85

$SE_{b_3}$ = 0.351
$SE_{\text{OLS}}$ = 0.190

$SE_b$ = 0.351 for feedback, and its uncertainty comes almost entirely from the classroom variance, which 40 classrooms estimate only roughly. The formula's denominator, $\widehat{\text{Var}}(SE_b^2)$ , is how much that squared standard error would vary from sample to sample: the less certain it is, the fewer df. Satterthwaite's $\nu$ comes out at 36.85, close to 40 classrooms minus the 3 classroom-level coefficients. An ordinary regression that ignored the classrooms would have claimed $SE$ = 0.19 and 796 df.
7. Compare Feedback's Effect with Its Uncertainty

As in any regression, a coefficient only means something relative to its standard error. t counts how many standard errors it sits from zero, and the Satterthwaite df say which t distribution to judge it against.

$t = \dfrac{b}{SE_b}$ = $\dfrac{-0.82}{0.351}$

= −2.34

$p$ = .025

$b$ is feedback's −0.82 points and $SE_b$ is 0.351. −0.82 ÷ 0.351 = −2.34 standard errors. With 36.85 df, the two-sided $p$ is .025.
8. Give a Range of Plausible Values for Feedback's Effect

Instead of a single best guess, give the range of population effects that fit the data, using the same t distribution.

$b \pm t^{*}_{\nu} \times SE_b$ = $-0.82 \pm 2.026 \times 0.351$

= [−1.53, −0.11] points

$t^{*}_{\nu}$ = 2.026 holds the middle 95% of a t distribution with 36.85 df. 2.026 × 0.351 = 0.71, so the range is −0.82 ± 0.71, from −1.53 to −0.11 points per point of feedback.
9. Put Feedback's Effect on a Common Scale

A point on a 1–7 observer rating isn't a natural unit. Scaling the effect by the spread of feedback across classrooms and the spread of anxiety across students gives a size that can be compared with other studies.

$\beta_j = b_j \times \dfrac{s_{x_j}}{s_y}$ = $-0.82 \times \dfrac{1.21}{6.58}$

= −0.15

$s_x$ = 1.21
$\beta_{\text{between}}$ = −0.41

$s_x$ is the 1.21-point SD of feedback across the 40 classrooms (the data box's 1.20 counts each student, so larger classes weigh more) and $s_y$ the 6.58-point SD of students' anxiety. A classroom 1 SD higher in feedback averages 0.82 × 1.21 = 0.99 points less anxiety, 0.99 ÷ 6.58 = 0.15 of a student SD. Against the 2.40-point SD of the classrooms' true averages ( $\sqrt{\tau_{00}} = \sqrt{5.78}$ from the empty model; the 40 observed averages spread a little more, because each also carries sampling error), the same 0.99 points is 0.41.

Run It in Software

The code below reads the running example's data from data.csv; save it in the same folder you run the code from.

SPSS (version 30, as documented by IBM)

First make the centered predictors. Data → Aggregate: break by classroom, and add the mean of prior_grade as a new variable class_mean_grade (keep the option to add it to the active dataset). Then Transform → Compute Variable three times, as in the syntax below: grade_within is each student's grade minus their classroom mean; class_grade_c and feedback_c subtract the average classroom's values (75.30 and 3.86: the mean of the 40 classroom values, not of the 800 students).
Analyze → Mixed Models → Linear. Move classroom into Subjects and click Continue.
Move math_anxiety into Dependent Variable and grade_within, class_grade_c and feedback_c into Covariate(s).
Fixed: add all three covariates as main effects. Random: check Include intercept, move classroom into Combinations, and leave Covariance Type on Variance Components.
Estimation: leave Method on Restricted Maximum Likelihood (REML). Statistics: check Parameter estimates and Tests for covariance parameters. Click OK.
For the intraclass correlation, run it once more with no covariates (the empty model) and divide the intercept variance by the sum of the two variances.

Or paste this syntax:

AGGREGATE
  /OUTFILE=* MODE=ADDVARIABLES
  /BREAK=classroom
  /class_mean_grade=MEAN(prior_grade).
COMPUTE grade_within = prior_grade - class_mean_grade.
* 75.30475 and 3.8625 are the means of the 40 classroom values.
COMPUTE class_grade_c = class_mean_grade - 75.30475.
COMPUTE feedback_c = growth_feedback - 3.8625.
EXECUTE.

* The empty model, for the intraclass correlation.
MIXED math_anxiety
  /FIXED=| SSTYPE(3)
  /METHOD=REML
  /PRINT=SOLUTION TESTCOV
  /RANDOM=INTERCEPT | SUBJECT(classroom) COVTYPE(VC).

* The full model.
MIXED math_anxiety WITH grade_within class_grade_c feedback_c
  /FIXED=grade_within class_grade_c feedback_c | SSTYPE(3)
  /METHOD=REML
  /PRINT=SOLUTION TESTCOV
  /RANDOM=INTERCEPT | SUBJECT(classroom) COVTYPE(VC).

Reading the Output

Start with the empty model's Estimates of Covariance Parameters: the intercept variance (5.78) over the total (5.78 + 37.56) is the intraclass correlation, .13. Then, in the full model, read Estimates of Fixed Effects: each row's Estimate and its 95% Confidence Interval, then t, df and Sig. SPSS's MIXED also gives fractional degrees of freedom for each fixed effect; the table below uses R's Satterthwaite values, and SPSS's should be very close.

Estimates of Fixed Effects (full model)

Parameter	Estimate	Std. Error	df	t	Sig.	95% CI Lower	Upper
Intercept	21.931397	.420040	36.996	52.213	<.001	21.080313	22.782482
grade_within	-.110444	.023191	759.169	-4.762	<.001	-.155970	-.064919
class_grade_c	.008970	.076527	36.933	.117	.907	-.146098	.164038
feedback_c	-.821072	.351308	36.852	-2.337	.025	-1.532987	-.109158

Estimates of Covariance Parameters (full model)

Parameter	Estimate
Residual	36.522552
Intercept [subject = classroom] Variance	5.190324

SPSS also prints each variance's Std. Error and a Wald Z test. The Wald test of a variance is a poor guide, since a variance can't go below zero; to test whether classrooms differ, compare the models with and without the random intercept with a likelihood ratio test.

Output laid out from this page's computed results, not captured from a live run; labels and layout may differ slightly in your version.

R

library(lmerTest)   # loads lme4 and adds Satterthwaite df and p-values

classes <- read.csv("data.csv")

# Center prior grade within each classroom, and enter the classroom mean at level 2
classes$class_mean_grade <- ave(classes$prior_grade, classes$classroom)
classes$grade_within <- classes$prior_grade - classes$class_mean_grade
one_row <- classes[!duplicated(classes$classroom), ]   # one row per classroom
classes$class_grade_c <- classes$class_mean_grade - mean(one_row$class_mean_grade)
classes$feedback_c <- classes$growth_feedback - mean(one_row$growth_feedback)

# The empty model: how much variance lies between classrooms?
empty <- lmer(math_anxiety ~ 1 + (1 | classroom), data = classes)
print(VarCorr(empty), comp = "Variance")
vc <- as.data.frame(VarCorr(empty))$vcov
vc[1] / sum(vc)   # ICC

# The full model, by REML
fit <- lmer(math_anxiety ~ grade_within + class_grade_c + feedback_c + (1 | classroom),
            data = classes)
summary(fit)

# 95% CI for feedback with the Satterthwaite df
contest1D(fit, c(0, 0, 0, 1), confint = TRUE)

 Groups    Name        Variance
 classroom (Intercept)  5.7805
 Residual              37.5600
[1] 0.1333738

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: math_anxiety ~ grade_within + class_grade_c + feedback_c + (1 |
    classroom)
   Data: classes

REML criterion at convergence: 5207.3

Scaled residuals:
     Min       1Q   Median       3Q      Max
-2.41250 -0.71721 -0.01332  0.65196  3.01582

Random effects:
 Groups    Name        Variance Std.Dev.
 classroom (Intercept)  5.19    2.278
 Residual              36.52    6.043
Number of obs: 800, groups:  classroom, 40

Fixed effects:
               Estimate Std. Error        df t value Pr(>|t|)
(Intercept)    21.93140    0.42004  36.99624  52.213  < 2e-16 ***
grade_within   -0.11044    0.02319 759.16929  -4.762 2.29e-06 ***
class_grade_c   0.00897    0.07653  36.93254   0.117    0.907
feedback_c     -0.82107    0.35131  36.85176  -2.337    0.025 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) grd_wt clss__
grade_withn  0.000
class_grd_c -0.001  0.000
feedback_c  -0.003  0.000 -0.049

    Estimate Std. Error       df   t value     lower      upper   Pr(>|t|)
1 -0.8210723  0.3513081 36.85176 -2.337186 -1.532987 -0.1091579 0.02497259

Read the Random effects block for the variance components, then each row of Fixed effects: Estimate, Std. Error, the Satterthwaite df, t value and Pr(>|t|). The last line gives feedback's 95% interval with the same df. Plain lme4::lmer prints no df or p; loading lmerTest adds them. confint(fit) gives profile-likelihood intervals instead, which are slightly different. To test fixed effects by a likelihood ratio test, compare models with anova(), which refits them by maximum likelihood first; performance::r2_nakagawa(fit) gives the marginal and conditional R².

Output from a real run: R 4.6.1 (lme4 2.0.6, lmerTest 3.2.1), run 2026-09-29.

Python

import pandas as pd
import statsmodels.formula.api as smf

classes = pd.read_csv("data.csv")

# Center prior grade within classrooms; classroom-level variables at the classrooms' average
classes["class_mean_grade"] = classes.groupby("classroom")["prior_grade"].transform("mean")
classes["grade_within"] = classes["prior_grade"] - classes["class_mean_grade"]
one_row = classes.drop_duplicates("classroom")
classes["class_grade_c"] = classes["class_mean_grade"] - one_row["class_mean_grade"].mean()
classes["feedback_c"] = classes["growth_feedback"] - one_row["growth_feedback"].mean()

empty = smf.mixedlm("math_anxiety ~ 1", classes, groups=classes["classroom"]).fit(reml=True)
tau, sigma2 = empty.cov_re.iloc[0, 0], empty.scale
print(f"ICC = {tau / (tau + sigma2):.3f}")

fit = smf.mixedlm("math_anxiety ~ grade_within + class_grade_c + feedback_c",
                  classes, groups=classes["classroom"]).fit(reml=True)
print(fit.summary())

ICC = 0.133
          Mixed Linear Model Regression Results
==========================================================
Model:            MixedLM Dependent Variable: math_anxiety
No. Observations: 800     Method:             REML
No. Groups:       40      Scale:              36.5226
Min. group size:  15      Log-Likelihood:     -2603.6543
Max. group size:  25      Converged:          Yes
Mean group size:  20.0
----------------------------------------------------------
                Coef.  Std.Err.   z    P>|z| [0.025 0.975]
----------------------------------------------------------
Intercept       21.931    0.420 52.208 0.000 21.108 22.755
grade_within    -0.110    0.023 -4.762 0.000 -0.156 -0.065
class_grade_c    0.009    0.077  0.117 0.907 -0.141  0.159
feedback_c      -0.821    0.351 -2.337 0.019 -1.510 -0.133
Group Var        5.190    0.277
==========================================================

The estimates and standard errors match R's. The tests don't: statsmodels judges each coefficient against the normal distribution (z), as if there were unlimited classrooms, so feedback's p is .019 and its interval [−1.51, −0.13], a little too narrow. With 40 classrooms, use R's Satterthwaite t (p = .025) or compute it from these numbers with 36.85 df. Scale is the student-level variance $\sigma^2$ and Group Var the classroom variance $\tau_{00}$ .

Output from a real run: Python 3.9.6 (statsmodels 0.14.6, pandas 2.3.3), run 2026-10-03.

Impact

Effect Size

A multilevel model has no single effect size. The fixed effect b gives each predictor's size in the outcome's own units, which is usually the most honest summary. A standardized coefficient puts it on a common scale, but it depends on which standard deviations you divide by, since the outcome now varies at two levels.

\beta = b \times \dfrac{s_x}{s_y} = -0.82 \times \dfrac{1.21}{6.58} = -0.15

A classroom one standard deviation higher in feedback (1.21 points on the 1–7 rating) averages 0.99 points less anxiety: 0.15 of the spread between students, or 0.41 of the spread between classrooms' true averages ( $\sqrt{\tau_{00}}$ , 2.40 points). Same difference, two sizes, because classrooms differ much less than students do.

"Variance explained" splits too ("explained" here means predicted statistically, not caused). Marginal R² is the share of the total variance the fixed effects predict (Nakagawa & Schielzeth, 2013); ours is .05. The drop in each variance component from the empty model gives a share per level: 10% of the classroom variance and 3% of the student variance. These level-by-level shares are rough. Adding a predictor can make one of these shares shrink, or even turn negative, because the variance components are estimated, not partitioned like sums of squares, so say which measure you report.

Whether a classroom-level difference of about a point on a 9–45 anxiety scale matters is a question about students and schools, not statistics.

There are no conventional benchmarks for a standardized coefficient in a multilevel model, and it depends on which standard deviations you divide by. Judge the size in the outcome's own units, and say which standard deviations you used.

Statistical Significance

Each fixed effect's t test asks: if this predictor's population coefficient were zero, holding the others constant, how often would its coefficient land at least this far from zero, in either direction? Below the conventional α = .05, the result is called statistically significant. The catch is the reference distribution. Variance components are estimated, so the t statistic doesn't follow an exact t distribution; Satterthwaite's approximation, lmerTest's default, picks the degrees of freedom from how uncertain the standard error is (Kuznetsova et al., 2017). The lme4 authors call such approximations "at best ad hoc," but they work far better than treating t as z when there are few groups.

If feedback were unrelated to classroom anxiety (holding prior grades constant), a coefficient at least this far from zero, in either direction, would turn up about 2.5% of the time: t(36.85) = −2.34, p = .025. Had we ignored the classrooms, ordinary regression would have reported p < .001, nearly 1 in 90,000, because it would have treated 800 students as 800 independent pieces of evidence about 40 teachers. Judged against a normal distribution instead of t (as statsmodels does), p would be .019.

Two more tests turn up in multilevel work. To compare models that differ in their fixed effects, use a likelihood ratio test on models refitted by maximum likelihood, not REML: for feedback, dropping it raises the deviance by 5.50, χ²(1) = 5.50, p = .019. To ask whether groups differ at all, compare the model with and without the random intercept: here χ²(1) = 49.21, p < .001, so the classrooms clearly differ. Because a variance can't be negative, the null value sits on the edge of what's possible, and the ordinary χ²(1) p is about twice too large; halve it (Self & Liang, 1987). Here that makes no difference.

p is not the probability that feedback is useless, and it says nothing about how much feedback matters. With 40 classrooms, a p of .025 also means the estimate is imprecise: the interval runs from 0.11 to 1.53 points.

Benchmarks p at a Glance

Very Strong

p < .001

Strong

.001 ≤ p < .01

Significant

.01 ≤ p < .05

Weak

.05 ≤ p < .10

Little

p ≥ .10

Strength of evidence against the null hypothesis.

Full benchmarks and caveats ↑

Power and Sample Size

Power is the chance that a study detects an effect that's really there. In a multilevel study it depends on the size of the effect, α, the intraclass correlation, and two sample sizes: the number of groups and the number of people per group. For a group-level predictor like feedback, the number of groups matters far more. Adding students to each classroom only sharpens each classroom's average, and once those averages are sharp, only more classrooms help (Snijders, 2005; Hedges & Rhoads, 2010). Few groups also cost accuracy, though less than you might fear: in Maas and Hox's (2005) simulations, the regression coefficients and variance components were estimated accurately, and only with 50 or fewer groups were the standard errors at the group level biased. Our 40 classrooms fall inside that range, so treat the classroom-level standard errors as somewhat approximate.

For the simplest case there's a formula. With one group-level predictor, a random intercept and equal groups, the test works like a t-test on the groups' averages, with J − 2 degrees of freedom for J groups: each group of m people counts for m ÷ [1 + (m − 1) × ICC] independent people, the design effect again. For anything more complex, such as random slopes, several correlated predictors or unequal groups, power is usually found by simulation: generate many datasets from the model you expect, fit the model to each, and count how often the effect comes out significant (Kumle et al., 2021). R's simr package does this for lme4 models.

Say a new study wants to detect an effect the size of ours: 0.82 points per point of feedback, or 0.99 points per standard deviation of feedback across classrooms (1.21). With our variance components (5.19 between classrooms, 36.52 within), anxiety's SD among classrooms with the same feedback is √(5.19 + 36.52) = 6.46 and the ICC is .12, so the standardized effect is 0.99 ÷ 6.46 = 0.15. With 40 classrooms of 20 students, that study would have only about 63% power. Doubling each class to 40 students raises it to 69%; even 1,000 students per class would give just 76%. Adding classrooms works: 60 classrooms of 20 give 81%, and 80 give 91%. For 80% power, plan on about 60 classrooms of 20, 52 of 40, or 75 of 10.

A simulation of 2,000 studies with these values, each fitted with lme4 and tested with lmerTest, agreed: 63% at 40 classrooms and 80% at 60. Our estimate is itself rough (its interval runs from 0.11 to 1.53 points), so when planning, try smaller effects too. And don't compute power at the effect you just found to explain your own result: that "observed power" only restates p and adds nothing (Hoenig & Heisey, 2001).

Researchers usually aim for at least 80% power. In multilevel studies, plan the number of groups first, for the smallest group-level effect that would matter, using an intraclass correlation from earlier studies; then decide how many people per group. The calculator below covers one group-level predictor, or two conditions assigned to whole groups, in a random-intercept model; for random slopes or several correlated predictors, simulate. Try it:

Other Info

Follow-Up and Post Hoc Tests

None needed. This test answers a single question, so there are no planned contrasts or post hoc tests to run.

Related Methods

General CaseMultilevel SEM: a multilevel regression is a multilevel SEM with one observed outcome and no latent variables. Multilevel SEM adds measurement models and separate paths within and between groups.
Related DesignMultiple Linear Regression: the one-level version. When the groups don't differ at all (a classroom variance of zero), the fixed effects and their standard errors match an ordinary regression's; the more groups differ, the more ordinary regression understates the uncertainty of group-level predictors.
Related DesignGeneralized Estimating Equations (GEE): fits the same clustered data without modeling the groups: it estimates average coefficients with cluster-robust standard errors and reports no variance components.
Related DesignMultilevel Logistic Regression: the same model for a yes/no outcome, with the group variance on the log-odds scale.
Related DesignMultilevel Growth-Curve Modeling: the same model with time as the level-1 predictor and people as the groups, so each person gets their own starting point and, with a random slope, their own rate of change.
Special CaseFactorial Repeated-Measures ANOVA: A mixed model with random intercepts for each person and each person-by-factor combination can reproduce the factorial repeated-measures ANOVA's F-tests with complete, balanced data and compatible variance and df assumptions; other covariance structures or estimation methods can differ. The ANOVA can't keep people with missing cells or be fit to individual trials; this model can, though keeping incomplete people doesn't by itself remove bias from dropout: valid results still depend on why data are missing and on the model's assumptions.
Special CaseIntraclass Correlation Coefficient (ICC): A reliability ICC is the targets' share of the total variance, the same ratio a random-intercept model estimates with ratings nested in targets. With every target rated by every rater the two agree; the multilevel model can also handle missing ratings.
Special CaseMixed Factorial ANOVA: With a between-subjects group factor, a within-subjects time factor, complete, balanced data and compatible covariance, estimation and df assumptions, a model with a random intercept for each person can reproduce a mixed factorial ANOVA's F tests; other covariance structures or df conventions can differ. Unlike the ANOVA, it keeps people who miss a time point, though that doesn't by itself remove bias from dropout.
Special CaseOne-Way Repeated-Measures ANOVA: With every person measured in every condition, a model with condition as a predictor and a random intercept for each person gives the same F as an uncorrected one-way repeated-measures ANOVA, because both assume equal variances and correlations across conditions. Unlike the ANOVA, a mixed model can relax that with a more flexible covariance structure and keep people with missing time points; its results then depend on the model and estimation method.

References

Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 38(4), 963–974. https://doi.org/10.2307/2529876
“A general family of models is discussed, which includes both growth models and repeated-measures models as special cases.” (Abstract)
Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical linear models: Applications and data analysis methods (2nd ed.). Sage.
Patterson, H. D., & Thompson, R. (1971). Recovery of inter-block information when block sizes are unequal. Biometrika, 58(3), 545–554. https://doi.org/10.1093/biomet/58.3.545
Enders, C. K., & Tofighi, D. (2007). Centering predictor variables in cross-sectional multilevel models: A new look at an old issue. Psychological Methods, 12(2), 121–138. https://doi.org/10.1037/1082-989X.12.2.121
“The article provides a number of practical recommendations designed to facilitate centering decisions in MLM applications.” (Abstract)
Bell, A., Jones, K., & Fairbrother, M. (2018). Understanding and misunderstanding group mean centering: A commentary on Kelley et al.'s dangerous practice. Quality & Quantity, 52(5), 2031–2036. https://doi.org/10.1007/s11135-017-0593-5
“They also fail to acknowledge the importance of including the group-mean variable in the model to achieve unbiased within effects.” (Group mean centering and within, between and contextual effects)
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. https://doi.org/10.18637/jss.v067.i01
“the default behavior is to refit the models with ML if fitted with REML = TRUE, which is necessary in order to get sensible answers when comparing models that differ in their fixed effects” (§5.2, p. 34)

“Thus approximate methods for computing the approximate degrees of freedom for t distributions, or the denominator degrees of freedom for F statistics (Satterthwaite 1946; Kenward and Roger 1997), are at best ad hoc solutions.” (§5.2, p. 34)
Kuznetsova, A., Brockhoff, P. B., & Christensen, R. H. B. (2017). lmerTest package: Tests in linear mixed effects models. Journal of Statistical Software, 82(13), 1–26. https://doi.org/10.18637/jss.v082.i13
“We have implemented the Satterthwaite's method for approximating degrees of freedom for the t and F tests.” (Abstract)
Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin, 2(6), 110–114. https://doi.org/10.2307/3002019
Nakagawa, S., & Schielzeth, H. (2013). A general and simple method for obtaining R² from generalized linear mixed-effects models. Methods in Ecology and Evolution, 4(2), 133–142. https://doi.org/10.1111/j.2041-210x.2012.00261.x
“most definitions of R² for mixed-effects have theoretical problems (e.g. decreased or negative R² values in larger models)” (Summary)
Maas, C. J. M., & Hox, J. J. (2005). Sufficient sample sizes for multilevel modeling. Methodology, 1(3), 86–92. https://doi.org/10.1027/1614-2241.1.3.86
“only a small sample size at level two (meaning a sample of 50 or less) leads to biased estimates of the second-level standard errors.” (Abstract)
Self, S. G., & Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605–610. https://doi.org/10.1080/01621459.1987.10478472
Hoenig, J. M., & Heisey, D. M. (2001). The abuse of power: The pervasive fallacy of power calculations for data analysis. The American Statistician, 55(1), 19–24. https://doi.org/10.1198/000313001300339897
Kumle, L., Võ, M. L.-H., & Draschkow, D. (2021). Estimating power in (generalized) linear mixed models: An open introduction and tutorial in R. Behavior Research Methods, 53(6), 2528–2543. https://doi.org/10.3758/s13428-021-01546-0
“A flexible and very intuitive alternative to analytic power solutions are simulation-based power analyses.” (Abstract)
Snijders, T. A. B. (2005). Power and sample size in multilevel linear models. In B. S. Everitt & D. C. Howell (Eds.), Encyclopedia of statistics in behavioral science (Vol. 3, pp. 1570–1573). Wiley. https://doi.org/10.1002/0470013192.bsa492
“For estimating a population mean (case 1) or the effect of a level-two variable (case 4), on the other hand, the multilevel design always is less efficient, and more seriously so as the cluster size and the intraclass correlation are larger.” (Case 4)
Hedges, L. V., & Rhoads, C. (2010). Statistical power analysis in education research (NCSER 2010-3006). National Center for Special Education Research, Institute of Education Sciences, U.S. Department of Education. https://ies.ed.gov/ncser/pubs/20103006/pdf/20103006.pdf
“beyond a point (which occurs when n is rather modest), obtaining a larger sample size by increasing n has little effect on power.” (p. 19)
Schielzeth, H., Dingemanse, N. J., Nakagawa, S., Westneat, D. F., Allegue, H., Teplitsky, C., Réale, D., Dochtermann, N. A., Garamszegi, L. Z., & Araya-Ajoy, Y. G. (2020). Robustness of linear mixed-effects models to violations of distributional assumptions. Methods in Ecology and Evolution, 11(9), 1141–1152. https://doi.org/10.1111/2041-210X.13434
“Model estimates were usually robust to violations of assumptions, with the exception of slight upward biases in estimates of random effect variance if the generating distribution was bimodal but was modelled by Gaussian error distributions.” (Abstract)
Raykov, T., & DiStefano, C. (2022). Design effect in multilevel settings: A commentary on a latent variable modeling procedure for its evaluation. Educational and Psychological Measurement, 82(5), 1020–1030. https://doi.org/10.1177/00131644211019447
“represents the ratio of the squared SE (i.e., the sampling variance) of the grand mean estimate when evaluated within an explicit two-level setting, to its squared SE (i.e., sampling variance) if this parameter is evaluated disregarding the clustering effects” (Design Effect in Two-Level Studies)
Muff, S., Nilsen, E. B., O'Hara, R. B., & Nater, C. R. (2022). Rewriting results sections in the language of evidence. Trends in Ecology & Evolution, 37(3), 203–210. https://doi.org/10.1016/j.tree.2021.10.009
“we suggest language of evidence that allows for a more nuanced approach to communicate scientific findings as a simple and intuitive alternative to statistical significance testing” (Abstract)
Williams, M. N., Gómez Grajales, C. A., & Kurkiewicz, D. (2013). Assumptions of multiple regression: Correcting two misconceptions. Practical Assessment, Research & Evaluation, 18(11), 1–14. https://doi.org/10.7275/55hn-wk47
“ordinary least squares estimates will be unbiased and consistent as long as the errors are independent, but will not be efficient” (p. 9)
Hayes, A. F., & Cai, L. (2007). Using heteroskedasticity-consistent standard error estimators in OLS regression: An introduction and software implementation. Behavior Research Methods, 39(4), 709–722. https://doi.org/10.3758/BF03192961