Search for a Multi-Regime Shift Configuration on a Phylogeny

Uses a greedy, single-pass search based on multivariate mvgls() fits. One-shift candidates are ranked once and then tested in that fixed order. The selected configuration is not guaranteed to be globally optimal. Optional per-shift IC weights compare the final model with refits that each omit one accepted shift.

The function fits the supplied response variables directly; it does not perform PCA. Models assume multi-rate Brownian motion with proportional VCV scaling across regimes.

Usage

searchOptimalConfiguration(
  baseline_tree,
  trait_data,
  formula = "trait_data ~ 1",
  min_descendant_tips = 10,
  num_cores = 2,
  ic_uncertainty_threshold = 1,
  shift_acceptance_threshold = 20,
  uncertaintyweights = FALSE,
  uncertaintyweights_par = FALSE,
  plot = FALSE,
  IC = "GIC",
  store_model_fit_history = TRUE,
  verbose = FALSE,
  ...,
  progress = TRUE,
  start_strategy = c("legacy_1.2.1", "native")
)

Arguments

baseline_tree: A rooted phylo (or SIMMAP/phylo) object representing the starting tree. It does not need to already be painted: existing SIMMAP states and within-edge segments are discarded before candidate generation, and the function internally paints a single baseline state at the root. Search candidates represent cladogenic shifts painted at descendant nodes. Tip labels must match trait_data.
trait_data: A matrix or data.frame of continuous trait values with row names matching baseline_tree$tip.label (same order). For the default formula = "trait_data ~ 1", trait_data is typically supplied as a numeric matrix, but a numeric response-only data.frame is also accepted. When using more general formulas (e.g., pGLS-style models), a data.frame with named columns can be used instead.
formula: Character string or formula object passed to mvgls. Defaults to "trait_data ~ 1", which fits an intercept-only model treating the supplied multivariate trait matrix as the response. This is the appropriate choice for most morphometric data where there are no predictor variables. For more general models, formula can reference subsets of trait_data explicitly, for example "trait_data[, 1:5] ~ 1" to treat columns 1-5 as a multivariate response, "trait_data[, 1:5] ~ trait_data[, 6]" to fit a multivariate pGLS with an indexed predictor, or cbind(y1, y2) ~ size + grp to fit a named-column pGLS with numeric or factor predictors. Indexed predictor terms must each select one column: use trait_data[, 3] + trait_data[, 4] rather than trait_data[, 3:4]. Named formulas such as cbind(y1, y2) ~ x1 + x2 are preferred for multiple predictors. An intercept is included by default; use cbind(y1, y2) ~ 0 + size or cbind(y1, y2) ~ size - 1 to omit it. With a factor predictor, cbind(y1, y2) ~ 0 + grp estimates a separate mean for each group rather than fixing group means at zero.
min_descendant_tips: Integer ($\ge$2). Minimum number of tips required for an internal node to be considered as a candidate shift (forwarded to generatePaintedTrees). Defaults to 10, the value evaluated and used in the focal analysis of Berv et al. (2026). Larger values reduce the number of candidate shifts by excluding very small clades. Smaller values trigger a runtime advisory and should be assessed for the dataset at hand. This filter is applied when the initial candidate pool is generated. It is not reapplied to each regime after shifts are accepted.
num_cores: Integer. Maximum number of simultaneous fits during candidate scoring and parallel IC-weight estimation; defaults to 2. The forward pass and serial IC weights run one fit at a time. See "Parallelization" for backend and numerical-thread behavior.
ic_uncertainty_threshold: Numeric ($\ge$0). Reserved for future development in post-search pruning and uncertainty analysis; currently not used by searchOptimalConfiguration().
shift_acceptance_threshold: Numeric ($\ge$0). Minimum IC improvement (current_best_ic - proposal_ic) required to accept a candidate during the forward search. Larger thresholds require more improvement at each decision. Because the search is greedy, the final number and locations of shifts need not change monotonically with the threshold. Defaults to 20, matching the focal $\Delta$GIC setting of Berv et al. (2026); their GIC and BIC simulations also evaluated $\Delta$IC = 10. Values below 10 trigger a runtime advisory. These values are empirical reference points; users should explore alternative thresholds for their own datasets. The threshold must be a single finite nonnegative numeric value.
uncertaintyweights: Logical. If TRUE, compute per-shift IC weights serially by refitting the final model with each shift removed in turn. Exactly one of uncertaintyweights or uncertaintyweights_par must be TRUE to trigger IC-weight calculations; setting both to TRUE will result in an error. When enabled, the per-shift weights are returned in the $ic_weights component of the result.
uncertaintyweights_par: Logical. As above, but compute per-shift IC weights in parallel using future.apply. Exactly one of uncertaintyweights or uncertaintyweights_par must be TRUE to trigger IC-weight calculations.
plot: Logical. If TRUE, draw/update a SIMMAP plot as the search proceeds (requires phytools).
IC: Character. Which information criterion to use, one of "GIC" or "BIC" (case-sensitive).
store_model_fit_history: Logical; TRUE by default. Retain each forward-proposal fit, its IC, acceptance decision, and any recorded error. Fits are written to temporary files during the search and loaded into the returned object afterward, which can require substantial storage and memory. Initial one-shift scoring fits and drop-one IC-weight refits are not included.
verbose: Logical. If TRUE, emit detailed candidate-generation, fitting, acceptance, rejection, and IC messages. By default, these details are emitted via message(). When plot = TRUE in an interactive RStudio session, they are written via cat() so they remain visible while plots are updating. This is independent of progress.
...: Additional arguments passed to mvgls (e.g., method, penalty, target, error, or REML). Setting error = TRUE asks mvgls() to estimate a nuisance measurement-error (intraspecific-variance) term. See the method guidance below.
progress: Logical. If TRUE (default), show three persistent CLI progress lines for candidate scoring, the greedy shift search, and optional IC-weight re-estimation. On dynamic terminals, the spinner redraws continuously while a model fit is running; counts, percentages, and ETA advance only after completed fits. To keep the main process responsive for these redraws, otherwise serial fits run one at a time in a background Future without imposing a one-thread numerical limit. Set to FALSE to disable these lines and the heartbeat path. Both settings preserve the same fitting random sequence for the same seed and search settings. Use progress = FALSE, verbose = FALSE to suppress progress and verbose messages; warnings and errors are still reported.
start_strategy: Starting-value policy for every baseline, candidate, forward-search, and IC-weight fit. The default, "legacy_1.2.1", computes data-dependent BM/BMM starts using the mvMORPH 1.2.1 development strategy used for the cached vignettes. This includes changes after the CRAN 1.2.1 release, particularly the random second-tip safeguard for singleton regimes. It passes these starting values to the installed mvgls(). Set "native" to use the installed mvMORPH initialization unchanged. Explicit start values or grid.search = FALSE in ... take precedence. EmpBayes requires "native". This compatibility policy preserves the older initialization, including its dependence on trait scale and random state. To reproduce an older sequential search, use the same seed, preceding searches, and serial/parallel settings. Turning progress on or off does not change the random sequence. This policy does not guarantee convergence or identical results across numerical environments or future mvMORPH versions.

Value

A named list with (at minimum):

user_input: captured call expressions, including resolved progress and start_strategy settings and flattened additional mvgls() arguments. This does not capture the seed, a data snapshot, or the full software environment.
initialization: requested and actual starting policy for the baseline fit (legacy_1.2.1, native, supplied, or native_no_grid), plus the installed mvMORPH version. Each fitted model also records its own initialization information in its bifrost_initialization attribute.
tree_no_uncertainty_transformed: SIMMAP tree from the final model on the transformed scale used internally by mvgls.
tree_no_uncertainty_untransformed: same topology with original edge lengths restored.
model_no_uncertainty: the final mvgls model object.
shift_nodes_no_uncertainty: integer vector of accepted shift nodes; empty (possibly NULL) when no shifts are accepted.
optimal_ic: final IC value; baseline_ic: baseline IC.
IC_used: "GIC" or "BIC"; num_candidates: count of candidate one-shift models evaluated; candidate_nodes: integer node identifiers for those candidates.
model_fit_history: if store_model_fit_history = TRUE, a list of per-iteration fits (loaded from temporary files written during the run) with proposal step, candidate node, regime ID, IC, status, acceptance decision, and any error, plus an ic_acceptance_matrix (IC value and acceptance flag per step).
VCVs: named list of regime-specific VCV matrices extracted from the final model (penalized-likelihood estimates if PL was used). For a baseline-only BM fit, the sole entry "0" is its fitted global covariance matrix.

Additional components appear conditionally:

ic_weights: a data frame with one row per accepted shift when uncertaintyweights or uncertaintyweights_par is TRUE. Each row compares the final model (ic_with_shift) with a refit omitting that shift (ic_without_shift). Here, delta_ic = ic_with_shift - ic_without_shift, so negative values favor retaining the shift. The two IC weights sum to one within each row; they are not normalized across shifts. evidence_ratio is ic_weight_withshift / ic_weight_withoutshift.
warnings: when present, a list of warning and error messages collected during forward-proposal evaluation. It is not a complete record of warnings from every fitting stage or of optimizer convergence codes.

Details

Input requirements.

Tree: baseline_tree should be a rooted phylo tree with branch lengths expressed on a common scale. An ultrametric tree is not required. The starting tree does not need to already be painted; searchOptimalConfiguration() paints a single baseline regime internally before building shifted candidates.
Trait data alignment: rownames(trait_data) must match baseline_tree$tip.label in both names and order; any tips without data should be pruned beforehand.
Data type: trait_data is typically a numeric matrix of continuous traits; numeric response-only data.frames are also supported for intercept-only searches, and named mixed-type data.frames are supported for richer formulas. High-dimensional settings (p $\ge$ n) are supported via penalized-likelihood mvgls() fits.

Search outline. Fit the single-regime baseline and each eligible one-shift candidate. Rank candidates once by their IC improvement over the baseline. Visit that fixed list once, fitting each proposal on the current accepted tree. Accept it when its improvement over the current accepted model meets shift_acceptance_threshold. Rejected candidates are not revisited after later acceptances. Optional IC weights compare the final model with a refit omitting each accepted shift, using aicw.

If no eligible non-root candidates remain, or no proposed shifts are accepted, the returned final model is the single-regime baseline and optimal_ic equals baseline_ic.

Search-setting diagnostics. min_descendant_tips filters candidate nodes using their original number of descendant tips. It is applied only when the initial candidate pool is built; it is not reapplied to the resulting regimes after nested shifts are accepted. Consequently, a resulting regime can contain fewer tips than the threshold. An eligible small clade is evaluated like any other candidate and must pass shift_acceptance_threshold. Under Brownian motion, high trait disparity concentrated on short branches implies a high evolutionary-rate estimate, but a runtime rule cannot determine whether that pattern is biological or unstable. The permissive-settings advisory is therefore a guardrail rather than a power calculation. Examine per-shift IC weights by enabling uncertaintyweights or uncertaintyweights_par, repeat searches across plausible settings, and use dataset-specific simulations when conclusions depend on individual shifts.

Parallelization. num_cores caps simultaneous model fits: candidate scoring and parallel IC-weight re-estimation may use all requested workers, while the greedy search and serial IC weights always run one fit at a time. With progress enabled, fits run in Future workers so the main process can poll and redraw the spinner without advancing completion. Where forking is supported on Unix outside RStudio, multicore is used; otherwise multisession is used. Numerical-thread limits depend on the number of fitting workers, not on progress. With multiple workers, thread-limit environment variables and RhpcBLASctl runtime controls request one numerical thread in each worker. Single-worker stages preserve the caller's thread-limit environment, allowing within-fit multithreading when supported by the numerical library. For consistent behavior across fresh R sessions, configure numerical-thread environment variables before starting R; runtime-only changes in the main process are not automatically copied to multisession workers. Stage workers are shut down and the previous thread-limit environment and Future plan are restored afterward. Workers signal conditions through progressr; only the main R process renders progress, verbose messages, and plots.

Plotting. If plot = TRUE, trees are rendered with plotSimmap(); shift IDs are labeled with nodelabels().

Regime VCVs. For a multi-regime model, the returned $VCVs contain regime-specific covariance estimates extracted via extractRegimeVCVs. When no shifts are accepted, the final model is single-regime "BM": its trees retain baseline label "0", and VCVs[["0"]] contains the global covariance matrix from model_no_uncertainty$sigma$Pinv, including regularization when used. This entry does not represent an inferred shift or an additional fit; the model's param field remains NA.

Fitting methods. Choose an mvgls method appropriate for the formula and data. "H&L" with the "RidgeArch" penalty is restricted to intercept-only models. For models with predictors, "LL" is one option when an unregularized covariance can be estimated; "PL-LOOCV" is a penalized alternative. When IC = "BIC" and method is omitted, bifrost supplies method = "LL". For IC = "GIC", an omitted method is left to the installed mvMORPH default. See mvgls for method and penalty choices.

Convergence and robustness

The search is greedy and visits candidates in an initial ranked order. The selected configuration is not guaranteed to minimize the information criterion over all possible shift configurations.

Numerical fits can depend on starting values and optimizer stopping criteria with either method = "H&L" or method = "LL". For conclusions that depend on individual shifts, repeat the search with the same data, tree, formula, and search settings using several recorded seeds when initialization is stochastic. Comparisons of supported start_strategy values can provide a separate initialization-sensitivity check. Record the software versions and parallel settings, and compare the accepted shifts, information-criterion values, and fitted rates.

Inspect optimizer codes and messages for the final model and saved forward proposals when available. A fit with a nonzero optimizer convergence code can still contribute an IC to candidate ranking or a forward acceptance decision. A nonzero code indicates that the optimizer reported a warning or failure; code 0 does not establish a global optimum. Inspect optimizer codes and messages rather than relying on warnings alone. Agreement across repeated runs is useful evidence of stability, but changing the seed may not change a deterministic initialization.

Compact example objects may omit the forward-proposal fits.

Separately, vary shift_acceptance_threshold, min_descendant_tips, or the information criterion to assess sensitivity to search settings. Per-shift IC weights are conditional on the selected fitted model; they do not assess convergence or alternative search paths.

References

Berv, J. S. et al. (2026). Rates of passerine body plan evolution in time and space. Nature Ecology & Evolution. doi:10.1038/s41559-026-03110-5.

Examples

library(ape)
library(phytools)
#> Loading required package: maps
#> 
#> Attaching package: ‘maps’
#> The following object is masked from ‘package:viridis’:
#> 
#>     unemp
library(mvMORPH)
#> Loading required package: corpcor
#> Loading required package: subplex
#> ##
#> ## mvMORPH package (1.2.3)
#> ## On est deja grands : ................. ##
#> ##               ,--.  ,--.               ##
#> ##              /   | /    \              ##
#> ##              `|  ||  ()  |             ##
#> ##               |  | \    /              ##
#> ##               `--'  `--'               ##
#> ##                _|____|_                ##
#> ##               |        |               ##
#> ##            +--------------+            ##
#> ##            |              |            ##
#> ##        +----------------------+        ##
#> ##        |                      |        ##
#> ##   + ------------------------------ +   ##
#> ##            __  __  ___  ___ ___ _  _   ##
#> ##   _ ____ _|  \/  |/ _ \| _ \ _ \ || |  ##
#> ##  | '  \ V / |\/| | (_) |   /  _/ __ |  ##
#> ##  |_|_|_\_/|_|  |_|\___/|_|_\_| |_||_|  ##
#> ##... et ce n'est que le debut :)         ##
#> ## Multivariate evolutionary models
#> ##
#> ## See the tutorials: browseVignettes("mvMORPH")
#> ##
#> ## To cite package 'mvMORPH': citation("mvMORPH")
#> ##
set.seed(1)

# Simulate a tree
tr <- pbtree(n = 50, scale = 1)

# Define baseline regime "0" and high-rate regime "1" on a subset of tips
states <- setNames(rep("0", Ntip(tr)), tr$tip.label)
high_clade_tips <- tr$tip.label[1:20]
states[high_clade_tips] <- "1"

# Make a SIMMAP tree for the BMM simulation
simmap <- phytools::make.simmap(tr, states, model = "ER", nsim = 1)
#> make.simmap is sampling character histories conditioned on
#> the transition matrix
#> 
#> Q =
#>            0          1
#> 0 -0.3147857  0.3147857
#> 1  0.3147857 -0.3147857
#> (estimated using likelihood);
#> and (mean) root node prior probabilities
#> pi =
#>   0   1 
#> 0.5 0.5 
#> Done.

# Simulate traits under a BMM model with ~10x higher rate in regime "1"
sigma <- list(
  "0" = diag(0.1, 2),
  "1" = diag(1.0, 2)
)
theta <- c(0, 0)

sim <- mvMORPH::mvSIM(
  tree  = simmap,
  nsim  = 1,
  model = "BMM",
  param = list(
    ntraits = 2,
    sigma   = sigma,
    theta   = theta
  )
)

# mvSIM returns either a matrix or a list of matrices depending on mvMORPH version
X <- if (is.list(sim)) sim[[1]] else sim
rownames(X) <- simmap$tip.label

# Run the search on the unpainted single-regime tree
res <- searchOptimalConfiguration(
  baseline_tree              = as.phylo(simmap),
  trait_data                 = X,
  formula                    = "trait_data ~ 1",
  min_descendant_tips        = 10,
  num_cores                  = 1,   # keep it simple / CRAN-safe
  shift_acceptance_threshold = 20,  # conservative GIC threshold
  IC                         = "GIC",
  plot                       = FALSE,
  store_model_fit_history    = TRUE,
  verbose                    = FALSE,
  progress                   = FALSE
)

res$shift_nodes_no_uncertainty
#> [1] 54
res$baseline_ic - res$optimal_ic  # total IC improvement
#> [1] 38.3827
res$model_no_uncertainty$opt[c("convergence", "message")]
#> $convergence
#> [1] 0
#> 
#> $message
#> [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
#> 
head(icTrajectory(res))           # forward proposals and decisions
#>   step        ic accepted   best_ic   delta_ic   status candidate_node
#> 1    0 34.741374       NA 34.741374         NA baseline             NA
#> 2    1 -3.641323     TRUE -3.641323 38.3826973 accepted             54
#> 3    2 -2.695849    FALSE -3.641323 -0.9454740 rejected             53
#> 4    3 -5.053236    FALSE -3.641323  1.4119126 rejected             90
#> 5    4 -3.395075    FALSE -3.641323 -0.2462477 rejected             52
#> 6    5 -2.277586    FALSE -3.641323 -1.3637371 rejected             79
#>   regime_id error
#> 1         0  <NA>
#> 2         1  <NA>
#> 3         2  <NA>
#> 4         3  <NA>
#> 5         4  <NA>
#> 6         5  <NA>

# \donttest{
# Intercept-only empirical-style search:
# H&L + measurement error
res_hl <- searchOptimalConfiguration(
  baseline_tree              = as.phylo(simmap),
  trait_data                 = X,
  formula                    = "trait_data ~ 1",
  min_descendant_tips        = 10,
  num_cores                  = 1,
  shift_acceptance_threshold = 20,
  IC                         = "GIC",
  plot                       = FALSE,
  method                     = "H&L",
  error                      = TRUE,
  store_model_fit_history    = FALSE,
  verbose                    = FALSE,
  progress                   = FALSE
)

# Formula-based search with a predictor:
# LL is one option for a model with predictors
dat <- data.frame(
  trait1    = X[, 1],
  trait2    = X[, 2],
  predictor = rnorm(nrow(X))
)
rownames(dat) <- simmap$tip.label

res_ll <- searchOptimalConfiguration(
  baseline_tree              = as.phylo(simmap),
  trait_data                 = dat,
  formula                    = "trait_data[, 1:2] ~ trait_data[, 3]",
  min_descendant_tips        = 10,
  num_cores                  = 1,
  shift_acceptance_threshold = 20,
  IC                         = "GIC",
  plot                       = FALSE,
  method                     = "LL",
  error                      = TRUE,
  store_model_fit_history    = FALSE,
  verbose                    = FALSE,
  progress                   = FALSE
)
# }