Avoid profile likelihood fits in Firth models

Author: wolski

NEWS.md

@@ -13,6 +13,8 @@
 - Hardened `plot_pca()`: errors early on duplicated sample names, an all-missing matrix, or too few samples instead of returning `NULL` (which broke `pca_plotly()`); joins scores to annotation with an explicit `by`; makes `prcomp(center = TRUE, scale. = FALSE)` explicit.
 - Hardened abundance heatmaps for sparse significant-feature subsets: when row or column distances are non-finite because of missing values, `plot_heatmap()` now falls back to the input order instead of returning a `ComplexHeatmap` object that fails during drawing.
 - `LFQDataPlotter$heatmap()` now shows only the `top_n` most variable features (default 1000), ranked by the prolfqua per-feature statistic (CV for untransformed data, sd for transformed, via `LFQDataStats`). Row clustering uses `stats::hclust`, which errors above 65536 features, so peptide-list / entrapment searches with tens of thousands of degenerate protein groups no longer crash the QC heatmap. Pass `top_n = NULL` (or `Inf`) to keep every feature.
+- `StrategyLogistf` now uses Wald confidence intervals instead of profiling every coefficient, preventing `firth_nested`
+  analyses from stalling for days on proteins with hundreds or thousands of peptides.
 
 # prolfqua 1.6.1
 

R/logistf.R

@@ -365,7 +365,7 @@ StrategyLogistf <- R6::R6Class(
       DFT <- x |>
         dplyr::group_by(dplyr::across(dplyr::all_of(predictor_vars))) |>
         dplyr::summarize(Freq = dplyr::n(), .groups = "drop")
-      tryCatch(logistf::logistf(self$formula, data = DFT, weights = Freq), error = .error_handler)
+      tryCatch(logistf::logistf(self$formula, data = DFT, weights = Freq, pl = FALSE), error = .error_handler)
     },
 
     #' @description Check if model is singular (NA coefficients or df < 2)

tests/testthat/test-Model.R

@@ -102,6 +102,17 @@ test_that("ModelFirth", {
   expect_s3_class(cv$plot, "ggplot")
 })
 
+test_that("StrategyLogistf avoids coefficient-wise profile likelihood fits", {
+  data <- data.frame(
+    bin_resp = c(1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L),
+    condition = rep(c("A", "B"), each = 4)
+  )
+
+  fit <- strategy_logistf("bin_resp ~ condition")$model_fun(data)
+
+  expect_setequal(fit$method.ci, "Wald")
+})
+
 test_that("build_model_glm_protein and build_model_glm_peptide return ModelFirth", {
   prot <- sim_lfq_data_protein_config(Nprot = 10, with_missing = TRUE, weight_missing = 0.5, seed = 7)
   prot_lfq <- LFQData$new(prot$data, prot$config)