Multi-level Models

Overview

The ml_lvm family implements two-level random-intercept latent variable models for clustered data — for example, students nested in schools, employees nested in teams, or repeated measures nested in persons (when temporal ordering is not of interest; for panel models with dynamics, see the DLVM1 family).

The model decomposes each observation into a within-cluster component (how units deviate from their cluster average) and a between-cluster component (how cluster averages differ from each other), each with its own measurement model and latent structure. As everywhere in psychonetrics, every latent and residual covariance matrix on either level can independently be parameterized as a Gaussian Graphical Model.

This page first covers the latent-variable family (ml_lvm and its network variants); its observed-variable counterpart — ml_varcov, ml_ggm and ml_corr, which model the within- and between-cluster (co)variance structures directly, without a latent layer — is described below. For multi-level time-series (ESM) data, see the multi-level VAR family (ml_var / ml_gvar).

The Two-Level Model

For unit $p$ in cluster $c$, the observed response vector is decomposed as:

$$\boldsymbol{y}_{pc} = \boldsymbol{\nu} + \boldsymbol{\Lambda}\left(\boldsymbol{\eta}^{(B)}_{c} + \boldsymbol{\eta}^{(W)}_{pc}\right) + \boldsymbol{\varepsilon}^{(B)}_{c} + \boldsymbol{\varepsilon}^{(W)}_{pc}$$

• $\boldsymbol{\Lambda}$ (lambda): factor loadings, shared across levels.
• $\boldsymbol{\eta}^{(W)}_{pc}$, $\boldsymbol{\eta}^{(B)}_{c}$: within- and between-cluster latent variables, with covariance structures controlled by within_latent and between_latent ("cov", "chol", "prec", or "ggm").
• $\boldsymbol{\varepsilon}^{(W)}$, $\boldsymbol{\varepsilon}^{(B)}$: residuals on each level, controlled by within_residual and between_residual.
• Structural regressions are possible on both levels via beta_within and beta_between.

Key matrices include sigma_zeta_within / sigma_zeta_between (or omega_zeta_within / omega_zeta_between for GGM parameterizations) and sigma_epsilon_within / sigma_epsilon_between. Estimation handles unbalanced cluster sizes naturally, and as of version 0.16 a fast sufficient-statistics two-level ML estimator is used by default for complete data.

Wrapper Functions

The cluster membership is supplied via the clusters argument (a column name in the data).

Example: Two-Level CFA on Simulated Data

We simulate data for 60 clusters of 10 units each, with two latent factors measured by six indicators on both levels, and recover the structure with ml_lvm(). The example is fully reproducible and runs in a few seconds:

library("psychonetrics")
library("dplyr")

# --- Simulate two-level data ---
set.seed(2)
nCluster <- 60; perCluster <- 10; nv <- 6
cl <- rep(seq_len(nCluster), each = perCluster)

# True loadings (2 factors, 3 indicators each):
loadings <- matrix(0, nv, 2)
loadings[1:3, 1] <- c(1, .8, .6)
loadings[4:6, 2] <- c(1, .8, .6)

# Between-cluster latents, within latents, and residuals on both levels:
eta_b <- matrix(rnorm(nCluster * 2, 0, sqrt(0.5)), nCluster, 2)
eta_w <- matrix(rnorm(nCluster * perCluster * 2), nCluster * perCluster, 2)
u_b   <- matrix(rnorm(nCluster * nv, 0, sqrt(0.3)), nCluster, nv)
Y <- (eta_b[cl, ] + eta_w) %*% t(loadings) + u_b[cl, ] +
     matrix(rnorm(nCluster * perCluster * nv, 0, sqrt(0.5)), ncol = nv)

Data <- as.data.frame(Y)
colnames(Data) <- paste0("y", 1:nv)
Data$cluster <- cl

# --- Fit the two-level CFA ---
Lambda <- matrix(0, nv, 2)
Lambda[1:3, 1] <- 1
Lambda[4:6, 2] <- 1

mod_ml <- ml_lvm(Data, lambda = Lambda, clusters = "cluster",
                 vars = paste0("y", 1:nv)) %>% runmodel

mod_ml %>% fit

# Latent covariances on both levels:
round(getmatrix(mod_ml, "sigma_zeta_between"), 2)
round(getmatrix(mod_ml, "sigma_zeta_within"), 2)

The model fits well (χ²(20) = 15.92, p = 0.72) and recovers the simulated latent variances on both levels (between ≈ 0.5 scaled by the marker loading, within ≈ 1).

Estimation: two-level ML vs FIML 0.16

Historically ml_lvm() estimated the model by full-information maximum likelihood (FIML) on the wide (one row per cluster) data layout. That is exact, but it scales poorly: each cluster contributes a likelihood term whose dimension grows with its size, so large clusters are expensive. Version 0.16 adds a sufficient-statistics two-level ML estimator (estimator = "ML"; McDonald & Goldstein, 1989; Muthen, 1990) that optimizes the numerically identical objective but works from the within- and between-cluster sufficient statistics, with full C++ standard errors and modification indices.

The estimator argument now defaults to "default", which picks "ML" for complete data whose largest cluster exceeds 5 units, and "FIML" otherwise. Explicit estimator = "FIML" reproduces previous results exactly.

# Same simulation as above, but with larger clusters (perCluster <- 18):
mod_ml   <- ml_lvm(Data, lambda = Lambda, clusters = "cluster",
                   vars = paste0("y", 1:nv), estimator = "ML")   %>% runmodel
mod_fiml <- ml_lvm(Data, lambda = Lambda, clusters = "cluster",
                   vars = paste0("y", 1:nv), estimator = "FIML") %>% runmodel

mod_ml   %>% fit
mod_fiml %>% fit

The two estimators reach the same fit; "ML" is about 3× faster here, and the gap widens sharply as clusters grow (the NEWS reports a 100-cluster model with cluster sizes 18–25 fitting in under a second versus ~90 seconds with FIML).

Network Variants

Replacing a covariance structure with a GGM turns the corresponding level into a network model. For example, ml_lnm() models the between-cluster latent structure as a partial-correlation network:

# Multi-level latent network model:
mod_mlnm <- ml_lnm(Data, lambda = Lambda, clusters = "cluster",
                   vars = paste0("y", 1:nv)) %>% runmodel

# Within-level latent network (partial correlations):
round(getmatrix(mod_mlnm, "omega_zeta_within"), 2)

# Prune non-significant network edges as usual:
mod_mlnm %>% prune(alpha = 0.05)

The same prune() / stepup() / compare() workflow applies as in all other psychonetrics families.

Observed-variable models: ml_varcov, ml_ggm and ml_corr

The ml_varcov family is the two-level analogue of the single-level varcov family: it models the within-cluster and between-cluster variance–covariance matrices of clustered data directly, with no latent layer. It is the observed-variable counterpart of ml_lvm — formally the special case with identity factor loadings and no residual variance — but it is implemented as its own framework, with its own native matrix names and its own simpler analytic derivatives. The headline use case is the multi-level Gaussian graphical model (ml_ggm), which estimates a within-cluster and a between-cluster partial-correlation network simultaneously.

The Two-Level Model

For unit $j$ in cluster $i$, the observed response vector is decomposed into independent within- and between-cluster components:

$$\boldsymbol{y}_{ij} = \boldsymbol{\mu} + \boldsymbol{\eta}^{(W)}_{ij} + \boldsymbol{\eta}^{(B)}_{i}, \qquad \boldsymbol{\eta}^{(W)}_{ij} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}_{W}), \quad \boldsymbol{\eta}^{(B)}_{i} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}_{B})$$

with $\boldsymbol{\eta}^{(W)}_{ij}$ and $\boldsymbol{\eta}^{(B)}_{i}$ mutually independent. Because all units in a cluster share the same between-cluster component, this induces a simple two-level covariance structure:

$$\text{var}(\boldsymbol{y}_{ij}) = \boldsymbol{\Sigma}_{W} + \boldsymbol{\Sigma}_{B}, \qquad \text{cov}(\boldsymbol{y}_{ij}, \boldsymbol{y}_{ij'}) = \boldsymbol{\Sigma}_{B} \quad (j \neq j')$$

• $\boldsymbol{\Sigma}_{W}$ (sigma_within): the within-cluster covariance, parameterized through the within argument.
• $\boldsymbol{\Sigma}_{B}$ (sigma_between): the between-cluster covariance, parameterized through the between argument.
• Each of within and between can independently be any of the five varcov parameterizations — "cov", "chol", "prec", "ggm", or "cor" — exactly like the single-level type argument.

For the "ggm" parameterization each level is written as a Gaussian graphical model,

$$\boldsymbol{\Sigma} = \boldsymbol{\Delta}(\boldsymbol{I} - \boldsymbol{\Omega})^{-1}\boldsymbol{\Delta},$$

with a partial-correlation network $\boldsymbol{\Omega}$ (omega_within / omega_between) and diagonal scaling $\boldsymbol{\Delta}$ (delta_within / delta_between). The remaining parameterizations expose their usual matrices with _within / _between suffixes (lowertri_*, kappa_*, rho_*, SD_*), and the mean structure is the single vector mu.

Both ml_lvm estimators are supported and behave identically to their ml_lvm equivalents: the fast two-level sufficient-statistics Gaussian ML estimator (estimator = "ML"; cost independent of cluster size, unbalanced clusters handled exactly) and wide-format full-information ML (estimator = "FIML"; the natural choice with missing data). The two maximize the same likelihood and give identical estimates on complete data. The analytic gradients of both were verified against numeric differentiation, against an independent C++ implementation, and against ml_lvm itself (identical fit function and gradient at identical parameter values, to machine precision).

Wrapper Functions

As in ml_lvm, the variables are supplied via vars and cluster membership via the clusters argument (a column name in the data). Multi-group models (groups) and cross-group equality constraints (equal) work exactly as elsewhere in psychonetrics.

Example: Multi-Level GGM on Simulated Data

We simulate 200 clusters of 8 units on four variables, giving the within-cluster covariance a single edge (between V1 and V2) and the between-cluster covariance a different edge (between V3 and V4), then recover both networks with ml_ggm(). The example is fully reproducible and runs in a couple of seconds:

library("psychonetrics"); library("dplyr")

# --- Simulate two-level data ---
set.seed(1)
nCluster <- 200; perCluster <- 8; p <- 4
cl <- rep(seq_len(nCluster), each = perCluster)

SigmaW <- diag(p); SigmaW[1,2] <- SigmaW[2,1] <- 0.4        # within-cluster edge V1-V2
SigmaB <- diag(0.5, p); SigmaB[3,4] <- SigmaB[4,3] <- 0.2   # between-cluster edge V3-V4

b <- MASS::mvrnorm(nCluster, rep(0, p), SigmaB)
Y <- b[cl, ] + MASS::mvrnorm(nCluster * perCluster, rep(0, p), SigmaW)
Data <- as.data.frame(Y); names(Data) <- paste0("V", 1:p); Data$cluster <- cl

# --- Fit the multi-level GGM ---
mod <- ml_ggm(Data, vars = paste0("V", 1:p), clusters = "cluster") %>% runmodel

# Within- and between-cluster partial-correlation networks:
getmatrix(mod, "omega_within")
getmatrix(mod, "omega_between")

# Prune non-significant edges as usual:
mod <- mod %>% prune(alpha = 0.05)

getmatrix() returns the two partial-correlation networks: here the within-cluster network recovers the V1–V2 edge and the between-cluster network the V3–V4 edge, matching the simulation. The full psychonetrics toolbox — fit indices, prune() / stepup() search, and compare() for confirmatory testing — applies exactly as for the single-level ggm().

Summary

• ml_lvm() fits two-level random-intercept latent variable models for clustered data; a fast sufficient-statistics two-level ML estimator (default) optimizes the same objective as FIML but much faster for large clusters
• Within-cluster and between-cluster structures each get their own latent and residual covariance matrices
• Any of these matrices can be parameterized as a network (GGM): ml_lnm(), ml_rnm(), ml_lrnm()
• Cluster membership is a single column passed via clusters
• For clustered data on observed variables with no latent layer, the ml_varcov / ml_ggm / ml_corr family models the within- and between-cluster (co)variance structures and networks directly
• For clustered longitudinal data with temporal dynamics, use the DLVM1 / panel family instead