Fix term matching in check_collinearity() for clm scale and nominal effects, more vglm families, and backticked multinomial names

NEWS.md

@@ -22,6 +22,21 @@
   values of the smaller model to the length of the larger one. Now all LOOIC
   weights are `NA`, with a warning.
 
+* `check_collinearity()` now uses only the location slopes of `ordinal::clm()`
+  models with scale effects. Before, the VIF of a location slope also used the
+  scale effect with the same name. For `clm()` models with scale or nominal
+  effects, a message names the effects that the VIFs leave out.
+
+* `check_collinearity()` now works for `VGAM::vglm()` models with the
+  `sratio()`, `cratio()`, `acat()`, and `poissonff()` families. Before, it
+  failed with "'list' object cannot be coerced to type 'integer'" or gave VIFs
+  for the wrong terms.
+
+* `check_collinearity()` now works for multinomial models from *nnet*,
+  *brglm2*, and *mclogit* with backticked variable names, such as
+  `` `my x` ``. Before, it returned `NULL` with "Could not extract model
+  terms".
+
 * `check_collinearity()` now returns correct VIFs for the zero-inflation
   component of *glmmTMB* models whose zero-inflation formula differs from the
   conditional formula. Before, the zero-inflation terms used the term blocks of

R/check_collinearity.R

@@ -41,12 +41,12 @@
 #'
 #' For multinomial models from `nnet::multinom()`, `mclogit::mblogit()`,
 #' `brglm2::brmultinom()` and `brglm2::bracl()`, and for `VGAM::vglm()` models
-#' with the multinomial or cumulative family, each term can have its own set
-#' of slopes for each outcome level except one.
+#' with the multinomial, cumulative, sratio, cratio, or acat family, each term
+#' can have its own set of slopes for each outcome level except one.
 #' These models get one VIF per term, computed from the slopes of that term
 #' across all outcome levels, so the VIF does not depend on the reference
-#' level. Parallel models, such as `bracl(parallel = TRUE)`, have one slope
-#' per model matrix column. The adjusted VIF is then `VIF^(1/(2*<n>))`, where
+#' level. Parallel models, such as `bracl(parallel = TRUE)`, and `vglm()`
+#' models with the poissonff family have one slope per model matrix column. The adjusted VIF is then `VIF^(1/(2*<n>))`, where
 #' `<n>` is the number of slopes of the term across all outcome levels.
 #'
 #' The confidence interval of a VIF is an interval for the population R2 of the
@@ -534,7 +534,7 @@
 }
 
 
 .check_collinearity <- function(x, component, ci = 0.95, verbose = TRUE) {
   v <- .safe(insight::get_varcov(x, component = component, verbose = FALSE))
 
   # fix class for fixest, which returns a "fixest_vcov" here
@@ -575,10 +575,23 @@
 
   # ordinal models have one threshold for each outcome level except one, in
   # place of one intercept. polr, svyolr and stan_polr put the thresholds after
-  # the slopes, clm and clm2 put them first, so keep the slopes by name
+  # the slopes, clm2 puts them first, so keep the slopes by name
   if (.is_threshold_model(x)) {
     slopes <- .threshold_slopes(x, component)
-    keep_idx <- colnames(v) %in% slopes$names
+    if (inherits(x, "clm") && is.logical(x$aliased$alpha)) {
+      # clm puts the thresholds first, including those of nominal effects,
+      # then the location slopes without aliased ones, then the scale effects.
+      # A scale effect can have the same name as a location slope, so keep
+      # the slopes by position
+      n_alpha <- sum(!x$aliased$alpha)
+      keep_idx <- seq_len(ncol(v)) %in% (n_alpha + seq_along(slopes$names))
+      # no message when there are too few terms for VIFs
+      if (isTRUE(verbose) && length(unique(stats::na.omit(slopes$assign))) > 1) {
+        .alert_clm_effects(x)
+      }
+    } else {
+      keep_idx <- colnames(v) %in% slopes$names
+    }
     term_assign <- slopes$assign[match(colnames(v)[keep_idx], slopes$names)]
     if (!length(term_assign) || anyNA(term_assign)) {
       if (verbose) {
@@ -606,9 +619,16 @@
     } else {
       mm <- .safe(insight::get_modelmatrix(x))
     }
-    coef_names <- .multinom_coef_names(x, colnames(v), colnames(mm))
+    # the column names of the variance-covariance matrix from insight have no
+    # backticks, but coefficient and model matrix names can keep them
+    mm_names <- gsub("`", "", colnames(mm), fixed = TRUE)
+    coef_names <- .multinom_coef_names(
+      x,
+      gsub("`", "", colnames(v), fixed = TRUE),
+      mm_names
+    )
     keep_idx <- !is.na(coef_names) & coef_names != "(Intercept)"
-    term_assign <- attr(mm, "assign")[match(coef_names[keep_idx], colnames(mm))]
+    term_assign <- attr(mm, "assign")[match(coef_names[keep_idx], mm_names)]
     # a coefficient without a term would give VIFs for the wrong terms
     if (!length(term_assign) || anyNA(term_assign) || anyNA(coef_names)) {
       if (verbose) {
@@ -688,6 +708,7 @@
 
   # we have rank-deficiency here. remove NA columns from assignment
   if (isTRUE(attr(v, "rank_deficient")) || anyNA(v)) {
+    na_cols <- NULL
     if (!is.null(attr(v, "na_columns_index"))) {
       # If this attribute exists, then NA values were already removed from v
       term_assign <- term_assign[-attr(v, "na_columns_index")]
@@ -759,7 +780,7 @@
     return(NULL)
   }
 
   R <- stats::cov2cor(v)
   detR <- det(R)
 
   result <- vector("numeric")
@@ -859,7 +880,7 @@
 }
 
 
 .vif_ci <- function(x, result, ci, n_coef) {
   # CIs for the population R2 of each term's auxiliary regression, which is
   # 1 - 1 / VIF. The F statistic of the auxiliary regression follows a
   # noncentral F distribution with k and n - k - 1 df, where k is the number
@@ -1158,9 +1179,9 @@
 # to, and the model matrix columns the fit dropped as aliased
 .threshold_slopes <- function(x, component) {
   if (inherits(x, c("clm", "clmm", "clm2"))) {
     beta <- x$beta
   } else {
     beta <- x$coefficients
   }
   # clm keeps aliased slopes as NA, polr and clmm drop them. Coefficient names
   # can keep backticks, which the column names of the variance-covariance
@@ -1190,20 +1211,49 @@
   used <- colSums(abs(mm), na.rm = TRUE) > 0
   mm_names <- setdiff(colnames(mm)[used], "(Intercept)")
 
+  # clm drops a location term that is also a nominal effect, which is not
+  # collinearity
+  aliased <- setdiff(mm_names, slope_names)
+  if (inherits(x, "clm")) {
+    nominal <- .safe(labels(stats::terms(x$formulas$nominal)))
+    aliased_terms <- labels(stats::terms(f))[
+      attr(mm, "assign")[match(aliased, colnames(mm))]
+    ]
+    aliased <- aliased[!aliased_terms %in% nominal]
+  }
+
   list(
     names = slope_names,
     assign = attr(mm, "assign")[match(slope_names, colnames(mm))],
-    aliased = setdiff(mm_names, slope_names)
+    aliased = aliased
   )
 }
 
+# clm models can have scale and nominal effects, which are not slopes of the
+# linear predictor, so the VIFs leave them out
+.alert_clm_effects <- function(x) {
+  left_out <- NULL
+  for (i in c("scale", "nominal")) {
+    effects <- .safe(labels(stats::terms(x$formulas[[i]])))
+    if (length(effects)) {
+      left_out <- c(left_out, sprintf("%s effects (%s)", i, toString(effects)))
+    }
+  }
+  if (length(left_out)) {
+    insight::format_alert(
+      "VIFs are computed for the location effects only.",
+      sprintf("Left out: %s.", paste(left_out, collapse = " and "))
+    )
+  }
+}
+
 .warn_rank_deficient <- function(na_cols = NULL) {
   if (length(na_cols) > 0) {
     insight::format_warning(
       "Model matrix is rank deficient. VIFs may not be sensible.",
       paste0(
         "The following coefficients have VIF = Inf / tolerance = 0: ",
         paste0(na_cols, collapse = ", ")
       )
     )
   } else {
@@ -1213,11 +1263,15 @@
   }
 }
 
-# VGAM models with the multinomial or cumulative family
+# VGAM models with the multinomial, cumulative, sratio, cratio, acat, or
+# poissonff family, whose model matrix has an "assign" list
 .is_vglm_categorical <- function(x) {
   inherits(x, "vglm") &&
     !inherits(x, "vgam") &&
-    isTRUE(.safe(x@family@vfamily[1]) %in% c("multinomial", "cumulative"))
+    isTRUE(
+      .safe(x@family@vfamily[1]) %in%
+        c("multinomial", "cumulative", "sratio", "cratio", "acat", "poissonff")
+    )
 }
 
 # coefficient names of multinomial models are prefixed with the outcome level,
@@ -1255,6 +1309,7 @@
   if (!is.matrix(coef_mat)) {
     return(coef_names)
   }
+  colnames(coef_mat) <- gsub("`", "", colnames(coef_mat), fixed = TRUE)
   prefixed <- outer(rownames(coef_mat), colnames(coef_mat), paste, sep = sep)
   # level and coefficient names can combine into the same name, for example
   # level "b" with "x1:x3" and level "b:x1" with "x3"

tests/testthat/test-check_collinearity.R

@@ -121,7 +121,7 @@
     .package = "insight"
   )
   testthat::local_mocked_bindings(
     model.matrix = function(...) stop("rebuild failed"),
     .package = "stats"
   )
   expect_null(
@@ -130,6 +130,23 @@
 })
 
 
+test_that("check_collinearity, rank-deficient flag without NA columns", {
+  m <- lm(mpg ~ wt + hp + disp, data = mtcars)
+  v <- stats::vcov(m)
+  expected <- check_collinearity(m)$VIF
+  # the flag is set, but the matrix has no NA and no index of NA columns
+  testthat::local_mocked_bindings(
+    get_varcov = function(...) structure(v, rank_deficient = TRUE),
+    .package = "insight"
+  )
+  expect_warning(
+    out <- check_collinearity(m, verbose = TRUE),
+    "Model matrix is rank deficient"
+  )
+  expect_equal(out$VIF, expected, tolerance = 1e-8)
+})
+
+
 test_that("check_collinearity", {
   skip_if_not_installed("glmmTMB")
   skip_if_not(getRversion() >= "4.0.0")
@@ -607,6 +624,90 @@
   expect_equal(out$VIF, c(1.11, 1.11), tolerance = 0.05)
 })
 
+# generalized VIF of each block of columns of `v`, from the correlation matrix
+# of all columns in `blocks`
+.gvif_blocks <- function(v, blocks) {
+  cols <- unlist(blocks)
+  cor_mat <- stats::cov2cor(v[cols, cols])
+  pos <- split(seq_along(cols), rep(seq_along(blocks), lengths(blocks)))
+  vapply(
+    pos,
+    function(idx) {
+      det(cor_mat[idx, idx, drop = FALSE]) *
+        det(cor_mat[-idx, -idx, drop = FALSE]) /
+        det(cor_mat)
+    },
+    numeric(1),
+    USE.NAMES = FALSE
+  )
+}
+
+.clm_data <- function(n = 400) {
+  set.seed(42)
+  d <- data.frame(x1 = rnorm(n), x3 = rnorm(n), x4 = rnorm(n))
+  d$x2 <- d$x1 + rnorm(n, sd = 0.7)
+  latent <- 0.5 * d$x1 + 0.4 * d$x2 - 0.3 * d$x3 + stats::rlogis(n)
+  d$y <- cut(
+    latent,
+    stats::quantile(latent, c(0, 0.25, 0.5, 0.75, 1)),
+    include.lowest = TRUE,
+    labels = c("a", "b", "c", "d"),
+    ordered_result = TRUE
+  )
+  d
+}
+
+test_that("check_collinearity, ordinal clm models with scale effects", {
+  skip_if_not_installed("ordinal")
+  d <- .clm_data()
+  m <- ordinal::clm(y ~ x1 + x2 + x3, scale = ~x1, data = d)
+  expect_message(
+    out <- check_collinearity(m),
+    "scale effects \\(x1\\)"
+  )
+  expect_identical(out$Term, c("x1", "x2", "x3"))
+  # the location slopes follow the thresholds; the scale effect of x1 has a
+  # column with the same name after them
+  slope_cols <- length(m$alpha) + seq_along(m$beta)
+  expect_identical(colnames(stats::vcov(m))[slope_cols], names(m$beta))
+  expect_identical(which(colnames(stats::vcov(m)) == "x1")[1], slope_cols[1])
+  expected <- .gvif_blocks(stats::vcov(m), as.list(slope_cols))
+  expect_equal(out$VIF, expected, tolerance = 1e-8)
+
+  # no message when one location term leaves too few terms for VIFs
+  m2 <- ordinal::clm(y ~ x1, scale = ~x2, data = d)
+  expect_no_message(
+    check_collinearity(m2),
+    message = "location effects only"
+  )
+})
+
+test_that("check_collinearity, ordinal clm models with nominal effects", {
+  skip_if_not_installed("ordinal")
+  d <- .clm_data()
+  m <- ordinal::clm(y ~ x1 + x2 + x4, nominal = ~x3, data = d)
+  expect_message(
+    out <- check_collinearity(m),
+    "nominal effects \\(x3\\)"
+  )
+  expect_identical(out$Term, c("x1", "x2", "x4"))
+  slope_cols <- length(m$alpha) + seq_along(m$beta)
+  expected <- .gvif_blocks(stats::vcov(m), as.list(slope_cols))
+  expect_equal(out$VIF, expected, tolerance = 1e-8)
+
+  # no message without scale or nominal effects
+  m2 <- ordinal::clm(y ~ x1 + x2 + x4, data = d)
+  expect_no_message(check_collinearity(m2))
+
+  # a location term that is also a nominal effect is dropped by clm, which is
+  # not rank deficiency
+  m3 <- ordinal::clm(y ~ x1 + x2 + x3, nominal = ~x3, data = d)
+  expect_no_warning(
+    out3 <- suppressMessages(check_collinearity(m3))
+  )
+  expect_identical(out3$Term, c("x1", "x2"))
+})
+
 test_that("check_collinearity, standard lm models with offset", {
   # Standard linear model with an offset
   m_lm_offset <- lm(mpg ~ wt + cyl + offset(disp), data = mtcars)
@@ -949,6 +1050,101 @@
   expect_equal(out2$SE_factor, out2$VIF^(1 / (2 * m_coef2)), tolerance = 1e-8)
 })
 
+test_that("check_collinearity, VGAM::vglm sratio, acat, cratio, and poissonff", {
+  skip_if_not_installed("VGAM")
+  d <- .multinom_data()
+  d$o <- cut(
+    0.5 * d$x1 - 0.3 * d$x3 + 0.2 * (d$f == "q") + stats::rlogis(nrow(d)),
+    c(-Inf, -0.5, 0.5, Inf),
+    labels = c("a", "b", "c"),
+    ordered_result = TRUE
+  )
+  d$count <- stats::rpois(nrow(d), exp(0.3 * d$x1 - 0.2 * d$x3))
+  fits <- list(
+    sratio = VGAM::vglm(o ~ x1 + x2 + f, family = VGAM::sratio(), data = d),
+    acat = VGAM::vglm(o ~ x1 + x2 + f, family = VGAM::acat(), data = d),
+    cratio = VGAM::vglm(o ~ x1 + x2 + f, family = VGAM::cratio(), data = d),
+    poissonff = VGAM::vglm(count ~ x1 + x2 + f, family = VGAM::poissonff(), data = d),
+    sratio_parallel = VGAM::vglm(
+      o ~ x1 + x2 + f,
+      family = VGAM::sratio(parallel = TRUE),
+      data = d
+    )
+  )
+  for (i in names(fits)) {
+    m <- fits[[i]]
+    # the "assign" list maps each term to its columns of the
+    # variance-covariance matrix
+    blocks <- attr(stats::model.matrix(m), "assign")
+    blocks <- blocks[names(blocks) != "(Intercept)"]
+    out <- check_collinearity(m)
+    expect_identical(out$Term, c("x1", "x2", "f"), info = i)
+    expected <- .gvif_blocks(stats::vcov(m), blocks)
+    expect_equal(out$VIF, expected, tolerance = 1e-8, info = i)
+    n_cols <- unname(lengths(blocks))
+    expect_equal(out$SE_factor, out$VIF^(1 / (2 * n_cols)), tolerance = 1e-8, info = i)
+  }
+
+  # the poissonff fit is the same model as a poisson glm
+  m_glm <- stats::glm(count ~ x1 + x2 + f, family = stats::poisson(), data = d)
+  expect_equal(
+    check_collinearity(fits$poissonff)$VIF,
+    check_collinearity(m_glm)$VIF,
+    tolerance = 1e-6
+  )
+})
+
+test_that("check_collinearity, multinomial models with backticked names", {
+  skip_if_not_installed("nnet")
+  skip_if_not_installed("brglm2")
+  d <- .multinom_data()
+  d$`my x` <- d$x1
+  d$my_x <- d$x1
+  fits <- list(
+    nnet = nnet::multinom(y ~ `my x` + x2 + f + `my x`:f, data = d, trace = FALSE),
+    brmultinom = brglm2::brmultinom(y ~ `my x` + x2 + f + `my x`:f, data = d),
+    bracl = brglm2::bracl(y ~ `my x` + x2 + f + `my x`:f, data = d),
+    bracl_parallel = brglm2::bracl(
+      y ~ `my x` + x2 + f + `my x`:f,
+      data = d,
+      parallel = TRUE
+    )
+  )
+  fits_renamed <- list(
+    nnet = nnet::multinom(y ~ my_x + x2 + f + my_x:f, data = d, trace = FALSE),
+    brmultinom = brglm2::brmultinom(y ~ my_x + x2 + f + my_x:f, data = d),
+    bracl = brglm2::bracl(y ~ my_x + x2 + f + my_x:f, data = d),
+    bracl_parallel = brglm2::bracl(
+      y ~ my_x + x2 + f + my_x:f,
+      data = d,
+      parallel = TRUE
+    )
+  )
+  for (i in names(fits)) {
+    out <- suppressMessages(check_collinearity(fits[[i]]))
+    out_renamed <- suppressMessages(check_collinearity(fits_renamed[[i]]))
+    expect_identical(out$Term, c("my x", "x2", "f", "my x:f"), info = i)
+    expect_equal(out$VIF, out_renamed$VIF, tolerance = 1e-8, info = i)
+  }
+})
+
+test_that("check_collinearity, mclogit::mblogit with backticked names", {
+  skip_if_not_installed("mclogit")
+  d <- .multinom_data()
+  d$`my x` <- d$x1
+  d$my_x <- d$x1
+  m <- suppressWarnings(
+    mclogit::mblogit(y ~ `my x` + x2 + f + `my x`:f, data = d, trace = FALSE)
+  )
+  m_renamed <- suppressWarnings(
+    mclogit::mblogit(y ~ my_x + x2 + f + my_x:f, data = d, trace = FALSE)
+  )
+  out <- suppressMessages(check_collinearity(m))
+  out_renamed <- suppressMessages(check_collinearity(m_renamed))
+  expect_identical(out$Term, c("my x", "x2", "f", "my x:f"))
+  expect_equal(out$VIF, out_renamed$VIF, tolerance = 1e-8)
+})
+
 test_that("check_collinearity, multinomial models without intercept", {
   skip_if_not_installed("nnet")
   skip_if_not_installed("mclogit")