Skip to contents

Empirically Calibrated Simulations for bifrost, Part 1: Performance

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

This is Part 1 of the empirically calibrated simulation workflow. Part 2 uses these performance scenarios to tune search settings and apply them to an empirical analysis.

As part of a bifrost workflow, it helps to know how the method behaves on datasets that resemble the one you actually care about. A practical way to answer that question is to calibrate simulations to the empirical dataset. Here, calibration means parameterizing simulations to match the empirical tree, trait dimensionality, fitted mean structure, and residual covariance scale.

We then compare search behavior under three scenarios:

  1. a null model with no shifts, to estimate expected false positives,
  2. the current generating model, where shifts are proportional covariance scalings, and
  3. a deliberately non-generating robustness scenario coupling reduced integration with faster marginal rates, or increased integration with slower marginal rates.

Empirically calibrated simulations can be used before or after an empirical search. Before the final search, they help choose search settings; after a search, they help evaluate whether those settings behave reasonably for the empirical tree, trait dimensionality, fitted mean structure, and residual covariance scale. In both cases, the workflow uses the same ingredients: a global calibration template, null and shifted simulation scenarios, false-positive and recovery summaries, and a way to compare candidate search settings. This is the same broad calibration strategy used in Berv et al. (2026). On large trees, repeated subtree resampling is helpful because it lets simulations vary across many local topologies and parameter draws without requiring a full-tree search on every replicate.

The package generator improves on the manuscript implementation for new simulation studies by drawing ancestral covariance structures around the complete fitted empirical covariance and by representing integration-rate shifts in both biological directions. This refinement provides a more direct robustness test and does not alter the manuscript’s empirical analyses or conclusions.

Start by loading bifrost.

# Load bifrost before constructing templates or running simulation studies.
library(bifrost)

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.

# Locate the packaged passerine tree and calibrated skeletal measurements.
tree_path <- pkg_file("extdata", "avian-skeleton", "passerine_bodyplan_tree.tre")
trait_path <- pkg_file("extdata", "avian-skeleton", "passerine_bodyplan_data.RDS")

# 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 comes from a single-regime global fit. That fit provides the empirical variance and covariance scale used to generate a fresh ancestral covariance matrix for each simulation replicate. Here the calibration model is a single-regime multivariate BM fit with mvMORPH::mvgls() using method = "LL", following the multivariate phylogenetic GLS framework described by Clavel, Aristide, and Morlon (2019) and Clavel and Morlon (2020).

# Fit the single-regime model that calibrates means and residual covariance.
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
)
#> Bifrost Simulation Template
#> ===========================
#> 
#> Tree
#>   Tips: 2057
#>   Response traits: 12
#>   Predictor columns: 1
#> 
#> Model
#>   Calibration formula: trait_data[, 1:12] ~ trait_data[, 13]
#>   Simulation search:   trait_data ~ 1
#>   Method: LL
#>   Error term: TRUE
#> 
#> Empirical Covariance Summaries
#>   Variance mean: 0.000215
#>   Variance sd:   0.000124
#>   Covariance mean: 0.000018
#>   Covariance sd:   0.000109

The template stores the aligned calibration data, the fitted global mvgls model, fitted means, and empirical variance/covariance summaries used to generate a fresh ancestral covariance matrix for each simulation replicate. Covariates such as body mass can shape those fitted means and residual covariance, but the simulation studies regenerate only the response block and evaluate intercept-only shift searches on those simulated responses. If no calibration covariates are needed, the same workflow reduces to the simpler trait_data ~ 1 template. In all cases, response columns must be raw numeric variables that can be regenerated in each replicate.

The examples below use simulation_generator = "empirical", the recommended option when the fitted full covariance structure should anchor new draws.

Simulation scenarios

The figure below places the three package scenarios side by side: an empirical ancestral covariance matrix, a proportional rescaling, and the integration-rate trade-off. The last scenario raises the eigenvalues of the ancestral correlation matrix to a power g, renormalizes to a correlation matrix, and multiplies the covariance matrix by 1 / g. Thus g < 1 reduces integration while increasing every marginal evolutionary variance; g > 1 increases integration while decreasing every marginal variance.

Schematic of the empirically calibrated null, proportional, and integration-rate simulation scenarios.

Figure 1. Empirically calibrated simulation scenarios. Panels A-C show the ancestral VCV, proportional scaling with s = 0.5, and a reduced-integration/faster-rate shift with g = 0.65. Panels D-F show their correlation matrices. Panels G-I summarize correlation effective dimensionality, covariance eigenvalue spectra, and mean marginal evolutionary variance. The proportional scenario changes scale without changing correlation structure; the integration-rate scenario changes covariance shape and rate together.

In the proportional panels, the covariance structure keeps the same shape and simply expands or contracts. In the integration-rate panels, weaker correlation structure is paired with larger marginal variance, so covariance shape and evolutionary rate change together. The shifted passerine examples below use these two simulation styles under repeated subtree resampling from the full tree, alongside a null no-shift baseline.

Calibration guardrails

Before launching a large study, keep four details in mind.

  • Simulation studies regenerate the response block only. Predictors shape the global calibration fit, but they are not re-fit inside each subtree-resampled search replicate.

  • Use the wrapper-level seed argument for reproducibility. Do not pass a separate seed inside simulation_options; the study wrappers manage replicate-level seeds.

  • The wrappers inherit template fit settings such as method and error, so search_options only needs to include deliberate overrides. Parallelize either the study wrapper or each embedded search, not both; the examples parallelize replicates and keep each search serial.

  • Check per_replicate$n_candidates. The value of min_descendant_tips changes the candidate-node universe rather than merely applying a different decision threshold to a fixed set of nodes. If it leaves no candidate shifts, that replicate is treated as unevaluable, and recovery specificity, false-positive rate, and balanced accuracy are reported as NA. Lower min_descendant_tips if too many replicates are unevaluable.

The next three sections show the individual study calls. The fixed-settings comparison then uses bundled preview results to show the same logic side by side.

1. Estimate False Positives

The null study asks the baseline error question: if the real dataset contained no shifts at all, how often would the search still infer one?

# Estimate false positives under an empirical-covariance null model.
fp_study <- runFalsePositiveSimulationStudy(
  template = bodyplan_template,
  n_replicates = 100,
  tree_tip_count = 250,
  simulation_options = list(
    simulation_generator = "empirical"
  ),
  search_options = list(
    formula = "trait_data ~ 1",
    min_descendant_tips = 10,
    shift_acceptance_threshold = 10,
    IC = "GIC",
    num_cores = 1,
    method = "LL",
    error = TRUE,
    plot = FALSE,
    store_model_fit_history = FALSE
  ),
  num_cores = 8,
  seed = 5
)
# Inspect the study object, representative replicates, and aggregate summary.
fp_study
head(fp_study$per_replicate)
fp_study$study_summary

Null simulations are useful for screening candidate choices of shift_acceptance_threshold, IC, and min_descendant_tips.

2. Estimate Recovery Under Proportional Scaling

Recovery under proportional scaling evaluates the scenario that matches the current bifrost BMM search assumptions: shifted covariance matrices are scalar multiples of the background covariance matrix.

# Measure recovery when simulated shifts match proportional BMM scaling.
recovery_prop <- runShiftRecoverySimulationStudy(
  template = bodyplan_template,
  n_replicates = 100,
  tree_tip_count = 250,
  simulation_options = list(
    simulation_generator = "empirical",
    num_shifts = 5,
    min_shift_tips = 10,
    max_shift_tips = 20,
    scale_mode = "proportional",
    scale_factor_range = c(0.1, 2.0),
    exclude_range = c(0.5, 1.5),
    buffer = 3
  ),
  search_options = list(
    formula = "trait_data ~ 1",
    min_descendant_tips = 10,
    shift_acceptance_threshold = 10,
    IC = "GIC",
    num_cores = 1,
    method = "LL",
    error = TRUE,
    uncertaintyweights_par = TRUE,
    plot = FALSE,
    store_model_fit_history = FALSE
  ),
  fuzzy_distance = 2,
  weighted = TRUE,
  num_cores = 8,
  seed = 5
)
# Compare strict, fuzzy, and distance-weighted recovery summaries.
recovery_prop
recovery_prop$study_summary
recovery_prop$evaluation$strict
recovery_prop$evaluation$fuzzy
recovery_prop$evaluation$weighted

3. Check the Integration-Rate Robustness Scenario

It is also useful to test a deliberately mismatched scenario in which trait integration and evolutionary rate evolve together. This asks how well the same search settings behave when the data depart from strict proportional VCV scaling. The moderate power bands below sample reduced integration from 0.5-0.8 and increased integration from 1.1-1.25; each power g scales marginal variances by exactly 1 / g.

# Stress-test the same search under integration-rate model misspecification.
recovery_integration <- runShiftRecoverySimulationStudy(
  template = bodyplan_template,
  n_replicates = 100,
  tree_tip_count = 250,
  simulation_options = list(
    simulation_generator = "empirical",
    num_shifts = 5,
    min_shift_tips = 10,
    max_shift_tips = 20,
    scale_mode = "correlation",
    integration_power_range = c(0.5, 1.25),
    integration_exclude_range = c(0.8, 1.1),
    buffer = 3
  ),
  search_options = list(
    formula = "trait_data ~ 1",
    min_descendant_tips = 10,
    shift_acceptance_threshold = 10,
    IC = "GIC",
    num_cores = 1,
    method = "LL",
    error = TRUE,
    uncertaintyweights_par = TRUE,
    plot = FALSE,
    store_model_fit_history = FALSE
  ),
  fuzzy_distance = 2,
  weighted = TRUE,
  num_cores = 8,
  seed = 5
)
# Review the overall study and the fuzzy recovery result used for reporting.
recovery_integration
recovery_integration$study_summary
recovery_integration$evaluation$fuzzy

This scenario is not the current generating model for bifrost, so the goal is not perfect recovery. It deliberately changes covariance shape while the fitted search model permits only proportional covariance shifts. Lower recovery than in the proportional scenario is therefore expected and measures robustness to biologically motivated model misspecification rather than failure under the model’s own assumptions.

4. Compare Fixed Settings Across Scenarios

After running each scenario on its own, it helps to put one fixed search configuration side by side across the null, proportional, and integration-rate cases. This worked example uses 100 repeated 250-tip subtree resamples from the passerine tree.

The comparison uses one search configuration across all three scenarios: shift_acceptance_threshold = 10 and min_descendant_tips = 10. In the shifted scenarios, each replicate adds five simulated shifts spanning clades with 10 to 20 descendant tips, matching the simulation-performance settings described in Berv et al. (2026). Every scenario uses the explicit empirical generator. The point here is not to optimize yet; it is to compare how null and shifted datasets behave under the same search settings.

The table below uses bundled summaries from the code shown here, so package builds do not rerun the longer simulation studies.

# Prevent nested BLAS/OpenMP parallelism during replicate-level parallel runs.
if (requireNamespace("RhpcBLASctl", quietly = TRUE)) {
  RhpcBLASctl::blas_set_num_threads(1)
  RhpcBLASctl::omp_set_num_threads(1)
}

# Reuse the fitted empirical template for all fixed-setting comparisons.
medium_case_template <- preview_template

# Define the proportional-shift scenario shared by GIC and BIC.
fixed_setting_proportional_options <- list(
  simulation_generator = "empirical",
  num_shifts = 5,
  min_shift_tips = 10,
  max_shift_tips = 20,
  scale_factor_range = c(0.1, 2.0),
  exclude_range = c(0.5, 1.5),
  buffer = 3
)

# Hold every search and simulation setting fixed except the IC family.
fixed_setting_args <- list(
  shift_acceptance_thresholds = 10,
  min_descendant_tips_values = 10,
  tree_tip_count = 250,
  null_replicates = 100,
  recovery_replicates = 100,
  null_simulation_options = list(
    simulation_generator = "empirical"
  ),
  proportional_simulation_options = fixed_setting_proportional_options,
  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 = TRUE,
    num_cores = 1,
    uncertaintyweights_par = TRUE,
    plot = FALSE,
    store_model_fit_history = FALSE,
    progress = FALSE
  ),
  fuzzy_distance = 2,
  weighted = TRUE,
  num_cores = 8,
  seed = 5,
  store_studies = FALSE
)

# Run the identical comparison once with GIC and once with BIC.
gic_fixed <- do.call(
  runSearchTuningGrid,
  c(list(template = medium_case_template, IC = "GIC"), fixed_setting_args)
)
bic_fixed <- do.call(
  runSearchTuningGrid,
  c(list(template = medium_case_template, IC = "BIC"), fixed_setting_args)
)

Null and shifted simulations answer different questions, so they are reported separately. The first table shows null false-positive behavior; the second shows recovery only where true shifts are known. This avoids treating an inapplicable recovery metric in a null simulation as a missing result.

Table 1. Null behavior for the fixed GIC and BIC settings across 100 repeated 250-tip subsamples. Mean FP and Any FP directly summarize false-positive control.
IC Mean FP Any FP Mean shifts
GIC 0.0017 0.100 0.10
BIC 0.0007 0.040 0.04
Table 2. Shift recovery for the same fixed settings and repeated subsamples. Fuzzy scoring credits near misses within two nodes of a true shift. Balanced accuracy is the mean of fuzzy recall and fuzzy specificity; F1 is retained as the complementary precision-recall summary.
IC Scenario Fuzzy rec. Fuzzy spec. Fuzzy F1 Fuzzy BA Mean shifts
GIC Proportional 0.846 0.993 0.883 0.920 4.58
GIC Integration-rate 0.667 0.982 0.717 0.824 4.26
BIC Proportional 0.768 0.996 0.848 0.882 4.06
BIC Integration-rate 0.588 0.992 0.702 0.790 3.31

Candidate nodes without a true shift generally outnumber true-shift nodes, creating class imbalance. Raw accuracy can therefore look strong even when sensitivity is poor. Fuzzy balanced accuracy gives recall and specificity equal influence, while fuzzy matching credits biologically close near misses instead of counting every small node displacement as a complete failure. Balanced accuracy complements rather than replaces F1: balanced accuracy exposes the sensitivity-specificity trade-off, whereas F1 focuses on precision and recall and remains useful when false discoveries among inferred shifts are the main concern.

For GIC, the cache reports fuzzy balanced accuracy of 0.920 under proportional shifts and 0.824 under the integration-rate trade-off. The corresponding BIC values are 0.882 and 0.790. These cache-derived comparisons describe performance under the fitted generating model and under deliberate covariance-shape misspecification without hard-coding a particular cache result into the narrative.

Performance Summary

These checks establish false-positive behavior under the null, recovery under the proportional generating model, and robustness to the integration-rate trade-off. Part 2 uses those scenarios to tune search settings and apply them to the empirical analysis.

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.
  • ggplot2, patchwork, viridis, and reshape2 draw the schematic simulation-scenario figure.
  • 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.