What's New

Multi-level VAR & GVAR for ESM data: ml_var() & ml_gvar() 0.16.7

July 2026 — experimental; see the multi-level VAR page

A new model family estimates the multi-level (random intercept) VAR for intensive longitudinal / experience-sampling data: ml_gvar() yields the three multi-level GVAR networks — the temporal network (beta / PDC), the contemporaneous network (omega_zeta_within) and the between-persons network (omega_zeta_between) — with confirmatory model tests, model search, and multi-group support. The input mirrors the mlVAR package (idvar, dayvar, beepvar).

mod <- ml_gvar(Data, vars = vars, idvar = "id", dayvar = "day", beepvar = "beep") %>%
  runmodel

mod %>% summary                              # incl. AIC / BIC
getmatrix(mod, "PDC")                        # temporal network
getmatrix(mod, "omega_zeta_within")          # contemporaneous network
getmatrix(mod, "omega_zeta_between")         # between-persons network

Two estimators are available. The two-level summary-statistics pseudo-ML (estimator = "ML") compresses each subject's series once into two-level sufficient statistics, so its cost is independent of the number of time points — the right choice for long ESM series. estimator = "FIML" instead routes to the panelvar framework for exact full-information ML on short (panel-like) occasion grids. The default estimator = "auto" picks between them at 10 occasions per subject and prints its choice. Estimates closely track both panelvar() and mlVAR() on the same data; see the statistical caveats for the pseudo-likelihood fine print.

Naming cleanup, summary() method & meta-model guards 0.16.7

July 2026

The plain names are now the primary names of several model families; the historical *1 names (the “1” referred to the lag order) remain fully supported aliases, so all existing code keeps working:

• dlvm() = dlvm1(), tsdlvm() = tsdlvm1(), meta_var() = meta_var1(), and ml_var() / ml_gvar() = ml_var1() / ml_gvar1().
• var1() / gvar() are unchanged: exporting var() would mask stats::var.
• within / between are now the preferred argument names for the contemporaneous and between-person structures in the latent-free frameworks panelvar() and ml_var(); the *_latent names remain accepted everywhere.

In addition, psychonetrics models gained a summary() method (the print output extended with log-likelihood, number of parameters, AIC, BIC, RMSEA, CFI and TLI — previously summary() silently fell through to the useless S4 default), and meta_var() / meta_gvar() now fail early with an informative error when more than 11 variables are used (the meta-level outcome dimension grows as p(p+1)/2 + p², which exceeds internal limits beyond that) and warn when the number of studies is smaller than the number of meta-level outcomes. For subjects of a single intensive longitudinal study, ml_var() / ml_gvar() is the recommended tool instead.

Multi-level variance-covariance & GGM: ml_varcov(), ml_ggm() & ml_corr() 0.16

July 2026 — see the multi-level page

A new framework brings the varcov family to two-level (clustered) data. ml_varcov() models the within-cluster and between-cluster variance–covariance matrices directly — with no latent layer — supporting the same two estimators as ml_lvm(): the fast sufficient-statistics two-level ML estimator (cost independent of cluster size) and wide-format FIML for missing data. Each level can independently take any of the five varcov parameterizations (cov, chol, prec, ggm, cor) through the within and between arguments. The headline wrapper is ml_ggm(), the multi-level Gaussian graphical model, which estimates a within- and a between-cluster partial-correlation network simultaneously; ml_corr() gives the correlation-matrix version.

# Multi-level GGM: within- and between-cluster networks from clustered data
mod <- ml_ggm(Data, vars = paste0("V", 1:4), clusters = "cluster") %>% runmodel

getmatrix(mod, "omega_within")    # within-cluster partial-correlation network
getmatrix(mod, "omega_between")   # between-cluster partial-correlation network

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

The framework has its own native matrix names (sigma_within, omega_within, sigma_between, omega_between, …) and simpler analytic derivatives, verified against numeric differentiation and an independent C++ implementation. The multi-level GGM is also available in the mlVAR package's mlGGM(), which uses a different estimator — nodewise multilevel regression (one mixed-effects model per variable); psychonetrics::ml_ggm instead estimates the full two-level covariance structure by maximum likelihood and additionally provides model fit indices and confirmatory testing.

Sampling weights (pseudo-ML) for complex survey data 0.16.5

July 2026 — see the lavaan interoperability page

varcov() and lvm() (and their ggm/precision/… wrappers) gain a sampling_weights argument for pseudo-maximum-likelihood estimation with survey weights. Point the argument at a column of data: the moments are then computed as weighted means/covariances and the model is fit with robust (MLR) Huber–White sandwich standard errors and a Yuan–Bentler–Mplus scaled chi-square, exactly as lavaan's sampling.weights= does.

• Matches lavaan's sampling.weights= to numerical precision — estimates, standard errors, the scaled chi-square and the robust fit indices — both single-group and multi-group (weights normalized so their total over all groups equals the sample size, lavaan's sampling.weights.normalization = "total").
• Supplying sampling_weights automatically selects the internal ML estimator with the MLR robust configuration and stores the raw data for the sandwich; constant weights reduce exactly to an unweighted MLR fit.
• Requires complete-data continuous input; not (yet) supported for ordinal data, missing-data FIML, or likelihood = "wishart". References: Asparouhov (2005); Rosseel (2012).

Model the PDC (standardized temporal network) directly 0.16.4

July 2026

The temporal effects of all lag-1 frameworks can now be parameterized directly in terms of partial directed correlations (PDC) — the standardized temporal network usually interpreted in graphical VAR analyses: temporal = "PDC" in var1()/gvar(), panelvar()/panelgvar() and tsdlvm1(), and temporal_latent = "PDC" / temporal_residual = "PDC" in dlvm1(). The PDC parameter matrix encodes from = row, to = column (matching getmatrix(..., "PDC")), with elements bounded in (−1, 1).

• Fit-equivalent to the raw parameterization for any zero pattern (the map is smooth, zero- and sign-preserving); the fitted PDC parameters equal computePDC() of the raw fit exactly at the optimum.
• Inference on the interpreted scale: Wald standard errors and CIs (CIplot(), arrows drawn from → to), modification indices, prune()/stepup() defaults, and cross-group equality constraints (groupequal("PDC")) — e.g. testing whether the standardized temporal network is invariant across groups while innovation variances differ.
• Both engines (R and C++); gradients verified against numeric derivatives for all innovation types (cov, chol, prec, ggm, cor) in every family, and PDC standard errors verified against delta-method standard errors from the raw parameterization.

dlvm1(): residual temporal effects (beta_epsilon) 0.16.3

July 2026 — see the panel data page

dlvm1() gained a beta_epsilon argument: lag-1 temporal effects at the residual level, epsilon_t = beta_epsilon epsilon_(t-1) + u_t. The default "zero" reproduces the previous white-noise residuals exactly; beta_epsilon = "diag" gives each indicator its own residual autoregression, modeling item-specific carry-over (autocorrelated uniquenesses) with a geometrically decaying structure — the classic remedy for the bias that arises when repeated items share unique variance over time (Cole, Ciesla & Steiger, 2007).

• With a non-zero beta_epsilon, sigma_epsilon_within (any type: cov, chol, prec, ggm, cor) parameterizes the residual innovation covariance; the stationary residual covariance is available as getmatrix(model, "sigma_epsilon_within_stationary").
• Together with sigma_epsilon_between (the perfectly stable item component), each indicator gets a STARTS-like decomposition: stable trait + autoregressive part + latent dynamics.
• The white-noise model is nested (beta_epsilon = 0), so modification indices on the beta_epsilon elements diagnose residual autocorrelation directly, and the likelihood-ratio test against the default model is the natural test.
• At least three waves are required (a warning is issued with two); four or more are recommended.
• Implemented in both the R and C++ engines; gradients and implied structures verified against independent implementations and Monte Carlo simulation.

panelvar() / panelgvar(): dedicated, faster model framework 0.16.2

July 2026 — see the panel data page

panelvar() and panelgvar() are now implemented as their own model framework rather than as a dlvm1() model with dummy matrices (identity factor loadings and zero residual variances). The model is unchanged — a within-person stationary VAR(1) process plus a between-person (random intercept) structure — and results are identical to the previous route: same estimates, standard errors, degrees of freedom and fit measures. Estimation now skips all factor-loading algebra in both the R and C++ engines, making runmodel() roughly 1.3–2.6× faster, with the gain growing in the number of variables and waves.

• The model matrices keep their names (beta, sigma_zeta_within, sigma_zeta_between and typed variants), so getmatrix(), prune(), stepup(), transmod(), CIplot(), esa() and modelsearch() keep working unchanged.
• The mean structure is a single mu matrix of observed stationary means (previously the observed intercepts appeared as nu with latent means mu_eta fixed to zero).
• Variable names default to the design-matrix rownames instead of Eta_1, Eta_2, ....
• Analytic gradients verified against numeric ones for all within/between types (cov, chol, prec, ggm, cor), including multi-group models with equality constraints, FIML, and designs with missing waves.
• panellvgvar() (the latent panel GVAR) is unaffected and still uses the dlvm1() framework.

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

June 2026 — experimental — 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 — experimental — 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 — experimental — 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 — experimental — 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 — experimental — 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.