check_collinearity() returns artificially inflated VIFs for ordinal::clmm models

[jmgirard]

I have noticed an issue where performance::check_collinearity() seems to miscalculate or heavily inflate Variance Inflation Factors (VIFs) when applied to Cumulative Link Mixed Models (clmm objects from the ordinal package).

When passing perfectly orthogonal, simulated predictors into a clmm, the resulting VIFs are flagged as highly collinear. Fitting a standard linear mixed model (lmer) to the exact same predictors correctly returns VIFs near 1.0.

I suspect the VIF calculation algorithm might be accidentally catching the multiple ordinal threshold estimates (the internal intercepts) from the clmm variance-covariance matrix, which artificially explodes the collinearity penalty for the fixed predictors.

Reproducible Example:

library(ordinal)
library(lme4)
#> Loading required package: Matrix
library(performance)

set.seed(999)
n <- 500

# 1. Simulate perfectly orthogonal predictors
x_continuous <- rnorm(n, mean = 0, sd = 1)
x_binary <- sample(c(-0.5, 0.5), size = n, replace = TRUE, prob = c(0.85, 0.15)) 
subject_id <- factor(rep(1:50, each = 10))

# 2. Generate an ordinal outcome with MANY categories
random_intercepts <- rnorm(50, 0, 1)
latent_y <- 2 * x_continuous + 3 * x_binary + random_intercepts[as.numeric(subject_id)] + rlogis(n)

# Cut into 15 categories to generate 14 distinct thresholds
y_ordinal <- cut(
  latent_y, 
  breaks = 15, 
  ordered_result = TRUE
)

dat <- data.frame(y_ordinal, x_continuous, x_binary, subject_id)

# 3. Fit models
mod_lmer <- lmer(as.numeric(y_ordinal) ~ x_continuous + x_binary + (1 | subject_id), data = dat)
mod_clmm <- clmm(y_ordinal ~ x_continuous + x_binary + (1 | subject_id), data = dat)

# 4. Compare Collinearity Checks
check_collinearity(mod_lmer)
#> # Check for Multicollinearity
#> 
#> Low Correlation
#> 
#>          Term  VIF  VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
#>  x_continuous 1.00 [1.00, Inf]     1.00      1.00     [0.00, 1.00]
#>      x_binary 1.00 [1.00, Inf]     1.00      1.00     [0.00, 1.00]

check_collinearity(mod_clmm)
#> # Check for Multicollinearity
#> 
#> Low Correlation
#> 
#>          Term  VIF   VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
#>  x_continuous 1.79 [1.60, 2.04]     1.34      0.56     [0.49, 0.63]
#>      x_binary 4.36 [3.77, 5.08]     2.09      0.23     [0.20, 0.27]

^{Created on 2026-04-21 with reprex v2.1.1}

Expected Behavior:

The VIF values for the fixed effects in the clmm object should roughly match the VIF values of the lmer object, as the fixed-effect design matrix is identical.

System Information:

R version 4.5.3 (2026-03-11 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] performance_0.16.0 lme4_2.0-1         Matrix_1.7-5      
[4] ordinal_2025.12-29

loaded via a namespace (and not attached):
 [1] numDeriv_2016.8-1.1 lattice_0.22-9      splines_4.5.3      
 [4] ucminf_1.2.3        cli_3.6.6           Rdpack_2.6.6       
 [7] nloptr_2.2.1        grid_4.5.3          reformulas_0.4.4   
[10] compiler_4.5.3      boot_1.3-32         rbibutils_2.4.1    
[13] rstudioapi_0.18.0   tools_4.5.3         nlme_3.1-169       
[16] minqa_1.2.8         Rcpp_1.1.1-1        rlang_1.2.0        
[19] MASS_7.3-65         insight_1.5.0      

[mattansb]

I believe the current implementation of the VIF computation is based on the vcov matrix and not the correlation matrix of the design matrix. So different models have different standard errors and would then produce different VIFs under this method.

We can see that for linear fixed effects models (lm) the vcov method is the same as what we would get from the model matrix covariance, but not so for non-linear (generalized) methods:

library(performance)

m1 <- lm(am ~ wt + cyl + gear, data = mtcars)
m2 <- glm(am ~ wt + cyl + gear, data = mtcars, family = binomial)
#> Warning: glm.fit: algorithm did not converge
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

# hack to make vcov (what check_collinearity uses) return the covariance matrix of the model matrix:
m3 <- m2
class(m3) <- c("lmlmlm", class(m3))
vcov.lmlmlm <- function(object, ...) {
  cov(model.matrix(object))
}

check_collinearity(m1)
#> # Check for Multicollinearity
#> 
#> Low Correlation
#> 
#>  Term  VIF   VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
#>    wt 2.98 [1.98, 4.98]     1.73      0.34     [0.20, 0.50]
#>   cyl 2.59 [1.76, 4.32]     1.61      0.39     [0.23, 0.57]
#>  gear 1.52 [1.18, 2.56]     1.23      0.66     [0.39, 0.85]
check_collinearity(m2)
#> # Check for Multicollinearity
#> 
#> Low Correlation
#> 
#>  Term  VIF   VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
#>    wt 3.03 [2.01, 5.08]     1.74      0.33     [0.20, 0.50]
#>   cyl 1.98 [1.42, 3.29]     1.41      0.50     [0.30, 0.70]
#>  gear 1.86 [1.35, 3.08]     1.36      0.54     [0.32, 0.74]
check_collinearity(m3)
#> # Check for Multicollinearity
#> 
#> Low Correlation
#> 
#>  Term  VIF   VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
#>    wt 2.98 [1.98, 4.98]     1.73      0.34     [0.20, 0.50]
#>   cyl 2.59 [1.76, 4.32]     1.61      0.39     [0.23, 0.57]
#>  gear 1.52 [1.18, 2.56]     1.23      0.66     [0.39, 0.85]

^{Created on 2026-04-22 with reprex v2.1.1}

Would it make sense to switch to using the model matrix covariance instead? I think so...

[jmgirard]

@mattansb raises an interesting broader question here about using vcov() versus cov(model.matrix).

But in case you all decide to stick with the current vcov() approach (to keep the Fox & Monette GVIF behavior), I dug into why clmm models specifically blow up here.

The issue is that check_collinearity() pulls the full variance-covariance matrix. For ordinal::clmm models, that matrix includes all the threshold estimates (which act as multiple intercepts) and the random effect variances. Leaving these structural parameters in the matrix artificially inflates the VIFs for the actual fixed predictors.

To fix this within the current architecture, we need to subset the matrix down to the true slope parameters. We just have to make sure we also subset the internal term_assign tracking vector at the same time so we don't break the matrix indexing later on.

In check_collinearity.R, we could replace the current missing-intercept check with something like this:

  # Filter to true slope parameters (handles multiple intercepts in ordinal models)
  if (inherits(x, c("clm", "clmm"))) {
    slope_names <- names(x$beta)
    keep_idx <- which(colnames(v) %in% slope_names)
  } else if (insight::has_intercept(x)) {
    # Standard behavior: drop the first column/row (the singular intercept)
    keep_idx <- seq_len(ncol(v))[-1]
  } else {
    keep_idx <- seq_len(ncol(v))
    if (isTRUE(verbose)) {
      insight::format_alert("Model without intercept. VIFs may not be sensible.")
    }
  }
  
  # Safely subset both the matrix and the assignment vector
  if (length(keep_idx) < ncol(v)) {
    v <- v[keep_idx, keep_idx, drop = FALSE]
    term_assign <- term_assign[keep_idx]
  }

When I tested this against the reproducible example in the original post, the calculated VIFs dropped from the artificially inflated 4.36 down to 1.12. This 1.12 correctly reflects the minor covariance induced between the parameter estimates within the ordinal likelihood space (unlike the raw design matrix, where they are perfectly orthogonal).

Let me know how you all want to handle the larger vcov vs model.matrix design choice, but I'm happy to put together a PR for this specific clmm matrix subsetting if you want to keep the current approach!

[bwiernik]

I think it makes sense to keep using vcov() and subset to only include fixed effects

[mattansb]

Let's keep the vcov behavior - happy to review your PR (: