check_collinearity() silently drops interaction terms for brmsfit models (no warning)

Disclaimer: this issue was written by Claude AI (I am not a programmer so I enlisted it's help)

Description

When calling check_collinearity() on a brmsfit model that includes interaction terms, the output silently omits every interaction term — no warning, no error, no indication anything is missing. Main effects and factor terms (e.g. grouping variables) are computed and reported correctly; only the interaction rows disappear.

Reproducible example

r
library(brms)
library(performance)

set.seed(1)
n <- 200
dat <- data.frame(
y = rbinom(n, 1, 0.5),
sex = rnorm(n),
x1 = rnorm(n),
x2 = rnorm(n),
grp = factor(sample(1:5, n, replace = TRUE))
)

fit <- brm(y ~ sex * (x1 + x2) + grp, data = dat, family = bernoulli(),
chains = 1, iter = 500, refresh = 0)

check_collinearity(fit)

formula(fit) and rownames(fixef(fit)) confirm the interaction terms (sex:x1, sex:x2) are genuinely present as population-level effects in the fitted model. But check_collinearity(fit) returns only rows for sex, x1, x2, grp — the two interaction terms are missing from the output entirely, with no message.

Root cause (traced through source)

In .check_collinearity(), .term_assignments() tries attr(insight::get_modelmatrix(x), "assign") first. For brmsfit, this doesn't return a usable assign attribute, so it falls back to .find_term_assignment().

That fallback builds its term lookup from insight::find_predictors(x), which returns unique variable names only (e.g. sex, x1, x2), not full formula term labels with interactions preserved (e.g. sex:x1). When it then does:

r
match(insight::clean_names(model_params), parms)

for a coefficient like "sex:x1", there's no matching entry in parms (built only from atomic variable names), so it returns NA for every interaction-term coefficient.

Downstream, .check_collinearity() loops over the real formula terms from stats::terms(f[[component]]) (which does include the interaction terms, e.g. 5 main terms + 2 interactions = 7 total). For term indices beyond the main effects, which(term_assign == term) finds nothing (since term_assign only ever contains values corresponding to the main-effect variable positions, or NA), so those terms get pushed into na_terms and silently removed from model_terms before the final output is built. This happens without triggering the interaction-VIF warning further down, since that check only fires when any(result > 10).

Suggested fix

Replace the fallback's term-assignment logic with an actual model.matrix()-derived assign attribute, which handles interactions correctly the same way it does for lm/glm:

r
.find_term_assignment_fixed <- function(x, component, verbose = TRUE) {
f <- insight::find_formula(x, verbose = FALSE)[[component]]
if (is.null(f)) return(NULL)

dat <- insight::get_data(x, verbose = FALSE)
mm <- tryCatch(
stats::model.matrix(stats::terms(f), data = dat),
error = function(e) NULL
)
if (is.null(mm)) return(NULL)

assign_vec <- attr(mm, "assign")
coef_names <- insight::clean_names(colnames(mm))
model_params <- insight::find_parameters(x)[[component]]

idx <- match(insight::clean_names(model_params), coef_names)
assign_vec[idx]
}

I verified this fix (via assignInNamespace()) against an independent glm()-based VIF check on the same design matrix — the patched check_collinearity(fit) output matches the glm()/car::vif() reference values closely, confirming the corrected term-assignment produces accurate GVIFs for all terms including interactions.

Impact

This is a silent failure — the function returns a plausible-looking, correctly-formatted table with no error or warning, so users have no indication interaction terms were dropped. For any brmsfit with interaction terms, reported "no multicollinearity" conclusions based on this output would be incomplete rather than wrong per se, but anyone specifically checking interaction-term VIFs (a common recommendation, since interactions often show inflated VIF) gets no information about them at all.

r

fit
Family: bernoulli
Links: mu = logit
Formula: y ~ sex * (x1 + x2) + grp
Data: dat (Number of observations: 200)
Draws: 1 chains, each with iter = 500; warmup = 250; thin = 1;
total post-warmup draws = 250

Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.32 0.33 -0.21 1.06 1.00 222 133
sex -0.12 0.17 -0.47 0.18 1.01 144 192
x1 -0.08 0.15 -0.35 0.23 1.00 290 163
x2 0.08 0.15 -0.21 0.33 1.00 595 165
grp2 -0.63 0.49 -1.68 0.22 1.00 187 153
grp3 -0.31 0.46 -1.22 0.51 1.00 197 167
grp4 -0.36 0.47 -1.22 0.49 1.00 261 184
grp5 -0.20 0.42 -1.02 0.49 1.00 215 221
sex:x1 -0.22 0.15 -0.49 0.07 1.00 182 150
sex:x2 0.03 0.14 -0.25 0.29 1.00 291 202

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

check_collinearity(fit)

Check for Multicollinearity

Low Correlation

Term VIF VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
sex 1.23 [1.10, 1.53] 1.11 0.81 [0.66, 0.91]
x1 1.12 [1.03, 1.46] 1.06 0.89 [0.69, 0.97]
x2 1.07 [1.01, 1.62] 1.03 0.94 [0.62, 0.99]
grp 1.20 [1.08, 1.49] 1.02 0.84 [0.67, 0.93]

Start with .check_collinearity(), .term_assignments(), and the .find_term_assignment() fallback described in the issue. Run the provided brmsfit reproducible example and inspect how interaction coefficients are assigned to formula terms. Done means check_collinearity() reports sex:x1 and sex:x2 instead of silently omitting them, with results matching the independent glm()/car::vif() reference.