Version 0.7.1 is the candidate for the first CRAN release of gllvmTMB.
It fits multivariate
response models using stacked long data or a wide-data formula, with
independent, full, or latent trait covariance and selected animal,
phylogenetic, spatial, and kernel correlation sources. This bounded 0.7.1
source adds no new response family, likelihood, integration engine,
random-slope capability, iSDM route, or broad predict(newdata = ) claim.
Availability on CRAN begins only after CRAN accepts and publishes the package.
Validation depth varies across these routes. In complete three-trait Gaussian
data, native-Laplace point estimation with indep() and dep() is covered at
sample sizes n = 60 and 240. Native-Laplace Gaussian rank-1 latent
evidence is partial: n = 240 passed, n = 60 was held, three admitted
correlation-stress cells failed, and three boundary/Psi challenges were held.
Native-Laplace Poisson-log rank-1 latent evidence is partial in the clean-data
n = 100 and 300 setting; a rare-level cell showed fixed-effect bias.
Native-Laplace NB2-log rank-1 latent evidence is mixed: n = 300 passed, but
n = 100 and missing-response cells failed with catastrophic dispersion
errors, while the rare-level cell was held. Evidence for other non-Gaussian
fixed-effect or extractor routes does not establish latent model performance.
The available covariance keywords do not imply equal validation for every
family and combination. Broad interval calibration is not established; one
narrow exception is the documented 0.94 coverage floor for a two-sided
Gaussian total-variance profile interval under its stated conditions. The
existing random-slope and MSPL routes remain experimental and are not expanded
by this release.
Parametric bootstrap refits are limited to two workers.
bootstrap_Sigma() and the predictor-informed latent-effect bootstrap
accept n_cores = 1 or 2; larger values return an error.
Explicitly unused optional grouping slots now warn (#1190). Supplying
unit_obs or cluster is useful only when a covariance keyword consumes
that column. If no keyword does, gllvmTMB() now warns and tells the caller
to omit the slot or use its column in the intended keyword. Omitted and
explicit-NULL slots remain silent; the warning does not change the model,
likelihood, or parameterisation.
extract_Sigma_B() and extract_Sigma_W() remain soft-deprecated
compatibility wrappers (#1194). They stay exported and retain their
historical return names. Their help now gives the direct migration to
extract_Sigma(fit, level = "unit") or extract_Sigma(fit, level = "unit_obs"); ordinary package summaries and residual helpers do not emit
the wrapper warning.
Variational approximation remains an opt-in experimental route (#1189). Native Laplace remains the default. This release makes no VA calibration, standard-error, confidence-interval, low-prevalence, or MSPL claim.
Existing random-slope documentation and existing MSPL material are retained. They are not new 0.7.1 capability claims: MSPL remains opt-in experimental, and this release neither expands nor promotes either topic.
extract_Sigma_B() / extract_Sigma_W() now warn on call (#1194).
The wrappers stay exported and still return the historical Sigma_B /
R_B and Sigma_W / R_W names. Callers of the old spelling now see
lifecycle::deprecate_soft() pointing at extract_Sigma(fit, level = "unit") / level = "unit_obs". Runtime warning only; no unexport.
The augmented-LHS parser guard now honours a non-default trait =
column (#1188). .assert_no_augmented_lhs() decided whether a bare
covstruct keyword (latent, unique, indep, dep, scalar,
spatial_dep) was written in the supported per-trait-intercept form
0 + trait | g by comparing the LHS symbol against the string
literal "trait", never consulting the resolved trait = argument.
A user with trait = "variable" writing the CORRECT long-format spec
latent(0 + variable | unit, d = K) was falsely rejected as an
unsupported augmented LHS, and was left with no way to write per-trait
intercepts in long format under their own column name. The guard now
accepts either the resolved trait column name or the literal
"trait", and its abort message names the user's own trait column.
⚠️ This blocked an external user from writing a correct model and led
to a wrong published analysis (the fallback spec collapsed 29 per-item
intercepts to one shared intercept). On the reporter's direct refit,
the d = 2 LV1-loading/prevalence association fell from R² = 0.742 with
the shared-intercept specification to R² = 0.125 with per-item
intercepts. Thanks to @iwogross for the report.
gllvmTMB()'s unit_obs and cluster arguments now default to
NULL instead of the concrete strings "site_species" / "species".
Both slots are optional -- most callers have no reason to supply them --
but a default that reads as a required column name invites callers to
manufacture columns just to satisfy the signature. An external systematic
map user (paper x item design, no sites, no species) wrote
site_species = paste(study_ID, variable, sep = "_") purely because
unit_obs defaulted to "site_species", and set cluster = "variable";
neither column had any meaning in the design or any effect on the fit.
⚠️ API/signature change; fitting behaviour preserved. The visible
defaults are now NULL, while omitted and explicit-NULL calls still
resolve internally to "site_species" / "species". Existing calls
that relied on or explicitly supplied those historical names fit
bit-identically. Also fixed: the "Column %s not found in data" error,
which never said which argument wanted the column, now names the argument
(trait, unit, or unit_obs) and the value that was looked for.
mesh= with no spatial term is no longer silently ignored (#1165).
Mesh validation used to run only when a spatial_*() term was present,
so gllvmTMB(mesh = <raw fmesher mesh>) with an ordinary non-spatial
formula fitted identically to the no-mesh call -- a clean converged
fit, and no signal that the spatial mechanism never ran. ⚠️
Behaviour change: a supplied mesh and no spatial term now warns,
naming both possible mistakes (the formula is missing a
spatial_*() term, or mesh was left over from a term that was
removed). A raw fmesher/INLA mesh is then rejected with the existing
make_mesh() error -- the check that already ran when a spatial term
was present. A valid make_mesh() object still fits as a
non-spatial model after the warning. Fits with a spatial term, and
fits with mesh = NULL, are unchanged.
predict(newdata = ) now also re-adds ordinary (1 | group) random
intercepts (#1138). This corrects an earlier entry, which recorded the
re_int tier as unreachable because "its group mapping is not a top-level
field on the fit". That was wrong -- fit$re_int carries groups,
n_groups and offsets, which is exactly what is needed. The tier is now
reconstructed and held to the same acceptance test as the others:
predict(newdata = training rows) reproduces report$eta exactly,
including with several (1 | g) terms, which is the case that
exercises the end-to-end packing via re_int_offsets -- an off-by-one
there would give a plausible value for the first term and a silently wrong
one for the second. An unseen level in a grouping factor still aborts, as
it did before: for an arbitrary grouping factor there is no defined
fixed-effects fallback, so refusing beats guessing.
predict(newdata = ) no longer requires the response column (#1154).
A prediction grid has no response by construction -- that is the point of
predicting on one -- but the newdata path built its design matrix with
model.matrix() on the two-sided model formula, which constructs a
model.frame() first and so evaluates the left-hand side. Predicting on a
grid failed with object '<response>' not found, and the workaround (add
a dummy column) was undiscoverable from the error and left the reader
wondering whether the placeholder affected the prediction. The design is
now built from the right-hand side only. Predictions are unchanged where
the column was present.
predict() and fitted() now return the family/source column on
mixed-family fits, and the effort-free scale is documented (#1133). On an
isdm_sources() fit the est column mixes scales by design -- Poisson
expected counts beside cloglog detection probabilities -- and the in-sample
path returned no column saying which row was which. The newdata path
always did (it returns all of newdata), so the default call, and
fitted() which wraps it, were the ones missing the label. The fit's
family_var column is now carried through; est remains the last column
and its values are unchanged. Fits with no family_var column --
every single-family fit -- return exactly what they always did.
Separately, ?predict.gllvmTMB_multi now documents that type = "response" includes the row's offset (so it is an expected count at that
effort), and that the effort-free / relative-intensity scale is obtained
by zeroing the offset variable in newdata -- exact, because the offset is
re-evaluated against newdata, which is why no new type = is needed.
vignette("response-families") gains seven compact worked examples,
one per family-specific dispersion trap (#1082). Gamma (phi_gamma is
a shape, not a dispersion), student (sigma_student is a scale, not the
response SD), Beta (phi_beta is a precision; boundary 0/1 values are
not Beta data), multi-trial binomial (cbind() vs weights = n_trials,
and a now-fixed simulate() gotcha for pre-2026-08-17 saved draws),
ordinal_probit() (cutpoints depend on which categories are actually
observed, per trait), the truncated count families (phi_truncnb2 is
its own vector, separate from phi_nbinom2), and delta_gamma()
(phi_gamma_delta is a CV, the inverse-square of its near-identical
neighbor phi_gamma) each get a minimal runnable fit and the documented
route to reading its dispersion. The residuals(type = "randomized_quantile") diagnostics comment is also corrected: exact
residuals cover 13 families now, not the 4 it previously named. These
are worked examples of point estimation and fitted diagnostics, not an
evidence-tier promotion.
predict(newdata = ) now re-adds three more random-effect tiers, and
propto is no longer skipped for every row when the first row has an
unseen species (#1138, #1132). diag_species, rr_W and diag_W are
reconstructed on the newdata path, each pinned by the same acceptance
test the SPDE tier had to meet: predict(newdata = training rows) must
reproduce report$eta exactly. That standard matters here because
q_sp is indexed (trait, species) -- the transpose of p_phy -- so a
swapped index would still yield a plausible non-zero contribution.
The rr_W/diag_W pair is keyed on the unit-observation column and is
added only when newdata carries it; otherwise the existing warning names
them as omitted. Also fixed: the propto re-add was gated on
!is.na(sp_id[1]) -- row one only -- so a single unseen species in the
first row silently dropped the tier for the whole frame; each row is now
guarded on its own.
Still deliberately not re-added, and still warned about: equalto (indexed
by observation, so it has no meaning for new rows), and diag_cluster2 /
the *_slope and phylo-diagonal blocks, whose reshape conventions are not
established (see getREsd()).
Boundary screening now sees two spatial/kernel Psi companions it
previously missed, and check_gllvmTMB()'s spatial psi row no longer
returns zero rows (#1119). ⚠️ Behaviour change: fits that previously
passed clean may now emit new boundary flags / WARN rows. Two real,
non-trivial quantities were absent from .gllvmTMB_boundary_flags()'s
screening list entirely: sd_spde_unique (the spatial *_unique() Psi
companion) and sd_kernel_diag (its multi-kernel analogue) -- a
collapsed spatial or kernel Psi went unflagged. Both are now screened;
sd_kernel_diag is filtered to only the tiers that actually carry a
fitted diag, since the current multi-kernel grammar cannot request one
yet (it is otherwise always an all-zero placeholder that would false-fire
on every 2+-named-kernel fit). Separately, check_gllvmTMB()'s own
spatial psi row used the dead literal "sd_spde" -- no such name is
ever REPORTed (the real name is sd_spde_unique) -- so it silently
produced zero rows (no PASS, no WARN) for every spatial fixture; it now
reports correctly. The dead literals "sd_phy" / "sd_spde" are also
removed from .gllvmTMB_boundary_flags()'s own list (they never matched
a REPORTed name there either; sd_phy_diag and sd_spde_b were already
screened under their real names).
Unnamed mixed-family family = list(...) lists no longer silently
swap which family fits which trait (#1120). .align_mixed_family_list()
left an unnamed family list unchanged, while the selector column's levels
for a character column were built via sort(unique(...)) -- alphabetical,
not the user's list order. A fit such as family = list(student(), gaussian()) against a plain character family column
rep(c("student", "gaussian"), each = n) silently fit the student data as
Gaussian and the Gaussian data as Student-t: zero warnings, a converged
fit, every downstream quantity wrong. This belongs to a recurring class in
this package -- the code chose a fallback and then reported success --
alongside fitted()'s and deviance()'s silent NULL (#1114, #1118);
the reassuring-message variants (a converged fit, a clean summary) are the
dangerous form, because nothing about the output signals that a guess was
made.
This is a behaviour change, but not a blanket one: an unnamed list is
now validated by computing BOTH the order the list was written in and the
order implied by each family object's own name (e.g. a "student" level
naturally means student()). When they agree, nothing was ever
ambiguous and the fit proceeds exactly as before -- most existing code,
including every other mixed-family test in this package's own suite,
already writes list order to match level order and is unaffected. When
they disagree, that disagreement IS the #1120 defect firing, and the fit
now errors (class gllvmTMB_mixed_family_unnamed_ambiguous) showing both
readings explicitly and telling the user to name the list
(list(student = student(), gaussian = gaussian())) or supply an
explicit factor. Selector columns with no name evidence at all (e.g.
arbitrary labels like "count"/"binary") still use list order, but the
resolved pairing is now reported once so it is auditable rather than
assumed. A level whose text matches more than one family object's name,
or a list where some levels have name evidence and others don't, is
refused as ambiguous rather than partially guessed. The same defect
pattern existed a second time, independently, in the multinomial
mixed-family trait-identification path (expand_multinomial_response(),
R/gllvmTMB.R, reachable before .align_mixed_family_list() runs) and is
fixed the same way.
predict(newdata = ) no longer silently drops the spatial field,
ignores re_form, or applies the wrong arm's inverse link (#1132).
Three defects in one function, all measured in Design 126 §3:
newdata branch re-added only the rr_B, diag_B and propto
tiers out of ~37, so a spatial fit's SPDE field was absent from every
newdata prediction -- at training locations too -- while the branch
reported that random effects had been added. Measured on a converged
spatial fit: dropped-piece sd 0.381 against a linear-predictor sd of
0.949, and reproduced independently on a plain gaussian() spatial fit
(sd 0.516 vs 0.786). The SPDE contribution is now re-added by rebuilding
the mesh projection with fmesher::fm_basis(), and any remaining active
tier this path cannot reconstruct now raises a warning naming it
rather than vanishing. newdata = NULL was always correct and is
unchanged.re_form was read only on the newdata path, and there only as the
literal ~ 0. On the package's default calling convention
predict(fit, re_form = ~ 0) therefore returned the full conditional
predictor, contradicting the documentation; NA -- a documented form --
and numeric 0 were ignored on both paths. All three forms are now
honoured on both paths, and an unrecognised value warns instead of
silently including the random effects. This changes the numbers
returned by existing re_form = ~ 0 calls, including
fitted(object, re_form = ~ 0), which forwards to the same path.type = "response" on newdata reduced the per-row family/link ids to
a per-trait modal id, which cannot represent an isdm_sources() fit
(the family varies by source within trait). Detection-arm rows were
returned through the count arm's inverse link, giving "probabilities" in
[0.253, 2.32] -- above 1, silently. The ids are now recovered per row
from the fit's family_var column, keeping each (family, link) pair
together. Single-family and by-trait fits are unaffected.lv_B score mean, the
rr_B re-add used the z_B innovation alone and dropped the
predictor-informed mean X_lv_B alpha_lv_B; it now uses the reported
U_B_total.BEHAVIOUR CHANGE: gllvmTMB_conditional_residual_saturated now also
warns for Gamma, Beta, and student-t fits (#1083). A diagonal random
effect indexed at the observed (unit x trait) resolution
(latent(0 + trait | unit, d = 1) with the default unique = TRUE)
perfectly interpolates every observation, driving the family's
dispersion/scale parameter toward a degenerate confound: sigma_eps -> 0
for gaussian/lognormal, Gamma's shape -> Inf, Beta's precision -> Inf,
student's sigma -> 0. residuals(fit, type = "randomized_quantile")
already warned about this for gaussian and (as of the previous entry
below) lognormal fits, but Gamma/Beta/student fits with the identical
structure were silent. A 15-seed sweep found Gamma's phi_gamma running
away past 1e6 (true 6) in 9/15 seeds and student's sigma_student
collapsing below 0.1 (true 0.4) in 6/15 seeds under this structure;
Poisson, swept the same way, showed no such collapse (it has no
continuous dispersion parameter to degenerate, so it is deliberately
not included, along with binomial, NB1/NB2, tweedie, and
beta-binomial -- all discrete-pmf families bounded by probability 1
rather than continuous densities that can diverge to infinity). Models
that previously fit silently under this structure will now warn; the fit
itself is unchanged, only the residuals-time diagnostic.
gllvmTMB_conditional_residual_saturated now also warns for lognormal
fits (#1083). The warning was gated to family_id == 0L (gaussian)
even though gaussian and lognormal share one literal sigma_eps and are
auto-suppressed identically at fit time
(any_sigma_eps <- any(family_id_vec %in% c(0L, 3L)),
R/fit-multi.R:5177) under a per-row diagonal random effect. The gate
now matches that fit-time decision.
deviance() no longer returns a silent NULL on a gllvmTMB_multi
fit (#1118). With no method registered, deviance() fell through to
stats:::deviance.default, which reaches for object$deviance and
returns NULL without error -- the same shape as the fitted() defect
fixed in #1114. deviance.gllvmTMB_multi() now returns
-2 * logLik(object), delegating through logLik() so it inherits that
method's existing MAP-point disclosure warning on ridged
(aghq_ridge-penalised) fits rather than re-deriving it.
Mixed-family fits with a dispersion family now have a valid Hessian and
SEs (#1117). Per-trait dispersion/shape parameter vectors
(log_phi_nbinom2, log_phi_nbinom1, log_phi_gamma,
log_phi_tweedie + logit_p_tweedie, log_phi_beta,
log_phi_betabinom, log_sigma_student + log_df_student,
log_phi_truncnb2, log_sigma_lognormal_delta, log_phi_gamma_delta)
were previously gated WHOLE-VECTOR: mapped off only when the family was
absent from every trait. In a mixed-family fit where the family was
present on some traits but not others, the other traits' entries were
free parameters the likelihood never reads, a mechanically singular
Hessian (pdHess = FALSE), no valid sdreport/Wald/profile intervals,
and phantom entries enumerated by profile_targets()/confint(). Each
vector is now pinned per trait to the traits that actually use its
family. Single-family fits are unaffected (byte-identical).
Convergence diagnostics on ridged fits now judge the objective the fit
actually optimised (#1092). aghq_ridge applies a loading penalty in R,
outside the TMB objective, so fit_health$max_gradient previously reported
the gradient of the unpenalised likelihood at the penalised optimum, a
number that equals the missing ridge term (|Lambda|/tau^2) rather than
~0, making perfectly converged ridged fits read as unconverged. The
reported gradient, stationary_by_gradient, converged, sanity_multi()
and the AGHQ stop reason now evaluate the penalised objective; the new
fit_health$gradient_is_penalised field discloses when this applies, and
the raw likelihood gradient remains available via fit$tmb_obj$gr().
Unridged fits are unchanged. Relatedly, the ridge's reach is now pinned to
the ordinary loading block theta_rr_B only: an accidental merge had
briefly made it also penalise spatial latent loadings
(theta_rr_spde_lv), contradicting the documented exemption for spatial
terms; a regression test now proves spatial loadings carry no silent
penalty.
Multinomial fits now warn at fit time when their contrast structure is
degenerate. The screen added this release (collapsed contrast variance,
rail-correlated contrasts, or a collapsed spatial range) is now surfaced
automatically, under the existing gllvmTMBcontrol(warn_runaway = )
switch and with its own once-per-session slot, so it cannot suppress or be
suppressed by the binomial runaway warning. Set warn_runaway = FALSE to
silence both; check_gllvmTMB() and gllvmTMB_diagnose() are unaffected.
Ordinal fits deliberately emit no such warning; that family's arms ship
disarmed (see the calibration note below), so the row reports statistics
without a verdict.
This development release adds an opt-in separation-screening and LA-MSPL lane for complete single-trial Bernoulli GLLVMs. Ordinary ML remains unchanged and is still the default.
check_gllvmTMB()'s loading_absolute_thresh default raised from 6 to
8 (issue #1098). Some binomial fits that previously reported WARN on
the binomial_prevalence_loading row now report PASS; nothing that
reported PASS before starts reporting WARN. The earlier default was
calibrated on a pool whose true loading scale never reached the regime
where this arm misfires. A pool built to cross sigma_lambda in c(0.7, 3.0) (928 healthy / 272 degenerate binomial-probit fits; 3.0 chosen
to hit issue #847's aghq_ridge ridge-failure regime, not argued for
realism) measured a 25% false-positive rate at the old threshold, all
of it attributable to this one arm; raising the threshold to 8 lowers
that rate to 15.52% while sensitivity on the same pool's degenerate fits
falls only from 100.00% to 99.63% (one additional missed fit out of
272). This is an improvement, not a fix: the arm remains measurably
regime-dependent: false-positive rate 3.85% at a mild true loading
scale versus 49.08% at the scale aghq_ridge = 2 is already known to
struggle at (issue #847), so a fixed constant cannot be correct across
every loading scale a fit may have. This is effect-size dependence, not
the response-scale dependence #851/#855 otherwise describes: probit
fixes the residual variance at 1, so there is no free response scale
here to rescale against, and the class's usual per-fit device does not
obviously transfer.
aghq_ridge = 2 reduces the problem (46.0% -> 13.5% false positives at
that larger scale) but does not remove it.screen_gllvmTMB() now catches exact response-side dependencies
among three or more traits (issue Ayumi-495 /
urbanisation_map#23). The pairwise duplicate/complement check only
sees two-way exact relations; a response_affine_rank row (in $design)
and a new $response_dependencies table now also flag one-hot/simplex
blocks, most commonly a review-type, geographic-scope, or
temporal-scope variable recorded as several dummy columns that sum to
exactly 1 on every row, via an affine-rank check on the augmented
response matrix. Certificates are recovered on a best-effort basis (not
exhaustive minimal-subset discovery); the new known_groups argument
lets you declare a set of trait names as a one-hot block or a
nesting/containment chain for an exact, deterministic check instead.
See the new "Exact response-side dependencies" section of
vignette("pre-fit-response-screening", package = "gllvmTMB").
New ridge_path(): a loading-ridge sensitivity sweep (issue
Ayumi-495 / urbanisation_map#23). Refits a model across a grid of
gllvmTMBcontrol(loading_ridge = tau) scales and reports, per trait,
how the largest loading and its communality move as the penalty
weakens. print() classifies each trait as "interior (stabilises as
penalty weakens)" or "penalty-determined (moves toward boundary)" using
a documented slope rule on the last two finite-tau grid points. This
is sensitivity evidence, not an identification certificate; tau = Inf
(plain ML) can legitimately fail to converge for a pathological trait,
and that failure is itself part of the diagnostic.
fitted() now works on the default gllvmTMB_multi fit (#25).
fitted was only registered as an S3 method for the gllvmTMB_julia and
gllvmTMB_va engine classes, so fitted(fit) on an ordinary fit silently
returned NULL rather than erroring or dispatching. fitted.gllvmTMB_multi()
is a thin wrapper over predict(object, newdata = NULL, type = ), returning
the same long data frame (default type = "response", matching the
fitted() convention; predict()'s own default remains "link").
Exact randomized-quantile residuals now cover nine more families.
residuals(fit, type = "randomized_quantile") and
predictive_check(fit, type = "rq_qq") compute exact family-CDF
residuals for binomial (logit/probit/cloglog), lognormal, Gamma, Beta,
betabinomial, Student-t, zero-truncated Poisson, zero-truncated NB2, and
ordinal-probit rows, in addition to the existing Gaussian, Poisson, NB1,
and NB2 support. Rows for tweedie, the delta/hurdle families
(delta_lognormal, delta_gamma), and multinomial remain
status = "unsupported_family": tweedie has no closed-form CDF without a
new dependency, the delta/hurdle families need a design decision for
splitting the point mass at zero, and multinomial's categories are
unordered so a randomized-quantile residual is undefined for it. See
?residuals.gllvmTMB_multi for the full scope statement.
A check_gllvmTMB() row for multinomial K-1 contrast degeneracy.
multinomial() (fid 16) fits as K-1 baseline-category contrast
pseudo-traits per response; check_gllvmTMB() now reports a
multinomial_contrast_degeneracy row that screens those contrasts
against each other for three failure modes generic loading diagnostics
cannot see: one contrast's loading energy collapsing to ~0 (absolutely,
or relative to its sibling contrasts), two contrasts of the same
response loading almost perfectly on the same axis (evaluated only where
the tier's rank is 2 or more, since rank-1 tiers reach that correlation
on every healthy fit by construction), and a spatial field's practical
range collapsing relative to the coordinate domain. This is a check-row
addition only -- no fit-time warning, no new export, no behaviour
change to any existing fit.
Three screens ship armed at calibrated defaults: multinomial_collapse_floor = 1e-10 (contrast variance collapse), multinomial_rail_thresh = 0.99 (two contrasts railed onto one axis), and
multinomial_range_collapse_thresh = 0.02 (spatial range collapse).
multinomial_collapse_rel_thresh stays disarmed (Inf) -- untested.
Measured against pre-registered labels (128 fits, 122 converged with a PD
Hessian): the variance-collapse screen detected 6/7 labeled collapses
plus 7/7 on a later, entirely out-of-sample cell; the contrast-rail
screen detected 8/8 labeled rails, plus four individually
railed fits hidden inside a cell whose aggregate gate had passed
(refitting confirms |rho| = 1.00000, against controls at 0.49 and
-0.15); the spatial-range screen detected 3/3. Zero false positives on 40 informative healthy fits, and
the denominator is stated honestly: fits with no loading tier emit no row
at all and cannot evidence specificity, so the rule-of-three bound is
about 7.5%, not a verified zero. The contrast-rail screen is deliberately silent on rank-1
tiers, where |rho| = 1 holds on every healthy fit by construction
(verified 0/20 out-of-sample).
The spatial-range screen was fixed before arming: it originally read the low-rank loading
matrix Lambda_spde, which the engine reports only for
spatial_latent()/spatial_dep() fits, so on spatial_indep() fits --
exactly the ones it was built for -- it produced no row at all. It now
branches on the engine route, taking the range from log_tau_spde on
the diagonal route: rows emitted 0/20 to 20/20, detection 0/3 to 3/3,
with 0/11 false positives on the same cell's healthy fits.
A check_gllvmTMB() row for ordinal-probit loading degeneracy (fixes
the coverage gap reported in #897, where a degenerate ordinal_probit()
fit had no detector at all -- 239/239 unflagged where the binomial screen
caught 272/272). A mechanism probe found no empirical support for
cutpoint-underflow saturation across 24 measured degenerate fits
(flat-row share exactly 0 throughout) -- that negative finding stands.
The probe's third measurement, originally read as positive evidence for
category-level quasi-complete separation (the same mechanism the
binomial row already screens for), was later shown, on a larger 315-fit
calibration, to fire on 86.3% of healthy fits and so does not
discriminate; category-level separation remains the residual hypothesis,
not a demonstrated one. The pathology's shape is a single trait's
loading column running away while sibling traits stay near truth.
check_gllvmTMB() now reports an
ordinal_liability_loading row with two arms modeled directly on the
binomial row: a trait's largest loading relative to the typical loading
among the other ordinal traits, and the largest loading on the link
(liability) scale, unit tiers only -- scale-free because the
probit-liability residual variance is exactly 1 by the Wright/
Falconer/Hadfield threshold convention. This is a check-row addition
only -- no fit-time warning, no new export, no behaviour change to
any existing fit.
Both arms ship DISARMED at Inf/Inf -- the row still computes and
reports its statistics, but neither threshold is armed by default.
Calibration ran 315 fits across four pre-registered design arms, scored
under the frozen pre-registration (sensitivity on the degenerate arm's
rel_frob > 10 fits; false positives measured across healthy +
transport + mixed combined), and found no threshold that clears
the pre-registered target of 90% sensitivity with zero false alarms --
and none can: the healthy pool reaches a loading magnitude of 216.9
while the degenerate arm starts at 13.5, so the classes are not
separable on this statistic at all. At binomial's own threshold of 6,
the absolute arm reaches 100% sensitivity but 39.2% false positives
and the relative arm 61.0% sensitivity at 28.6% false positives --
worse than the 25% false-positive rate #897 reports for binomial
itself, so borrowing binomial's threshold would have shipped a bigger
problem than the one the issue complains about. The false alarms
concentrate in designs mixing very different per-trait loading scales:
an absolute liability-scale threshold cannot transport across
heterogeneous trait scales, so a future screen needs a scale-invariant
statistic rather than a better constant. Honest limits: no evidence at
n = 1600 (that arm was dropped for run time). What the campaign
establishes: link saturation is refuted as the mechanism (solid);
category-level separation, the residual hypothesis, is NOT demonstrated
because the originally cited evidence does not discriminate; and
the threshold question is answered negatively with a stated path forward.
Neither categorical screen changes what fitting itself does: gllvmTMB()
warns exactly as before, and both rows appear only when you call
check_gllvmTMB() on the fit.
Liability-scale phylogenetic heritability for categorical families
(this corrects a real defect: on a phylogeny-only ordinal_probit() or
multinomial() fit, extract_phylo_signal() previously reported
H2 = 1.0 for every trait and every contrast. That number was the
species-level-latent proportion and was arithmetically correct as such,
but as a heritability it was silent nonsense -- the fixed liability
residual never entered the denominator, so a phylogenetic signal that
should read around 0.3-0.4 read as 1. If you have reported a categorical
H2 from this function, re-run it with link_residual = "auto".)
extract_phylo_signal() gains a link_residual argument. The default
("none") keeps the historical species-level-latent denominator
unchanged, so no existing non-categorical result moves;
link_residual = "auto" adds each trait's fixed
distribution-specific latent residual to the denominator, returning the
conventional liability-scale phylogenetic heritability of Mizuno et al.
(2025, J. Evol. Biol. 38:1699-1715, eq 4/18/19):
ordinal_probit() reports per-trait H2 = V_a / (V_a + 1) and
multinomial() reports per-contrast
H2 = V_a(k) / (V_a(k) + pi^2/3). Multinomial contrast heritabilities
are baseline-referenced (the softmax link residual couples contrasts
through the shared baseline category) and are reported one row per
contrast, never collapsed to a scalar. An advisory now fires when a
categorical fit is summarised with the default denominator, and
ci = TRUE with link_residual = "auto" refuses with a typed error
rather than returning uncalibrated intervals. Verified against an
MCMCglmm family = "ordinal" comparator on a shared phylogenetic
fixture.
A Species Distribution Models article collection. The pkgdown site
gains a dedicated navbar menu ordering the SDM material as a curriculum:
the joint species distribution model guide, a new presence-only opener,
and the three integrated-source articles. The new article, Joint
ecological intensity from opportunistic records
(vignettes/articles/gbif-joint-intensity.Rmd), fits a joint
relative-intensity model to GBIF-style records alone through the ordinary
gllvmTMB() call (Poisson counts, an effort offset, a named recording-bias
covariate, and one latent() factor). It estimates relative ecological
intensity only, no abundance, absolute occurrence, or detectability
claim, and states the design assumption its bias term rests on.
A fifth SDM article: repeated survey visits. What do repeated survey
visits add to an integrated model?
(vignettes/articles/integrated-repeated-visits.Rmd) shows visit-varying
conditions entering the survey branch on the complementary-log-log scale,
demonstrates that within-cell replication under varying conditions is what
separates a poor visit from low intensity, and states the boundary with
separate-detection occupancy models. The two-source article is renamed
Integrating opportunistic records with a designed survey to keep the two
subjects distinct.
A sixth SDM article: rare species and separation. When a rare
species breaks the JSDM (vignettes/articles/rare-species-jsdm.Rmd) is
the species-facing door to the separation workflow: a heat specialist
recorded only at the warmest sites, the screen_gllvmTMB() certificate,
what ordinary ML does at that boundary, the finite opt-in experimental
MSPL point with its fences, and a live demonstration that the loading
ridge does not repair separation.
A four-lens editorial audit of the SDM collection. Every article under
the Species Distribution Models menu was audited for readability, claim
discipline, statistical correctness, and figure integrity. Fixes include:
the two-source and design articles' displayed equations now carry the
per-species reporting level rho_s their code always fitted; the stale
"more than two sources: not available" claims are corrected to point at
isdm_sources(); health-check outcomes are stated after every
check_gllvmTMB() chunk; pd_hessian and Fisher-z are glossed at first
use; and figures gain alt text, equal-axis identity lines, and
non-overlapping labels.
LA-MSPL worked example (development site only). Rare items and runaway estimates in Paper × Items evidence synthesis is available on the development site. It is not included in the CRAN 0.7.1 package. The article grounds its example in an evidence-map-shaped corpus and demonstrates, on one dataset, that runaway loadings and fixed-design separation are different diseases with different matched remedies: the opt-in loading ridge for the first and opt-in MSPL for the second. It also shows that the ridge does not repair separation. The article states the boundaries for probit and logit ridge regimes, AGHQ's large-sample evidence, and MSPL's inference surface.
Complete-Bernoulli GLLVMs can opt in to LA-MSPL point estimation. Set
estimator = "mspl" to fit the experimental maximum softly penalised
Laplace objective for one ordinary latent(d = 1:2) block with no free Psi,
standalone spatial_indep(), or standalone spatial_latent(d = 1:2), with
a common logit, probit, or complementary-log-log link. estimator = "ml"
remains the unchanged default, and a separation warning never switches
estimators automatically. The admitted surface also requires a full-rank
resolved fixed-effect design and all-zero offsets. The implementation follows
the construction of Sterzinger and Kosmidis (2023,
doi:10.1007/s11222-023-10217-3) for its fixed-effect
Jeffreys component and extends it with rotation-invariant radial penalties
for estimated GLLVM loadings and reference-scale spatial coordinates.
At this development checkpoint, the compiled objective and point-estimation
API are implemented, while frozen ordinary and spatial simulation promotion
remains in progress. Likelihood
comparison, standard errors, confidence intervals, profiles, and hypothesis
tests fail closed. Missing responses, grouped binomial data, mixed families,
rank above two, extra covariance tiers, VA/AGHQ/Julia MSPL, and general
inference remain planned rather than advertised.
The loading MAP penalty now has the integration-neutral spelling
loading_ridge. It is an alias for aghq_ridge, not a new penalty, and
preserves unpenalised Laplace fits unless named explicitly. Loading ridge and
LA-MSPL are distinct estimators; supplying both is an error.
screen_gllvmTMB() can opt in to fixed-design separation certificates.
Use screen_control(separation = "fixed") to classify maximal
coefficient-connected Bernoulli fixed-effect blocks as overlap, complete
separation, or quasi-complete separation. The detector supports the logit,
probit, and complementary log-log links, finite known offsets, shared
coefficients, and structural-zero Xcoef_fixed constraints. It requires the
optional detectseparation package. The default remains
separation = "none" and does not load that package. This is a fixed-design
diagnostic only: it does not certify finite latent loadings or covariance
parameters, and it never selects a penalized estimator automatically.
screen_gllvmTMB()'s known_groups argument no longer misses a
partial-order nesting declaration or a certified one-hot block in its
unresolved-dependency count (reported by @Ayumi-495,
urbanisation_map#23). Two related defects in the known_groups
feature shipped in #1123: (1) a group of three or more traits whose
containment relations form a partial order rather than a single total
chain -- for example a broad realm indicator ("water") containing two
mutually incomparable narrower ones ("freshwater", "marine") --
previously fell through the chain-only check (which tested only the
declared order and its reverse) and silently PASSed
known_group_checked instead of FAILing known_nesting; the check now
tests every pairwise containment among the declared members and names
every relation it finds. (2) $response_dependencies' unresolved
affine-dependency count never subtracted the one-hot certificates a
user declared via known_groups, so a correctly certified dependency
could still be reported as "N further exact affine dependencies ...
could not be automatically resolved". The count is now the RANK of the
pooled certified null-vector span (automatic certificates + declared
one-hot blocks), not a subtraction of row counts, so declaring the same
dependency under two names -- or via two overlapping groups -- cannot
under-report a genuinely unresolved dependency. A nesting/containment
relation is an inequality, not an exact affine relation, and still
contributes nothing toward unresolved. Also fixed in the same pass: a
typo'd trait name in known_groups was silently accepted whenever the
affine/known-group check itself was infeasible for the data (no unit
column, non-Bernoulli rows, duplicate unit-trait rows, too few complete
units); trait names are now validated before that feasibility check, so
a typo always aborts.
screen_gllvmTMB()'s known_groups one-hot test was whole-group
only and is now a bounded exhaustive subset search (#1154, reported by
@Ayumi-495). known_groups = list(g = c("A", "B", "C", "D")), where
{A, B, C} is a genuine one-hot block and D is an unrelated trait
declared in the same group, previously reported known_group_checked
PASS -- the same "declaring a slightly-too-large group weakens the
verdict" shape as the nesting defect fixed just above, since the test
was all(abs(rowSums(Yg) - 1) < tol) over the whole declared group.
A declared group of k <= 12 members now gets an exhaustive search
over every subset of size 2 or more (at most 2^12 - 13 = 4083
row-sum checks, microseconds), and every MINIMAL one-hot subset found
is reported as its own known_one_hot_subset certificate -- A + B + C = 1 once, never also as a redundant larger superset -- and
contributes to the unresolved affine-dependency count exactly like a
whole-group known_one_hot certificate. k > 12 is never silently
skipped: the whole-group one-hot and pairwise nesting checks still
run, but an explicit known_group_subset_not_attempted row records
that the subset search itself was not attempted at that size, rather
than folding into a silent PASS. A subset that sums to 1 only because
every one of its members is constant is not certified (the same guard
the whole-group check already had). Nesting was not affected: its
pairwise scan already covers every pair within a declared group
regardless of unrelated members, so it had no equivalent blindness.
fit_health$boundary_flags no longer flags the auto-Psi skip block's
mapped-off sd_B placeholders as near_zero_sd_B (#25). The same
skip_psi_b_t pinning described in the near_zero_psi_unit fix below
(single-trial Bernoulli traits, and every multinomial contrast
pseudo-trait, pinned to sd_B = 1e-6 and mapped off) was also read raw
by .gllvmTMB_boundary_flags(), which had no diag_B_skip filter of its
own, so it unconditionally flagged the pinned placeholder. A default
auto-skip fit and its explicit latent(..., unique = FALSE) mirror could
therefore disagree on boundary_flags for an otherwise identical model.
Both readers now share one helper, .gllvmTMB_estimable_components(),
which drops the mapped-off entries before either screen runs; a genuine
collapse among the remaining, estimable traits is still caught. This is a
behaviour change: fits with a mapped-off Psi trait that previously
carried a spurious near_zero_sd_B flag no longer do.
simulate() now draws from each row's true family instead of silently
substituting a Gaussian-on-link-scale number for nine families, and
multi-trial binomial() responses are drawn at the correct trial count.
.draw_y_per_family() previously fell through to a single shared branch
for any family it did not explicitly recognise, so tweedie, Beta,
betabinomial, Student-t, truncated Poisson, truncated NB2,
delta_lognormal, delta_gamma, and ordinal-probit rows were redrawn as
plain Gaussian noise on the link scale, a plausible-looking but wrong
number, not an error. simulate(), predictive_check()'s simulation-based
plot types (stat_grouped, dens_overlay), and anything built on
.gllvmTMB_predictive_draws() (bootstrap/coverage helpers included) now
draw correctly for 16 of these: every family except tweedie, which has no
exact draw without a new dependency and now returns NA with a per-call
warning rather than the previous wrong-distribution substitute. Multinomial
rows draw a single grouped categorical outcome per observation rather than
an independent draw per contrast row. Binomial rows previously ignored
trial count and drew Bernoulli regardless of cbind(success, failure) /
weights = n_trials; they now draw rbinom() at the row's actual
n_trials.
check_gllvmTMB()'s near_zero_psi_unit screen no longer flags traits
the auto-Psi skip block deliberately pinned off. R/fit-multi.R's
skip_psi_b_t block maps a trait's between-unit Psi off (single-trial
Bernoulli, and every multinomial contrast pseudo-trait) by pinning
theta_diag_B at log(1e-6), but src/gllvmTMB.cpp still REPORTs
sd_B for every trait including the pinned ones, so the pinned 1e-6
entry always cleared both the absolute (psi_thresh = 1e-4) and relative
(psi_rel_thresh) collapse thresholds. This was a structural false
positive that predates the Design 123 arc above: check_gllvmTMB() WARNed
near_zero_psi_unit on every fit with a single-trial-Bernoulli or
multinomial trait sharing a latent() term with a free partner trait,
regardless of whether the free trait's Psi was healthy. The screen now
drops pinned entries (via tmb_data$diag_B_skip) before evaluating the
unit-level psi row; a genuine collapse among the remaining free traits is
still caught.
multinomial() structured-term admission is now fail-closed (Slice 0,
Design 108/123). Several deferred keywords previously desugared
(R/brms-sugar.R) onto the same internal engine flag as an admitted
keyword and silently reached an untested categorical path instead of
erroring: dep() at the unit tier, phylo_dep(), phylo_indep() /
phylo_unique() (standalone), animal_latent(), single-name
kernel_*(), and phylo_scalar() / animal_scalar(). Every one of these
now aborts, as the documentation always said they should. If you fitted
a multinomial() trait combined with any of the keywords above, that fit
ran on an unvalidated structured-term path and should be re-checked;
the currently admitted set is unchanged: fixed effects, an ordinary shared
latent(0 + trait | unit, d = k) ordination, and intercept-only
phylo_latent() (default unique = FALSE; it emits no Psi companion at
all). A mi() predictor term combined with a multinomial trait, previously
invisible to the admission scan due to definition order, now also aborts.
A follow-up adversarial review found the fence itself was not yet
load-bearing everywhere: augmented (intercept + slope) latent() /
phylo_latent() random regressions, phylo_latent(unique = TRUE) (a free
phylogenetic Psi, never admitted, but the fence's first pass wrongly
labelled it admitted), and meta_V() / equalto() (no established route
on a categorical-contrast pseudo-trait) are now all explicitly blocked,
with a shared classed condition (gllvmTMB_multinomial_structured_not_admitted)
on every path.
The multinomial() structured-term surface now admits a bounded set of among-category and grouping
structures, each gated on a signed, pre-registered recovery campaign
rather than construction alone. Every admission below is enforced by
the same fail-closed classifier (R/multinomial-fence.R); anything not
named below still aborts typed
(gllvmTMB_multinomial_structured_not_admitted). The full per-cell table
and supporting evidence are available in the
package repository.
Admitted: the phylogenetic/relatedness surface -- intercept-only
phylo_latent()/animal_latent()/single-name kernel_latent()
(loadings-only), phylo_dep()/animal_dep()/kernel_dep() (full
unstructured V, the IDENTICAL parameterisation as
phylo_latent(d = K - 1)), and phylo_indep()/animal_indep()/
kernel_indep() (diagonal V); the spatial (SPDE) surface --
spatial_latent(), spatial_indep(), and spatial_dep() (verified
identical to spatial_latent(d = n_traits)); a generic (1 | group)
random intercept (baseline-vs-rest semantics, sigma_re
reference-category-specific); and the non-phylogenetic cluster/
cluster2 diagonal tier, indep(0 + trait | g).
The honest evidence, not softened: on the phylogenetic surface, the
ONE-CATEGORICAL-DRAW-PER-SPECIES recovery gate FAILED for both the
loadings-only route (FAM-20C: rail rate 8/20, exceeding the 6/20
threshold, identically for animal_latent()/kernel_latent() by proven
engine identity to phylo_latent()) and the mode-axis route (FAM-20D:
phylo_dep() rails 8/20; phylo_indep()'s corrected diagonal-truth
rerun shows larger contrast variances recover fine, median ratio 0.78,
17/20 in band, but smaller ones collapse, median ratio 0.24, 9/20 --
7 of those 20 seeds collapse the smaller contrast variance to
numerical zero (≤1e-9), each with convergence = 0 AND a PD Hessian,
and no runtime detector currently flags it (R/diagnose.R's degeneracy
gate is family_id == 1-only, issue #897's class) -- and the planted-zero
check FAILS -- a full-V phylo_latent() refit on
diagonal-truth data rails to median |rho| = 1.0). A pre-registered
replication rescue PASSED: five categorical draws per species (n_sp = 300, n_rep = 5) recovers V with rail rate 4/20, median rho 0.680
among the 16 NON-RAILED seeds (0.696 among all 20 conv+PD seeds --
both bands, true 0.6), SD ratios 0.89/0.85 -- one categorical draw per
species does not identify V; five draws per species does. This rescue
transfers exactly to phylo_dep() (the identical parameterisation) but
has NOT been tested for the diagonal-V mode (phylo_indep()). The
spatial kappa/tau gate, by contrast, PASSED all three cells (median
practical-range ratios 1.75 / 1.12 / 1.75 among the 14 conv+PD seeds per
cell, band 0.33-3.0; rails 0, 3, 0 of those 14 against the frozen
threshold, restated in its pre-registered form: >6/20; 6/20 seeds were
non-PD per cell, excluded from the ratio band). Per-seed dispersion is
wider than the median suggests: 4 of the 14 PD seeds fall outside
[0.33, 3.0] for latent()/dep() (up to a 4.56 ratio), and
spatial_indep() collapses the range in 6 of all 20 seeds (ratios
down to 7e-5): 3 of those 6 are PD -- the rails counted under the frozen
criteria, at ratios 2.3e-4 to 3.4e-4 -- and the other 3 sit among the
6/20 non-PD exclusions. A collapsed field can pass the Hessian check,
the same pattern as the phylogenetic surface's zero-collapse above. The group-intercept gate ((1 | group)
ONLY) also PASSED: 20/20 converged with a PD Hessian, median
sigma_re ratio 0.947, range [0.60, 1.51]. This PASSED verdict does
NOT extend to the cluster/cluster2 diagonal tier -- the s4 campaign
fit re_int exclusively; cluster/cluster2 are admitted with
construction-level evidence only (the fit constructs, extract_Sigma()
returns a well-formed diagonal), and their recovery axis remains OPEN
(a correction to an earlier draft of this entry, which wrongly implied
the same PASSED verdict covered both).
Refused, not merely deferred: phylo_scalar()/animal_scalar()/
kernel_scalar()/spatial_scalar() and common = TRUE on the cluster/
cluster2 tier -- a single shared level across the K-1 contrasts has no
interpretable null on the (I+J) contrast geometry. Null-DGP evidence:
on V_true = 0 data, phylo_indep() correctly recovers near-zero
variance in 5/5 seeds, but phylo_dep()'s rho_hat rails toward +-1
in 4/5 seeds despite a PD Hessian.
Behaviour changes: a (1 | group) or cluster/cluster2 indep() term
whose grouping factor covers exactly one categorical observation per
level now aborts typed (gllvmTMB_multinomial_olre_not_admitted) -- it
is an observation-level random effect in disguise, unidentifiable
because the softmax latent scale is fixed. meta_V()/equalto()
(known-sampling-covariance) remains fail-closed for multinomial()
traits -- confirmed Gaussian-only, no established route on a
categorical-contrast pseudo-trait. cluster2 co-admission alongside
cluster was a maintainer decision (2026-08-16): use_diag_species and
use_diag_cluster2 are literally identical engine math on two different
grouping columns.
This release focuses on multivariate stacked-trait models fitted through the R/TMB engine. Models are fitted by Laplace approximation by default; an opt-in scalar variational research route is also available. The optional Julia bridge remains experimental and is not required for the main workflow.
Integrated models now take any number of named sources, declared with
isdm_sources() (experimental). A portal stream, digitised literature
records, checklists, and a structured survey can each be a named source with
its own observation law, all sharing one ecological process:
family = isdm_sources(gbif = poisson(), literature = poisson(), survey = binomial("cloglog")), with an isdm_source column in the data naming each
row's source. Two laws are admitted: a Poisson count stream and a
complementary-log-log detection stream, because both observe a thinning of
the same shared intensity; that argument holds arm by arm, so it does not
weaken as sources are added. Everything the two-source route refused stays
refused: logit or probit detection, dispersion-carrying families, weights,
multi-trial rows, and any trait not observed by every declared source. The
existing two-source form (list(gbif = ..., survey_pa = ...)) keeps working
unchanged and gives identical fits; it is now the two-source case of the
same rule. A worked example, Integrating three data sources at once:
portal, atlas, and survey, fits a three-source model end to end through the
declared route. Everything reported remains relative intensity, and the
interface remains experimental.
Integrated two-source models can now be fitted through gllvmTMB()
itself (experimental). Opportunistic presence-only records and a
structured detection/non-detection survey of the same species can enter one
likelihood, sharing an ecological linear predictor while each source keeps
its own observation model. There is no separate function for this: supply
family = list(gbif = poisson(), survey_pa = binomial("cloglog")) with
attr(family, "family_var") <- "isdm_family", a source column, and a
matching isdm_family column. This is the one case where the response
family may vary within a trait rather than only between traits, and it is
admitted because the cloglog link makes a detection observation consistent
with the same underlying intensity that generates the counts. An offset()
is admitted on that arm as a known change-of-support term. Anything short of
that exact contract still gets the ordinary one-family-per-trait error, and
Poisson-log mixed with binomial-logit or -probit remains refused.
In: the route is reachable, documented, and contract-tested. Two articles cover it: Integrating opportunistic records with a repeated survey fits and renders through the public route, and How big does an integrated survey design need to be? reads a known-truth design curve to ask whether a design is large enough before you fit at all. Partial: everything reported is relative intensity; presence-only data cannot identify absolute abundance, occupancy, or detectability, and none are estimated. Source-specific spatial structure is only weakly identified on small designs; treat a portal-only field as a nuisance adjustment unless the design is large. Estimator accuracy at the advertised scope is not certified, and the spatial arm in particular rests on development experience rather than a cleared recovery gate. Not included: more than two sources, calibrated intervals for this route, and weighted joint likelihoods. Expect the interface to change.
The exact-aa confirmation completed 44,800 production fits and returned
aggregate HOLD. Gaussian indep(), Gaussian dep(),
Poisson-log rank-1 latent(unique = TRUE), and Binomial(10)-logit rank-1
latent(unique = TRUE) passed their prespecified small/large core
point-estimation gates. Gaussian latent and NB2 latent remained
characterization-only. This does not establish a package-wide dependable
core: the full eligible cell set was incomplete, every campaign contained
failed cells, and the observable terminal-status detector identified only
19 of 80 catastrophic truth errors. The result supports narrow tested-regime
point-estimation statements only; it does not certify intervals, structured
sources, slopes, mixed families, alternative integration engines, or
reliable silent-failure detection. The package remains version 0.6.0 while
this evidence and the remaining pre-0.7 issue backlog are reconciled.
A fit without standard errors no longer returns a silent all-NA answer.
When a model is fitted with gllvmTMBcontrol(se = FALSE), there is no
sd_report, so a Wald interval has nothing to be built from. confint()
used to return a matrix of NA bounds with no error and no warning: a
non-answer that looks like an answer, and that flows onward into tables and
plots with nothing marking it.
confint(method = "wald") now raises a typed error
(gllvmTMB_confint_no_sdreport) naming both remedies:
fit <- standard_errors(fit), or method = "profile", which does not need
standard errors. This covers both Wald routes: the fixed-effects path and
the variance-component target path (parm = "sigma_eps" and friends). The
extractor-style consumers getREsd(), getLV(se = TRUE) and
predict(se.fit = TRUE) already behaved this way.
The other way to have no usable standard errors is also covered. A fit
whose Hessian is not positive-definite has an sd_report, but its standard
errors come back non-finite, and those used to print as a wall of bare NaN.
summary() now says so, and confint() aborts
(gllvmTMB_confint_nonfinite_se). The advice deliberately differs: this is a
property of the fit, so standard_errors() cannot help and the message points
at gllvmTMB_diagnose() and method = "profile" instead.
A single NA is left alone. A coefficient fixed via Xcoef_fixed has no
standard error, and NA is the right answer there; only an entirely
non-finite set is treated as the pathology.
summary() deliberately still works. Fitting fast and reading point
estimates is a legitimate workflow, so summary() prints as before, but it
now says why the Std.Err column is empty instead of leaving a column of
bare NAs to be read as a computed result. extract_cutpoints() reports the
same way for its tau_se column.
This is a deliberate behaviour change, made rather than staged: the package is pre-1.0 and experimental, so no released code depends on the old return, and anything that did would have been depending on a wrong answer.
gllvmTMBcontrol(integration = "va") now admits all 18 scalar family/link cells and defaults to seven-node Gauss-Hermite evaluation. Each cell passed an independent arithmetic, compiled, and light-fit gate before H = 7 and automatic GH routing were promoted; analytic cells retain exact fast paths, explicit JJ remains available only for binomial-logit comparisons, and multinomial or other non-scalar architectures remain excluded. The preregistered 36,000-fit confirmation is complete: only the Poisson-log q = 5 cell passed the overall point route, while 24 of 36 family-by-rank cells failed and 11 were inconclusive. Against the campaign's cell-specific criteria, fixed-effect VA-Wald outcomes were 20 pass / 16 fail, while latent posterior-SD outcomes were 15 pass / 20 fail / 1 inconclusive. These are campaign outcomes, not public calibration certificates: vcov() and confint() remain labelled calibrated = FALSE, getLV(se = TRUE) remains a variational posterior SD rather than a frequentist standard error, and Laplace remains the package default. Both campaign stages ran on Totoro, so this is not cross-platform confirmation.
vcov() and coef() now work on multi-trait fits. coef(fit) returns
the named fixed-effect estimates and vcov(fit) their covariance, taken from
the fit's sdreport(). Both were previously registered only for the
variational gllvmTMB_va class, where coef() deliberately refuses and
vcov() is now restricted to uncalibrated fixed-effect VA-Wald inference,
so calling either on an ordinary fit raised "no applicable method", despite
the documentation saying otherwise.
coef() works on any fit, including one made with se = FALSE: point
estimates do not depend on sdreport(). vcov() does, and raises the same
typed conditions confint() does when it is missing or non-finite, so a
caller that handles one handles the other. Rows and columns for coefficients
held fixed via Xcoef_fixed are NA; a parameter that was not estimated has
no sampling covariance.
standard_errors() computes standard errors after fitting. Fitting with
gllvmTMBcontrol(se = FALSE) skips the TMB sdreport() and is meaningfully
faster, but it used to be a one-way door: the only route to standard errors
afterwards was to fit the model again. fit <- standard_errors(fit) now
computes them on demand from the fitted object.
The result is the fit-time calculation deferred, not a different one: the
same single sdreport() call on the same converged parameter vector,
verified bit-exact (tolerance = 0) against a fit made with se = TRUE.
Every existing consumer (summary(), getLV(), getREsd(),
and confint(method = "wald")) then works on the returned object.
Note the R semantics: the fit is returned, not modified in place. Assign the result, or the standard errors are discarded.
Same-session only. A TMB ADFun holds external pointers that do not
survive saveRDS() or a new R session, so this cannot revive a saved fit.
That limitation is shared by every part of the package that reuses the fitted
TMB object; what is new is that this function says so with a clear, typed
error instead of failing obscurely. This adds no new inference and changes no
likelihood, parameterisation, or honesty caveat; Wald standard errors carry
exactly the caveats they carried before.
Spatial mesh, CRS, and range-plot helpers were substantially rewritten
within gllvmTMB. The retained implementation descends from earlier
GPL-3 sdmTMB helpers and keeps that attribution in inst/COPYRIGHTS.
make_mesh() returns a gllvmTMBmesh built through
fmesher's public mesh, finite-element, and basis APIs, while valid legacy
sdmTMBmesh objects receive a temporary lifecycle warning and conversion.
plot_anisotropy() and plot_anisotropy2() now show the fitted isotropic
practical range, sqrt(8) / kappa, for native gllvmTMB spatial fits. Equal
axes are explicitly the model assumption H = I; gllvmTMB does not estimate
directional anisotropy. Delta and spatiotemporal states remain unsupported.
This changes no TMB likelihood or spatial covariance parameterisation.
The mesh/plotting helpers themselves are covered by focused tests; broader
spatial-family evidence remains partial; directional anisotropy,
delta/spatiotemporal fields, barriers, and new spatial likelihoods are
rejected here.
offset() now works, for count responses. offset(log(trap_nights))
in a Poisson or negative-binomial model is the standard way to model a rate
rather than a raw count, and until now gllvmTMB rejected it. Both closest
comparators offer offsets, so this was a gap rather than a rough edge, and
for counts there is no workaround, since an effort adjustment cannot be
folded into the response the way a Gaussian one can be centred by hand.
It is deliberately restricted to count families: poisson(),
nbinom1(), nbinom2(), truncated_poisson(), truncated_nbinom2().
An offset is a multiplicative rate adjustment on the log link; under
gaussian() it would be an unexplained mean shift and under binomial()
a fixed shift in log-odds, neither of which is what the term is for. In a
stacked model the check is per trait, so the refusal says which one:
offsets are supported for count families (poisson, nbinom) only;
trait `t2` uses `gaussian`.
A mixed-family design allows the offset to vary by trait. A single-family package can only recycle one offset across every response.
A zero offset is allowed everywhere and does nothing, because zero on
the log scale is a multiplier of one. That is how a mixed-family model
gives an offset to its count traits and not the rest; set the offset
column to 0 on the other rows. In wide format, offset(w) applies one
unit-level column to every trait, while offset(e1, e2) gives one column
per trait in traits() order.
The offset also reaches simulate(), bootstrap_Sigma(),
coverage_study(), and predict(newdata = ), all of which rebuild the
linear predictor themselves and would otherwise have used a model you did
not fit. predict(newdata = ) needs the offset variable present in
newdata and says so if it is missing. engine = "julia" has no offset in
its predictor and refuses the term rather than dropping it.
Offsets inside an lv = ~ ... sub-formula, or inside an impute / mi()
covariate formula, are still not supported.
A loading penalty is available for fits that run away, via
gllvmTMBcontrol(aghq_ridge = tau). Binomial fits at small sample sizes can
drive one trait's loading to an absurd value while reporting every
conventional sign of health: convergence = 0 and a positive-definite
Hessian, because quasi-complete separation makes that solution the genuine
maximum of the likelihood. A Gaussian ridge on the loadings adds curvature
where the likelihood is flat and removes the runaway: measured at 47% of
fits down to 0% at n = 100, and on one reproduction fit it takes the largest
implied loading norm from 979.1 to 3.35.
On the default Laplace route the penalty is opt-in and never applied unless
you name it, so no existing Laplace fit changes. AGHQ is itself opt-in; once
AGHQ is enabled, its default is the penalised tau = 2 MAP route unless
aghq_ridge = Inf requests an unpenalised fit. tau is the prior standard
deviation on each free loading, and
its scale is set by the model rather than tuned: the latent scores are
standard normal by identification, so a loading is the trait's latent standard
deviation in link units, making tau = 2 weakly informative. Its influence
also vanishes as the sample grows, because a fixed penalty contributes a
constant against a log-likelihood growing with n. The penalty is
rotation-invariant.
Two costs, stated plainly. A penalised fit is a maximum-a-posteriori point,
not a maximum-likelihood estimate, so logLik(), AIC() and BIC() no
longer describe it; set aghq_ridge = Inf and refit every model being
compared if you need likelihood-based comparison. And the penalty currently
covers the unit-tier loadings only. check_gllvmTMB() now names this remedy
when it reports a runaway loading.
aghq_ridge = "auto" adds a narrower experimental, opt-in scale-aware
route for pure single-trial Bernoulli models with one ordinary unit-tier
latent() block. It first runs an unpenalised 9-node multi-start AGHQ pilot,
then uses
tau = min(6, max(1, ||Lambda_pilot||_F / sqrt(p q))). A plain Laplace fit is
never used as the scale yardstick. If that pilot or the scale-aware final fit
is unusable, the function returns an independently started tau = 2 fit and
records the reason in fit$aghq$ridge_auto; it does not silently claim that
auto-selection succeeded. Both pilot and returned AGHQ fits use the fixed
9-node multi-start rule evaluated for failure/runaway avoidance in this
Bernoulli grid; this is not interval calibration or general estimator
certification. Conflicting aghq or aghq_multistart controls
are replaced with a warning.
This scope is supported as a runaway/failure-avoidance capability, not as
a general accuracy improvement. In the 600-replicate evaluation,
the auto fit was usable in 135 replicates. A transparent auto-when-usable,
otherwise-tau = 2 policy did not increase per-cell failure or runaway rates
and stayed within the preregistered +0.02 loading-error non-inferiority margin,
but its six-cell macro-mean loading-error difference was +0.00282 (slightly
worse).
Caps 5, 6, and 8 were indistinguishable for fresh valid pilots; cap 6 is the
smallest candidate supported by the stored wider-range rescore. Other
families, covariance tiers, and likelihood-based model comparison remain
outside this claim. The measured grid used logit, p = 6, q = 2, and
n = 100, 400, or 1600; other links and dimensions are explicit
extrapolations, not covered accuracy claims. The package default stays
tau = 2.
Adaptive quadrature (gllvmTMBcontrol(aghq = k)) now tries two starting
points and keeps the better fit. It previously ran from a single start: the
Laplace optimum, on the grounds that without a penalty there is nothing to
choose between two starts. That holds for two starting points; it does not
hold for two converged fits, which can simply be compared. It matters
because the Laplace optimum is sometimes itself a runaway, and quadrature then
inherits it: measured over 40 binomial fits at n = 100, sixteen ran away
catastrophically, and on all sixteen the second start reached a strictly
better objective, by 1.1 to 12.9 in negative log-likelihood, and a far more
plausible loading matrix. Catastrophic fits fell from 16 in 40 to 1 in 40.
The second start is data-driven and uses no knowledge of the truth; the cost
is one extra adaptation run. aghq_multistart = FALSE restores the previous
single-start behaviour exactly. Note that this switch was documented but
did not previously take effect, because it was not among the arguments
gllvmTMBcontrol() accepted.
Quadrature fits can now report convergence at realistic sample sizes, and
every stop says why. The convergence test compared the gradient against a
fixed threshold. Because a likelihood's gradient grows with the amount of
data, that threshold became unreachable as the sample grew: no fit at n = 400
or n = 1600 could be certified as converged, in any family tried, including
cases where the quadrature had landed on precisely the point the Laplace fit
reported as converged. The test now also accepts a relative gradient
(aghq_grad_tol_rel). This changes the verdict, not the estimate: fits are
identical either way.
Every stopping condition now reports its gradient, including the one that
previously reported none, so a genuine stall can be distinguished from a near
miss. Read fit$aghq$converged for the verdict: not
fit$opt$convergence, which on this path records the optimiser's per-pass
iteration cap and therefore reports a limit even on a healthy fit. fit$aghq
also now carries grad_max, grad_rel, and how many starts were run.
getLV() gains an se = TRUE argument that returns the standard error of
every unit-level (or within-unit) latent score, alongside the scores, as
list(scores, se). The default se = FALSE is unchanged (a bare matrix).
SEs come from the fitted model's TMB sdreport() random-effect block and
were cross-checked against an independent joint-precision-inversion route
to machine precision. Requires rotate = "none"; predictor-informed
latent(..., lv = ~ x) fits and engine = "julia" bridge fits are not
yet supported and raise an informative error.
gllvmTMB() now accepts both canonical long data and wide data through a
traits(...) left-hand side. The two forms use the same fitting engine.
Random slopes accept the || uncorrelated coupling: mode(1 + x || g)
fits the intercept and slope with no intercept-slope covariance (equivalently
mode(1 | g) + mode(0 + x | g)), alongside the correlated single-bar
mode(1 + x | g). Available for phylo_indep/animal_indep/kernel_indep
(per-trait diagonal), phylo_dep/animal_dep/kernel_dep (block
Sigma_int (+) Sigma_slope), and the source-tier phylo_latent/animal_latent/
spatial_latent (which were already the uncorrelated form).
Dense-kernel random slopes: kernel_indep(1 + x | g, K = K) and
kernel_dep(1 + x | g, K = K) (and their || forms) fit a random regression on
a supplied dense kernel K, byte-equivalent to the phylogenetic path with
vcv = K.
Random-slope routes are runtime-permitted for lognormal and Student-t responses. This is fit admission only: direct route-specific recovery and inference evidence are not yet covered, so this release makes no scientific-validation claim for those combinations.
New multinomial() response family for an unordered categorical response
with three or more categories (baseline-category logit / softmax). It recovers the
per-category intercepts and slopes as contrasts against a reference category, and
predict(type = "response") returns per-category probabilities. Use
multinomial(baseline = ...) to choose the reference category. The validation
boundary is explicit: fixed-effect point recovery is validated: no
detectable bias (|bias| ≤ 0.02 per coefficient across a 500-seed calibration,
with recovery asserted on a 20-seed aggregate at K = 3 and K = 4 rather
than any single fit), and that is a statement about the estimates, not about
their intervals. The two covariance routes: a single phylo_latent() term,
and the narrow ordinary shared-latent() cross-family route, are
only partially validated, meaning
they fit and report but their recovery has not been certified. The latter reports the nominal
trait as its K - 1 baseline-contrast block rather than inventing one scalar
categorical correlation; it permits one multinomial trait per fit and rejects
unsupported tiers before TMB construction. Multiple multinomial traits,
augmented slopes, explicit multinomial unique()/indep(), and unlisted
source tiers remain blocked. A two-category response is
binomial(link = "logit"). For
ordered categories use ordinal_probit(). See the Unordered categories with
multinomial() article for a worked diet-guild example.
phylo_latent() on a multinomial() trait (partially validated) reports the
(K-1) x (K-1) among-category phylogenetic covariance V (how the category
liabilities coevolve) via
extract_Sigma(fit, level = "phy", part = "shared", link_residual = "none").
The default total/"auto" extraction instead reports V plus the fixed softmax
residual (pi^2/6)(I + J). Two honest caveats:
recovery of V is data-hungry (it needs per-species replication or large N; a
single categorical draw per species is weakly informative, so one-per-species point
estimates are high-variance and can reach the +/-1 boundary), and V is on the
baseline-contrast scale, so a diagonal V is not independence -- the null contrast
covariance is (I + J)-structured (equal diagonal, equal off-diagonal; the
observation-scale link residual is applied as (pi^2/6)(I + J) -- the softmax
analog of binomial's pi^2/3). Treat this phylogenetic V route as
recovery-oriented and data-hungry, not universally validated.
For the admitted cross-family nominal route (partially validated), ordinary
latent() keeps its default diagonal companion but the current engine maps
off multinomial-contrast Psi. That variance is not identified with one
categorical draw per unit; replication can identify it in principle, but the
current conservative implementation still suppresses it. Explicitly adding
unique() or indep() for those contrasts remains fail-closed. Point
extraction and target-specific
Wald/bootstrap interval plumbing exist, but their repeated-sampling
calibration is not covered. Nonlinear profile intervals are withdrawn.
extract_cross_correlations() now restricts level to the ordinary unit tier
for every method. Previously only method = "profile" enforced this, so
level = "unit_obs", "phy", or "spatial" combined with method = "point",
"wald", or "bootstrap" were reachable. Those combinations now raise a typed
error. This is a deliberate reduction rather than a regression: the estimand for
a source-tier cross-family correlation was never validated on those paths, and
returning an uncalibrated number was worse than refusing. Use
extract_Sigma() for source-tier covariance.
Known limitation: random-slope covariance is not calibrated when each
cluster carries little information. This is a limitation of the data
regime, not of one keyword: it applies to any random-slope covariance fitted
on single-trial binary responses with few observations per grouping level,
and it affects both the current phylo_indep() / animal_indep() /
spatial_indep() slope forms and the soft-deprecated *_unique() forms.
Measured on a phylogenetic slope fit with a logit link (60 species, 4
replicates, 3 traits, 12 single-Bernoulli observations per species), the
whole 2x2 slope covariance is over-estimated, not just its slope entry:
| target | true | relative error | |---|---|---| | intercept variance | 0.40 | 0.82 | | slope variance | 0.30 | 0.78 | | intercept-slope correlation | 0.50 | 0.367 (absolute) |
The bias does not shrink with more clusters: it persists across 60, 120
and 240 species, on fits that are otherwise healthy (converged,
positive-definite Hessian, valid sdreport). The cause is too little
information per cluster: with a handful of single-trial binary observations
per species, the sampling variance of each species' estimated slope is
comparable to the true between-species variance itself, so roughly half the
spread across species is sampling noise. The identical design recovers cleanly
under a Gaussian response, which is what rules out an engine problem.
Do not read a random-slope variance or correlation from sparse binary data
as calibrated. The remedy is more information per grouping level: more
replicates per species, or multi-trial cbind(successes, failures) data
instead of single 0/1 draws, rather than more species. Note also that the
binomial slope routes are covered by a structural contract only: those
tests check that the model fits and reports the right shapes, and
deliberately do not certify variance recovery or interval calibration. The
corresponding recovery test is deliberately skipped rather than passed by
retuning its data-generating truth.
The reader-facing covariance grammar crosses five correlation sources
(none, animal, phylo, spatial, and kernel) with three taught modes:
independent, dependent, and latent. The one-shared-variance ("scalar") case is
the parsimony modifier common = TRUE on any indep term.
The retained pkgdown guides now focus on runnable, numerically inspected workflows: morphometrics, Gaussian latent-rank selection, binary joint species distribution modelling, behavioural covariance, reaction norms, phylogenetic covariance, missing responses, response screening, fit diagnostics, profile routes, and formula/reference concepts.
predictive_check(), diagnostic residuals, diagnostic_table(),
check_gllvmTMB(), and gllvmTMB_diagnose() provide complementary fitted-model
checks. They diagnose a fitted response distribution and numerical health; they
do not prove latent rank or interval calibration.
extract_Sigma_table(), plot_Sigma_table(), and plot_correlations() provide
report-oriented covariance and correlation displays. Correlation extraction is
point-only by default; interval routes are explicit and carry an uncalibrated
status rather than a coverage certificate.
The default number of bootstrap replicates is now 999, raised from 200.
bootstrap_Sigma() is the only exported function whose signature changes.
The same raise was applied to the internal bootstrap paths behind
extract_lv_effects(), extract_communality(), extract_repeatability(),
the loading intervals, and the phylogenetic-signal intervals, so every
bootstrap interval in the package now uses one replicate count: those
extractors take longer and return slightly different bounds without any
change to their own arguments. Calls that relied on the old default return
slightly different interval bounds and take roughly five times longer. Pass
n_boot = 200 to bootstrap_Sigma() to restore the previous behaviour; for
exploratory work that remains a reasonable time-for-precision trade.
Two distinct things bound this argument, and only the lower one is a
correctness constraint. A percentile interval built from B draws is bounded
by its widest possible realisation, whose coverage cannot exceed
(B - 1) / (B + 1) whatever the data are, so below roughly 2 / (1 - conf)
the interval cannot reach the requested level. Above that threshold what
remains is Monte Carlo error in the endpoints, which is a precision question.
The old default of 200 already cleared the ceiling comfortably (0.990 at
conf = 0.95), so this change buys endpoint precision for intervals that get
reported; it does not fix a correctness bug. 999 rather than 1000 so that
(1 - conf) / 2 * (B + 1) is a whole number at conf = 0.95, letting the
bounds land on order statistics instead of being interpolated between them.
bootstrap_Sigma() refuses an n_boot below the arithmetic floor, warns
below the default, and returns $coverage_ceiling so a simulation campaign
can assert coverage_ceiling >= conf on its own configuration before trusting
its own numbers. The guard and its tests are covered; bootstrap interval
calibration for non-Gaussian and mixed-family fits remains uncertified.
Ordinary latent() now represents
Sigma = Lambda Lambda^T + Psi by default. Use
latent(..., unique = FALSE) for the earlier loadings-only subset.
Source-specific and kernel latent terms remain loadings-only by default; pass
unique = TRUE when their intended covariance includes the diagonal companion.
The one-shared-variance ("scalar") covariance is now the parsimony modifier
common = TRUE on any indep term: indep(..., common = TRUE),
phylo_indep(..., common = TRUE), animal_indep(..., common = TRUE),
spatial_indep(..., common = TRUE), and kernel_indep(..., common = TRUE) fit
one variance shared across all traits (intercept-only). The covariance grid is
taught as three modes -- independent, dependent, latent -- with common = as
the scalar sub-case, rather than a separate fourth mode.
Current phylo_indep(), animal_indep(), and spatial_indep()
intercept-and-slope terms fit one independent 2 x 2 (intercept, slope)
block per trait: within-trait correlation is estimated for |, fixed to
zero for ||, and cross-trait covariance is zero. Current *_dep() routes
instead use a full 2T x 2T augmented covariance. The soft-deprecated
phylo_unique(), animal_unique(), and spatial_unique() slope forms retain
their legacy shared 2 x 2 channels; they are not aliases for the current
*_indep() shape. Admission is decided separately for each response family
and each random-effect route, so a combination that fits is not thereby
validated: some routes are admitted with recovery evidence, others are
permitted at fit time only. This release does not publish a per-route
coverage table, so the documentation will not tell you which of the two a
given combination is. Treat a successful fit as evidence that the model is
admissible, not that its variance components or intervals have been
validated. Where a route is known to be weak this changelog says so
explicitly; see the random-slope limitation above.
fit$fit_health separates optimiser success, raw and objective-scaled
gradients, Hessian health, and sdreport() availability. Its converged field
is conservative: optimiser success, a finite objective, and a small raw maximum
gradient are all required. Hessian health remains a separate inference check.
extract_correlations() now returns point estimates by default. Fisher-z/Wald
bounds are heuristic sensitivity summaries; bootstrap routes are
target-specific and are not labelled as universally calibrated. The nonlinear
penalty-profile prototype is no longer a public route.
Reader-facing pages no longer expose internal validation identifiers, development phases, agent roles, or capability bookkeeping.
Standardized loading inference now uses the model-implied total variance,
rho[t,k] = Lambda[t,k] / sqrt(Sigma_total[t,t]), rather than an entrywise
loading-plus-scalar approximation. loading_ci(method = "wald_asym"),
suggest_lambda_constraint(method = "wald_retention"), and
suggest_lambda_constraint(method = "varimax_threshold") now account for
every latent axis in the trait denominator; the Wald routes propagate the
full joint fixed-parameter covariance, including fitted variance components
and parameter-dependent link residuals. Loading CI, flagging, plotting,
bootstrap, and confint(..., parm = "Lambda") outputs now label their
loading_scale; raw symmetric Wald inference remains the default.
Deterministic algebra and routing are covered, but standardized interval coverage is not yet calibrated. Wald retention treats the fitted varimax rotation as fixed, and standardized profile/bootstrap intervals are not currently available. Per-axis loading intervals therefore remain rotation-conditional sensitivity summaries rather than the primary rotation-invariant evidence.
Ordination no longer collapses when the response is on a large scale. Starting values for the latent structure were built as though the response had a standard deviation of about 1: a hardcoded loading start, a matching variance start, and latent scores beginning at exactly zero. Standardising the latent scores is precisely what pushes the response scale into the loadings, so where that assumption did not hold the fit could collapse: loadings, the covariance, the fixed effects, and the correlations and communality that most users actually report all came back wrong, with every convergence signal green. Nothing warned you.
All three starts now follow the data. The loadings and the variance term are placed on the scale of the working residuals, and the latent scores are seeded from the data instead of from zero, scaled to unit variance, since the scores are standardised by definition and it is the starting direction that was missing, not the magnitude.
For an ordinary single-tier latent() model on 4 traits, multiplying the
response by 100 now reproduces every expected transformation to within about
one part in 100,000. At a factor of 5000 it holds to about 1% in the worst of
eight simulated datasets, looser, but well inside the 2% we accept, and the
quantities most people report (the correlations and the communality) hold
there too. Two other implementations of the same model do not hold that law on
the same data: their worst cases are roughly 100% and 200% out, meaning a
loading that should have doubled did not move. If you have worked around this
by rescaling your response by hand, you no longer need to.
This covers Gaussian responses, and the ordinary latent structure on your
unit grouping: latent() and the variance term it carries. Everything else
deliberately keeps its previous starting values: any model with a
non-Gaussian response, the phylogenetic, spatial, kernel and random-slope
latent terms, and the second grouping of a nested two-tier fit (the
within-unit level, e.g. species-within-site). Starting values for the latent
scores are also left alone when you give latent() a predictor with
lv = ~ ..., because there the scores have a fitted mean rather than a free
one.
The reason for drawing the line there rather than wider: "multiply the response by 100" is only a meaningful thing to do to an unbounded continuous response, and that is the only case the collapse was measured on. A count or a presence/absence on an awkward scale may well have the same problem, but moving a starting value that has not been measured trades a known problem for an unmeasured one, which is the whole reason the old default was dangerous. If you fit any of the excluded forms on a response far from unit scale, rescaling it by hand is still worth doing.
One case is not fully resolved. In a nested two-tier fit the residual scale error is around 2%, and that remainder is a property of the likelihood surface rather than of the starting values: a fit can report convergence while sitting some distance from the optimum, because a wide region of parameter space is nearly flat there. That is tracked separately and no starting value will address it.
check_gllvmTMB() no longer passes a binomial fit whose loading has run
away. The loading row could only fire when the trait's marginal prevalence
was also extreme (at or beyond 0.9). But quasi-complete separation is a
property of the fitted linear predictor, not of the marginal rate, so it runs
a loading away while prevalence stays entirely ordinary, and the row was
keyed on a quantity the pathology does not move. Across 3,944 simulated
binomial fits the worst-affected trait's prevalence never left 0.20 to 0.807,
and its distance from 0.5 was essentially uncorrelated with the size of the
blow-up, while loadings reached 24,000 times the typical trait's. On one
Bernoulli fit the row reported PASS with a loading 6,980 times typical and
every fitted probability saturated, while the implied covariance was wrong by
a factor of 156,645. A loading at or beyond loading_runaway_thresh (new
argument, default 25) now reports on its own, naming the improper solution
(Heywood case) and pointing at the loading penalty in
gllvmTMBcontrol(aghq_ridge = ). The existing loading_relative_thresh of 8
keeps its prevalence conjunct, because healthy fits with a sparse loading
structure reach that level routinely, so nothing that was flagged before
stops being flagged, including a genuinely near-constant trait, which the row
still reports as it always did.
A second, complementary criterion is added alongside it. A relative
criterion cannot see a loading matrix that is inflated as a whole, because
scaling every loading leaves every ratio unchanged, so
loading_absolute_thresh (new argument, default 6) reports a loading that is
simply too large on the link scale. That threshold is meaningful because the
latent scores are standard normal by identification, making a binomial
loading the trait's latent standard deviation in link units; a value of 6
already implies a fitted probability indistinguishable from 0 or 1 across an
ordinary swing of the axis. Measured over the same 3,944 fits: no healthy fit
exceeded 3.99, none was flagged, and it reported 97.3% of degenerate fits,
catching 14 that the relative criterion missed. Being a link-scale quantity
it does not transport to families whose response scale is arbitrary, which is
why this row remains binomial-only.
Two limits on the calibration are stated plainly. It was measured on single-family binomial fits at the true latent rank, where it reports 96.3% of fits whose implied covariance is wrong by a factor of five or more, with no healthy fit reaching the threshold (the largest was 12.1). When the fitted rank is larger than the truth, several traits can inflate together, which lifts the very yardstick the ratio is measured against; in one such run the row missed 3 of 8 degenerate fits, and they were the three worst. A scale diagnostic, not a ratio, is the right instrument for that case and none is wired yet.
check_gllvmTMB() now reports a unique variance that has collapsed only
relative to its siblings. A Heywood case in a Gaussian or Poisson fit
usually appears as a per-trait unique variance driven to the boundary, not as
a runaway loading, and because psi is estimated on the log scale, the
boundary is an interior point of the transformed space, so pdHess stays
positive definite and nothing else objects. Across 360 fits with a
deliberately over-specified latent rank, 58% drove a unique standard
deviation below a tenth of its true value while reporting convergence = 0
and pdHess = TRUE: one reached 6e-50. The covariance itself was still
recovered to within 7%, so this is a failure of the
Lambda Lambda' + Psi decomposition rather than of the fitted covariance,
and no recovery-based check can see it.
psi_rel_thresh is raised from 0.001 to 0.01, which reports 96.2% of
those fits rather than 73.7%. The measured false-positive rate is zero
both on 151 healthy fits and on 359 healthy fits whose true unique variances
differ by up to a factor of 1000, the case that decides whether the number
transports, since a small ratio is then correct rather than pathological.
Looser values do not transport: 0.1 reaches full sensitivity but flags 19% of
those healthy fits. Some fits that previously passed will now warn; on
this evidence they are fits with a genuinely boundary-pinned component.
The typical loading size is now taken over the traits being screened. Previously it pooled every trait in the fit regardless of family, so in a mixed-family model a trait on a large response scale could set the yardstick for a binomial one, masking a genuine runaway, or manufacturing a spurious one when the other family's loadings were small.
Profile confidence intervals no longer lose their bounds at higher
confidence levels. The profile search used a fixed deviance budget that did
not depend on level, and that budget was smaller than the threshold a
profile must cross for any level above roughly 0.955. Above it the search
could stop short and the bound was reported as infinite. On a four-trait
Gaussian fit, level = 0.99 returned an infinite bound for four of ten
targets where level = 0.95 returned all ten finite. The budget is now sized
from the requested level, with headroom; merely reaching the threshold is
not enough, because the bound is located by interpolating across it. The
default level = 0.95 also gains margin it did not previously have.
An unbounded interval is no longer reported where the bound is simply
unknown. A profile that stops without crossing its threshold can mean two
different things, and only one was reported. If the profile has flattened out,
no finite bound exists and an infinite bound is the honest answer. If the
profile was still climbing when the search stopped, the bound is unknown;
calling it infinite asserts an unbounded parameter on the strength of having
stopped looking. The two are now distinguished from the shape of the profile
itself, and the second returns NA.
Interval bounds are more accurate. Bounds are now interpolated on the
signed-square-root (zeta) scale rather than the deviance scale. For a
quadratic log-likelihood zeta is exactly linear in the parameter, so the
interpolation is exact there; on the deviance scale it carried curvature error
that grew with the spacing of the profile grid. Reported bounds may shift
slightly, and are closer to what a finer grid would give.
The internal coverage-study helper no longer counts an interval it could not compute as a success. A missing bound was already treated as non-coverage, but an infinite bound fell through that guard, was scored as covering, and stayed in the denominator. A method returning an unbounded interval on every replicate would have reported perfect coverage. Missing and infinite bounds are now both excluded and reported as excluded. The helper is internal and unexported; the correction is recorded because any coverage figure produced with an earlier version of it is affected.
Near-zero variance components are now detected relative to their
siblings, not only against an absolute threshold. A boundary-pinned
(Heywood) component, one trait's unique variance collapsing to zero, could
previously pass every check the package had: check_gllvmTMB() reported
near_zero_psi_unit … PASS and fit_health$boundary_flags stayed empty for a
component whose variance was six orders of magnitude below the others. The
absolute threshold is expressed on the standard-deviation scale, so the old
1e-4 demanded a variance below 1e-8 before flagging anything. pdHess
could not catch it either: psi is estimated on the log scale, so a collapsed
component is an interior point of the transformed parameter space and the
Hessian stays positive definite there. check_gllvmTMB() gains
psi_rel_thresh (default 0.001) and reports the offending ratio in the
check message; .gllvmTMB_boundary_flags() applies the same relative test to
every variance block it already covered. Fits that were flagged before are
still flagged.
A diagonal covariance term is no longer duplicated when the unit and
cluster columns are the same grouping factor. This removes a flat variance
split and restores coherent covariance extraction and Wald infrastructure.
extract_phylo_signal() now uses the declared species-level denominator in
crossed site-by-species designs instead of silently returning one for every
trait when non-phylogenetic species variance was stored at the cluster tier.
phylo_dep() is treated as a full covariance parameterisation, not as a set of
exchangeable latent axes, in rotation and weak-axis diagnostics.
Missing or undefined link-scale residual variances now propagate as NA
instead of being replaced by zero or another finite fallback. In particular,
Student-t variance is undefined when its degrees of freedom are at most two.
Several optional Julia-bridge shape, dispatch, missing-cell, and confidence- interval error paths now fail explicitly instead of silently returning malformed output. The bridge remains experimental.
The sigma_d2 argument to loading_ci(),
suggest_lambda_constraint(), and suggest_lambda_constraints() is
deprecated and ignored. Standardized loading denominators are now derived
from each fitted model's total trait variance.
gllvmTMBcontrol(start_method = list(method = "res")) is soft-deprecated
and warns once per session. It still fits; prefer the default starts.
It was retired on measurement. Across 89 simulated fits (Gaussian, Poisson
and negative-binomial; d = 1 to 3; three and five traits), the residual
start was never materially better than the default start (its three best
margins were 0.07, 0.29 and 0.66 log-likelihood units, the scale of landing on
a slightly different point of the same optimum), was materially worse eight
times, once by 14.6 units, and was exactly neutral at d >= 2
(objectives agreed to seven decimals in 34 of 34 fits). Every failure reported
convergence == 0 and a positive-definite Hessian on both sides, so the worse
fit was silent, and restarts were not a guard: in one such fit all five
n_init = 5 restarts returned the same worse optimum.
The damage concentrates at d = 1 with three traits, the exactly identified
corner where a per-trait variance can collapse to zero. Seeding from the
residual covariance commits the optimiser to that matrix's leading direction,
which there is not the best-likelihood factor. This is not a residual-noise
problem: a start built from noise-free random-effect estimates lands in the
same wrong place, so it is not fixable by a better residual.
The formula parser continues to accept unique() as compatibility syntax;
source-specific *_unique() functions remain exported soft-deprecated
aliases. Use indep() / *_indep() in new standalone diagonal formulas.
The scalar family -- scalar(), phylo_scalar(), animal_scalar(),
spatial_scalar(), kernel_scalar() -- is soft-deprecated compatibility
syntax that emits a one-time warning and keeps working. Use
indep(..., common = TRUE) / *_indep(..., common = TRUE), which fits the
same model.
gllvmTMB_wide() remains available for migration, but new wide examples use
gllvmTMB(traits(...) ~ ...).
meta_known_V() remains a deprecated alias of meta_V().
profile_ci_total_variance() with
interval_status = "certified-0.94". That label marks a documented 0.94
coverage floor under its preregistered conditions, not nominal 95% coverage
or a guarantee for an individual interval. Other covariance routes have
cell-specific evidence only. Intervals are route output, not a coverage
guarantee.check_identifiability() and coverage_study() prototypes
have been withdrawn from the exported surface. Their fitted-model simulation
designs did not establish unknown generating rank or retain every attempted
replicate in the coverage denominator. They remain internal until redesigned
around predeclared known-data-generating targets and complete failure accounting.meta_V() remains an important development target, but a dedicated public
article will wait until its supported estimands, diagnostics, and validation
evidence form a complete reader path.Earlier development release establishing the stacked-trait R/TMB engine, the long-format API, initial covariance keywords, simulation helpers, and extractor infrastructure. The 0.7.1 notes above describe the current taught syntax and reader-facing scope.