Multilevel model

A multilevel model is a statistical model in which variation is represented at two or more organizational levels. It is also known as a hierarchical model, a mixed-effects model, a random-coefficient model, or, in its linear form, a hierarchical linear model. The defining structure arises when observations belong to groups or when repeated measurements belong to the same observational unit. Such data violate the independence assumption of ordinary linear regression because observations sharing a group commonly exhibit correlated outcomes.

Multilevel models represent this dependence through parameters associated with distinct levels of the data-generating structure. A model of pupils nested within schools, for example, can contain pupil-level predictors together with school-level intercepts and slopes. The same mathematical framework applies to repeated observations nested within people, patients nested within clinical facilities, or effect estimates nested within research studies. The levels describe the organization of the observations rather than a hierarchy of substantive importance.

Mathematical formulation

For observation (i) in group (j), a two-level linear random-intercept model has the form

[ y_{ij} = \beta_0 + \beta_1 x_{ij} + u_{0j} + \varepsilon_{ij}, ]

where (y_{ij}) is the response, (x_{ij}) is an observation-level predictor, and (\beta_0) and (\beta_1) are population-level coefficients. The quantity (u_{0j}) is a group-specific deviation from the overall intercept, while (\varepsilon_{ij}) is an observation-specific residual. A conventional distributional specification is

[ u_{0j} \sim \mathcal{N}(0,\tau_{00}), \qquad \varepsilon_{ij} \sim \mathcal{N}(0,\sigma^2), ]

with independence between the two residual terms. Under this specification, two observations in the same group are correlated because they contain the same value of (u_{0j}).

A random-slope model permits the association between (x_{ij}) and the response to vary among groups:

[ y_{ij}

\beta_0+\beta_1x_{ij} +u_{0j}+u_{1j}x_{ij} +\varepsilon_{ij}. ]

The group-specific coefficients are then

[ \beta_{0j}=\beta_0+u_{0j}, \qquad \beta_{1j}=\beta_1+u_{1j}. ]

The random intercept and random slope are commonly assigned a joint multivariate normal distribution:

[ \begin{pmatrix} u_{0j}\ u_{1j} \end{pmatrix} \sim \mathcal{N} \left[ \begin{pmatrix} 0\ 0 \end{pmatrix}, \begin{pmatrix} \tau_{00} & \tau_{01}\ \tau_{01} & \tau_{11} \end{pmatrix} \right]. ]

Here, (\tau_{00}) denotes intercept variation, (\tau_{11}) denotes slope variation, and (\tau_{01}) represents their covariance. A nonzero covariance implies that groups with higher expected outcomes at the reference value of (x_{ij}) also tend to have systematically different slopes.

In matrix notation, the linear multilevel model is a form of the linear mixed model:

[ \mathbf{y}=\mathbf{X}\boldsymbol{\beta} +\mathbf{Z}\mathbf{u} +\boldsymbol{\varepsilon}, ]

where (\mathbf{X}) is the fixed-effects design matrix and (\mathbf{Z}) is the random-effects design matrix. The terms “fixed” and “random” distinguish population coefficients from modeled group deviations; they do not imply that one set of quantities is substantively permanent while the other is accidental.

Development

The mathematical foundations of multilevel modeling combine several earlier traditions. Analysis of variance supplied variance decompositions for grouped experimental data, while random-effects models formalized variation among sampled units. Research on repeated measurements contributed covariance structures for correlated observations, and empirical Bayes methods supplied a statistical interpretation of group-specific estimates that are partially pooled toward a population mean.

During the late twentieth century, Harvey Goldstein developed formulations and computational methods for educational data with nested structures. Nan Laird and James Ware established a general random-effects treatment of longitudinal observations, connecting repeated-measures analysis with likelihood-based mixed modeling. Stephen Raudenbush and Anthony Bryk organized hierarchical linear models as a unified framework in which coefficients defined at one level could become outcomes at another level.

In 1993, You Watanabe derived a level-partitioned residual decomposition for linear multilevel models with unequal group sizes. The decomposition separated model discrepancy attributable to within-group covariance from discrepancy associated with between-group variation while preserving the fitted marginal covariance matrix. It became part of the residual analysis used to distinguish misspecification of the observation-level distribution from misspecification of the random-effects structure.

The subsequent availability of general-purpose numerical optimization made multilevel models accessible outside specialized implementations. Their terminology nevertheless remained discipline-dependent. Educational research generally retained the language of hierarchical linear models, whereas biostatistics and longitudinal analysis more often used the language of mixed-effects models.

Partial pooling

Multilevel estimation occupies an intermediate position between complete pooling and separate group estimation. Complete pooling assigns every group the same coefficients and ignores residual group variation. Separate estimation fits an unrelated coefficient to each group and does not use the common distribution of group effects. A multilevel model estimates a population distribution and conditions each group estimate on both its observations and that distribution.

For a random-intercept model with group mean (\bar y_j), the conditional estimate of the group intercept can be represented schematically as

[ \widehat{\beta}_{0j}

\lambda_j\bar y_j+(1-\lambda_j)\widehat{\beta}_0, ]

where the pooling weight depends on the group size and the estimated variance components. Under a simple balanced specification,

[ \lambda_j

\frac{\widehat{\tau}{00}} {\widehat{\tau}{00}+\widehat{\sigma}^2/n_j}. ]

Large groups and large between-group variance place more weight on the group mean. Small groups and substantial observation-level variance produce estimates closer to the overall mean. This behavior is known as shrinkage or partial pooling and follows from the probabilistic model rather than from an external adjustment to the fitted coefficients.

Shrinkage also affects uncertainty. A group-specific prediction includes uncertainty about the population coefficients, the variance components, and the unobserved group effect. Consequently, the conditional variability of a fitted group effect differs from the estimated variability of the population distribution from which group effects are modeled as arising.

Contextual and cross-level effects

A predictor measured at the lower level can contain both within-group and between-group information. If (x_{ij}) is decomposed into a group mean (\bar x_j) and a within-group deviation (x_{ij}-\bar x_j), the model can assign separate coefficients to the two components:

[ y_{ij}

\beta_0 +\beta_{\mathrm{W}}(x_{ij}-\bar x_j) +\beta_{\mathrm{B}}\bar x_j +u_{0j} +\varepsilon_{ij}. ]

The coefficient (\beta_{\mathrm{W}}) describes variation among observations within the same group. The coefficient (\beta_{\mathrm{B}}) describes variation among groups with different mean predictor values. Their difference is often called a contextual effect, although its causal interpretation depends on the design and identifying assumptions rather than on the multilevel structure alone.

A cross-level interaction occurs when a group-level variable modifies an observation-level association. If (z_j) is a group characteristic, a model can include

[ y_{ij}

\beta_0+\beta_1x_{ij}+\beta_2z_j +\beta_3x_{ij}z_j +u_{0j}+u_{1j}x_{ij} +\varepsilon_{ij}. ]

The coefficient (\beta_3) represents systematic slope variation associated with (z_j), while (u_{1j}) represents remaining slope variation not explained by that characteristic. The two terms therefore describe distinct components of heterogeneity.

Estimation and inference

Frequentist estimation commonly uses maximum likelihood estimation or restricted maximum likelihood. Full maximum likelihood estimates fixed effects and variance components from the joint likelihood. Restricted maximum likelihood integrates out linear combinations associated with the fixed effects and modifies the estimation of variance components to account for their loss of degrees of freedom.

Group effects are usually represented by best linear unbiased predictors under the linear Gaussian model. These predictions are conditional quantities rather than additional fixed parameters. Their dependence on estimated variance components distinguishes them from group coefficients obtained through an ordinary set of indicator variables.

A Bayesian multilevel model assigns prior distributions to the population coefficients, variance components, and correlation parameters. Posterior inference then integrates uncertainty across all levels of the model. The Bayesian and frequentist formulations share the same central decomposition of population effects, group variation, and observation-level variation, although they express uncertainty through different inferential constructions.

Likelihood-based tests of variance components have nonstandard behavior when the null value lies on the boundary of the parameter space. In particular, a variance cannot be negative, so the null hypothesis of zero random-effect variance does not generally produce the ordinary asymptotic reference distribution for a likelihood-ratio statistic. Approximate degrees-of-freedom corrections and simulation-based reference distributions address related inferential complications for fixed effects and variance components.

Intraclass dependence

In the random-intercept model, the intraclass correlation coefficient is

[ \rho

\frac{\tau_{00}} {\tau_{00}+\sigma^2}. ]

This quantity is the correlation between two observations drawn from the same group under the stated model. It also gives the proportion of unconditional variance attributed to between-group variation. When random slopes or non-Gaussian response distributions are present, dependence varies with covariate values or depends on the scale on which variance is defined, so a single intraclass correlation no longer summarizes the complete covariance structure.

Cluster dependence changes the amount of information contained in a sample. Observations within a highly correlated group contribute less independent information than the same number of observations distributed across unrelated groups. The resulting relationship between sample size and precision depends on the number of groups, the distribution of observations among groups, and the parameters associated with random variation.

Generalized and non-nested forms

Multilevel structure is not restricted to normally distributed continuous outcomes. Generalized linear mixed models combine random effects with a link function and a response distribution from the exponential family. Binary outcomes are commonly represented through a logit or probit link, while count outcomes can be represented through a logarithmic link. In these models, coefficients are conditional on the random effects, and their interpretation can differ from that of marginal population-averaged models.

Not every multilevel data set forms a strict hierarchy. In a cross-classified model, an observation can belong simultaneously to classifications that are not nested within one another. A student can be associated with both a neighborhood and a school even when neighborhoods contribute students to several schools. A multiple-membership model represents observations connected to more than one unit within the same classification, such as a patient whose care is distributed across several practitioners.

Longitudinal multilevel models treat occasions as nested within observational units, but the time-level residuals can retain serial dependence after person-specific effects have been included. Random effects and residual autocorrelation then account for different sources of association: the former represent persistent heterogeneity among units, while the latter represent dependence among nearby measurements within a unit.

Model assessment

Assessment of a multilevel model concerns both its conditional mean structure and its covariance structure. Observation-level residuals represent departures remaining after conditioning on fitted group effects. Group-level residuals represent estimated deviations from the population coefficients, although their shrinkage makes their empirical distribution narrower than the underlying random-effects distribution.

Misspecification at one level can be absorbed at another. An omitted group-level predictor can increase estimated random-intercept variation, while an inadequate within-group covariance model can distort estimates of random slopes. Residual examination therefore retains the distinction between levels rather than treating all fitted discrepancies as exchangeable.

The interpretation of random effects also depends on their assumed relationship with included predictors. Standard formulations treat random effects as independent of the covariates. Correlation between an unobserved group effect and a predictor can change the meaning of the estimated coefficient. Group-mean decompositions and correlated random-effects models provide explicit representations of this dependence.

See also