diff --git a/NEWS.md b/NEWS.md index 49aa04b67..8aa38ff7f 100644 --- a/NEWS.md +++ b/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 diff --git a/R/check_collinearity.R b/R/check_collinearity.R index aae771456..2ed67d16e 100644 --- a/R/check_collinearity.R +++ b/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*))`, 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*))`, where #' `` 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 @@ -575,10 +575,23 @@ check_collinearity.zerocount <- function( # 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 @@ check_collinearity.zerocount <- function( } 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 @@ check_collinearity.zerocount <- function( # 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")] @@ -1190,13 +1211,42 @@ check_collinearity.zerocount <- function( 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( @@ -1213,11 +1263,15 @@ check_collinearity.zerocount <- function( } } -# 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 @@ check_collinearity.zerocount <- function( 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" diff --git a/man/check_collinearity.Rd b/man/check_collinearity.Rd index 03d84684f..2373cc530 100644 --- a/man/check_collinearity.Rd +++ b/man/check_collinearity.Rd @@ -61,12 +61,12 @@ the square root of the VIF for terms with one coefficient. For multinomial models from \code{nnet::multinom()}, \code{mclogit::mblogit()}, \code{brglm2::brmultinom()} and \code{brglm2::bracl()}, and for \code{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 \code{bracl(parallel = TRUE)}, have one slope -per model matrix column. The adjusted VIF is then \verb{VIF^(1/(2*))}, where +level. Parallel models, such as \code{bracl(parallel = TRUE)}, and \code{vglm()} +models with the poissonff family have one slope per model matrix column. The adjusted VIF is then \verb{VIF^(1/(2*))}, where \verb{} 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 diff --git a/tests/testthat/test-check_collinearity.R b/tests/testthat/test-check_collinearity.R index 4cd2f6d37..431228498 100644 --- a/tests/testthat/test-check_collinearity.R +++ b/tests/testthat/test-check_collinearity.R @@ -130,6 +130,23 @@ test_that("check_collinearity, fallback returns NULL if matrix rebuild fails", { }) +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 @@ test_that("check_collinearity, ordinal clm models with offset", { 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 @@ test_that("check_collinearity, VGAM::vglm cumulative", { 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")