What's New

Robust maximum likelihood (MLM / MLMV / MLMVS / MLR) 0.16

June 2026 — see the new robust & scaled estimation page

psychonetrics now offers the full family of robust maximum-likelihood estimators, selectable through the model constructors or with setestimator(). The point estimates are ordinary maximum likelihood, but the inference is corrected for non-normality:

• "MLM", "MLMV", "MLMVS" — Browne (1984) robust sandwich standard errors with the Satorra–Bentler, scaled-and-shifted, and mean-and-variance-adjusted scaled chi-square, respectively (complete data).
• "MLR" — Huber–White sandwich standard errors with the Yuan–Bentler–Mplus scaled chi-square; supports missing data via FIML with the Savalei (2010) FIML-corrected robust fit indices.

All estimators also report robust RMSEA / CFI / TLI. Everything was validated against lavaan.

library("psychonetrics")
library("lavaan")
library("dplyr")
data(HolzingerSwineford1939)

Lambda <- matrix(0, 9, 3)
Lambda[1:3, 1] <- 1; Lambda[4:6, 2] <- 1; Lambda[7:9, 3] <- 1

# Robust ML (MLR): sandwich SEs + scaled chi-square + robust fit indices
mod_mlr <- lvm(HolzingerSwineford1939, lambda = Lambda,
               vars = paste0("x", 1:9),
               latents = c("visual", "textual", "speed"),
               identification = "variance",
               estimator = "MLR") %>% runmodel

mod_mlr %>% fit

These match lavaan::cfa(…, estimator = "MLR") to three decimals (scaled χ² = 87.13, scaling factor 0.979, robust CFI 0.930, robust TLI 0.895). See the robust estimation page for MLM and the setestimator() workflow.

Fast multilevel SEM: two-level ML for ml_lvm() 0.16

June 2026 — see the multi-level page

ml_lvm() gains a sufficient-statistics two-level ML estimator (estimator = "ML"; McDonald & Goldstein, 1989; Muthen, 1990). It optimizes the numerically identical objective as the wide-format FIML estimator but is dramatically faster for models with large clusters, with full C++ standard errors and modification indices. The two are the same exact likelihood computed two ways (not an approximation, and not the pseudo-balanced “MUML”); the "ML"/"FIML" labels denote the computational route — summary statistics versus row by row, as elsewhere in psychonetrics — and both can handle within-cluster missing data (MCAR/MAR). The estimator argument now defaults to "default", which selects "ML" for complete data with the largest cluster > 5 units and "FIML" otherwise; explicit estimator = "FIML" reproduces previous results exactly.

# 60 clusters of 18 units, 2 factors x 3 indicators on each level
mod_ml   <- ml_lvm(Data, lambda = Lambda, clusters = "cluster",
                   vars = paste0("y", 1:6), estimator = "ML")   %>% runmodel
mod_fiml <- ml_lvm(Data, lambda = Lambda, clusters = "cluster",
                   vars = paste0("y", 1:6), estimator = "FIML") %>% runmodel

The two estimators give the same fit; "ML" is here about 3× faster, and the gap grows with cluster size. (CFI/RMSEA use the number of clusters J for the two-level estimator, matching lavaan's two-level SEM, so they can differ slightly from the FIML values.)

Interoperability with lavaan: fromlavaan() & tolavaan() 0.16

June 2026 — see the lavaan interoperability page

Two new converters bridge psychonetrics and lavaan. fromlavaan() turns a fitted lavaan object (or model syntax + data) into an equivalent psychonetrics lvm, reproducing estimates, standard errors, chi-square, df and log-likelihood for standard ML fits (multi-group, equality constraints, std.lv and FIML are supported). tolavaan() goes the other way, converting a psychonetrics lvm into a fitted lavaan object or parameter table — so you can specify in one framework and prune, search, or report in the other.

library("psychonetrics"); library("lavaan"); library("dplyr")
data(HolzingerSwineford1939)

HS.model <- 'visual  =~ x1 + x2 + x3
             textual =~ x4 + x5 + x6
             speed   =~ x7 + x8 + x9'
fit_lav <- cfa(HS.model, data = HolzingerSwineford1939, meanstructure = TRUE)

# lavaan  ->  psychonetrics:
mod_pn <- fit_lav %>% fromlavaan %>% runmodel

# ... and back again:
back   <- tolavaan(mod_pn)

Estimates, standard errors, chi-square and log-likelihood round-trip exactly. Unsupported lavaan features raise informative errors rather than failing silently.

Standardized solutions, the Wishart likelihood & fixed exogenous covariates 0.16

June 2026 — see the LVM and varcov pages

Three smaller additions to the lvm and varcov families round out the release.

Completely standardized solutions

parameters(x, standardized = TRUE) now returns the completely standardized (std.all) solution with delta-method standard errors (new std and se_std columns), matching lavaan::standardizedSolution():

mod_cfa <- lvm(HolzingerSwineford1939, lambda = Lambda,
               vars = paste0("x", 1:9),
               latents = c("visual", "textual", "speed"),
               identification = "variance") %>% runmodel

parameters(mod_cfa, standardized = TRUE)   # adds std + se_std columns

The std / se_std columns equal lavaan's est.std / se exactly.

The Wishart likelihood

likelihood = "wishart" for lvm() and varcov() uses the unbiased (n−1) sample covariance with an (n−1)-scaled chi-square and standard errors, matching lavaan(likelihood = "wishart"). The default likelihood = "normal" (the n denominator) is unchanged.

mod_w <- lvm(HolzingerSwineford1939, lambda = Lambda, vars = paste0("x", 1:9),
             latents = c("visual", "textual", "speed"),
             identification = "variance",
             likelihood = "wishart") %>% runmodel
mod_w %>% fit    # chisq = 85.022, df = 24  (matches lavaan likelihood = "wishart")

Fixed exogenous covariates

fixed_x fixes declared exogenous covariates' means and mutual (co)variances to their sample values (conditioning on x), matching lavaan(fixed.x = TRUE). In varcov() it names exogenous observed variables; in lvm() it names the single-indicator latents carrying those covariates.

# Same covariance model, with x3 and x4 treated as fixed exogenous covariates:
m0 <- varcov(HolzingerSwineford1939, vars = c("x1","x2","x3","x4"),
             type = "cov", sigma = S) %>% runmodel
m1 <- varcov(HolzingerSwineford1939, vars = c("x1","x2","x3","x4"),
             type = "cov", sigma = S, fixed_x = c("x3","x4")) %>% runmodel

Estimates, standard errors, chi-square, df and the conditional log-likelihood match lavaan(fixed.x = TRUE); incremental fit indices use psychonetrics' unconditional baseline.

Scaled chi-square difference tests 0.16

June 2026 — see the robust estimation page

For nested models fitted with a robust or scaled estimator, compare() now reports a Satorra–Bentler-family scaled chi-square difference test (new Chisq_diff_scaled and p_value_scaled columns). A new scaled.test.method argument selects "satorra.bentler.2001" (default), "satorra.bentler.2010", or "satorra.2000" (scaled-and-shifted; appropriate for MLMV/WLSMV). Output for non-robust estimators is unchanged. This reproduces lavaan's lavTestLRT().

# Two nested MLR models (configural vs equal loadings across schools):
mod_config <- lvm(HolzingerSwineford1939, lambda = Lambda, vars = paste0("x", 1:9),
                  latents = c("visual","textual","speed"),
                  identification = "variance",
                  groups = "school", estimator = "MLR") %>% runmodel
mod_weak   <- mod_config %>% groupequal("lambda") %>% runmodel

compare(configural = mod_config, weak = mod_weak)

The scaled difference (6.57, p = 0.36) matches lavaan::lavTestLRT() (6.567, p = 0.363).

Stability & performance release 0.16

June 2026

Version 0.16 implements the findings of a large internal code audit: bug fixes across all model families, robustness improvements, substantial speedups, and a new automated test suite. No user action is needed — existing scripts keep working, but many run faster and more reliably.

Robustness

• Panel-model divergence fixed. The optimizer penalty for improper (non-positive-definite) implied covariance matrices was re-enabled; previously the default optimizer could diverge to absurd estimates on panelvar() / panelgvar() and SEM models.
• Row-wise FIML (fullFIML = TRUE) works again for all standard model families.
• Standard errors: inversion problems are no longer silently suppressed, and negative variances now yield NA standard errors with a clear warning.
• Many further fixes: factorscores() with structural regressions, multiple-comparison adjustment in prune(), getmatrix(diag = FALSE), equal = "omega_zeta" in latent network models, multi-group EPC columns, group-stratified bootstrap(), and more.

Performance

• About 3× faster model construction (e.g., a 30-node GGM) through sparse, cached algebra helper matrices.
• Up to several hundred times faster equality-constrained modification indices on large multi-group models (used by runmodel() and addMIs()), with identical results.
• About 1.3–1.6× faster Ising estimation through a reworked state-enumeration loop.

Quality assurance

• A new automated test suite (tinytest) covers the examples on this website plus targeted regression tests, with reference values cross-checked against lavaan.
• Large documentation cleanup and removal of unused dependencies.

Meta-analysis: covariance input & Olkin–Siotani V methods 0.16

June 2026 — see the new meta-analysis page

meta_varcov() now genuinely supports covariance matrices as input (previously covariances were silently modeled as correlations), and two new Vmethod options — "OS_individual" and "OS_pooled" — implement the Olkin–Siotani (1976) sampling covariance of sample correlations, validated against metaSEM::asyCov. These are recommended when the between-study heterogeneity (random-effects) variances are of substantive interest. The default Vmethod is unchanged, so published analyses reproduce exactly.

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

# cors: a list of per-study correlation matrices; nobs: their sample sizes
# (see the meta-analysis page for a full simulated example)
mod_meta <- meta_ggm(cors = cors, nobs = nobs,
                     Vmethod = "OS_pooled") %>% runmodel

# Pooled partial-correlation network:
getmatrix(mod_meta, "omega_y")

# Prune to a sparse pooled network:
mod_meta %>% prune(alpha = 0.05, matrices = "omega_y")

Also in 0.16: meta_varcov() / meta_ggm() with estimator = "ML" and bivariate (single-correlation) meta-analyses no longer crash, and meta models give an informative error when two variables are never observed together in any study.

Blume-Capel & multi-state Ising models 0.16

See the updated Ising & Blume-Capel page

The Ising model has been generalized in two directions. First, Ising() now accepts any number of ordered response options (e.g., c(-1, 0, 1) or a Likert scale), not just binary responses — the parameterization is unchanged, and binary results are identical to before. Second, the new BlumeCapel() model adds a per-node quadratic potential δ that controls the tendency toward extreme versus neutral responses:

# Data with three ordered responses (-1, 0, 1):
mod_ising <- Ising(Data, responses = c(-1L, 0L, 1L)) %>% runmodel

# Blume-Capel: adds a quadratic potential (delta) per node:
mod_bc <- BlumeCapel(Data, responses = c(-1L, 0L, 1L)) %>% runmodel

# The Blume-Capel model fits these data much better:
compare(ising = mod_ising, blume_capel = mod_bc)
#        model      AIC      BIC
#        ising  8656.65  8726.92
#  blume_capel  8407.63  8501.32

Because the exact maximum likelihood enumerates all length(responses)^nNode response patterns, multi-state models are limited to small networks; the new maxStates argument controls this limit. In 0.16 the partition function additionally uses log-sum-exp accumulation, so extreme parameter values no longer produce Inf/NaN fits.

RI-CLPM: random-intercept cross-lagged panel models 0.16

See the panel-data page for a full walkthrough

The random intercept cross-lagged panel model framework has been finalized. ri_clpm() accepts the same flexible data input as the dlvm1 family (wide format with a design matrix, or long format via idvar), builds an independence baseline so CFI/TLI are available, and supports standardization options. The new ri_clpm_search() runs a sequential stationarity search: it adds stationarity of intercepts, temporal effects, contemporaneous relations, and innovation variances one at a time — and finally the panel-VAR model — keeping each constraint only if it does not worsen fit:

# Data: wide-format panel data; design: variables x waves matrix of column names
mod <- ri_clpm(Data, vars = design) %>% runmodel

# Sequential stationarity search:
search <- ri_clpm_search(mod, criterion = "BIC")
#            model       constraint DF      AIC      BIC Chisq p_diff decision
#             base           (none) 21 18346.63 18637.43 24.78           start
#       intercepts       intercepts 30 18332.83 18585.71 28.99  0.897   accept
#         temporal         temporal 48 18315.12 18492.14 47.28  0.437   accept
#  contemporaneous  contemporaneous 54 18307.42 18459.15 51.58  0.636   accept
#       innovation       innovation 60 18300.87 18427.31 57.03  0.487   accept
#         panelVAR wave1_endogenous 66 18305.19 18406.34 73.34  0.012   accept
#
# Selected model: panelVAR

In 0.16 the stationarity constraints were further hardened: they now select parameters structurally per wave, so pruned or partially fixed models constrain correctly instead of silently corrupting the constraint set.

Multi-group tooling 0.16

Several smaller releases improved multi-group workflows:

• equalityScoreTest() 0.16: a score (Lagrange multiplier) test for cross-group equality constraints that reproduces lavaan::lavTestScore to numerical precision. MIs(type = "free") now also reports joint multivariate score statistics.
• partialprune() 0.16: better specificity at low sample sizes (release_model = "pruned"), joint score tests for ranking candidate releases, and large speedups — an N = 15, two-group Ising partialprune() dropped from ~700s to ~43s.
• allequal() 0.16: constrain all free elements of a matrix to a single shared parameter within each group — e.g., all thresholds equal, or an equal-edge-weight network.
• Saturated-model fix 0.16: the internal saturated reference model no longer inherits cross-group equality from equal = ..., fixing inflated chi-square p-values and negative degrees of freedom.
• Multi-group ordinal data fix 0.16: per-group polychoric correlations are now computed from each group's own rows (previously all groups silently received pooled statistics). Robust sandwich standard errors for ULS/DWLS estimators followed in 0.16, matching lavaan's se = "robust".

Penalized (regularized) estimation 0.16

Experimental — see the tutorial for details

The 0.16 release introduced penalized maximum likelihood (estimator = "PML") and penalized FIML ("PFIML") with LASSO, ridge, and elastic-net penalties. The penalty parameter is selected automatically via an EBIC-based grid search with warm starts (find_penalized_lambda()), and refit() re-estimates the selected sparse structure without the penalty for valid inference:

library("psychonetrics")
library("dplyr")
library("psych")
data(bfi)

# LASSO-penalized GGM with automatic EBIC lambda selection:
mod_pml <- ggm(na.omit(bfi[, 1:10]), estimator = "PML") %>% runmodel
mod_pml@optim$penalty$lambda     # selected penalty: 0.088 (15 of 45 edges set to zero)

# Refit the selected structure without penalty for valid inference:
mod_refit <- mod_pml %>% refit

Recent fixes: the elastic-net penalty (0 < penalty_alpha < 1) no longer applies the ridge component twice (0.16; LASSO results were unaffected), and the default beta_min thresholding now matches the documentation.

Multi-level & meta-analytic frameworks matured 0.16

See the new multi-level and meta-analysis pages

The multi-level latent variable model ml_lvm() — a two-level random-intercept model decomposing within-cluster and between-cluster structures — received important correctness fixes in this series: the latent intercepts are now properly identified (0.16; standard errors and degrees of freedom were affected, estimates were not), the analytic Jacobian includes the previously omitted mean derivative with respect to the between-level regressions (0.16), and variable-name matching and incomplete-design handling were hardened (0.16).

# Two-level CFA: 60 clusters, 2 factors, 6 indicators
mod_ml <- ml_lvm(Data, lambda = Lambda, clusters = "cluster") %>% runmodel
mod_ml %>% fit    # verified: X2(20) = 15.92, p = 0.72, CFI = 1.00

# Network variants: ml_lnm() (between-cluster latent network),
# ml_rnm(), ml_lrnm()

On the meta-analysis side, the 0.14.20–0.14.27 releases introduced meta_lvm() (single-stage meta-analytic CFA/SEM extending the MAGNA framework) and meta_var1() / meta_gvar() (meta-analytic temporal and contemporaneous networks pooled across studies), both with fast C++ pipelines.