Rare items and runaway estimates in Paper × Items evidence synthesis

Advanced evidence-synthesis case study. This is not the beginner route through gllvmTMB. It deliberately studies difficult binary-data failures and uses opt-in remedies, including the experimental MSPL estimator. Start with Get started and Choose a response family for ordinary modelling workflows.

Evidence synthesis teams increasingly work with a papers-by-items table: each row a published study, each column a binary indicator — did the paper report an effect size, share its data, preregister, use blinded scoring? The scientific questions are ordination questions: which practices travel together across papers, and does a methods gradient explain them? That is a job for a joint latent-variable model of binary responses — the same machinery as a joint species distribution model, with papers in place of sites and reporting items in place of species.

This article fits that model, then walks into the two distinct ways such fits break, on purpose, and shows the matched remedy for each:

1. Runaway loadings. With a strong latent structure and only a few dozen papers, plain maximum likelihood can report convergence while the latent loadings drift to absurd magnitudes. The matched remedy is an opt-in loading ridge (gllvmTMBcontrol(loading_ridge = 2)).
2. A separated rare item. When a rare item is recorded only by papers at one end of a measured gradient, its fixed-effect estimate has no finite maximum. The matched remedy is the opt-in, experimental MSPL estimator (estimator = "mspl"), after the separation is certified by screen_gllvmTMB().

Each remedy treats only its own disease: they do not substitute for one another — the middle of the article demonstrates that too — and they are never combined in one fit. Ordinary estimator = "ml" with no penalty remains the default throughout the package; nothing here changes any default.

A papers-by-items corpus

Sixty papers scored on eight binary items, with one measured covariate: a design_score summarising each paper’s methodological rigour. Two latent axes generate residual dependence between items — think of them as unmeasured “methods culture” gradients across labs. The loadings are deliberately strong, because that is the regime where the first failure mode lives.

library(gllvmTMB)

set.seed(202)
n_paper <- 60
items <- c("defines_urban", "reports_effect_size", "shares_data",
           "preregistered", "blinded_scoring", "reports_power",
           "maps_species", "escalation_plan")
design <- as.numeric(scale(seq_len(n_paper)))
z   <- matrix(rnorm(n_paper * 2), n_paper, 2)     # two latent methods axes
Lam <- matrix(rnorm(16, 0, 3), 8, 2)              # strong true loadings
alpha <- c(0.3, -0.2, -0.5, -0.8, 0.0, -0.4, 0.6, -2.6)
beta  <- c(0.4,  0.6,  0.5,  0.7, 0.2,  0.5, 0.1,  1.2)

eta <- sweep(z %*% t(Lam), 2, alpha, "+") + outer(design, beta)
Y <- matrix(rbinom(n_paper * 8, 1, pnorm(eta)), n_paper, 8,
            dimnames = list(NULL, items))
papers <- data.frame(paper = factor(seq_len(n_paper)),
                     design_score = design, Y)
round(colMeans(Y), 2)
#>       defines_urban reports_effect_size         shares_data       preregistered 
#>                0.43                0.53                0.42                0.30 
#>     blinded_scoring       reports_power        maps_species     escalation_plan 
#>                0.50                0.60                0.60                0.27

Every item is assessed for every paper, and each cell is a single 0/1 record — a zero means “assessed and absent”, not “not reported”. That completeness is a real requirement: missing paper-by-item cells are not supported on the routes this article uses.

The model is the wide-format call an evidence table naturally invites — probit link, two latent dimensions, loadings-only residual structure (unique = FALSE, because a free per-item variance is not identifiable from single 0/1 trials):

form <- traits(defines_urban, reports_effect_size, shares_data,
               preregistered, blinded_scoring, reports_power,
               maps_species, escalation_plan) ~ 1 + design_score +
  latent(1 | paper, d = 2, unique = FALSE)

Failure mode 1: convergence with runaway loadings

Fit it exactly as you would any other model:

fit_ml <- gllvmTMB(
  form, data = papers, unit = "paper",
  family = binomial("probit"), estimator = "ml",
  control = gllvmTMBcontrol(n_init = 1, se = FALSE, warn_runaway = FALSE)
)
fit_ml$opt$convergence
#> [1] 0
max(abs(fit_ml$report$Lambda))
#> [1] 78.86288

The optimizer reports 0 — it found a point it judges stationary. But look at the loadings. On the probit scale a loading is a latent standard deviation in link units: the true values in this simulation top out near 4, and probabilities saturate within a few units of zero. An estimated magnitude in the tens is not a strong effect; it is the likelihood drifting along a flat direction while the fitted probabilities for those cells pin at 0 and 1. Convergence codes do not catch this — the flat direction looks stationary.

The matched remedy: an opt-in loading ridge

gllvmTMBcontrol(loading_ridge = 2) adds a Gaussian penalty on the loadings — a MAP-style shrinkage toward zero with scale tau = 2 — and nothing else:

fit_ridge <- gllvmTMB(
  form, data = papers, unit = "paper",
  family = binomial("probit"), estimator = "ml",
  control = gllvmTMBcontrol(n_init = 1, se = FALSE, warn_runaway = FALSE,
                            loading_ridge = 2)
)
fit_ridge$opt$convergence
#> [1] 0
max(abs(fit_ridge$report$Lambda))
#> [1] 3.095273

The loadings return to the scale of the truth, and convergence is clean. This is not a lucky draw: in this package’s internal calibration on the binomial-probit grid, the same penalty took the runaway rate at a comparable design size from roughly one fit in two to almost none. Two honest boundaries travel with that record. First, the measured failure corner of the ridge itself needs the logit link and few items together — on logit with small item counts and large paper counts the penalised fit can still run away, so the protective record quoted here is the probit record. Second, shrinkage is not free: the penalty biases every loading toward zero a little to prevent a few from escaping to infinity. The ridge is an opt-in stabiliser with a measured track record in this regime, not a default and not a certificate.

The same penalty composes with the package’s opt-in adaptive quadrature engine (gllvmTMBcontrol(aghq = 3, aghq_ridge = 2) — aghq_ridge is the same penalty under its original name). AGHQ replaces the Laplace approximation of the latent integral with an adaptive quadrature; its measured accuracy gains are for binomial models at unit counts far larger than sixty papers, so at this corpus size it is the ridge, not the quadrature, doing the visible work.

Failure mode 2: a separated rare item

Now the second, entirely different pathology. Suppose the escalation-plan item had been recorded only by the handful of papers at the very top of the design gradient — a realistic pattern for a practice that only the most rigorous papers adopt:

papers2 <- papers
papers2$escalation_plan <- as.integer(design > quantile(design, 0.93))
mean(papers2$escalation_plan)
#> [1] 0.08333333

A straight line on design_score can now split every presence from every absence. That is complete separation, and it is a property of the fixed design, not of the latent structure: the maximum-likelihood slope for this item is infinite, no matter how the latent side is handled.

First certify it

screen_gllvmTMB() runs the formal separation test item by item (it uses the optional detectseparation package):

scr <- screen_gllvmTMB(
  traits(defines_urban, reports_effect_size, shares_data,
         preregistered, blinded_scoring, reports_power,
         maps_species, escalation_plan) ~ 1 + design_score,
  data = papers2, unit = "paper",
  family = binomial("probit"),
  control = screen_control(separation = "fixed")
)
tab <- screen_table(scr, "separation")
tab[, c("traits", "status", "severity", "separated", "infinite_terms")]

#> Install the optional detectseparation package to run this certificate.

Seven items pass with overlap; the escalation item is flagged complete, and the infinite_terms column names the exact model terms whose estimates are infinite. The screen only diagnoses; it never changes the fit.

What each route does with it

Fit the separated corpus three ways and read the escalation item’s slope:

routes <- list(
  ml    = gllvmTMBcontrol(n_init = 1, se = FALSE, warn_runaway = FALSE),
  ridge = gllvmTMBcontrol(n_init = 1, se = FALSE, warn_runaway = FALSE,
                          loading_ridge = 2)
)
slopes <- sapply(names(routes), function(r) {
  f <- suppressWarnings(gllvmTMB(form, data = papers2, unit = "paper",
    family = binomial("probit"), estimator = "ml", control = routes[[r]]))
  co <- coef(f)
  unname(co[grep("escalation_plan:design_score", names(co))])
})
fit_mspl <- suppressWarnings(gllvmTMB(
  form, data = papers2, unit = "paper",
  family = binomial("probit"), estimator = "mspl",
  control = gllvmTMBcontrol(n_init = 1, se = FALSE)
))
co <- coef(fit_mspl)
slopes <- c(slopes,
            mspl = unname(co[grep("escalation_plan:design_score", names(co))]))
round(slopes, 1)
#>    ml ridge  mspl 
#>  53.8 276.2  16.3

The pattern is the whole lesson. Plain ML runs away, as the certificate predicted. The loading ridge does not help — it penalises loadings, and this pathology lives in the fixed effects; the two failure modes are genuinely different diseases. MSPL returns a finite slope, because its penalty — a Jeffreys-prior term on the fixed design, following Sterzinger & Kosmidis (2023) — targets exactly this boundary. (Supplying loading_ridge together with estimator = "mspl" is refused with an error: that hybrid has no supporting theory.)

The finite MSPL point is a penalised estimate, and its fences are strict by design: no logLik(), no AIC/BIC or likelihood-ratio comparison against the ML fit, and no standard errors or intervals while repeated-sampling calibration remains in progress — each of those fails with a typed error rather than returning an uncalibrated number. What you get is a finite, reportable description of the item’s association with the gradient, instead of a silently deleted item or an unusable coefficient in the hundreds.

Which remedy, when

Symptom	Diagnosis	Matched remedy
convergence = 0 but loadings at absurd link-scale magnitudes	Flat likelihood direction in the latent structure	gllvmTMBcontrol(loading_ridge = 2) — opt-in shrinkage, probit-grid track record
screen_gllvmTMB() flags an item complete / quasi_complete	Fixed-design separation; the ML slope is infinite	estimator = "mspl" — opt-in, experimental penalised point estimate with explicit fences
Both symptoms in the same fit	Two diseases at once — but still one fit	Still one outer penalty: treat the separation with estimator = "mspl". The ridge+MSPL combination is refused, so a co-occurring loading runaway has no packaged remedy — simplify the latent structure or collect more information

Report, in either case: the certificate table, which route was fitted and why, the penalty scale (summary(fit)$estimation stores c_n for MSPL), and the convergence and gradient diagnostics — not just the finite number.

What each remedy actually penalises

Neither remedy is a general-purpose fix for a rare or unusual binary trait, but the difference between them is not “loadings versus separation” — MSPL is a family of soft penalty atoms, and more than one of them can be active at once. For an ordinary-structure fit it adds:

Atom	Acts on	The pathology it addresses
Jeffreys	fixed effects, including intercepts	separation (a covariate slope, or a completely separated intercept-only item)
V_loading (radial)	the loadings Lambda	runaway latent loadings — the same target as the loading ridge
Hirose, sum_j S_jj / psi_j	the unique variances	a variance collapsing to the boundary (a Heywood case)
V_covariance	spatial range and scale	range or amplitude drifting to an unidentified extreme

Two practical consequences follow. First, the loading ridge and MSPL overlap: both act on Lambda, so they are alternatives rather than complements, and the package refuses to stack them. Second, an intercept-only model is not outside MSPL’s reach — its loading atom applies with or without a covariate, and a completely separated item gives the Jeffreys atom a target through the intercept.

The remaining difference is how the penalty is scaled, and it decides which tool answers which question. MSPL’s scale is fixed by the data — c_n is on the order of sqrt(p / N) and vanishes as the sample grows, which makes it a finite-sample correction that steps out of the way asymptotically. The ridge’s tau is a number you choose, and it never stops acting.

So the ridge is the weaker estimator and the better instrument: because you can weaken it deliberately, you can watch what an estimate does as the penalty recedes (see ridge_path()). MSPL has no such knob, so it cannot produce that evidence. Conversely, if a loading is not identified by the data at all — it keeps moving toward the boundary as the penalty weakens — then a penalty designed to vanish with information will not rescue it; it will return a finite, reasonable-looking number for a quantity the data cannot support. Diagnose first, and treat a rare or skewed intercept-only item as a marginal-support question (Screen binary traits before fitting: prevalence, minority count) before reaching for either penalty.

In this package MSPL remains experimental and opt-in, and its interval inference is not calibrated — that, not any absence of a target, is the reason not to reach for it in an applied analysis today.

Neither remedy, in turn, makes an exactly deterministic response block a valid set of traits. If several binary columns sum to exactly 1 on every row — a review-type, geographic-scope, or temporal-scope dummy set is a common source — that is one categorical variable recorded across several columns, not several independent binary traits. Loading shrinkage does not change that structure, and neither does a finite MSPL slope: the fix is to recode the block as multinomial() or to fix a deliberate reference level, not to penalise the fit around it. screen_gllvmTMB()’s response-dependency check (see Screen binary traits before fitting) flags this before either remedy is reached.

See also

• Joint species distribution model — the same latent machinery on sites × species occurrence data.
• Screen binary traits before fitting — the wider pre-fit screening workflow around screen_gllvmTMB().
• Can I trust this fit? — the general diagnostics this article’s convergence caveats point into.

References

• Kosmidis I, Firth D (2021). Jeffreys-prior penalty, finiteness and shrinkage in binomial-response generalized linear models. Biometrika 108:71–82. https://doi.org/10.1093/biomet/asaa052
• Sterzinger P, Kosmidis I (2023). Maximum softly-penalized likelihood for mixed effects logistic regression. Statistics and Computing 33:53. https://doi.org/10.1007/s11222-023-10217-3
• Sterzinger P, Kosmidis I, Moustaki I (2026). Maximum softly penalized likelihood in factor analysis. Psychometrika 91:494–507. https://doi.org/10.1017/psy.2026.10092 (adjacent theory for the factor extension; not proof for this implementation)
• Warton DI, Blanchet FG, O’Hara RB, Ovaskainen O, Taskinen S, Walker SC, Hui FKC (2015). So many variables: joint modeling in community ecology. Trends in Ecology & Evolution 30:766–779.
• van der Veen B, Hui FKC, Hovstad KA, O’Hara RB (2023). Concurrent ordination: simultaneous unconstrained and constrained latent variable modelling. Methods in Ecology and Evolution 14:683–695. https://doi.org/10.1111/2041-210X.14035