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 result from the default heuristic screen therefore does not prove that separation cannot occur; request the formal certificate when the fixed-design question matters.
Set screen_control(separation = "fixed") to run the
optional formal fixed-design certificate. It classifies
coefficient-connected blocks as overlap,
complete, or quasi_complete for logit, probit,
and complementary-log-log links after fixed-effect maps, ties, finite
offsets, and structural-zero constraints are resolved. It requires the
optional detectseparation package and never changes the
fitted estimator. The certificate does not inspect latent loadings or
covariance parameters. Finite known offsets are accepted for this
diagnostic certificate, but gllvmTMB() currently rejects
nonzero offsets for binomial fits; diagnostic support does not imply
fitting support.
When a scientifically necessary binary trait is separated, the related LA-MSPL worked example shows the opt-in finite point estimator. Filtering remains a scientific/data-quality decision, not the screen’s automatic response.
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 8 | WARN 1 | PASS 13
#> Inspect 9 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
#> 1 trait WARN inspect
#> 2 trait FAIL exclude_from_latent_block
#> 3 pair FAIL collapse_or_recode
#> 4 pair FAIL collapse_or_recode
#> 5 pair FAIL collapse_or_recode
#> 6 design FAIL inspect
#> 7 response_dependency FAIL collapse_or_recode
#> 8 response_dependency FAIL collapse_or_recode
#> 9 response_dependency FAIL collapse_or_recode
#> trait
#> 1 region_north_america
#> 2 level_ecosystem
#> 3 outcome_richness / outcome_diversity
#> 4 outcome_richness / outcome_no_change
#> 5 outcome_diversity / outcome_no_change
#> 6 <NA>
#> 7 level_ecosystem
#> 8 outcome_richness, outcome_diversity
#> 9 outcome_richness, 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 rows
#> 6 response matrix augmented with an intercept is rank deficient by 3 on 40 complete rows: at least one exact affine dependency exists among the screened traits
#> 7 level_ecosystem = 0: this affine dependency is already reported at trait/pair level (a constant trait) and was projected out before the one-hot certificate search
#> 8 outcome_diversity = outcome_richness: this affine dependency is already reported at trait/pair level (a duplicate pair) and was projected out before the one-hot certificate search
#> 9 outcome_richness + outcome_no_change = 1: this affine dependency is already reported at trait/pair level (a complement pair) and was projected out before the one-hot certificate searchThe 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
#> 4 response_affine_rank FAIL 3 6
#> 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
#> 4 response matrix augmented with an intercept is rank deficient by 3 on 40 complete rows: at least one exact affine dependency exists among the screened 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.
Exact response-side dependencies
The pairwise check above only sees two-way exact relations: duplicates and complements. It cannot see a higher-order exact relation among three or more traits, and systematic-map data commonly contain exactly that: a review-type, geographic-scope, or temporal-scope variable recorded as several dummy columns rather than one categorical column, so those columns sum to exactly 1 on every review.
screen_gllvmTMB() checks for this with an affine-rank
check on the response matrix (augmented with an intercept column),
restricted to reviews with every screened trait observed. A rank
deficiency means at least one exact affine dependency exists among the
traits; where it decomposes into a clean one-hot block, a human-readable
certificate is reported. This is a best-effort search, not exhaustive
minimal-subset discovery: with more than one such block at once it can
fail to decompose them cleanly (see ?screen_gllvmTMB), so
declaring known structure with known_groups is the more
reliable route whenever it is already known.
scope <- data.frame(
review = factor(seq_len(30)),
scope_global = rep(c(1, 0, 0), 10),
scope_regional = rep(c(0, 1, 0), 10),
scope_local = rep(c(0, 0, 1), 10),
outcome_richness = rep(c(0, 1), 15)
)scope_global, scope_regional, and
scope_local are the dummy coding of one geographic-scope
variable: exactly one of the three is 1 for every review.
scr_scope <- suppressWarnings(screen_gllvmTMB(
traits(scope_global, scope_regional, scope_local, outcome_richness) ~
1 + latent(1 | review, d = 2),
data = scope,
unit = "review",
family = binomial("probit")
))
design_scope <- screen_table(scr_scope, "design")
design_scope[
design_scope$component == "response_affine_rank",
c("component", "status", "value", "threshold", "message")
]
#> component status value threshold
#> 4 response_affine_rank FAIL 4 5
#> message
#> 4 response matrix augmented with an intercept is rank deficient by 1 on 30 complete rows: at least one exact affine dependency exists among the screened traits
screen_table(scr_scope, "response_dependencies")[, c(
"scope", "type", "status", "traits", "certificate"
)]
#> scope type status traits
#> 1 affine one_hot_block FAIL scope_global, scope_regional, scope_local
#> certificate
#> 1 scope_global + scope_regional + scope_local = 1The design table’s response_affine_rank row
reports the rank deficiency; the response_dependencies
table adds the certificate,
"scope_global + scope_regional + scope_local = 1", and
names exactly those three traits — outcome_richness is
uninvolved.
When the grouping is already known, known_groups gives
the same certificate as an exact, deterministic check rather than a
best-effort search, and also covers a declared nesting/containment chain
(for example a broad realm indicator together with a narrower nested
realm):
scr_known <- screen_gllvmTMB(
traits(scope_global, scope_regional, scope_local, outcome_richness) ~
1 + latent(1 | review, d = 2),
data = scope,
unit = "review",
family = binomial("probit"),
known_groups = list(
geographic_scope = c("scope_global", "scope_regional", "scope_local")
)
)
deps_known <- screen_table(scr_known, "response_dependencies")
deps_known[
deps_known$scope == "known_group",
c("group", "type", "status", "certificate")
]
#> group type status
#> 2 geographic_scope known_one_hot FAIL
#> certificate
#> 2 scope_global + scope_regional + scope_local = 1Neither check removes or recodes the block automatically. An exact
one-hot block is one categorical variable: recode it as multinomial() or fix a
deliberate reference level, rather than fitting it as several
unconstrained binary traits in the same latent block.
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 = screen_control(separation = "fixed")
)
screen_table(scr_candidate, "recommendations")[, c(
"scope", "status", "action", "trait", "evidence"
)]
#> scope status action trait
#> 1 trait WARN inspect region_north_america
#> 2 separation NOT_CHECKED install_or_disable_separation outcome_richness
#> 3 separation NOT_CHECKED install_or_disable_separation region_north_america
#> evidence
#> 1 the minority outcome has very few observed trials
#> 2 install detectseparation (>= 0.4.0), or set separation = "none"
#> 3 install detectseparation (>= 0.4.0), or set separation = "none"
screen_table(scr_candidate, "separation")[, c(
"traits", "status", "severity", "n_columns_active", "action"
)]
#> traits status severity n_columns_active
#> 1 outcome_richness NOT_CHECKED dependency_missing 1
#> 2 region_north_america NOT_CHECKED dependency_missing 1
#> action
#> 1 install_or_disable_separation
#> 2 install_or_disable_separationThis rerun contains only the deliberately retained rarity warning from the heuristic checks and now also applies the formal fixed-design certificate to the final candidate formula. After fixed-design overlap, use ML unless another scientific reason calls for regularisation. When complete or quasi-complete separation affects a scientifically necessary trait and the model lies inside the admitted complete-Bernoulli surface, consider opt-in MSPL. Do not use MSPL to rescue a constant response. If the model instead contains missing response cells, grouped binomial trials, mixed families, or another unsupported structure, redesign or filter transparently or await the corresponding extension. 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. The default heuristic screen does not replace a
formal separation analysis; request separation = "fixed"
when that formal fixed-design certificate is needed.
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.
- Kosmidis I, Firth D (2021). Jeffreys-prior penalty, finiteness and shrinkage in binomial-response generalized linear models. Biometrika 108:71–82. https://doi.org/10.1093/biomet/asaa052.
- Sterzinger P, Kosmidis I (2023). Maximum softly-penalized likelihood for mixed effects logistic regression. Statistics and Computing 33:53. https://doi.org/10.1007/s11222-023-10217-3.
- 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.