Screen candidate responses before fitting a gllvmTMB model
Source:R/screen-gllvmTMB.R
screen_gllvmTMB.Rdscreen_gllvmTMB() is a formula-aware pre-fit screen. It summarises
response and formula conditions before fitting a stacked-trait GLLVM.
The screen is advisory: it does not fit the model, remove traits, choose
a latent rank, prove identifiability, or guarantee convergence.
Binary/binomial screening is implemented. Non-binary modules, optional
comparator checks, and high-dimensional benchmarks are not currently part
of this function.
Usage
screen_gllvmTMB(
formula,
data,
family,
unit = NULL,
trait = "trait",
weights = NULL,
missing = miss_control(),
control = screen_control(),
Xcoef_fixed = NULL,
known_groups = NULL
)Arguments
- formula
A
gllvmTMB()formula. Wide data can usetraits()on the left-hand side; long data should use a response column on the left-hand side plustrait =.- data
A data frame.
- family
A response family. Version 1 screens single-family
binomial()fits withlogit,probit, orclogloglinks.- unit
Optional unit/grouping column. If omitted, the first covariance-structure grouping column is used when available.
- trait
Trait column for long data. Wide
traits()calls create this column internally.- weights
Optional weights vector, with the same binomial trial-count semantics as
gllvmTMB()for long-format binomial fits.- missing
Missing-data control; defaults to
miss_control().- control
A
screen_control()object.- Xcoef_fixed
Optional named structural-zero fixed-effect constraints, normalized against the expanded fixed-effect columns as in
gllvmTMB(). This argument is used by the opt-in fixed-design separation screen.- known_groups
Optional named list of character vectors of trait names, each with at least two distinct names (duplicated names within one group are an error). Each element declares a set of traits the user believes forms an exact one-hot (simplex) block or one or more exact pairwise containment (nesting) relations – for example the dummy columns of one categorical review-type or geographic-scope variable, or a broad-realm indicator together with one or more narrower nested realms. Each declared group is checked against the observed complete rows with an exact deterministic test: row sums for a one-hot block, or, for nesting, EVERY pairwise containment among the declared members (so it does not matter which order – or how many mutually incomparable narrower members – the group happens to be written in). When the relations found form a single total chain, the certificate reports that chain; when they form a partial order instead (for example one broad member containing two mutually incomparable narrower ones), every relation found is named and the group still fails as nesting. A group whose members are all constant on the complete rows is reported separately (not as nesting: a containment test is vacuously true against a constant and would not be genuine evidence); a member that is constant while other members of the same group are not is excluded from every pairwise comparison for the same reason. The one-hot test is not whole-group only: when the WHOLE declared group does not sum to 1 (for example a genuine one-hot block declared together with one unrelated trait), a bounded exhaustive search checks every subset of size 2 or more, smallest first, and reports each MINIMAL one-hot subset found as its own certificate –
A + B + C = 1is reported once even inside a larger declared group, never also as a redundant superset. This search is attempted only when the declared group has at most 12 members (2^12 - 13 = 4083subset checks, microseconds even for many rows; declared groups are small by construction – one categorical variable's dummy columns or one nesting relation). A group larger than that bound is never silently reported PASS: the whole-group one-hot and pairwise nesting checks still run, but an explicitknown_group_subset_not_attemptedrow records that the subset search itself was skipped for size, so a "no subset found" result can be told apart from "not looked for". A subset that sums to 1 only because every one of its members is constant is not reported (the same guard as the whole-group case). Each result is its own row in theresponse_dependenciestable, independent of the automatic affine-rank screen below – except that a declared one-hot block, including a certified subset, also counts toward that screen'sunresolvedaffine-dependency count (a nesting relation never does: it is an inequality, not an exact affine relation among the columns). Screening is Bernoulli-only, matching the automatic screen.
Value
A gllvmTMB_screen object. Use screen_table() to extract
report-ready tables.
Details
The first implemented module covers binomial traits. It distinguishes
Bernoulli responses from multi-trial binomial responses, reports the
relevant denominators, and flags constants, sparse minority outcomes,
duplicate or near-duplicate binary traits, rank-deficient fixed-effect
designs, and grouping/rank conditions that should be inspected before
interpretation. In systematic maps, a trait may be a content item or
indicator. In a binary JSDM, a trait may be a species presence-absence
response. Rare or constant species are flagged for inspection and possible
sensitivity analysis, not automatically removed.
The screen also checks the response side for exact higher-order
dependencies that the pairwise duplicate/complement check cannot see –
most commonly a one-hot/simplex block, where several trait columns are
really the dummy coding of one categorical variable and sum to exactly 1
on every complete row (a review-type, geographic-scope, or temporal-scope
set is a typical source). It works by an affine-rank check on the
augmented response matrix M = cbind(1, Y) over complete Bernoulli rows
(rows where every screened trait is observed), using base R's qr() with
tolerance sqrt(.Machine$double.eps) * max(dim(M)) (the usual
dimension-scaled machine-precision rule, so the check tightens
automatically on a larger response matrix rather than using one fixed
constant regardless of size). The $design table's
response_affine_rank row reports the matrix's rank (value) and its
column count (threshold, = ncol(M)); a rank deficiency –
threshold - value – means at least one exact affine dependency exists
among the traits. The $response_dependencies table then adds, where
cheaply recoverable from the matrix's null space, a human-readable
certificate naming the traits and the relation (for example
"A + B + C = 1"). This is a best-effort report, not exhaustive
minimal-subset discovery: with more than one simultaneous exact
dependency the automatic search can fail to decompose the null space into
individually clean certificates, in which case the remaining deficiency
is reported without a certificate (type = "unresolved"), never silently
dropped. Declaring known structure via known_groups sidesteps that
limitation with an exact, deterministic, per-group check. An exact
one-hot block is not itself a bug to silently filter – it is one
categorical variable, and belongs in multinomial() or a deliberate
reference-level coding, not in several unconstrained binary traits.
With screen_control(separation = "fixed"), the screen uses the optional
detectseparation package to classify the observed fixed design. This
includes finite known offsets, whose presence is recorded in the returned
separation table. The offset does not change the fixed-design recession
directions used by the certificate. This is diagnostic support only:
gllvmTMB() currently rejects nonzero offsets for binomial responses.
This certificate covers fixed-effect geometry only. It does not establish finiteness of latent loadings or covariance parameters, validate the full marginal GLLVM likelihood, or select a penalized estimator.
References
Albert A, Anderson JA (1984). On the existence of maximum likelihood estimates in logistic regression models. Biometrika 71:1–10. doi:10.1093/biomet/71.1.1.
Peduzzi P, Concato J, Kemper E, Holford TR, Feinstein AR (1996). A simulation study of the number of events per variable in logistic regression analysis. Journal of Clinical Epidemiology 49:1373–1379. doi:10.1016/S0895-4356(96)00236-3.
Vittinghoff E, McCulloch CE (2007). Relaxing the rule of ten events per variable in logistic and Cox regression. American Journal of Epidemiology 165:710–718. doi:10.1093/aje/kwk052.
Chalmers RP (2012). mirt: A multidimensional item response theory package for the R environment. Journal of Statistical Software 48(6):1–29. doi:10.18637/jss.v048.i06.
Kosmidis I, Schumacher D, Schwendinger F (2026). detectseparation: Detect and Check for Separation and Infinite Maximum Likelihood Estimates. R package version 0.4.0. doi:10.32614/CRAN.package.detectseparation.
Examples
df <- data.frame(
study = factor(seq_len(12)),
a = c(rep(1, 10), 0, 0),
b = c(rep(1, 10), 0, 0),
c = rep(c(0, 1), 6)
)
# The wrapper keeps the rendered example focused on the screen tables;
# inspect warnings in interactive work.
scr <- suppressWarnings(screen_gllvmTMB(
traits(a, b, c) ~ 1 + latent(1 | study, d = 2),
data = df,
unit = "study",
family = binomial()
))
screen_table(scr, "traits")
#> trait response_mode status severity n_obs n_valid n_units n_success n_failure
#> 1 a bernoulli WARN strong 12 12 12 10 2
#> 2 b bernoulli WARN strong 12 12 12 10 2
#> 3 c bernoulli WARN moderate 12 12 12 6 6
#> total_trials prevalence minority_rate minority_count info_fraction invalid_n
#> 1 12 0.8333333 0.1666667 2 0.5555556 0
#> 2 12 0.8333333 0.1666667 2 0.5555556 0
#> 3 12 0.5000000 0.5000000 6 1.0000000 0
#> action message
#> 1 inspect the minority outcome has very few observed trials
#> 2 inspect the minority outcome has very few observed trials
#> 3 inspect the minority outcome has few observed trials
screen_table(scr, "recommendations")
#> scope status action trait
#> 1 trait WARN inspect a
#> 2 trait WARN inspect b
#> 3 trait WARN inspect c
#> 4 pair FAIL collapse_or_recode a / b
#> 5 pair WARN inspect a / c
#> 6 pair WARN inspect b / c
#> 7 design FAIL inspect <NA>
#> 8 response_dependency FAIL collapse_or_recode a, b
#> evidence
#> 1 the minority outcome has very few observed trials
#> 2 the minority outcome has very few observed trials
#> 3 the minority outcome has few observed trials
#> 4 the two traits are exact duplicates on paired rows
#> 5 the two traits are near-duplicates on paired rows
#> 6 the two traits are near-duplicates on paired rows
#> 7 response matrix augmented with an intercept is rank deficient by 1 on 12 complete rows: at least one exact affine dependency exists among the screened traits
#> 8 b = a: this affine dependency is already reported at trait/pair level (a duplicate pair) and was projected out before the one-hot certificate search
#> model_implication
#> 1 inspect coding, missingness, and latent-response inclusion before fitting
#> 2 inspect coding, missingness, and latent-response inclusion before fitting
#> 3 inspect coding, missingness, and latent-response inclusion before fitting
#> 4 inspect whether both traits should enter the first latent-response block
#> 5 inspect whether both traits should enter the first latent-response block
#> 6 inspect whether both traits should enter the first latent-response block
#> 7 inspect the formula before fitting
#> 8 this is one categorical variable coded as several binary columns; keep it out of the latent block as-is, and either recode it as multinomial() or fix a deliberate reference level