Getting Started with MOBID

Introduction

MOBID (Multi-Objective Biomarker Identification and Discovery) helps researchers build compact biomarker signatures when multiple cohorts are available and the number of measured features greatly exceeds the number of samples. The package delegates the heavy lift of multi-objective search to rmoo (NSGA-II/NSGA-III) while exposing a friendly interface for selecting metrics, restricting the feature pool, and summarising the resulting Pareto front.

This vignette walks through three steps:

Prepare multi-cohort expression data.
Assemble objectives from the metric registry.
Run optimize_panel() to discover candidate panels and inspect the results.

1. Simulate a Multi-cohort Dataset

For demonstration we simulate four cohorts with simulate_expression_data(). The metadata element records which genes carry the planted signal.

cohorts <- simulate_expression_data(p = 250, n = 100, k = 4, seed = 20240220)

lapply(cohorts$x_list, dim)
#> $x1
#> [1] 100 250
#> 
#> $x2
#> [1] 100 250
#> 
#> $x3
#> [1] 100 250
#> 
#> $x4
#> [1] 100 250
colnames(cohorts$x_list[[1]])[1:6]
#> [1] "gene_0001" "gene_0002" "gene_0003" "gene_0004" "gene_0005" "gene_0006"

Each element of x_list is a log-scale expression matrix with the same feature ordering across cohorts. The matching y_list contains binary outcome labels stored as factors with levels "No" and "Yes".

2. Configure Objectives

Metrics are registered in an extensible registry. You can inspect the available metrics with metric_registry().

names(metric_registry())
#>  [1] "auc"                        "balanced_accuracy"         
#>  [3] "cohort_auc_gap"             "cohort_auc_var"            
#>  [5] "cohort_leakage"             "f1"                        
#>  [7] "max_cohort_brier"           "min_cohort_auc"            
#>  [9] "npv"                        "num_features"              
#> [11] "pauc"                       "precision"                 
#> [13] "sensitivity"                "sensitivity_at_specificity"
#> [15] "specificity"                "specificity_at_sensitivity"

To emphasise clinical utility we combine sensitivity, specificity, and a panel size penalty. We also tighten the classification probability cutoff for sensitivity.

objectives <- define_objectives(
  metrics = c("sensitivity", "specificity", "num_features"),
  params = list(sensitivity = list(cutoff_prob = 0.6))
)
objectives[["sensitivity"]]$direction
#> [1] "maximize"

When multiple cohorts are available you can layer on transfer-oriented objectives that operate on cohort-level sensitivities, calibration, and feature shifts. These metrics automatically receive cohort labels during optimisation.

transfer_objectives <- define_objectives(
  metrics = c(
    "min_cohort_auc",
    "cohort_auc_gap",
    "cohort_auc_var",
    "max_cohort_brier"
  )
)

A Rule-Out Recipe

Rule-out diagnostics pin the operating point at high sensitivity and then push specificity as far as it will go. There is no dedicated helper for this — it is an ordinary define_objectives() call in which specificity_at_sensitivity’s target_sensitivity and pauc’s sens_floor both name the same operating point.

ruleout_objectives <- define_objectives(
  metrics = c("specificity_at_sensitivity", "pauc", "min_cohort_auc",
              "num_features"),
  params = list(
    specificity_at_sensitivity = list(target_sensitivity = 0.90),
    pauc = list(sens_floor = 0.90, partial_auc_correct = TRUE)
  )
)

No hard constraint is needed: because both metrics are already evaluated at the 0.90 sensitivity floor, candidates that cannot reach it score poorly on their own. Add a min_metric_constraint("sensitivity", threshold = 0.9) only if you want infeasible candidates discarded outright rather than ranked last.

Enforce Minimum Standards

Some clinical settings demand a hard floor on performance before optimising another objective. Supply constraint descriptors—such as min_metric_constraint()—via the constraints argument to optimize_panel(). The optimiser will discard any candidate failing the constraint.

constraints <- list(
  min_metric_constraint("sensitivity", threshold = 0.9,
                        params = list(cutoff_prob = 0.6))
)

3. Run NSGA-II Optimisation

optimize_panel() accepts single matrices or lists of cohorts. We restrict the feature pool to keep the search lightweight for this vignette.

set.seed(9001)
result <- optimize_panel(
  x = cohorts$x_list,
  y = cohorts$y_list,
  objectives = objectives,
  max_features = 5,
  feature_pool = colnames(cohorts$x_list[[1]])[1:50],
  nsga_control = list(popsize = 64, generations = 60)
)

The result is a BiomarkerPanelResult object that stores the selected features, objective scores, and metadata about the optimisation run.

slotNames(result)
result@features
result@metrics
head(result@objectives)

Held-Out Partitioning

The multi-cohort call above scores every candidate by 5-fold cross-validation within the training data by default. Pass fitness_mode together with train_ratio (and val_ratio) to carve off a stratified partition per cohort and score candidates against it instead: fitness_mode = "loco" scores each candidate on left-out cohorts, "within_cohort_val" and "within_cohort_rotating" score on a held-back validation share, and "cv" (the overall default) k-fold cross-validates within the training partition. Whenever train_ratio + val_ratio < 1, the remaining share is stored on result@control$heldout_x / heldout_y / heldout_cohort for calibrate_panel() to use afterwards.

Visualise Trade-offs

We can visualise Pareto-efficient solutions with plot_pareto_front(). Pick any pair of objectives and extract the corresponding columns.

library(ggplot2)
pareto_df <- subset(result@objectives, objective %in% c("sensitivity", "specificity"))
pareto_df <- reshape(
  pareto_df,
  idvar = "solution_id",
  timevar = "objective",
  direction = "wide"
)
plot_pareto_front(
  data.frame(
    objective_x = pareto_df$value.sensitivity,
    objective_y = pareto_df$value.specificity,
    label = paste0("Solution ", pareto_df$solution_id)
  ),
  xlab = "Sensitivity",
  ylab = "Specificity"
)

Next Steps

Explore additional metrics by registering custom ones with register_metric().
Run BiocCheck::BiocCheck() and R CMD check before submitting to Bioconductor.

Package Authors