feat(engine): default fitting backend = mgcv::bam + gam escape (v0.1.1 Feature 4)

Implements Feature 4 of the locked v0.1.1 plan
(PLAN_v011_features.md §2.4) — the one non-byte-identical change
in the release. Documented loudly in NEWS.md and @details.

Engine dispatch
* New `engine = c("bam", "gam")` argument, default "bam" on both
  janusplot() and janusplot_data(). bam is mgcv's "big additive
  model" backend — block-Lanczos / fREML / lower memory.
  ~3-10x wall-time speedup at janusplot's scale (k = 15-25 vars,
  600+ pairwise fits per call). No new dependency.
* `.fit_one_gam()` dispatch helper used at every fitting site
  (initial fit + refit-on-flagged-cells under auto_refit_k =
  TRUE) so the engine choice is consistent within a single call.
* `bam` objects inherit from `gam` (class is c("bam","gam","glm",
  "lm")), so downstream code paths (k.check, predict, derivative
  LP-matrix arithmetic, shape-metric extractors) are unaffected.

Method default
* `method` argument default is now NULL, resolved per-engine at
  fit time: "fREML" for bam (mgcv's recommended at scale), "REML"
  for gam (v0.1.0 behaviour). Users who passed method explicitly
  in v0.1.0 see no behaviour change.
* User-supplied method (e.g. "GCV.Cp") overrides the per-engine
  default — tested in test-engine.R.

bam-only knobs
* `discrete = FALSE` default. Opt-in to mgcv's covariate-
  discretisation. Further ~2-5x speedup at sub-pixel prediction
  shift cost.
* `nthreads = 1L` default. Intra-fit threading. Default 1 to
  avoid CPU oversubscription when combined with parallel = TRUE
  (which fans out across pair-fits via future.apply).
* Both ignored silently on engine = "gam".

Provenance
* `engine` + `method` columns appear in
  janusplot(..., with_data = TRUE)$data and on every
  janusplot_data()$pairs[[i]] entry. Paper figures can document
  exactly which backend produced their EDFs.

Numerical-shift documentation
* NEWS.md leads with a NUMERICAL-SHIFT NOTICE block ABOVE the
  per-feature entries — explicit callout of the ~1-3% EDF /
  asymmetry index shift between v0.1.0 (REML) and v0.1.1
  (fREML). Recovery via engine = "gam" stated prominently.
* @details on engine arg repeats the same.
* @PARAM method documents the new NULL default + per-engine
  resolution.

Backward-compat escape
* Every v0.1.0 user can call janusplot(..., engine = "gam") and
  get byte-identical fit results vs v0.1.0.
* The vdiffr visual-regression suite is now PINNED to
  engine = "gam" precisely so the snapshots track the v0.1.0
  visual contract; suite acts as a regression gate on the
  backward-compat escape itself. Same pinning applied in
  test-derivatives.R for the d1 derivative snapshots.

Tests: new tests/testthat/test-engine.R (10 expectations):
  - default engine surfaces "bam" + "fREML" in with_data
  - engine = "gam" surfaces "gam" + "REML"
  - both engines agree to <10% relative shift on linear DGP
  - both detect non-linearity on quadratic DGP (forward dir)
  - discrete = TRUE on bam renders, silently ignored on gam
  - explicit method overrides per-engine default
  - invalid engine errors via arg_match
  - janusplot_data carries engine + method per pair

devtools::test(): 473 expectations pass (was 449 after Feature 3).
R CMD check --no-manual: 0 errors / 0 warnings / 2 environmental
NOTEs (HTML Tidy + 4-component dev version). lintr 0 lints.

This completes the four-feature v0.1.1 plan. Next: vignette
expansion + CRAN-comments draft + version bump to 0.1.1.

Author: max578

NEWS.md

@@ -1,5 +1,57 @@
 # janusplot (development version)
 
+### v0.1.1 NUMERICAL-SHIFT NOTICE (engine default change)
+
+**v0.1.1 changes the default GAM fitting backend from `mgcv::gam`
+to `mgcv::bam`** (Feature 4 below). `bam` uses fREML estimation
+instead of REML, which differs by ~1-3% in effective degrees of
+freedom on identical data. **Existing v0.1.0 users will see
+slightly different EDFs, asymmetry indices, and per-cell colour
+fills when they upgrade.** This is the single non-byte-identical
+change in v0.1.1 — every other v0.1.1 feature is strictly additive.
+
+**Recovery.** Set `engine = "gam"` on `janusplot()` or
+`janusplot_data()` to reproduce v0.1.0 numerical output verbatim.
+The package's vdiffr visual-regression suite pins
+`engine = "gam"` for exactly this purpose — every old snapshot
+remains valid under the backward-compat escape.
+
+### v0.1.1 Feature 4 — Default fitting engine = `bam` (with `gam` escape)
+
+* **New `engine = c("bam", "gam")` argument**, default `"bam"`.
+  At janusplot's scale (k = 15-25 vars, 600+ fits per call)
+  `mgcv::bam`'s block-Lanczos solve + fREML estimation delivers
+  ~3-10x wall-time speedup vs `mgcv::gam` without any new
+  dependency (`mgcv` already exports both).
+* **fREML vs REML.** `bam` defaults to fREML (fast REML); `gam`
+  defaults to REML. The two methods optimise the same penalty
+  target via different paths — EDFs differ by ~1-3% on identical
+  data. v0.1.1 surfaces this as the one numerical break with
+  v0.1.0, documented prominently above.
+* **`method` argument default is now `NULL`**, resolved
+  per-engine: `"fREML"` for bam, `"REML"` for gam. Users who
+  passed `method` explicitly in v0.1.0 see no behaviour change.
+* **`discrete = FALSE` (bam-only)** — opt-in to mgcv's
+  covariate-discretisation optimisation. Further ~2-5x speedup
+  at sub-pixel prediction shift cost.
+* **`nthreads = 1L` (bam-only)** — intra-fit threading. Default
+  1 to avoid oversubscription when combined with `parallel =
+  TRUE` (which fans out across pair-fits already).
+* **`engine` + `method` provenance** — both columns now appear
+  in `janusplot(..., with_data = TRUE)$data` and on every
+  `janusplot_data()$pairs[[i]]` entry. Useful for paper figures
+  whose methodology section needs to document which backend
+  produced the EDFs they report.
+* **`bam` inherits from `gam`** — `class(bam_fit)` is
+  `c("bam", "gam", "glm", "lm")`, so `mgcv::k.check()`,
+  `predict.gam()`, derivative LP-matrix arithmetic, and every
+  shape-metric extractor work without modification. Engine is
+  plumbing, not redesign.
+* **vdiffr suite pinned to `engine = "gam"`** so every legacy
+  visual snapshot continues to validate under the v0.1.0
+  backward-compat escape. Acts as a regression gate on the
+  `engine = "gam"` recovery path.
+
 ### v0.1.1 Feature 3 — Axes rendering modes + figure-to-file output
 
 * **`axes`** rendering knob with four modes:

R/internal-fit.R

@@ -412,14 +412,57 @@
 # Fit a single pairwise GAM and summarise
 # ---------------------------------------------------------------
 
+# ---------------------------------------------------------------
+# Fitting-engine dispatch (Feature 4). `bam` is mgcv's "big additive
+# model" — same formula language as `gam`, but uses fREML (fast REML)
+# by default, block-Lanczos / discrete-method optimisations for the
+# basis-coefficient solve, lower memory, and a built-in `nthreads`
+# argument. bam objects inherit from gam (`class(b)` ==
+# c("bam", "gam", "glm", "lm")) so every downstream code path
+# (predict, summary, k.check, derivative LP-matrix arithmetic)
+# works without modification — engine is plumbing, not redesign.
+#
+# Default method-per-engine: `fREML` for bam (mgcv's recommended at
+# scale), `REML` for gam (v0.1.0 behaviour). A user-supplied `method`
+# overrides both.
+# ---------------------------------------------------------------
+
+.engine_default_method <- function(engine) {
+  if (identical(engine, "bam")) "fREML" else "REML"
+}
+
+.fit_one_gam <- function(fml, dat, engine, method, discrete, nthreads, ...) {
+  method <- if (is.null(method) || identical(method, "auto") ||
+                identical(method, "default")) {
+    .engine_default_method(engine)
+  } else {
+    method
+  }
+  if (identical(engine, "bam")) {
+    mgcv::bam(
+      formula  = fml,
+      data     = dat,
+      method   = method,
+      discrete = isTRUE(discrete),
+      nthreads = as.integer(nthreads %||% 1L),
+      ...
+    )
+  } else {
+    mgcv::gam(formula = fml, data = dat, method = method, ...)
+  }
+}
+
 .fit_pair <- function(x_name, y_name, data_full, adjust, method, k, bs,
                      na_action, n_grid = 100L,
                      derivatives = integer(),
                      derivative_ci = "pointwise",
                      derivative_ci_nsim = 1000L,
                      k_check_thresholds = .default_k_thresholds(),
                      auto_refit_k = FALSE,
-                     k_max_iter = 2L, ...) {
+                     k_max_iter = 2L,
+                     engine = "bam",
+                     discrete = FALSE,
+                     nthreads = 1L, ...) {
   k_val <- .resolve_k(k, x_name)
   if (na_action == "pairwise") {
     dat <- .complete_pair(data_full, x_name, y_name, adjust)
@@ -434,8 +477,16 @@
 
   # Initial fit at user-requested k (mgcv resolves k=-1 internally).
   fml <- .build_formula(y_name, x_name, k_val, bs, adjust)
+  resolved_method <- if (is.null(method) || identical(method, "auto") ||
+                         identical(method, "default")) {
+    .engine_default_method(engine)
+  } else {
+    method
+  }
   fit <- tryCatch(
-    mgcv::gam(fml, data = dat, method = method, ...),
+    .fit_one_gam(fml, dat,
+                 engine = engine, method = resolved_method,
+                 discrete = discrete, nthreads = nthreads, ...),
     error = function(e) e
   )
   if (inherits(fit, "error")) {
@@ -471,7 +522,9 @@
       k_new <- min(2 * k_cur, k_cap)
       fml_new <- .build_formula(y_name, x_name, k_new, bs, adjust)
       new_fit <- tryCatch(
-        mgcv::gam(fml_new, data = dat, method = method, ...),
+        .fit_one_gam(fml_new, dat,
+                     engine = engine, method = resolved_method,
+                     discrete = discrete, nthreads = nthreads, ...),
         error = function(e) e
       )
       if (inherits(new_fit, "error")) break
@@ -563,6 +616,8 @@
     dev_exp  = dev_exp,
     n_used   = n_used,
     error    = NA_character_,
+    engine   = engine,
+    method   = resolved_method,
     k_check  = list(
       k_prime         = k_diag$k_prime,
       k_index         = k_diag$k_index,
@@ -589,6 +644,7 @@
     raw = data.frame(),
     edf = NA_real_, pvalue = NA_real_, dev_exp = NA_real_,
     n_used = n_used, error = NA_character_,
+    engine = NA_character_, method = NA_character_,
     k_check = list(
       k_prime        = NA_real_,
       k_index        = NA_real_,
@@ -618,7 +674,10 @@
                            derivative_ci_nsim = 1000L,
                            k_check_thresholds = .default_k_thresholds(),
                            auto_refit_k = FALSE,
-                           k_max_iter = 2L, ...) {
+                           k_max_iter = 2L,
+                           engine = "bam",
+                           discrete = FALSE,
+                           nthreads = 1L, ...) {
   .check_parallel_plan(parallel)
   if (na_action == "complete") {
     data <- data[stats::complete.cases(
@@ -640,7 +699,9 @@
       derivative_ci_nsim = derivative_ci_nsim,
       k_check_thresholds = k_check_thresholds,
       auto_refit_k = auto_refit_k,
-      k_max_iter = k_max_iter, ...
+      k_max_iter = k_max_iter,
+      engine = engine, discrete = discrete,
+      nthreads = nthreads, ...
     )
   }
 

R/janusplot.R

@@ -25,7 +25,37 @@
 #'   `gam(y ~ s(x) + s(age) + s(site, bs = "re"))` for each pair.
 #'   Default `NULL` fits unadjusted pairwise smooths.
 #' @param method Smoothing-parameter estimation method passed to
-#'   [mgcv::gam()]. Default `"REML"` per mgcv recommendation.
+#'   the chosen fitting backend. Default `NULL` resolves
+#'   per-engine: `"fREML"` for `engine = "bam"` (mgcv's
+#'   recommended at scale), `"REML"` for `engine = "gam"`
+#'   (the v0.1.0 behaviour). Pass any valid mgcv method string
+#'   to override.
+#' @param engine One of `"bam"` (default, **new in v0.1.1**) or
+#'   `"gam"`. Selects mgcv's fitting backend:
+#'   * `"bam"` — [mgcv::bam()]. Block-Lanczos solve + fREML
+#'     estimation + lower memory. ~3-10x speedup at janusplot's
+#'     scale (k = 15-25 vars, 600+ pairwise fits per call). The
+#'     **default**, and the one non-byte-identical change in
+#'     v0.1.1: fREML differs from REML by ~1-3% in EDF on
+#'     identical data, so the asymmetry index may shift by
+#'     similar amounts vs v0.1.0 output. Recoverable verbatim
+#'     via `engine = "gam"`.
+#'   * `"gam"` — [mgcv::gam()]. The v0.1.0 backend. Use for
+#'     backward-compat reproduction, very small n (< 200) where
+#'     bam's setup overhead exceeds its solve gain, or
+#'     methodologically sensitive contexts that require REML
+#'     rather than fREML.
+#' @param discrete Logical. `bam`-only opt-in to mgcv's
+#'   covariate-discretisation optimisation. Further ~2-5x
+#'   speedup at the cost of small (sub-pixel at typical
+#'   janusplot resolution) prediction-shift. Default `FALSE`.
+#'   Ignored when `engine = "gam"`.
+#' @param nthreads Integer. `bam`-only intra-fit threading.
+#'   Default `1L` to avoid oversubscription when combined with
+#'   `parallel = TRUE` (`future.apply` fans out pair fits
+#'   across cores; nthreads > 1 within each pair would
+#'   double-book CPUs). Raise above 1 only when `parallel =
+#'   FALSE`. Ignored when `engine = "gam"`.
 #' @param k Integer, or named list mapping variable names to integers.
 #'   Basis dimension for `s()`. Default `-1L` (mgcv's automatic choice).
 #' @param bs Basis type for `s()`. Default `"tp"` (thin plate).
@@ -354,9 +384,12 @@ janusplot <- function(
     data,
     vars = NULL,
     adjust = NULL,
-    method = "REML",
+    method = NULL,
     k = -1L,
     bs = "tp",
+    engine = c("bam", "gam"),
+    discrete = FALSE,
+    nthreads = 1L,
     order = c("original", "hclust", "alphabetical"),
     show_data = TRUE,
     show_ci = TRUE,
@@ -404,6 +437,15 @@ janusplot <- function(
   na_action     <- rlang::arg_match(na_action)
   compact       <- rlang::arg_match(compact)
   axes          <- rlang::arg_match(axes)
+  engine        <- rlang::arg_match(engine)
+  if (!is.logical(discrete) || length(discrete) != 1L || is.na(discrete)) {
+    cli::cli_abort("{.arg discrete} must be TRUE or FALSE.")
+  }
+  if (!is.numeric(nthreads) || length(nthreads) != 1L ||
+      !is.finite(nthreads) || nthreads < 1L) {
+    cli::cli_abort("{.arg nthreads} must be a single positive integer.")
+  }
+  nthreads <- as.integer(nthreads)
   glyph_style   <- rlang::arg_match(glyph_style)
   labels        <- rlang::arg_match(labels)
   diagonal      <- rlang::arg_match(diagonal)
@@ -607,7 +649,8 @@ janusplot <- function(
     derivative_ci_nsim = derivative_ci_nsim,
     k_check_thresholds = k_thresholds,
     auto_refit_k = auto_refit_k,
-    k_max_iter = k_max_iter, ...
+    k_max_iter = k_max_iter,
+    engine = engine, discrete = discrete, nthreads = nthreads, ...
   )
 
   .summarise_k_check(fits, k_thresholds, auto_refit_k)              # nolint: object_usage_linter.
@@ -826,6 +869,8 @@ janusplot <- function(
       k_final            = kc$k_final,
       k_iterations       = as.integer(kc$k_iterations %||% 0L),
       k_at_cap           = isTRUE(kc$k_at_cap),
+      engine             = f$engine %||% NA_character_,
+      method             = f$method %||% NA_character_,
       stringsAsFactors   = FALSE
     )
   })
@@ -1020,7 +1065,7 @@ janusplot_data <- function(
     data,
     vars = NULL,
     adjust = NULL,
-    method = "REML",
+    method = NULL,
     k = -1L,
     bs = "tp",
     na_action = c("pairwise", "complete"),
@@ -1034,9 +1079,21 @@ janusplot_data <- function(
     k_check_thresholds = NULL,
     auto_refit_k       = FALSE,
     k_max_iter         = 2L,
+    engine             = c("bam", "gam"),
+    discrete           = FALSE,
+    nthreads           = 1L,
     ...) {
   na_action     <- rlang::arg_match(na_action)
   derivative_ci <- rlang::arg_match(derivative_ci)
+  engine        <- rlang::arg_match(engine)
+  if (!is.logical(discrete) || length(discrete) != 1L || is.na(discrete)) {
+    cli::cli_abort("{.arg discrete} must be TRUE or FALSE.")
+  }
+  if (!is.numeric(nthreads) || length(nthreads) != 1L ||
+      !is.finite(nthreads) || nthreads < 1L) {
+    cli::cli_abort("{.arg nthreads} must be a single positive integer.")
+  }
+  nthreads <- as.integer(nthreads)
   .validate_inputs(data, vars, adjust, na_action)
   vars <- .resolve_vars(data, vars)
   k_thresholds <- k_check_thresholds %||% .default_k_thresholds()   # nolint: object_usage_linter.
@@ -1099,7 +1156,8 @@ janusplot_data <- function(
     derivative_ci_nsim = derivative_ci_nsim,
     k_check_thresholds = k_thresholds,
     auto_refit_k = auto_refit_k,
-    k_max_iter = k_max_iter, ...
+    k_max_iter = k_max_iter,
+    engine = engine, discrete = discrete, nthreads = nthreads, ...
   )
 
   .summarise_k_check(fits, k_thresholds, auto_refit_k)              # nolint: object_usage_linter.
@@ -1163,7 +1221,9 @@ janusplot_data <- function(
         deriv_yx           = f_yx$deriv,
         deriv_xy           = f_xy$deriv,
         k_check_yx         = f_yx$k_check,
-        k_check_xy         = f_xy$k_check
+        k_check_xy         = f_xy$k_check,
+        engine             = f_yx$engine %||% f_xy$engine %||% NA_character_,
+        method             = f_yx$method %||% f_xy$method %||% NA_character_
       )
     }
   }