Skip to contents

Empirically Calibrated Simulations for bifrost, Part 2: Tuning and Application

Source: vignettes/simulation-study-part-2.Rmd
simulation-study-part-2.Rmd

Development note: This simulation and tuning workflow is under active development and remains somewhat experimental. The interface and recommendations may change as it is evaluated more broadly.

How should we choose search settings for an empirical dataset? More permissive settings may recover more evolutionary shifts, but they may also infer shifts where none occurred. Simulations help assess this tradeoff using datasets whose true shifts are known.

In Part 1, we evaluated an acceptance threshold of 10 and minimum clade size of 10 under null, proportional-shift, and integration-rate scenarios. Here we compare alternative settings on the same datasets, choose settings separately for GIC and BIC, and apply them to an empirical search. The simulations use the passerine body-plan data from Berv et al. (2026) for calibration; this tuning exercise is a new worked example, not an analysis reported in that study.

This vignette can be followed independently of Part 1. Saved simulation summaries let you reproduce the tables and parameter selection without rerunning the simulations. The longer simulation and empirical-search examples are optional and remain switched off unless you explicitly enable them.

Start by loading bifrost.

# Load bifrost before constructing templates or tuning search settings.
library(bifrost)

This part requests the passerine tree (passerine-tree, approximately 81 KB), trait matrix (passerine-traits, approximately 212 KB), and simulation preview tables (simulation-preview-tables, approximately 6 KB). First use downloads and verifies each file; later calls reuse the verified cache. The tree and trait files are the same ones used in the passerine articles. To download an updated copy of the saved summaries, use:

# Refresh the preview tables only when an updated copy is wanted.
bifrost_example_file("simulation-preview-tables", refresh = TRUE)

Build an empirically calibrated template

Here we use the empirical calibration structure from Berv et al. (2026), which models 12 skeletal traits jointly while keeping log body mass as a phylogenetic covariate.

# Resolve the passerine tree and calibrated skeletal measurements.
tree_path <- bifrost_example_file("passerine-tree")
trait_path <- bifrost_example_file("passerine-traits")

# Align the trait rows to the tree's tip order.
bird_tree <- ape::read.tree(tree_path)
bodyplan_data <- readRDS(trait_path)
bodyplan_data <- as.matrix(bodyplan_data[bird_tree$tip.label, , drop = FALSE])

# Keep body mass as the predictor and the skeletal variables as responses.
skeletal_cols <- setdiff(colnames(bodyplan_data), "vertnet_mass")
bodyplan_data <- bodyplan_data[, c(skeletal_cols, "vertnet_mass"), drop = FALSE]

# Express the multivariate response and body-mass covariate by column index.
formula_str <- sprintf(
  "trait_data[, 1:%d] ~ trait_data[, %d]",
  length(skeletal_cols),
  ncol(bodyplan_data)
)

# Confirm the aligned sample size and response dimensionality.
c(
  species = nrow(bodyplan_data),
  skeletal_traits = length(skeletal_cols),
  total_columns = ncol(bodyplan_data)
)
#>         species skeletal_traits   total_columns 
#>            2057              12              13

The template uses a single-regime multivariate BM fit with mvMORPH::mvgls() and method = "LL", following Clavel, Aristide, and Morlon (2019) and Clavel and Morlon (2020). Its covariance summaries provide the empirical scale for a fresh ancestral covariance matrix in each replicate.

# Opt in before fitting the calibration model and running new simulations.
run_simulations <- FALSE
if (run_simulations) {
  set.seed(5, kind = "Mersenne-Twister")
  bodyplan_template <- createSimulationTemplate(
    baseline_tree = bird_tree,
    trait_data = bodyplan_data,
    formula = formula_str,
    response_columns = seq_along(skeletal_cols),
    predictor_columns = ncol(bodyplan_data),
    method = "LL",
    error = TRUE
  )
}

You can skip this fit when working with the saved results. To run new simulations interactively or in Colab, change run_simulations to TRUE in the block above before executing it. The calibration fit uses error = TRUE, as in the empirical analysis.

# Load the compact cached grids for the quick, no-simulation path.
preview_tables <- readRDS(bifrost_example_file("simulation-preview-tables"))
gic_grid <- preview_tables$tuning_grids$gic
bic_grid <- preview_tables$tuning_grids$bic

1. Choose Settings from a Search Grid

Rather than evaluate one configuration, we compare a grid of values for shift_acceptance_threshold and min_descendant_tips. Larger acceptance thresholds require stronger support for a shift; larger minimum clade sizes restrict the search to shifts subtending more tips. We compare settings within GIC and within BIC, because the two criteria penalize model complexity differently.

Simulation design

For this example, we compare three acceptance thresholds (10, 20, 30) and two minimum clade sizes (10, 20): six combinations for each IC. Each simulation draws a 250-tip subtree from the passerine phylogeny and generates 12 response traits using the empirical template. We attempt 500 datasets under each of the three scenarios.

Each simulated dataset is analyzed under all six combinations with both GIC and BIC. This is a paired design: settings are compared using the same trees, trait values, and true shift locations, so differences in their results do not arise simply because they received different simulated datasets. The two calls below use matching simulation inputs and the same seed to preserve this pairing.

In the two scenarios with shifts, we place five shifts in clades containing 10–40 tips. The 10-tip search cutoff can test all of these true shift locations; the 20-tip cutoff excludes those in 10–19-tip clades. This lets us examine the tradeoff between finer phylogenetic resolution and a more restricted candidate pool. Simulated searches use error = FALSE, while the empirical calibration fit uses error = TRUE.

The code below shows how to run this design with the package functions. The saved results were produced with the separate reproduction scripts, which add checkpoints for longer runs. Both use master seed 5; the saved summaries also record the derived seeds for each scenario.

Replicate-level results are also available for exploring variation and paired comparisons; the accompanying notes explain the fields and aggregation.

# Define one shared optional grid and run it only after explicit opt-in.
full_grid_args <- list(
  shift_acceptance_thresholds = c(10, 20, 30),
  min_descendant_tips_values = c(10, 20),
  tree_tip_count = 250,
  null_replicates = 500,
  recovery_replicates = 500,
  null_simulation_options = list(
    simulation_generator = "empirical"
  ),
  proportional_simulation_options = list(
    simulation_generator = "empirical",
    num_shifts = 5,
    min_shift_tips = 10,
    max_shift_tips = 40,
    scale_factor_range = c(0.1, 2.0),
    exclude_range = c(0.5, 1.5),
    buffer = 3
  ),
  correlation_simulation_options = list(
    simulation_generator = "empirical",
    integration_power_range = c(0.5, 1.25),
    integration_exclude_range = c(0.8, 1.1)
  ),
  base_search_options = list(
    formula = "trait_data ~ 1",
    method = "LL",
    error = FALSE,
    ic_uncertainty_threshold = 10,
    num_cores = 1,
    plot = FALSE,
    store_model_fit_history = FALSE
  ),
  weighted = TRUE,
  num_cores = 8,
  seed = 5,
  store_studies = FALSE
)

if (run_simulations) {
  gic_grid <- NULL
  bic_grid <- NULL
  gic_grid <- do.call(
    runSearchTuningGrid,
    c(list(template = bodyplan_template, IC = "GIC"), full_grid_args)
  )
  bic_grid <- do.call(
    runSearchTuningGrid,
    c(list(template = bodyplan_template, IC = "BIC"), full_grid_args)
  )
}

Choose settings from the results

We first screen settings for false positives in the null simulations, where any inferred shift is a false positive. The tables show two summaries: Null FP (%) is the mean percentage of eligible candidate nodes assigned a false shift, and Null any FP (%) is the percentage of null datasets with at least one inferred shift. These answer different questions: a small percentage of falsely identified nodes can still mean that many datasets contain a false finding.

Here we allow a candidate-level false-positive rate of at most 10% and require that no more than 5% of null datasets contain any inferred shift. We also require sufficient candidate availability (min_evaluable_fraction = 0.50). This last check concerns whether performance can be evaluated on eligible nodes; it is not a search-completion check. These are criteria for choosing settings in this example, not guarantees for other datasets.

Among settings that pass these checks, we choose the highest Score: the equally weighted mean of fuzzy balanced accuracy under the proportional and integration-rate scenarios. Fuzzy evaluation allows nearby inferred shifts to match true shifts, as explained in Part 1. Exact score ties are resolved first by lower null false-positive summaries, then by larger acceptance thresholds and minimum clade sizes. With allow_infeasible = FALSE, the selector stops rather than recommending settings if none passes the checks.

# Select directly from unrounded compact grids under one shared policy.
tuning_policy <- list(
  max_false_positive_rate = 0.10,
  max_any_false_positive = 0.05,
  min_evaluable_fraction = 0.50,
  primary_metric = "fuzzy_balanced_accuracy",
  scenario_weights = c(proportional = 0.50, correlation = 0.50),
  tie_break = "conservative",
  allow_infeasible = FALSE
)
gic_tuned <- do.call(
  selectTunedSearchParameters,
  c(list(tuning_grid = gic_grid), tuning_policy)
)
bic_tuned <- do.call(
  selectTunedSearchParameters,
  c(list(tuning_grid = bic_grid), tuning_policy)
)
stopifnot(!gic_tuned$used_all_settings)
stopifnot(!bic_tuned$used_all_settings)

Interpreting the results

All 500 datasets per scenario were generated successfully, and all twelve searches completed on each dataset: 1,500 datasets and 18,000 searches in total. Of the 2,500 planted shifts in each shifted scenario, 1,411 proportional shifts and 1,467 integration-rate shifts were in 10–19-tip clades; the remainder were in 20–40-tip clades.

Each row below represents one threshold and minimum-clade-size combination. Prop. BA and Int.-rate BA report fuzzy balanced accuracy for the two shifted scenarios. Selected identifies the chosen setting, Eligible identifies another setting that passed the checks, and Excluded identifies one that did not.

GIC tuning grid

Table 1. GIC search performance across acceptance thresholds and minimum clade sizes. All settings were evaluated on the same simulated datasets.
Threshold Min clade Null FP (%) Null any FP (%) Prop. BA Int.-rate BA Score Status
10 10 0.116 6.4 0.917 0.833 0.875 Excluded
20 10 0.000 0.0 0.853 0.770 0.812 Selected
30 10 0.000 0.0 0.793 0.741 0.767 Eligible
10 20 0.120 3.8 0.770 0.723 0.747 Eligible
20 20 0.000 0.0 0.727 0.687 0.707 Eligible
30 20 0.000 0.0 0.692 0.671 0.681 Eligible

BIC tuning grid

Table 2. BIC search performance on the same simulated datasets, using the same selection criteria as GIC.
Threshold Min clade Null FP (%) Null any FP (%) Prop. BA Int.-rate BA Score Status
10 10 0.007 0.4 0.888 0.799 0.844 Selected
20 10 0.000 0.0 0.818 0.751 0.785 Eligible
30 10 0.000 0.0 0.764 0.727 0.746 Eligible
10 20 0.000 0.0 0.749 0.707 0.728 Eligible
20 20 0.000 0.0 0.709 0.680 0.694 Eligible
30 20 0.000 0.0 0.677 0.661 0.669 Eligible

Selected tuning settings

Table 3. Search settings selected separately for GIC and BIC using the false-positive and candidate-availability checks above.
IC Threshold Min clade Null FP (%) Null any FP (%) Prop. BA Int.-rate BA Score Status
GIC 20 10 0.000 0.0 0.853 0.770 0.812 Selected
BIC 10 10 0.007 0.4 0.888 0.799 0.844 Selected

For GIC, the selected acceptance threshold is 20 and the minimum clade size is 10. The combination of threshold 10 and minimum clade size 10 has a higher Score, but it infers shifts in 6.4% of null datasets, exceeding our 5% limit. BIC selects threshold 10 and minimum clade size 10, with shifts inferred in 0.4% of null datasets.

The selected GIC settings match those used by Berv et al. (2026), providing additional support for that choice under the simulation scenarios and false-positive limits considered here.

No shifts were inferred in the null datasets at acceptance thresholds of 20 or 30. Raising the threshold reduced recovery in the shifted scenarios, however, and zero observed false positives do not guarantee zero risk on another dataset.

Both selected settings use a 10-tip cutoff. At these settings, strict recall, which requires an exact node match, was 0.572 and 0.438 for GIC, and 0.615 and 0.475 for BIC (proportional and integration-rate scenarios, respectively).

The alternative 20-tip cutoff cannot recover the exact locations of shifts in smaller clades. Those shifts remain in the recovery denominator, although fuzzy evaluation can credit a nearby eligible node as a match. The minimum clade size should therefore reflect both the groups of biological interest and the data needed to estimate evolutionary parameters reliably.

The following example passes the selected GIC settings to searchOptimalConfiguration(), restoring the empirical formula with body mass as a covariate and error = TRUE. This demonstrates how to transfer the settings; the 250-tip simulations do not establish optimal settings for the full passerine phylogeny. Set run_empirical_search to TRUE only when you are ready to run the full search:

# Opt in before applying the tuned controls to the empirical dataset.
run_empirical_search <- FALSE
if (run_empirical_search) {
  empirical_result <- do.call(
    searchOptimalConfiguration,
    utils::modifyList(
      gic_tuned$recommended_search_options,
      list(
        formula = formula_str,
        error = TRUE,
        baseline_tree = bird_tree,
        trait_data = bodyplan_data,
        num_cores = 20,
        store_model_fit_history = TRUE,
        verbose = TRUE
      )
    )
  )
}

Intercept-Only Data

For datasets without a predictor, such as the GPA-aligned landmarks in the jaw-shape vignette, the calibration model uses the full trait matrix as its response. The same workflow applies, but the template formula is simpler. In this illustration, replace fish_tree and fish_data with your aligned tree and trait matrix before enabling the example:

# Opt in after replacing the illustrative objects with aligned data.
run_intercept_only_example <- FALSE
if (run_intercept_only_example) {
  jaw_template <- createSimulationTemplate(
    baseline_tree = fish_tree,
    trait_data = fish_data,
    formula = "trait_data ~ 1",
    method = "H&L",
    error = TRUE
  )
}

The downstream null, proportional, and integration-rate studies are then identical in structure; only the template and search options change.

Practical Takeaways

The aim of tuning is to find settings that recover shifts while limiting false positives on datasets resembling your own. Fit an empirical template, compare search settings under null and shifted scenarios, and choose settings within the IC you intend to use. Then apply those controls with the appropriate empirical model formula and fitting options.

When adapting this workflow, consider how tree size, trait dimensionality, covariance structure, and the size and magnitude of true shifts may affect performance.

For the definitions and performance interpretation of the three simulation scenarios, return to Part 1.

References

If you use these data or reproduce this workflow, the most relevant citations are:

  • Berv, Jacob S., Charlotte M. Probst, Santiago Claramunt, J. Ryan Shipley, Matt Friedman, Stephen A. Smith, David F. Fouhey, and Brian C. Weeks. 2026. “Rates of passerine body plan evolution in time and space.” Nature Ecology & Evolution. https://doi.org/10.1038/s41559-026-03110-5
  • Berv, Jacob S., Charlotte M. Probst, Santiago Claramunt, J. Ryan Shipley, Matt Friedman, Stephen A. Smith, David F. Fouhey, and Brian C. Weeks. 2026. “Supplementary data archive for Rates of passerine body plan evolution in time and space” (v1.0.0) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.19198393
  • Clavel, Julien, Laurent Aristide, and Helene Morlon. 2019. “A penalized likelihood framework for high-dimensional phylogenetic comparative methods and an application to New-World monkeys brain evolution.” Systematic Biology 68(1):93-116. https://doi.org/10.1093/sysbio/syy045
  • Clavel, Julien, and Helene Morlon. 2020. “Reliable phylogenetic regressions for multivariate comparative data: Illustration with the MANOVA and application to the effect of diet on mandible morphology in phyllostomid bats.” Systematic Biology 69(5):927-943. https://doi.org/10.1093/sysbio/syaa010
  • Clavel, Julien, Gilles Escarguel, and Gildas Merceron. 2015. “mvMORPH: an R package for fitting multivariate evolutionary models to morphometric data.” Methods in Ecology and Evolution 6(11):1311-1319. https://doi.org/10.1111/2041-210X.12420
  • Smith, Stephen A., Nicholas Walker-Hale, and Caroline T. Parins-Fukuchi. 2023. “Compositional shifts associated with major evolutionary transitions in plants.” New Phytologist 239(6):2404-2415. https://doi.org/10.1111/nph.19099

Software Used in This Vignette

  • bifrost for empirical simulation templates, null and shifted datasets, repeated simulation studies, recovery evaluation, and tuning-grid summaries.
  • mvMORPH/mvgls machinery underlies the global calibration fit and empirical shift searches.
  • ape and phytools underlie the phylogenetic tree, SIMMAP, and subtree operations used by the simulation helpers.
  • knitr and rmarkdown render the vignette and preview tables.

AI Assistance

This vignette was developed with assistance from OpenAI tools for drafting, editing, and figure refinement; all scientific content, interpretation, and final decisions were reviewed by the authors.