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]
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
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 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]