diff --git a/DESCRIPTION b/DESCRIPTION
index 4191ecb..4c84366 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -35,8 +35,6 @@ Suggests:
     xgboost
 VignetteBuilder: knitr
 Config/testthat/edition: 3
-Remotes:
-    vip=url::https://cran.r-project.org/src/contrib/Archive/vip/vip_0.4.6.tar.gz
 Imports:
     arrow,
     BiocParallel,
@@ -72,7 +70,6 @@ Imports:
     tidyr,
     tune,
     utils,
-    vip,
     withr,
     workflows,
     workflowsets,
diff --git a/NAMESPACE b/NAMESPACE
index e4fb286..6533bbb 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,9 +1,27 @@
 # Generated by roxygen2: do not edit by hand
 
+export(.getFeatureTypes)
 export(.getTargetVarName)
+export(.register_parquet_views)
 export(applyBenjaminiHochberg)
 export(buildLRModel)
+export(buildPerfPq)
+export(buildPerfPqCrossCountry)
+export(buildPerfPqCrossDrug)
+export(buildPerfPqCrossYear)
+export(buildPerfPqLOOCountry)
+export(buildPerfPqLOODrug)
+export(buildPerfPqLOOYear)
+export(buildPerfPqMDR)
+export(buildPerfPqYearCountry)
+export(buildPredPqMDR)
 export(buildRecipe)
+export(buildTopFeatsPq)
+export(buildTopFeatsPqLOOCountry)
+export(buildTopFeatsPqLOODrug)
+export(buildTopFeatsPqLOOYear)
+export(buildTopFeatsPqMDR)
+export(buildTopFeatsPqYearCountry)
 export(buildTuningGrid)
 export(buildWflow)
 export(calculateEvalMets)
@@ -18,6 +36,7 @@ export(generateMLInputs)
 export(getConfusionMatrix)
 export(getNumFeat)
 export(loadMLInputTibble)
+export(parse_ml_filename)
 export(plotBaselineComparison)
 export(plotCM)
 export(plotCrossDrug)
@@ -150,8 +169,6 @@ importFrom(tune,extract_fit_parsnip)
 importFrom(tune,finalize_workflow)
 importFrom(tune,select_best)
 importFrom(tune,tune_grid)
-importFrom(vip,vi)
-importFrom(vip,vip)
 importFrom(workflows,add_model)
 importFrom(workflows,add_recipe)
 importFrom(workflows,workflow)
diff --git a/R/core_ml.R b/R/core_ml.R
index 7e398f2..588b54e 100644
--- a/R/core_ml.R
+++ b/R/core_ml.R
@@ -43,7 +43,6 @@
 #' @importFrom tune finalize_workflow
 #' @importFrom tune select_best
 #' @importFrom tune tune_grid
-#' @importFrom vip vi
 #' @importFrom workflows add_model
 #' @importFrom workflows add_recipe
 #' @importFrom workflows workflow
@@ -751,6 +750,31 @@ calculateEvalMets <- function(test_data_plus_predictions) {
   return(c(f1, auprc, bal_acc, mcc, nmcc, log2_apop))
 }
 
+#' Variable importance for a binomial glmnet fit
+#'
+#' Internal replacement for `vip::vi()`, dropped as a dependency after it left
+#' CRAN. Returns `|coefficient|` as `Importance` and its direction as `Sign`,
+#' sorted by decreasing importance.
+#'
+#' Importance is taken at glmnet's minimum lambda, not the tuned `penalty`.
+#' That is what `vip::vi()` did.
+#'
+#' @param fit A fitted workflow backed by a binomial glmnet engine.
+#' @return A tibble with `Variable`, `Importance`, and `Sign` columns.
+#' @keywords internal
+.viGlmnet <- function(fit) {
+  glmnet_fit <- parsnip::extract_fit_engine(fit)
+  coefs <- stats::coef(glmnet_fit, s = min(glmnet_fit$lambda))[, 1, drop = TRUE]
+  coefs <- coefs[names(coefs) != "(Intercept)"]
+
+  tibble::tibble(
+    Variable = names(coefs),
+    Importance = abs(unname(coefs)),
+    Sign = ifelse(coefs > 0, "POS", "NEG")
+  ) |>
+    dplyr::arrange(dplyr::desc(Importance))
+}
+
 #' extractTopFeats()
 #'
 #' Returns the top features that an ML model found to be important for
@@ -793,30 +817,8 @@ extractTopFeats <- function(
     stop("Please specify either `n_top_feats` or `prop_vi_top_feats`.")
   }
 
-  feats_arranged <- fit |>
-    workflowsets::extract_fit_parsnip() |>
-    vip::vi() |>
-    dplyr::arrange(dplyr::desc(Importance))
-
-  if (!is.na(n_top_feats)) {
-    top_feats_and_VIs <- feats_arranged |> dplyr::slice(seq_len(n_top_feats))
-  } else if (any(!is.na(prop_vi_top_feats))) {
-    cum_vi_lower <- prop_vi_top_feats[1] * sum(feats_arranged$Importance)
-    cum_vi_upper <- prop_vi_top_feats[2] * sum(feats_arranged$Importance)
-
-    top_feats_and_VIs <- feats_arranged |>
-      dplyr::mutate(cum_imp = cumsum(Importance)) |>
-      dplyr::filter(cum_imp < cum_vi_upper & cum_imp > cum_vi_lower)
-  }
-
-  top_feat_tibble <- tibble::tibble(
-    Variable = dplyr::pull(top_feats_and_VIs, Variable),
-    Importance = dplyr::pull(top_feats_and_VIs, Importance),
-    Sign = dplyr::pull(top_feats_and_VIs, Sign)
-  )
-
-  # Take a different approach if using multi-class (the previous code would give
-  # a less meaningful result).
+  # Multi-class: glmnet returns one coefficient matrix per class, so there is no
+  # single ranking to slice. Handled before `.viGlmnet()`, which expects one.
   if (is(fit$fit$actions$model$spec, "multinom_reg")) {
     warning(
       "Extracting top features from a multi-class model. ",
@@ -844,7 +846,28 @@ extractTopFeats <- function(
         values_fill = 0 # Fill missing values with 0.
       ) |>
       dplyr::filter(dplyr::if_any(-Variable, ~ . != 0))
+
+    return(top_feat_tibble)
+  }
+
+  feats_arranged <- .viGlmnet(fit)
+
+  if (!is.na(n_top_feats)) {
+    top_feats_and_VIs <- feats_arranged |> dplyr::slice(seq_len(n_top_feats))
+  } else if (any(!is.na(prop_vi_top_feats))) {
+    cum_vi_lower <- prop_vi_top_feats[1] * sum(feats_arranged$Importance)
+    cum_vi_upper <- prop_vi_top_feats[2] * sum(feats_arranged$Importance)
+
+    top_feats_and_VIs <- feats_arranged |>
+      dplyr::mutate(cum_imp = cumsum(Importance)) |>
+      dplyr::filter(cum_imp < cum_vi_upper & cum_imp > cum_vi_lower)
   }
 
+  top_feat_tibble <- tibble::tibble(
+    Variable = dplyr::pull(top_feats_and_VIs, Variable),
+    Importance = dplyr::pull(top_feats_and_VIs, Importance),
+    Sign = dplyr::pull(top_feats_and_VIs, Sign)
+  )
+
   return(top_feat_tibble)
 }
diff --git a/R/plot_ml.R b/R/plot_ml.R
index 412c0ef..2c6eb62 100644
--- a/R/plot_ml.R
+++ b/R/plot_ml.R
@@ -18,7 +18,6 @@
 #' @importFrom ggplot2 ylim
 #' @importFrom graphics barplot
 #' @importFrom tune extract_fit_parsnip
-#' @importFrom vip vip
 #' @importFrom yardstick pr_curve
 NULL
 
diff --git a/R/run_ml_pipeline.R b/R/run_ml_pipeline.R
index db0fc64..f021db0 100644
--- a/R/run_ml_pipeline.R
+++ b/R/run_ml_pipeline.R
@@ -328,7 +328,7 @@ runMLPipeline <- function(
     prop_vi_top_feats = prop_vi_top_feats
   )
 
-  # Account for the case where `vip::vi()` assigned a variable importance of
+  # Account for the case where a feature was assigned a variable importance of
   # zero, thereby not including the full `n_top_feats` requested. Do this only
   # if `n_top_feats` was specified instead of `prop_vi_top_feats`. Features will
   # be randomly assigned as top hits. This is important for downstream
diff --git a/tests/testthat/test-core-ml.R b/tests/testthat/test-core-ml.R
index 4220ee0..694f11e 100644
--- a/tests/testthat/test-core-ml.R
+++ b/tests/testthat/test-core-ml.R
@@ -126,3 +126,53 @@ test_that("extractTopFeats rejects conflicting / missing selectors", {
   expect_setequal(colnames(top), c("Variable", "Importance", "Sign"))
   expect_lte(nrow(top), 2)
 })
+
+# A multi-class (MDR-style) fit, where glmnet returns one coefficient matrix per
+# class. This is the case that must not reach `.viGlmnet()`.
+make_multiclass_fit <- function() {
+  set.seed(42)
+  n <- 30
+  mc <- tibble::tibble(
+    genome_id = paste0("g", seq_len(n)),
+    resistant_classes = factor(rep(c("A", "B", "C"), each = n / 3))
+  )
+  for (j in 1:8) mc[[paste0("feat_", j)]] <- as.integer(rbinom(n, 1, 0.5))
+
+  mod <- parsnip::multinom_reg(penalty = 0.01, mixture = 0.5) |>
+    parsnip::set_engine("glmnet")
+  buildWflow(mod, buildRecipe(mc)) |> parsnip::fit(data = mc)
+}
+
+test_that("extractTopFeats handles multi-class fits with per-class columns", {
+  skip_if_missing_deps()
+  fit <- make_multiclass_fit()
+
+  expect_warning(
+    top <- extractTopFeats(fit, n_top_feats = 5),
+    "multi-class"
+  )
+
+  # One column per class, not the binary Variable/Importance/Sign triple.
+  expect_setequal(colnames(top), c("Variable", "A", "B", "C"))
+  expect_gt(nrow(top), 0)
+})
+
+test_that(".viGlmnet reproduces vip::vi() for binomial glmnet fits", {
+  skip_if_missing_deps()
+  fx <- make_pipeline_fixture()
+  mod <- parsnip::logistic_reg(penalty = 0.01, mixture = 0) |>
+    parsnip::set_engine("glmnet")
+  fit <- buildWflow(mod, buildRecipe(fx)) |> parsnip::fit(data = fx)
+
+  vi <- .viGlmnet(fit)
+  expect_setequal(colnames(vi), c("Variable", "Importance", "Sign"))
+  expect_true(all(vi$Sign %in% c("POS", "NEG")))
+  # Downstream slicing relies on the decreasing sort `vip::vi()` provided.
+  expect_equal(vi$Importance, sort(vi$Importance, decreasing = TRUE))
+
+  # Importance is |coefficient| at the minimum lambda, not the tuned penalty.
+  gf <- parsnip::extract_fit_engine(fit)
+  expected <- stats::coef(gf, s = min(gf$lambda))[, 1, drop = TRUE]
+  expected <- expected[names(expected) != "(Intercept)"]
+  expect_equal(vi$Importance, sort(abs(unname(expected)), decreasing = TRUE))
+})
diff --git a/tests/testthat/test-pipeline.R b/tests/testthat/test-pipeline.R
index 8333942..fbcfe1d 100644
--- a/tests/testthat/test-pipeline.R
+++ b/tests/testthat/test-pipeline.R
@@ -10,7 +10,6 @@ skip_pipeline_if_missing <- function() {
   skip_if_not_installed("tune")
   skip_if_not_installed("workflows")
   skip_if_not_installed("yardstick")
-  skip_if_not_installed("vip")
 }
 
 minimal_grid_args <- list(