Multilevel SEM: Interpretation and APA Reporting

The Basics

When to Use It

Use it when people are nested in groups (employees in teams, students in classrooms, clients in therapy groups) and your question has a structure: paths between several variables, possibly at both levels, possibly with latent variables. Multilevel SEM splits every person-level variable into a within level part and a between level part (the group's latent group mean), fits a path model to each level, and estimates them together, so paths at the two levels can be compared and tested (Muthén, 1994). For example:

Are employees who work longer hours than their teammates more exhausted, and are teams with a more supportive climate less exhausted?
Do students who feel more competent than their classmates enjoy math more, and do classrooms with more supportive teachers enjoy it more on average?
Does a therapy group's cohesion predict members' average symptom change, over and above each member's own attendance?

Vocabulary You Need

Definition Within Level

In data where people are nested in groups (employees in teams, students in classrooms), the level that looks inside each group: each person's score minus their own group's true mean. Relationships at this level compare people with their own teammates, so anything shared by the whole group, such as its manager or its workload norms, is held constant automatically. Also called level 1. A within-level relationship can be quite different from the between-level one for the same two variables. Glossary →

Definition Between Level

In data where people are nested in groups, the level that compares whole groups: each group's true mean on each variable, plus any variables measured once per group, such as a team's climate rating. Relationships at this level are about groups: that teams with a better climate are less exhausted describes their members on average, not each member, and a group average of a person-level variable can relate to the outcome quite differently than individuals' own scores do. Also called level 2. With only as many cases as there are groups, it's the level where sample size usually runs short. Glossary →

Definition Latent Group Mean

The mean a group would have if every possible member were measured, treated as an unobserved (latent) variable. The average of the members actually sampled estimates it with error: with 8 employees per team, part of the spread in team averages is sampling noise, about the within-group variance divided by 8. Multilevel SEM subtracts that noise, so relationships between groups are estimated from the latent means. Using observed averages instead, as a multilevel regression with group-mean predictors does, blurs between-group relationships toward the within-group ones, and the fewer members per group and the lower the ICC, the more (Lüdtke et al., 2008). Latent means fit when members are a sample of the group, or when their answers reflect a shared group property. When every member is measured and the group score is just their average, the observed mean can be the better choice (Lüdtke et al., 2008). Glossary →

Definition Contextual Effect (

b_{B} - b_{W}

)

When a person-level predictor, such as hours worked, relates to an outcome at both levels, the between-group slope (how group means go together) minus the within-group slope (how people's scores relative to their own group go together). If it's 0, the group average adds nothing beyond each person's own score: two people working the same hours are expected to be equally exhausted whatever their team's average hours (other team characteristics, such as climate, can still differ). A positive contextual effect means being in a group that scores high goes with a higher outcome over and above one's own score, as when employees in long-hours teams report more strain than their own hours predict. Estimated from latent group means, it avoids the shrinkage toward zero that observed group averages can bring when members are a sample of the group (Lüdtke et al., 2008). It rests on the number of groups, so it's usually estimated imprecisely. Glossary →

Assumptions

1. Enough Independent Clusters You have enough separate people or groups, not just enough observations, for the estimates you need.

What It Means: Anything estimated from differences between clusters rests on the number of clusters, not the number of observations inside them. In a GEE or any analysis with robust (sandwich) standard errors, the standard errors are learned from the clusters themselves; with a few dozen clusters or fewer they tend to be too small, so p-values and confidence intervals are too optimistic. In a multilevel (mixed) model, the standard errors come from the model instead, but the random-effect variances and any predictor measured at the cluster level are still estimated from differences between clusters, so with few clusters they are imprecise, the variances can be biased (Moineddin et al., 2007, found this for yes/no outcomes) and their standard errors come out too small.
Applies To: Design
How to Check: Count the independent clusters (people when each person gives repeated answers, groups when people are nested in groups), not the rows. Suggested minimums vary by method, from about 30 for sandwich standard errors (Halekoh et al., 2006) to about 50 for a multilevel model's cluster-level standard errors (Maas & Hox, 2005). With robust standard errors and only a few dozen clusters, use a small-sample correction such as Mancl and DeRouen's bias-corrected estimator; in their simulations with 10 to 40 clusters, the jackknife, like the uncorrected sandwich, could give false-positive rates two to three times the intended level. With a mixed model, add intervals that don't lean on large-sample normality, such as profile-likelihood or parametric-bootstrap intervals, and for a yes/no or count outcome check whether the results hold with more quadrature points.
How Much It Matters: A lot, at the between level. The climate path and the contextual effect rest on 60 teams, not 480 employees. In our simulation, 60 teams gave the climate path 95% power but a contextual effect of 0.02 points per hour only 10%. With 50 groups or fewer, second-level standard errors, those of the between-level variances in particular, come out too small (Maas & Hox, 2005), and in multilevel SEM how many groups the between part needs depends on the estimator (Hox et al., 2010). The within level, with hundreds of employees, is far less affected.

If It's Violated

Design Change Add groups rather than members. Between-level estimates depend mainly on how many groups there are; in our simulation, going from 60 to 200 teams of 8 only raised the power to detect the contextual effect from 10% to 23%.
Alternative Test Multilevel Regression (HLM). With few groups, a simpler multilevel regression (observed team means, REML, small-sample df corrections) is more stable, at the cost of the latent means.
Proceed with Caveat Estimate only the between-level paths the question needs, report their intervals, and treat their tests as approximate. Watch for a between-level variance near or below zero, or a model that won't converge: signs the groups can't support it. Read Wald intervals for variances as rough, since they're symmetric and a variance can't go below zero.

2. 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: Some. The model handles employees who share a team; it assumes the teams are independent of one another. Teams that share a department or a manager are related, so the between-level standard errors come out too small.

If It's Violated

Design Change If teams sit in departments, add departments as a third level, or sample one team per department.
Correction Correct the between-level standard errors for the higher-level clustering if a third level can't be modeled.

3. Correctly Specified Model Your model includes the paths and correlations that really matter.

What It Means: The model's pattern is right in the ways that count: every path, loading or covariance that matters is included, and every one the model fixes at zero really is near zero. In a factor model, the items load on the factors you assigned them to and any residuals that share something (similar wording, a shared method) are allowed to correlate; in a path model, the arrows left out of the diagram really are absent. Every estimate, and every fit index, is conditional on that pattern. An omitted path, cross-loading or residual correlation doesn't just cost fit; the model makes up for it by bending the parameters it does have.
Applies To: Covariance structure
How to Check: Specify the model from theory (and, for a factor model, item content) before seeing the data. Then look beyond the overall fit: the residual correlations for pairs the model misses, the modification indices for where it misses most, and paths, loadings or factor correlations that look implausible. Compare the model with plausible alternatives specified in advance. Fit statistics can show a model is wrong; they can't show it's right.
How Much It Matters: A lot, and the data can't warn you here. Both levels are saturated: the model estimates as many numbers as the data supply, so it reproduces the data exactly, whatever is wrong with it, and can't be tested. The risk is what isn't in it: a team-level cause of both climate and exhaustion (a demanding client, a poor manager) would be credited to climate, and a curved hours–exhaustion relationship would be summarized by a straight slope. The levels must be right too: a person-level variable left out of the team level loses its between-team part, and its effect there goes unmodeled. The model also assumes the within-team hours slope is the same in every team. If it varies, the single b_W is an average, and its standard error can be too small. Check with a multilevel regression that adds a random slope for hours; lavaan's two-level SEM can't fit random slopes, but Mplus can.

If It's Violated

Design Change Measure plausible team-level common causes of climate and exhaustion (workload targets, manager, department) and add them as between-level predictors.
Design Change Specify competing models in advance, such as with and without a contextual effect, and compare them, rather than adjusting one model after seeing the data.
Proceed with Caveat Report that the model is saturated, and don't cite its fit indexes as support.

4. Multivariate Normality The variables are roughly bell-shaped together (within each group, when there are groups), or the sample is big enough that it matters less.

What It Means: The variables jointly follow a multivariate normal distribution: each is roughly normal, and every weighted combination of them is too, so the scatterplot of any two is an oval cloud. When the design has groups, as in MANOVA, it applies within each group; in a multivariate regression, to the outcomes' residuals; in factor analysis and SEM, it applies to the observed variables that maximum likelihood is fitted to (for outcomes, after allowing for any predictors the model takes as given). P-values and standard errors rely on it; estimates are usually less affected, and robust standard errors or a scaled test statistic can take its place.
Applies To: Joint distribution
How to Check: Look at a histogram of each variable (or each outcome's residuals, in a model with predictors) and a scatterplot of each pair, within each group if there are groups, for strong skew, curves or stray points. Multivariate outliers show up as large squared Mahalanobis distances from the means (each group's own, if there are groups). Don't rely on a normality test as a gatekeeper: with small samples it misses real problems, and with large ones it flags trivial ones.
How Much It Matters: Some. Maximum likelihood assumes the within parts and the latent team means are multivariate normal. Skewed exhaustion scores or a few extreme teams make standard errors too small and likelihood-ratio tests too liberal; the estimates themselves hold up better.

If It's Violated

Robust Estimation Use MLR (estimator = "MLR" in lavaan), which corrects the standard errors and test statistics for non-normality. With few groups, Hox et al. (2010) don't recommend the robust methods.
Transformation Transform a strongly skewed variable, such as hours with a long tail, and interpret on the new scale.

Limitations

This Is Not Your Test If…

Your people aren't nested in groups Path Analysis →
You have one outcome and just want slopes adjusted for clustering Multilevel Regression (HLM) →
Your outcome is yes/no Multilevel Logistic Regression →
You expect the within-group slope to differ from group to group Multilevel Regression (HLM) →
Your repeated measures are nested in people and you want change over time Multilevel Growth-Curve Modeling →
You only need how much of the variation lies between groups Intraclass Correlation Coefficient (ICC) →

What It Doesn't Tell You

What happens to any one employee. The climate path says that employees in better-climate teams are, on average, less exhausted, among teams working the same average hours. It doesn't say every employee in such a team is, or that moving someone there would help (nothing was manipulated). For a person-level variable such as hours, reading the team-level slope as the individual one is the ecological fallacy (Robinson, 1950).
Why. Nothing here was manipulated. Longer hours may exhaust people, exhausted people may work less efficiently and stay later, or a demanding role may drive both. Climate wasn't randomized either, and a team-level cause of both climate and exhaustion, such as a poor manager, would look like a climate effect.
That the model is right. Both levels are saturated, so the model reproduces the data perfectly whatever its flaws. Nothing on this page tests it as a whole.
That there's no contextual effect. The test found none (p = .718), but its interval runs from −0.050 to 0.072 points per hour, and 60 teams had only about a 10% chance of detecting a contextual effect of 0.02.
Who else. These are simulated teams, with 8 employees sampled from each. Real teams vary in size and may be nested in departments, which would add a level.

Common Mistakes

Don'tRead a team-level slope of a person-level variable, such as team-average hours, as the slope for individual employees.

InsteadUse the within-team path for how an employee's own hours relate to their exhaustion, and say which level each result describes.
Don'tCount employees when judging whether the team-level part is big enough.

InsteadCount teams. Between-level estimates and tests rest on the number of groups.
Don'tReport CFI, RMSEA or SRMR as evidence of good fit for a saturated model or level.

InsteadSay the model has 0 degrees of freedom and fits perfectly by construction; test specific constraints by comparing nested models.
Don'tChoose latent or observed team averages without asking how the team score arises.

InsteadUse latent means when members are a sample of a larger group, or when each member's answer reflects a shared team property. When every member is measured and the team score is just their average, the observed mean can be the better choice (Lüdtke et al., 2008).
Don'tConclude there's no contextual effect because its test isn't significant.

InsteadReport its confidence interval and the number of groups. With 60 teams, a sizable contextual effect can go undetected.
Don'tRead a bigger χ² as a stronger effect.

InsteadReport an effect size suited to the design, such as Cramér's V, an odds ratio or a model's coefficients. χ² grows with sample size even when the effect doesn't.
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_W = 0.038, 95% CI [0.022, 0.055], b* = .22, p < .001; b_B,climate = −0.41, 95% CI [−0.59, −0.22], b* = −.61, z = −4.27, p < .001; contextual effect of hours = 0.011, 95% CI [−0.050, 0.072], χ²(1) = 0.13, p = .718

Table 1. Two-Level Path Model: Within- and Between-Team Estimates

Example

	Estimate	$SE$	95% CI	$b^*$	$p$
Within teams (level 1)
Hours → exhaustion	0.038	0.008	[0.022, 0.055]	.22	< .001
Residual variance of exhaustion	0.751	0.052	[0.649, 0.852]
Between teams (level 2)
Climate → exhaustion	−0.405	0.095	[−0.591, −0.219]	−.61	< .001
Team-average hours → exhaustion	0.050	0.029	[−0.006, 0.106]	.30	.082
Contextual effect of hours (team − within)	0.011	0.031	[−0.050, 0.072]
Hours ↔ climate (covariance)	−0.850	0.304	[−1.445, −0.254]	−.47	.005
Variance of team-mean hours	7.247	1.927	[3.470, 11.024]
Variance of climate	0.456	0.083	[0.293, 0.620]
Residual variance of team-mean exhaustion	0.075	0.032	[0.013, 0.137]

Note. N = 480 employees in 60 teams. Maximum likelihood; standard errors from the observed information (the curvature of the likelihood at its peak). b* = standardized at the path's own level (for the covariance, a correlation). Hours and exhaustion were split into within-team parts and latent team means; climate was measured once per team. The contextual effect is the team-level hours path minus the within-team one; its standard error uses the two paths' covariance. R² for exhaustion: .049 within teams, .627 between teams.

Table 2. How Each Variable's Variance Splits between the Levels

Example

	Within-Team Variance	Between-Team Variance	ICC	Design Effect
Weekly hours	26.468	7.247	.21	2.50
Emotional exhaustion	0.790	0.202	.20	2.43
Team climate		0.456

Note. Estimates from the unrestricted two-level model, which the full model reproduces exactly. Between-team variances are for the latent team means (the variance of the observed team averages, dividing by the number of teams as maximum likelihood does, minus within-team variance ÷ 8). ICC = between ÷ (between + within); design effect = 1 + 7 × ICC. Climate was measured once per team, so it has only a between-team variance.

Table 3. Testing the Contextual Effect: Two Nested Models

Example

	$\ln L$	Parameters	$df$	$\chi^2$ vs Full	$p$	Climate Path	$SE$
Full: team-level hours path free	−2188.814	11	0			−0.405	0.095
No contextual effect: team-level hours path = within-team path	−2188.879	10	1	0.13	.718	−0.424	0.080

Note. Both models are fitted by maximum likelihood to the same data. The full model frees the team-level hours path; the second constrains it to equal the within-team path (no contextual effect). The likelihood-ratio χ² is twice the difference in log-likelihood, on 1 degree of freedom. It compares the two models; it isn't a test of either model's fit. The full model has 0 degrees of freedom: it reproduces the unrestricted two-level model exactly, so there is nothing left to test.

Full Text (APA Style)

Example

A two-level path model (maximum likelihood, lavaan 0.7-2; Rosseel, 2012) split weekly hours and emotional exhaustion into within-team parts and latent team means; about a fifth of the variance in each lay between teams (ICC = .21 for hours, .20 for exhaustion; Table 2). Within teams, employees who worked more hours than their teammates were more exhausted, b_W = 0.038 points per hour, 95% CI [0.022, 0.055], b* = .22, z = 4.68, p < .001. Between teams, holding team-average hours constant, teams with a more supportive climate had lower mean exhaustion, b_B = −0.41 points per climate point, 95% CI [−0.59, −0.22], b* = −.61, z = −4.27, p < .001 (Table 1). The team-level hours slope (0.050) was close to the within-team slope: the contextual effect was 0.011 points per hour, 95% CI [−0.050, 0.072], and constraining it to zero did not significantly worsen fit, χ²(1) = 0.13, p = .718 (Table 3). With 60 teams, that interval is too wide to rule out a sizable contextual effect. Both levels of the model are saturated, so it has no test of overall fit.

Say how many groups and how many people per group, and give each variable's ICC, before any paths.
Label every path with its level (within or between), and say which variables exist at which level.
Name the software, the estimator and how the between parts were formed (latent group means here).
Report a contextual effect as an estimate with its interval, not only as a significance test, and say how many groups it rests on.
If a level is saturated, say so, and don't quote fit indexes for it.

Reading Each Piece

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

ICC = .20	About a fifth of the differences in exhaustion are differences between teams, and the same is true of hours (.21). That's plenty to matter: analyzed as 480 independent people, the sampling variance of a quantity such as the overall mean exhaustion would be understated by a factor of about 2.4, and its standard error by a factor of about 1.6.
b_W = 0.038	Within a team, an employee who works 10 hours a week more than their teammates is about 0.38 points more exhausted on the 0-to-6 scale.
b_B,climate = −0.41	Among teams working the same average hours, a team rated one point higher on supportive climate (say 4 instead of 3) is about 0.4 points less exhausted on average. That's an average over each team's members; it doesn't say every member is less exhausted.
95% CI [−0.59, −0.22]	The team-level climate path is plausibly anywhere from −0.59 to −0.22 points per climate point. Sixty teams pin it down only roughly, but the whole range is well below zero.
95% CI [0.022, 0.055]	The within-team path is plausibly anywhere from 0.022 to 0.055 points per hour: 0.2 to 0.55 points for 10 extra hours. Clearly above zero.
Contextual effect = 0.011	Beyond what each employee's own hours predict, working on a team that averages an extra hour goes with about 0.01 more points of exhaustion: close to nothing.
95% CI [−0.050, 0.072]	The contextual effect is plausibly anywhere from −0.050 to 0.072 points per hour. That range includes zero, but also a team-level slope nearly three times the within-team one. Sixty teams can't tell these apart.
b*_W = .22	One within-team SD of hours (5.1 hours) goes with about a fifth of a within-team SD of exhaustion. Hours predict 5% of the within-team variance in exhaustion (R² = .049).
b*_B,climate = −.61	One SD of climate (0.68 points) goes with 0.61 SD of the latent team means. The SD of true team exhaustion is small (0.45 points), so this large-sounding standardized path is about a quarter of a point. Climate and team-average hours together predict 63% of the between-team variance (R² = .627).
z = −4.27	The climate path is 4.27 standard errors below zero.
χ²(1) = 0.13	Forcing the team-level and within-team hours slopes to be equal barely changes the fit: the data are consistent with no contextual effect.
p < .001	If climate were unrelated to team exhaustion among teams with the same average hours, a path at least this far from zero, in either direction, would almost never turn up. The within-team hours path, at z = 4.68, also gives p < .001. Benchmarks ↑
p = .718	If the team-level and within-team hours slopes were equal in the population, a gap in fit at least this big would turn up about 72% of the time. That's no evidence of one, but not evidence against a moderate one either: in our simulation, 60 teams detected a contextual effect of 0.02 points per hour only 10% of the time. Benchmarks ↑

Calculations

The Running Example

Example

The team-exhaustion study: 480 employees, 8 from each of 60 teams, with weekly hours and emotional exhaustion measured on each employee and a supportive-climate rating on each team.

	Null Hypothesis	Alternative
Within-team hours path	$b_W = 0$	$b_W \neq 0$
Team climate path	$b_{B,\text{climate}} = 0$	$b_{B,\text{climate}} \neq 0$
Contextual effect of hours	$b_{B,\text{hours}} - b_W = 0$	$b_{B,\text{hours}} - b_W \neq 0$

$b_W$ is how much more exhausted, on average, an employee is per weekly hour they work above their own team's average, in the population of employees and teams like these. $b_{B,\text{climate}}$ is how much lower, on average, a team's true mean exhaustion is per point higher on climate, among teams with the same average hours. $b_{B,\text{hours}}$ is how much higher a team's true mean exhaustion is per hour of team-average hours, among teams with the same climate; the contextual effect is how far that team-level slope is from the within-team one (Lüdtke et al., 2008).

What Your Data Should Look Like

One row per person, with a column saying which group they belong to. Person-level measures (hours, exhaustion) vary within groups; group-level measures (a team's climate rating) repeat the same value on every row of the group. The software splits the person-level measures into within and between parts itself, so don't center them or add group averages. Count the groups as well as the people: the between-level part of the model rests on the number of groups, and a few dozen is a practical minimum.

Example

employee	team	climate	hours	exhaustion
1	1	3.2	45	4.4
2	1	3.2	37	3.6
3	1	3.2	45	4.0
4	1	3.2	33	3.6
5	1	3.2	48	4.0
6	1	3.2	39	1.8
… 474 more rows

	M	SD	1	2
1. Emotional exhaustion	2.55	1.00
2. Weekly hours	41.36	5.81	.30
3. Team climate	3.53	0.68	−.34	−.22

N = 480.

Correlations are Pearson r.

Simulated data.

Download all 480 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 team-exhaustion study's numbers. Collapse all to read just the headings: they tell the whole story on their own.

1. Measure How Much Employees Differ from Their Own Teammates

Start inside the teams. Take each employee's distance from their own team's average and pool those differences across all 60 teams. This is the within-team variance: differences between teammates, with everything the team shares taken out.

$\sigma^2_W = \dfrac{\sum_j \sum_i (y_{ij} - \bar{y}_j)^2}{N - G}$ = $\sigma^2_W = \dfrac{\sum_j \sum_i (y_{ij} - \bar{y}_j)^2}{480 - 60} = 0.790$

= 0.790

$y_{ij}$ is employee $i$ 's exhaustion in team $j$ and $\bar{y}_j$ that team's average. The squared distances, added over all 480 employees and divided by $N - G$ = 480 − 60 = 420, give 0.790 squared points: teammates typically differ by about √0.790 = 0.89 of a point.
2. Take the Sampling Noise Out of the Team Averages

The 60 team averages also spread out, but part of that spread is noise: each average rests on only 8 people. Subtracting the noise the within-team variance predicts leaves the variance of the teams' true means, the latent group means the model works with. That fits here because the 8 employees are a sample of each team's members. If you measure every member and want the team's actual average, the observed average can be the better choice (Lüdtke et al., 2008).

$\sigma^2_B = s^2_{\bar{y}} - \dfrac{\sigma^2_W}{n}$ = $0.301 - \dfrac{0.790}{8} = 0.301 - 0.099$

= 0.202

$s^2_{\bar{y}}$ = 0.301

$s^2_{\bar{y}}$ is how much the 60 observed team averages vary, 0.301, dividing by the 60 teams as maximum likelihood (ML, which picks the values that make the observed data most probable) does. (Dividing by 59, as a calculator's variance does, gives slightly larger values; that's what REML, restricted maximum likelihood, reports.) With $n$ = 8 per team, sampling noise alone contributes $\sigma^2_W / n$ = 0.790 ÷ 8 = 0.099. That leaves 0.202 for real differences between teams, a third less than the raw averages suggest.
3. Find the Share of Exhaustion That Lies between Teams

Put the two pieces together. The intraclass correlation says what share of all the variation in exhaustion is about which team someone is in. It also says how much the clustering shrinks the sample: the design effect.

$\text{ICC} = \dfrac{\sigma^2_B}{\sigma^2_B + \sigma^2_W}$ = $\dfrac{0.202}{0.202 + 0.790} = \dfrac{0.202}{0.992}$

= .20

$\text{ICC}_{\text{hours}}$ = .21
$\text{deff}$ = 2.43

$\sigma^2_B$ is the 0.202 between teams and $\sigma^2_W$ the 0.790 within them. 0.202 ÷ 0.992 = .20. With 8 per team, the design effect is 1 + 7 × .20 = 2.43, so for estimating the overall mean exhaustion our 480 employees carry about as much information as 198 independent ones. How much clustering costs depends on the estimate: within-team paths use every employee, while team-level paths rest on the 60 teams.
4. Estimate How Extra Hours Go with Exhaustion within a Team

Now the within-team path. It compares each employee with their own teammates: someone who works more hours than their team's average, are they more exhausted than their team's average? Anything the whole team shares, including its climate and workload norms, is held constant by comparing teammates. Differences between the employees themselves, such as their roles, are not.

$y_{ij} = y^{W}_{ij} + y^{B}_{j}, \quad y^{W}_{ij} = b_W\, x^{W}_{ij} + e_{ij}, \quad y^{B}_{j} = \alpha + b_{B,x}\, x^{B}_{j} + b_{B,z}\, z_j + u_j$ = $b_W = \dfrac{\text{cov}_W(x, y)}{\sigma^2_{x,W}} = \dfrac{1.017}{26.468} = 0.0384, \quad 0.0384 \pm 1.96 \times 0.0082$

= 0.038 points per hour

$95\%\ \text{CI}$ = [0.022, 0.055]
$b^*_W$ = .22
$R^2_W$ = .049
$\text{cov}_W(x, y)$ = 1.017

In the model, $x$ is hours, $y$ exhaustion and $z$ climate; $e$ and $u$ are what each level's paths leave over. $b_W$ links the within parts, and with one predictor it's an ordinary regression slope: the within-team covariance of hours and exhaustion, 1.017, divided by the within-team variance of hours, 26.468 (Table 2). So each weekly hour above the team's average goes with 0.0384 more points of exhaustion. Its standard error is 0.0082, so 1.96 × 0.0082 = 0.016 and the interval runs from 0.022 to 0.055. Ten extra hours go with about 0.38 points, over a third of a within-team SD.
5. Estimate How Team Climate Goes with Team-Average Exhaustion

The between-team paths work on the latent team means. Climate was rated once per team, so it exists only at this level. The model also lets team-average hours predict team exhaustion, so the climate path compares teams with the same average hours: do teams with a better climate have lower true mean exhaustion? The arithmetic is an ordinary two-predictor regression, run on the latent team means' variances and covariances instead of raw scores. That's the one new idea.

$y_{ij} = y^{W}_{ij} + y^{B}_{j}, \quad y^{W}_{ij} = b_W\, x^{W}_{ij} + e_{ij}, \quad y^{B}_{j} = \alpha + b_{B,x}\, x^{B}_{j} + b_{B,z}\, z_j + u_j$ = $b_{B,\text{climate}} = \dfrac{\text{cov}_B(z,y)\,\sigma^2_{x,B} - \text{cov}_B(x,y)\,\text{cov}_B(x,z)}{\sigma^2_{x,B}\,\sigma^2_{z,B} - \text{cov}_B(x,z)^2} = \dfrac{(-0.227)(7.247) - (0.704)(-0.850)}{(7.247)(0.456) - (-0.850)^2} = \dfrac{-1.047}{2.582} = -0.405, \quad -0.405 \pm 1.96 \times 0.095$

= −0.41 points per climate point

$95\%\ \text{CI}$ = [−0.59, −0.22]
$b^*_B$ = −.61
$R^2_B$ = .627
$\text{cov}_B(x, y)$ = 0.704
$\text{cov}_B(z, y)$ = −0.227

Between teams, the latent team means of hours ( $x$ ) and exhaustion ( $y$ ) covary by 0.704, climate ( $z$ ) and exhaustion by −0.227, and hours and climate by −0.850; the variances of team-mean hours and climate are 7.247 and 0.456 (Table 1). Solving the regression gives $b_{B,\text{climate}}$ , the climate path among teams with the same average hours: each point of climate goes with 0.41 fewer points of true mean exhaustion. The same arithmetic gives the team-level hours slope, 0.0496, used in step 7. Its standard error, 0.095, rests on 60 teams, so the interval is wide: −0.405 ± 0.186, from −0.59 to −0.22.
6. Ask How Surprising the Climate Path Would Be If Climate Were Unrelated to Exhaustion

Divide the path by its standard error to see how many standard errors it is from zero, and compare that with the normal curve. Maximum likelihood standard errors are large-sample approximations, and here the sample that counts is 60 teams.

$z = \dfrac{b}{SE_b}$ = $z = \dfrac{-0.4052}{0.0950} = -4.27, \quad p = 2 \times P(Z \ge 4.27)$

< .001

$z$ = −4.27
$SE_b$ = 0.095

$b$ is the −0.4052 climate path and $SE_b$ its 0.0950 standard error, so the path is 4.27 standard errors below zero. Both tails of the normal curve beyond 4.27 hold about .00002 of it, 1 in 50,000.
7. Compare the Team-Level Hours Slope with the Individual One

Hours relate to exhaustion at both levels. The team-level slope asks whether teams that work longer on average are more exhausted than you'd expect from each member's own hours; the difference between it and the within-team slope is the contextual effect (Lüdtke et al., 2008). Reading either level's slope as the other's is the ecological fallacy.

$b_{B,x} - b_W$ = $0.0496 - 0.0384 = 0.011$

= 0.011 points per hour

$b_{B,\text{hours}}$ = 0.050
$b_{\text{observed means}}$ = 0.045
$95\%\ \text{CI}$ = [−0.050, 0.072]

$b_{B,\text{hours}}$ is the team-level slope, 0.0496 points per hour of team-average hours (holding climate constant), and $b_W$ the within-team 0.0384. The gap is 0.011 points per hour. Computed from the noisy observed team averages instead of latent means, the team-level slope is 0.045, pulled toward the within-team value.
8. Count What the Model Has Left to Test

A model can only be tested against the data if it reproduces fewer variances, covariances and means than the data supply. Ours estimates every one of them, at both levels, so it fits perfectly by construction and has no test of fit.

$df = (\text{within moments} + \text{between moments} + \text{means}) - q$ = $11 - 11 = 0$

= 0

$\text{moments}$ = 11
$q$ = 11

The data supply 11 pieces of information: 2 within-team variances and covariances (lavaan treats the within-team spread of hours as given, as a regression treats its predictor, so it's neither estimated nor counted), 6 between-team ones and 3 means. The model estimates 11 parameters. 11 − 11 = 0, so χ² is 0 and fit indexes such as RMSEA, CFI and SRMR say nothing.
9. Test Whether the Team-Level Hours Slope Differs from the Individual One

A question the model can test is whether freeing the team-level hours slope was needed. Fit the model again with that slope constrained to equal the within-team one, then compare the two log-likelihoods. A log-likelihood measures how well a model's estimates reproduce the data (higher, closer to 0, is better). This likelihood-ratio test compares two models; it isn't evidence that either one fits.

$\chi^2 = 2\,(\ln L_{\text{full}} - \ln L_{\text{restricted}})$ = $2 \times \big(-2188.814 - (-2188.879)\big) = 2 \times 0.065$

= 0.13

$\ln L_{\text{full}}$ = −2188.814
$\ln L_{\text{restricted}}$ = −2188.879
$P(\chi^2_1 \ge 0.13)$ = 0.718
$b_{B,\text{climate}}^{\text{restricted}}$ = −0.42

$\ln L_{\text{full}}$ is the full model's log-likelihood, −2188.814, and $\ln L_{\text{restricted}}$ the constrained model's, −2188.879. Twice the 0.065 gap is 0.13, on 1 degree of freedom (one constraint).

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)

SPSS Statistics has no multilevel SEM, and neither does its SEM add-on, Amos. What SPSS can do is the first half of the analysis, splitting exhaustion's variance into within-team and between-team parts, with a random-intercept model:

Analyze → Mixed Models → Linear.
Move team into Subjects, and click Continue.
Move exhaustion into Dependent Variable.
Random: tick Include intercept, move team into Combinations, and click Continue.
Estimation: choose Maximum Likelihood (ML). Click OK.

For the two-level path model itself, use R's lavaan (below) or Mplus.

Or paste this syntax:

MIXED exhaustion
  /FIXED=| SSTYPE(3)
  /METHOD=ML
  /PRINT=SOLUTION TESTCOV
  /RANDOM=INTERCEPT | SUBJECT(team) COVTYPE(VC).

Reading the Output

Read the Estimates of Covariance Parameters table. With our data, the Residual row (the within-team variance) should be 0.790 and the Intercept [subject = team] Variance row (the between-team variance) 0.202, the values in Table 2. Divide the second by their sum for the ICC: 0.202 ÷ 0.992 = .20. Repeat with hours as the dependent variable for its ICC.

With equal team sizes, ML's between-team variance is the latent one from Napkin Math step 2, not the variance of the observed team averages. SPSS's default, REML, gives a slightly larger between-team variance (it corresponds to dividing the team averages' variance by 59 rather than 60). These values were computed from the data, not produced by running SPSS.

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(lavaan)
teams <- read.csv("data.csv")

model <- '
  level: 1
    exhaustion ~ bw*hours
  level: 2
    exhaustion ~ bb*hours + climate
    hours ~~ climate
  contextual := bb - bw
'
fit <- sem(model, data = teams, cluster = "team")
summary(fit, standardized = TRUE, rsquare = TRUE)
lavInspect(fit, "icc")

# No contextual effect: the team-level hours slope equals the within-team one
model0 <- '
  level: 1
    exhaustion ~ b*hours
  level: 2
    exhaustion ~ b*hours + climate
    hours ~~ climate
'
fit0 <- sem(model0, data = teams, cluster = "team")
lavTestLRT(fit0, fit)

lavaan 0.7-2 ended normally after 41 iterations

  Estimator                                         ML
  Number of model parameters                        11
  Number of observations                           480
  Number of clusters [team]                         60

Model Test User Model:
  Test statistic                                 0.000
  Degrees of freedom                                 0

Parameter Estimates:
  Standard errors                             Standard
  Information                                 Observed

Level 1 [within]:
Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  exhaustion ~
    hours     (bw)    0.038    0.008    4.676    0.000    0.038    0.222
Variances:
   .exhaustion        0.751    0.052   14.491    0.000    0.751    0.951
R-Square:
    exhaustion        0.049

Level 2 [team]:
Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  exhaustion ~
    hours     (bb)    0.050    0.029    1.738    0.082    0.050    0.297
    climate          -0.405    0.095   -4.267    0.000   -0.405   -0.609
Covariances:
  hours ~~
    climate          -0.850    0.304   -2.797    0.005   -0.850   -0.467
Intercepts:
   .exhaustion        1.931    1.397    1.383    0.167    1.931    4.293
    hours            41.362    0.419   98.617    0.000   41.362   15.365
    climate           3.530    0.087   40.473    0.000    3.530    5.225
Variances:
   .exhaustion        0.075    0.032    2.382    0.017    0.075    0.373
    hours             7.247    1.927    3.760    0.000    7.247    1.000
    climate           0.456    0.083    5.477    0.000    0.456    1.000
R-Square:
    exhaustion        0.627

Defined Parameters:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    contextual        0.011    0.031    0.361    0.718    0.011    0.074

exhaustion      hours
     0.204      0.215

Chi-Squared Difference Test

     Df    AIC    BIC  Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
fit   0 4399.6 4445.5 0.0000
fit0  1 4397.8 4439.5 0.1303    0.13027     0       1     0.7182

Output trimmed, not edited: a few setting lines, blank lines and a warning are left out. Read the paths first: Level 1 [within] holds the within-team path (0.038 per hour); Level 2 [team] holds the climate path (−0.405) and the team-level hours path (0.050). contextual is their difference from the within-team path, with its own standard error. Std.all standardizes each path with its own level's standard deviations. Skip the 0.074 under Std.all for contextual: it's the team-level standardized path minus the within-team one, two numbers on different scales, so it doesn't mean anything. The model's test statistic is 0 on 0 degrees of freedom: it reproduces the data exactly, so don't read fit indexes for it. The Chi-Squared Difference Test compares it with the constrained model: χ²(1) = 0.13, p = .718. The warning left out says level-specific fit measures aren't available, because the constraint crosses levels.

Two syntax points. hours ~~ climate makes lavaan estimate the team-level variances of hours and climate instead of treating them as fixed, which gives the full-likelihood results on this page; without it, the paths are the same and their z-values differ in the third decimal. And an employee-level variable named only under level: 1 is treated as within-only, so any variable you want split into both parts must also appear under level: 2, as hours does here. For non-normal data, add estimator = "MLR".

Output from a real run: R 4.6.1 (lavaan 0.7-2), run 2026-09-29.

Python

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

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

# Variance split and ICC for exhaustion: a random-intercept model with no predictors (ML)
null = smf.mixedlm("exhaustion ~ 1", teams, groups="team").fit(reml=False)
between, within = null.cov_re.iloc[0, 0], null.scale
print(f"between {between:.4f}  within {within:.4f}  ICC {between / (between + within):.4f}")

# Hours (centered at the grand mean) plus the observed team mean: the hours_team slope is the contextual effect
teams["hours_c"] = teams["hours"] - teams["hours"].mean()
teams["hours_team"] = teams.groupby("team")["hours_c"].transform("mean")
m = smf.mixedlm("exhaustion ~ hours_c + hours_team + climate", teams, groups="team").fit(reml=False, method=["powell"])
print("converged:", m.converged)
print(m.summary().tables[1].round(3))

between 0.2023  within 0.7896  ICC 0.2040
converged: True
             Coef. Std.Err.       z  P>|z|  [0.025  0.975]
Intercept    4.010    0.306  13.113  0.000   3.411   4.610
hours_c      0.038    0.008    4.676  0.000   0.022   0.055
hours_team   0.007    0.020    0.361  0.718  -0.031   0.045
climate     -0.413    0.085   -4.839  0.000  -0.580  -0.246
team Var     0.076    0.038

Python has no mainstream package for two-level SEM, so this is the nearest thing: a multilevel regression. The variance split and ICC match the page exactly, and so does the within-team hours slope. Entering an employee's own hours and their team's average hours together makes the hours_team coefficient the contextual effect directly: 0.007 points per hour, p = .718. It uses observed team averages rather than latent means, and here it comes out smaller than the SEM's 0.011, and climate's path differs slightly (−0.413 against −0.405). With statsmodels' default optimizer this model failed to converge; method=["powell"] fixed it, so always check converged.

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

Impact

Effect Size

In a multilevel model, "how big" has a level attached. Each path has an unstandardized size, in the outcome's own units, and a standardized one that uses that level's standard deviations. Because the within-team and between-team standard deviations differ, the two levels' standardized paths aren't on one scale; neither are their R² values, which are shares of each level's own variance.

b^*_W = b_W \times \dfrac{\sigma_{x,W}}{\sigma_{y,W}} = 0.038 \times \dfrac{5.14}{0.89} = .22, \qquad b^*_{B} = b_{B,\text{climate}} \times \dfrac{\sigma_{z,B}}{\sigma_{y,B}} = -0.405 \times \dfrac{0.68}{0.45} = -.61

Within teams, 10 extra hours go with 0.38 more points of exhaustion (.22 SD per SD), and 5% of the within-team variance in exhaustion lines up with hours. Between teams, one climate point goes with 0.41 fewer points of team-average exhaustion among teams working the same hours (−.61 SD per SD), and climate and team hours together predict 63% of the between-team variance.

The between-level standardized path looks much bigger, but it's standardized by the small spread of true team means (SD 0.45 points). In the scale's own units, a colleague working 10 hours more than their teammates is about as much more exhausted (0.38 points) as a team one climate point worse (0.41 points). There are no agreed benchmarks for either; judge them in points of exhaustion. The contextual effect is a size too: 0.011 points per hour of team-average hours, with an interval (−0.050 to 0.072) wide enough to include nothing and something several times larger. The clustering ICC is another: 20% of the variation in exhaustion lies between teams, enough that, for a quantity such as the overall mean exhaustion across all employees, ignoring the teams would treat 480 employees as worth about 2.4 times as many independent people as they are. Within-team paths don't pay that price: they draw on every employee.

There are no agreed benchmarks for standardized paths in multilevel models. Borrowing correlation benchmarks is especially shaky at the between level, where the SDs are those of the latent team means, which are small. Judge the unstandardized path in the outcome's units, and compare with similar studies.

Statistical Significance

This page has two kinds of p. For a path, p asks: if the path were 0 in the population (and the model right), how often would an estimate land at least this many standard errors from 0, in either direction? For the model comparison, p asks: if the team-level and within-team hours slopes were equal in the population, how often would constraining them to be equal cost at least this much log-likelihood?

The climate path is 4.27 standard errors below 0 and the within-team hours path 4.68 above it, both p < .001. Constraining the team-level hours slope to equal the within-team one gives χ²(1) = 0.13, p = .718, so the data give no evidence of a contextual effect. They don't show there isn't one: the interval for it runs from −0.050 to 0.072 points per hour.

There's no p for the model as a whole. With 0 degrees of freedom, its χ² is 0 by construction, and fit indexes such as RMSEA, CFI and SRMR carry no information. When a model does leave degrees of freedom at a level, lavaan also reports level-specific fit, with the other level saturated, because overall fit can hide misfit at the group level (Ryu & West, 2009). Between-level p-values rest on the number of groups, and maximum likelihood's standard errors are large-sample approximations, so treat them as rougher than their decimals suggest when groups are few (Maas & Hox, 2005; Hox et al., 2010). p is not the probability that a path is zero, and a non-significant contextual-effect test is not evidence that the effect is zero.

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 SEM, power differs by level. Within-level paths draw on every person and usually have plenty. Between-level paths and contextual effects draw on the groups, so the number of groups matters far more than the number of people in each. In simulations of multilevel regression, second-level standard errors, those of the group-level variances in particular, came out too small with 50 groups or fewer (Maas & Hox, 2005); for multilevel SEM, how many groups the between part needs also depends on the estimation method (Hox et al., 2010).

We simulated studies like ours from the model that generated the data: a climate path of −0.40 points per climate point, a contextual effect of 0.02 points per hour, teams of 8. The climate path came out significant 74% of the time with 30 teams, 95% with our 60, and nearly always with 100. The contextual effect was detected 8% of the time with 30 teams, 10% with 60, 14% with 100 and 23% with 200. A contextual effect that small would need well over a thousand teams of 8 (roughly 1,100 for 80% power, since the standard error shrinks only with the square root of the number of teams). With realistic numbers of groups, plan for the smallest contextual effect you could actually detect, or report it as an estimate with its interval.

Plan the number of groups first, for the smallest group-level path or contextual effect that would matter. There's no calculator for this test; simulation is the standard approach, generating many data sets from the model you expect and refitting it to each. We refitted the full two-level model, with latent team means, to every simulated data set using our own maximum-likelihood code in Python, since no mainstream Python package fits it. In R you can loop lavaan over simulated data sets; Mplus has a Monte Carlo facility built in.

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

Related DesignMultilevel Regression (HLM): A multilevel regression answers some of the same questions with one outcome and observed group means; with a person-level predictor and its group mean, it estimates the same contextual effect, shrunk toward zero when groups are small or the ICC is low, because observed means pull the group-level slope toward the within-group one. Multilevel SEM adds latent group means, several outcomes and paths between them at each level, latent variables, and tests of constraints across the whole model (Preacher et al., 2010).
Related DesignConfirmatory Factor Analysis (CFA): The same two-level machinery fits a factor model at each level: a multilevel CFA asks whether a scale's items hang together the same way within groups and between them.
Related DesignIntraclass Correlation Coefficient (ICC): The ICCs here are the clustering ICC, ICC(1) from the one-way random-effects model: the share of a variable's variance between teams. The reliability ICCs on that page use the same kind of variance ratio for raters, with different benchmarks.
Special CasePath Analysis: Path analysis with one level. Multilevel SEM fits a path model twice, once to the within-group parts and once to the latent group means, and estimates both together.

References

Muthén, B. O. (1994). Multilevel covariance structure analysis. Sociological Methods & Research, 22(3), 376–398. https://doi.org/10.1177/0049124194022003006
“This article gives an introduction to some new techniques for multilevel covariance structure modeling with latent variables.” (Abstract)
Lüdtke, O., Marsh, H. W., Robitzsch, A., Trautwein, U., Asparouhov, T., & Muthén, B. (2008). The multilevel latent covariate model: A new, more reliable approach to group-level effects in contextual studies. Psychological Methods, 13(3), 203–229. https://doi.org/10.1037/a0012869
“this MMC approach can result in substantially biased estimates of contextual effects and can substantially underestimate the associated standard errors” (Abstract)

“group-level (L2) characteristics are often measured by aggregating individual-level (L1) characteristics within each group so as to assess contextual effects” (Abstract)
Preacher, K. J., Zyphur, M. J., & Zhang, Z. (2010). A general multilevel SEM framework for assessing multilevel mediation. Psychological Methods, 15(3), 209–233. https://doi.org/10.1037/a0020141
“these MLM approaches do not accommodate mediation pathways with Level-2 outcomes and may produce conflated estimates of between- and within-level components of indirect effects” (Abstract)
Robinson, W. S. (1950). Ecological correlations and the behavior of individuals. American Sociological Review, 15(3), 351–357. https://doi.org/10.2307/2087176
“the only reasonable assumption is that an ecological correlation is almost certainly not equal to its corresponding individual correlation” (Final section (read in the 2009 reprint, International Journal of Epidemiology, 38(2), 337–341))
Ryu, E., & West, S. G. (2009). Level-specific evaluation of model fit in multilevel structural equation modeling. Structural Equation Modeling, 16(4), 583–601. https://doi.org/10.1080/10705510903203466
“A simulation study showed that the standard approach failed to detect the lack of fit at the group level.” (Abstract)
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)
Hox, J. J., Maas, C. J. M., & Brinkhuis, M. J. S. (2010). The effect of estimation method and sample size in multilevel structural equation modeling. Statistica Neerlandica, 64(2), 157–170. https://doi.org/10.1111/j.1467-9574.2009.00445.x
“The within‐group part of the model posed no problems.” (Abstract)
Rosseel, Y. (2012). lavaan: An R package for structural equation modeling. Journal of Statistical Software, 48(2), 1–36. https://doi.org/10.18637/jss.v048.i02
Wright, S. (1934). The method of path coefficients. The Annals of Mathematical Statistics, 5(3), 161–215. https://doi.org/10.1214/aoms/1177732676
“as a flexible means of relating the correlation coefficients between variables in a multiple system to the functional relations among them” (Abstract)
Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical linear models: Applications and data analysis methods (2nd ed.). Sage.
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)
Halekoh, U., Højsgaard, S., & Yan, J. (2006). The R package geepack for generalized estimating equations. Journal of Statistical Software, 15(2), 1–11. https://doi.org/10.18637/jss.v015.i02
“For a small number of clusters (K ≤ 30) the sandwich variance estimator exhibits bias and Paik (1988) recommended using the jackknife variance estimator.” (p. 6)
Mancl, L. A., & DeRouen, T. A. (2001). A covariance estimator for GEE with improved small-sample properties. Biometrics, 57(1), 126–134. https://doi.org/10.1111/j.0006-341X.2001.00126.x
“The bias-corrected covariance estimator gave tests with sizes close to the nominal level even when the number of subjects was 10 and cluster sizes were unequal, whereas the robust and jackknife covariance estimators gave tests with sizes that could be 2–3 times the nominal level.” (Summary)
Moineddin, R., Matheson, F. I., & Glazier, R. H. (2007). A simulation study of sample size for multilevel logistic regression models. BMC Medical Research Methodology, 7, 34. https://doi.org/10.1186/1471-2288-7-34
“The estimates of the variance covariance components are slightly biased even with 100 groups and group size of 50.” (Abstract)
MacCallum, R. C., Roznowski, M., & Necowitz, L. B. (1992). Model modifications in covariance structure analysis: The problem of capitalization on chance. Psychological Bulletin, 111(3), 490–504. https://doi.org/10.1037/0033-2909.111.3.490
“The use of alternative a priori models is recommended as a preferred strategy.” (Abstract)