Feature Request: Allow custom evaluation functions in the `stats` argument

[n8thangreen]

Hi projpred team,

Thank you for developing and maintaining such a fantastic package. I am using projpred for clinical risk prediction modelling and was hoping to evaluate the submodels using decision-theoretic metrics, specifically Vickers' Net Benefit.

Currently, the stats argument in evaluation functions (like summary.vsel() and plot.vsel()) is restricted to a hardcoded list of strings (e.g., "elpd", "brier", "auc"). This works perfectly for standard metrics, but makes it difficult to evaluate submodels using custom loss functions or utility metrics that require external parameters (such as a willingness-to-pay decision threshold).

There is growing clinical demand for this type of decision-theoretic variable selection. In the frequentist space, packages like NBvarsel have popularised selecting variables based on Net Benefit. However, they rely on computationally intensive stepwise selection and provide only single-point estimates. Allowing custom evaluation functions in projpred would effectively enable users to perform a Bayesian, projection-predictive equivalent—yielding the full posterior uncertainty required for Health Technology Assessment (HTA) while retaining the speed and mathematical elegance of the L1-penalised search path.

Proposed Solution
I would like to propose an enhancement that allows the stats argument to accept a custom function (or a function factory) in addition to the standard character strings.

For example, a user could define a metric that locks in a specific threshold:

# Example function factory for Vickers' Net Benefit
vickers_nb <- function(threshold = 0.1) {
  f <- function(y_true, y_pred) {
    n <- length(y_true)
    tp <- sum(y_pred >= threshold & y_true == 1)
    fp <- sum(y_pred >= threshold & y_true == 0)
    
    (tp / n) - (fp / n) * (threshold / (1 - threshold))
  }
  attr(f, "metric_name") <- sprintf("Net Benefit (pt = %s)", threshold)
  return(f)
}

# Proposed API behaviour
cv_varsel(fit, stats = vickers_nb(threshold = 0.1))

Implementation Details
Architecturally, I imagine this would require:

1. Adding an is.function(stats) control flow to bypass the current string matching.
2. Routing y_true and the projected probabilities (y_pred) to the user's function during the evaluation loop.
3. Checking for a metric_name attribute on the function to properly label the y-axis in plot.vsel().

This approach is metric-agnostic, meaning the core package would not need to maintain specific health economic metrics natively, but would provide the flexibility for users to evaluate submodels using any custom utility function.

I would be very happy to work on this and submit a pull request if the team is open to the architectural concept. Please let me know your thoughts or if you foresee any major internal roadblocks with this design.