Screen binary traits before fitting
Source:vignettes/articles/pre-fit-response-screening.Rmd
pre-fit-response-screening.Rmdscreen_gllvmTMB() is a pre-fit check for candidate
responses. It answers a small question before you fit a stacked-trait
GLLVM:
Are these traits, items, or indicators risky to put into this exact model?
It does not decide which responses to use, remove indicators, choose
a latent rank, prove identifiability, or guarantee convergence. It uses
the same wide-to-long stacking route as gllvmTMB(), then
reports data and formula conditions that are easy to miss before
fitting. It is not a complete dry run of every parser or engine
check.
The current screen covers one binomial() family using
the logit, probit, or cloglog link. Trait-level support works for
Bernoulli and multi-trial binomial responses; pairwise
duplicate/complement screening is Bernoulli-only. Other families return
NOT_CHECKED rather than pretending they were screened.
The screen has three separate jobs:
- Marginal response support: invalid values, unusable denominators, constants, rare outcomes, and strong imbalance.
- Response-pair redundancy: exact or near duplicates, complements, and unusually strong pairwise association on jointly observed units.
- Formula shape: fixed-effect rank, grouping-level support, and whether the requested latent dimension is smaller than the raw number of traits.
Formal complete or quasi-complete separation depends on the
response–predictor design. Rare or constant outcomes can motivate a
formal separation check, but duplicate or complementary response columns
are redundancy conditions, not classical separation. A clean
screen_gllvmTMB() result therefore does not prove that
separation cannot occur.
Setup
The example uses a systematic-map style data set. Each row is a
review. Each binary column records whether a review has a topic, region,
taxonomic group, or outcome indicator. In this setting the word
trait means “response column”; scientifically, it may be a
content item or coding indicator.
For a binary joint species distribution model (JSDM), read
review as the site, plot, survey, or other sampling unit
and read each indicator column as one species’ presence-absence
response. The screen asks whether those response columns are well
supported for the latent block you are about to fit.
n <- 40
reviews <- data.frame(
review = factor(seq_len(n)),
outcome_richness = rep(c(0, 1), length.out = n),
outcome_diversity = rep(c(0, 1), length.out = n),
outcome_no_change = 1 - rep(c(0, 1), length.out = n),
region_north_america = c(rep(1, 37), rep(0, 3)),
level_ecosystem = 0
)This candidate set contains three common problems.
outcome_diversity duplicates
outcome_richness. outcome_no_change is its
exact complement. region_north_america is near-universal,
with only three minority outcomes. level_ecosystem is
constant.
Wide-data screen
Most systematic-map data start wide: one row per review and one
column per indicator. Use the same traits(...) left-hand
side that you would use in a wide gllvmTMB() call.
Here latent(1 | review, d = 2) requests two ordinary
unit-tier latent variables for the review or site units. The screen
warns when d is not smaller than the raw number of traits.
It reports pair redundancy separately and does not estimate an effective
rank after collapsing duplicate responses. This page does not screen
phylogenetic, spatial, or meta-analytic covariance tiers.
The thresholds are advisory and configurable. This example stores the defaults explicitly so the active rules remain visible.
ctl <- screen_control()
data.frame(
rule = c(
"minority count: stronger warning",
"minority count: warning",
"prevalence: stronger information flag",
"prevalence: information flag",
"pair discordance: stronger warning",
"pair discordance: warning",
"pair Hamming rate: warning",
"absolute phi/Jaccard: stronger warning",
"absolute phi/Jaccard: warning"
),
default = c(
paste0("< ", ctl$rare_strong_n),
paste0("< ", ctl$rare_warn_n),
paste(ctl$prevalence_strong, collapse = " to "),
paste(ctl$prevalence_warn, collapse = " to "),
paste0("< ", ctl$discordant_strong_n),
paste0("< ", ctl$discordant_warn_n),
paste0("<= ", ctl$hamming_rate_warn),
paste0(">= ", ctl$phi_strong),
paste0(">= ", ctl$phi_warn)
)
)
#> rule default
#> 1 minority count: stronger warning < 5
#> 2 minority count: warning < 10
#> 3 prevalence: stronger information flag 0.02 to 0.98
#> 4 prevalence: information flag 0.05 to 0.95
#> 5 pair discordance: stronger warning < 5
#> 6 pair discordance: warning < 10
#> 7 pair Hamming rate: warning <= 0.01
#> 8 absolute phi/Jaccard: stronger warning >= 0.95
#> 9 absolute phi/Jaccard: warning >= 0.9Prevalence outside the displayed interval is flagged; it is not a deletion cutoff. Minority counts take priority, which prevents a large data set with a low prevalence from receiving the same warning as a tiny data set with only a few minority outcomes.
scr <- screen_gllvmTMB(
traits(
outcome_richness,
outcome_diversity,
outcome_no_change,
region_north_america,
level_ecosystem
) ~ 1 + latent(1 | review, d = 2),
data = reviews,
unit = "review",
family = binomial("probit"),
control = ctl
)
#> Warning: ! Ordinary `latent()` now includes a per-trait Psi by default (Sigma = Lambda
#> Lambda^T + Psi).
#> ℹ This changed in gllvmTMB 0.2.0; earlier `latent()` was loadings-only (Lambda
#> Lambda^T).
#> → Pass `latent(..., unique = FALSE)` for the old rotation-invariant
#> loadings-only fit.
scr
#> gllvmTMB pre-fit response screen
#> FAIL 4 | WARN 1 | PASS 13
#> Inspect 5 recommendations.
#> status action trait
#> WARN inspect region_north_america
#> FAIL exclude_from_latent_block level_ecosystem
#> FAIL collapse_or_recode outcome_richness / outcome_diversity
#> FAIL collapse_or_recode outcome_richness / outcome_no_change
#> FAIL collapse_or_recode outcome_diversity / outcome_no_change
#> evidence
#> the minority outcome has very few observed trials
#> all observed outcomes are on one side
#> the two traits are exact duplicates on paired rows
#> the two traits are exact complements on paired rows
#> the two traits are exact complements on paired rowsThe summary counts checks, not traits: FAIL 4 combines
one trait-level and three pair-level failures, WARN 1 is
one trait warning, and PASS 13 includes trait, pair, and
formula-shape checks. Statuses are scope-specific. A pair-level
PASS cannot cancel a trait-level FAIL, and a
trait-level PASS cannot cancel an exact redundant-pair
finding.
-
PASSmeans no condition was flagged at the current thresholds; it is not a guarantee of a safe fitted model. -
WARNandINFOare advisory.severitygrades the evidence within that status, soWARNwithseverity = "strong"is not contradictory. -
FAILmeans invalid binomial input, no usable denominator, a constant response that is uninformative for this latent block, or exact pair redundancy on the observed overlap. -
NOT_CHECKEDmeans the requested screen is unavailable for that row or family.
The trait table reports the denominators behind each warning. A 5%
prevalence with 20 Bernoulli trials is not the same evidence as a 5%
prevalence with 100000 trials. n_valid counts usable rows,
while total_trials includes the multi-trial
denominator.
trait_screen <- screen_table(scr, "traits")
trait_screen[, c(
"trait", "status", "severity", "n_valid", "total_trials",
"prevalence", "minority_count", "action"
)]
#> trait status severity n_valid total_trials prevalence
#> 1 outcome_richness PASS none 40 40 0.500
#> 2 outcome_diversity PASS none 40 40 0.500
#> 3 outcome_no_change PASS none 40 40 0.500
#> 4 region_north_america WARN strong 40 40 0.925
#> 5 level_ecosystem FAIL constant 40 40 0.000
#> minority_count action
#> 1 20 keep
#> 2 20 keep
#> 3 20 keep
#> 4 3 inspect
#> 5 0 exclude_from_latent_blockThe pair table catches exact duplicates, exact complements, and
near-duplicate Bernoulli indicators on jointly observed review rows.
n_pair is the overlap denominator;
discordant_n counts different outcomes, and
hamming_rate is that count divided by n_pair.
“Exact” therefore means exact among the paired rows shown, not
necessarily across rows where one response is missing.
pairs <- screen_table(scr, "pairs")
pairs[, c(
"trait_i", "trait_j", "status", "severity",
"n_pair", "discordant_n", "hamming_rate", "action"
)]
#> trait_i trait_j status severity n_pair
#> 1 outcome_richness outcome_diversity FAIL duplicate 40
#> 2 outcome_richness outcome_no_change FAIL complement 40
#> 3 outcome_richness region_north_america PASS none 40
#> 4 outcome_richness level_ecosystem PASS none 40
#> 5 outcome_diversity outcome_no_change FAIL complement 40
#> 6 outcome_diversity region_north_america PASS none 40
#> 7 outcome_diversity level_ecosystem PASS none 40
#> 8 outcome_no_change region_north_america PASS none 40
#> 9 outcome_no_change level_ecosystem PASS none 40
#> 10 region_north_america level_ecosystem PASS none 40
#> discordant_n hamming_rate action
#> 1 0 0.000 collapse_or_recode
#> 2 40 1.000 collapse_or_recode
#> 3 21 0.525 keep
#> 4 20 0.500 keep
#> 5 40 1.000 collapse_or_recode
#> 6 21 0.525 keep
#> 7 20 0.500 keep
#> 8 19 0.475 keep
#> 9 20 0.500 keep
#> 10 37 0.925 keepRecommendations are written as actions to inspect, not automatic decisions. For example, a near-universal indicator may be scientifically central. In that case, report why it stays in the first latent block and run a sensitivity model. If it is only a coding artefact, recode it, collapse it with a related item, or keep it descriptive instead of making it define a latent variable.
recommendations <- screen_table(scr, "recommendations")
recommendations[, c("scope", "status", "action", "trait", "evidence")]
#> scope status action trait
#> 1 trait WARN inspect region_north_america
#> 2 trait FAIL exclude_from_latent_block level_ecosystem
#> 3 pair FAIL collapse_or_recode outcome_richness / outcome_diversity
#> 4 pair FAIL collapse_or_recode outcome_richness / outcome_no_change
#> 5 pair FAIL collapse_or_recode outcome_diversity / outcome_no_change
#> evidence
#> 1 the minority outcome has very few observed trials
#> 2 all observed outcomes are on one side
#> 3 the two traits are exact duplicates on paired rows
#> 4 the two traits are exact complements on paired rows
#> 5 the two traits are exact complements on paired rowsThe formula and grouping checks are separate tables. Here all 40 units have both outcome sides across the screened block, and the requested rank is smaller than the raw trait count.
screen_table(scr, "design")[, c(
"component", "status", "value", "threshold", "message"
)]
#> component status value threshold
#> 1 fixed_effect_rank PASS 5 5
#> 2 unit_levels PASS 40 2
#> 3 latent_rank_1 PASS 2 5
#> message
#> 1 fixed-effect design has full column rank
#> 2 unit grouping has at least two levels
#> 3 requested latent rank is smaller than the number of traits
as.data.frame(table(
status = screen_table(scr, "units")$status
))
#> status Freq
#> 1 PASS 40Long-data screen
The same screen works after the indicators have already been stacked
long. The long form is useful when your workflow stores one row per
(review, indicator) cell.
long <- stats::reshape(
reviews,
varying = list(names(reviews)[names(reviews) != "review"]),
v.names = "present",
timevar = "trait",
times = names(reviews)[names(reviews) != "review"],
idvar = "review",
direction = "long"
)
long$trait <- factor(long$trait, levels = names(reviews)[names(reviews) != "review"])
row.names(long) <- NULL
scr_long <- screen_gllvmTMB(
present ~ 0 + trait + latent(0 + trait | review, d = 2),
data = long,
unit = "review",
trait = "trait",
family = binomial("probit"),
control = ctl
)
long_traits <- screen_table(scr_long, "traits")
long_traits[, c(
"trait", "status", "n_valid", "prevalence", "minority_count"
)]
#> trait status n_valid prevalence minority_count
#> 1 outcome_richness PASS 40 0.500 20
#> 2 outcome_diversity PASS 40 0.500 20
#> 3 outcome_no_change PASS 40 0.500 20
#> 4 region_north_america WARN 40 0.925 3
#> 5 level_ecosystem FAIL 40 0.000 0
identical(
trait_screen[, c("trait", "status", "n_success", "n_failure")],
long_traits[, c("trait", "status", "n_success", "n_failure")]
)
#> [1] TRUEThe printed comparison verifies that the wide and canonical long calls produce the same trait-support summaries in this example. It does not claim that every future table is identical under every formula. Use whichever data shape matches your workflow.
Apply the screen before fitting
The constant response cannot inform this latent block, and the duplicate and complementary outcomes should not enter the same block unchanged. The next candidate keeps one representative response, removes the constant, and retains the near-universal response because we decide it is scientifically important. That retained warning must be documented and revisited in a sensitivity fit.
candidate <- reviews[c(
"review",
"outcome_richness",
"region_north_america"
)]
scr_candidate <- screen_gllvmTMB(
traits(outcome_richness, region_north_america) ~
1 + latent(1 | review, d = 1),
data = candidate,
unit = "review",
family = binomial("probit"),
control = ctl
)
screen_table(scr_candidate, "recommendations")[, c(
"scope", "status", "action", "trait", "evidence"
)]
#> scope status action trait
#> 1 trait WARN inspect region_north_america
#> evidence
#> 1 the minority outcome has very few observed trialsThis rerun contains only the deliberately retained rarity warning. If the rare response is not scientifically essential, keep it descriptive or fit a sensitivity model without it. The screen supports that decision; it does not make the decision.
What to do with warnings
Read FAIL according to its scope. Invalid values or
denominators require data repair. A constant response is valid binomial
data but has no minority support for the proposed latent block; exclude
it from that block or collect informative variation. Exact duplicate or
complementary pairs should be selected, combined, or recoded before
fitting the same latent block unchanged. Always check the reported
n_pair first, especially when responses have different
missingness.
WARN means limited marginal support or unusually strong
pairwise redundancy or association. Inspect the coding, missingness,
overlap denominator, and scientific purpose. Do not lower the latent
rank and assume the problem is gone; a near-constant binary item can
dominate a one-axis model too.
INFO means the screen saw an imbalanced condition but
enough denominator support that a hard warning would overstate the
evidence. This distinction is why the table reports minority counts.
After fitting, still run check_gllvmTMB(). Pre-fit
screening and fitted-model diagnostics answer different questions.
screen_gllvmTMB() catches candidate indicator risks before
fitting; check_gllvmTMB() catches large fitted loadings,
weak axes, optimizer status, and uncertainty-health symptoms after
fitting. The full sequence is shown in Can I trust this fit?.
Evidence base
The defaults are evidence-informed heuristics, not universal deletion rules. Complete and quasi-complete separation is a well-studied maximum-likelihood failure mode (Albert and Anderson 1984; Mansournia et al. 2018), but its presence depends on the response–predictor design. Minority-count rules such as “10 events per variable” are context-dependent heuristics rather than hard laws (Peduzzi et al. 1996; Vittinghoff and McCulloch 2007). This screen reports conditions that deserve inspection; it does not replace a formal separation analysis.
Rare species are not automatically discarded in JSDMs. Shared community structure can make sparse responses scientifically valuable, while very sparse or constant responses can still make species-specific parameters hard to estimate (Ovaskainen and Soininen 2011; Warton et al. 2015; Norberg et al. 2019). Documented prevalence thresholds are analysis choices, not general rules. Retain an important rare species only with a stated reason and a sensitivity analysis.
Psychometric item analysis uses the same principle: item support,
redundancy, response frequencies, and interpretability are evidence for
review rather than automatic deletion. mirt::itemstats() is
a useful precedent for inspecting items before fitting an IRT model
(Chalmers 2012). screen_gllvmTMB() likewise reports
response risks without silently changing the formula.
Post-fit loading salience, rotation, and latent-rank interpretation belong after fitting. They are not substitutes for this pre-fit response screen.
Next decision
Fit only after every WARN row has a documented keep,
recode, combine, or omit decision. Then run Can I trust this fit?; a clean pre-fit
screen does not make the fitted optimiser, Hessian, or weak-axis checks
pass automatically. If a scientifically essential response remains
constant or has too little denominator support, report that limitation
and defer its species-specific effect rather than forcing it into the
model.
References
- Albert A, Anderson JA (1984). On the existence of maximum likelihood estimates in logistic regression models. Biometrika 71:1–10.
- 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.
- 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.
- Mansournia MA, Geroldinger A, Greenland S, Heinze G (2018). Separation in logistic regression: causes, consequences, and control. American Journal of Epidemiology 187:864–870.
- Chalmers RP (2012). mirt: A multidimensional item response theory package for the R environment. Journal of Statistical Software 48(6):1–29.
- Chalmers RP.
mirt::itemstats()documentation: generic item summary statistics without prior IRT fitting. https://rdrr.io/cran/mirt/man/itemstats.html - Norberg A, Abrego N, Blanchet FG, Adler FR, Anderson BJ, Anttila J, Araujo MB, Dallas T, Dunson D, Elith J, Foster SD, Fox R, Franklin J, Godsoe W, Guisan A, O’Hara B, Hill NA, Holt RD, Hui FKC, Husby M, Kalas JA, Lehikoinen A, Luoto M, Mod HK, Newell G, Renner I, Roslin T, Soininen J, Thuiller W, Vanhatalo J, Warton D, White M, Zimmermann NE, Gravel D, Ovaskainen O (2019). A comprehensive evaluation of predictive performance of 33 species distribution models at species and community levels. Ecological Monographs 89:e01370. https://doi.org/10.1002/ecm.1370.
- Ovaskainen O, Soininen J (2011). Making more out of sparse data: hierarchical modeling of species communities. Ecology 92:289–295. https://doi.org/10.1890/10-1251.1.
- Warton DI, Blanchet FG, O’Hara RB, Ovaskainen O, Taskinen S, Walker SC, Hui FKC (2015). So many variables: joint modeling in community ecology. Trends in Ecology & Evolution 30:766–779. https://doi.org/10.1016/j.tree.2015.09.007.