hOUwie performance and stability improvements
In OUwie: Analysis of Evolutionary Rates in an OU Framework

knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE,
  fig.align = "center",
  out.width = "90%"
)

performance <- data.frame(
  scenario = c(
    "Moderate tree",
    "More histories",
    "1,000-tip tree",
    "16 states / 74 parameters",
    "500 tips / 8 states"
  ),
  tips = c(96, 96, 1000, 250, 500),
  histories = c(25, 100, 25, 25, 50),
  states = c(4, 4, 4, 16, 8),
  parameters = c(8, 8, 8, 74, 22),
  previous_seconds = c(0.335, 1.283, 4.713, 1.908, 5.513),
  updated_seconds = c(0.131, 0.463, 1.372, 0.365, 1.267),
  stringsAsFactors = FALSE
)
performance$speedup <- performance$previous_seconds /
  performance$updated_seconds

performance_table <- data.frame(
  Scenario = performance$scenario,
  Tips = performance$tips,
  Histories = performance$histories,
  States = performance$states,
  Parameters = performance$parameters,
  `Previous version (s)` = sprintf("%.3f", performance$previous_seconds),
  `Updated version (s)` = sprintf("%.3f", performance$updated_seconds),
  Speedup = sprintf("%.2fx", performance$speedup),
  check.names = FALSE
)

hOUwie now runs between r sprintf("%.2f", min(performance$speedup)) and r sprintf("%.2f", max(performance$speedup)) times faster across the benchmark set. The largest improvements appear on large trees and models with many states and parameters.

How much faster is it?

knitr::kable(
  performance_table,
  align = c("l", rep("r", 7)),
  caption = "Median elapsed time across repeated runs. Lower times are better."
)

bar_colors <- ifelse(performance$speedup >= 4, "#2b8cbe", "#7bccc4")
old_par <- par(mar = c(8, 4.2, 1, 0.5))
bars <- barplot(
  performance$speedup,
  names.arg = performance$scenario,
  las = 2,
  ylim = c(0, max(performance$speedup) * 1.18),
  ylab = "Speedup (times faster)",
  col = bar_colors,
  border = NA
)
abline(h = 1, lty = 2, col = "grey50")
text(
  bars,
  performance$speedup,
  labels = sprintf("%.2fx", performance$speedup),
  pos = 3,
  cex = 0.9
)
par(old_par)

The moderate 96-tip analysis is 2.6 times faster. Increasing the number of histories from 25 to 100 raises that to 2.8 times. On the larger and more parameter-rich cases, the updated implementation is 3.4 to 5.2 times faster. A 1,000-tip run falls from 4.713 seconds to 1.372 seconds, while the 16-state, 74-parameter case falls from 1.908 seconds to 0.365 seconds.

Every run produced the requested number of stochastic histories. The median size of the returned history objects changed by less than one percent in every case.

Where the speedup comes from

Most of the improvement comes from preparing the tree structure once and reusing it. A phylogenetic tree has the same topology throughout a fit. The updated code prepares its paths, descendants, edge positions, and state relationships up front, then keeps compact lookup tables for the calculations that follow.

Tree traversals and tip-to-root paths are prepared once instead of rebuilt for every stochastic history.
Repeated searches through edge and node vectors are replaced with direct indexing.
The likelihood calculation works with a compact representation of each history. Full painted trees are still assembled for the returned result, but they no longer need to be the working format for every intermediate step.
Matrix and cache bookkeeping is lighter, especially when the number of states and parameters grows.

Together, these changes remove a large amount of repeated work from the inner likelihood loop. The 16-state analysis benefits most: the updated implementation spends more of its time on the statistical calculation and less on repeatedly organizing the same tree.

A steadier optimization problem

hOUwie estimates a likelihood by averaging over sampled evolutionary histories. That sampling introduces a little numerical movement from one evaluation to the next. An optimizer can mistake that movement for a meaningful improvement and spend time following noise.

The updated implementation compares candidate parameter values with a shared underlying random stream within each optimizer start. The histories still respond to the candidate parameters; they are not frozen. The shared stream simply makes one candidate's likelihood more directly comparable with the next. This behavior is enabled by default and can be controlled with the common_random_numbers argument to hOUwie().

On an easy example, both approaches found the same answer. On a more difficult hidden-state model, the steadier comparison produced better and more consistent fits on fresh evaluations:

| Measure | Ordinary sampling | Shared random stream | |:--|--:|--:| | Mean fresh-evaluation log likelihood | -27.9244 | -27.2323 | | Median fresh-evaluation log likelihood | -28.0759 | -27.1515 | | Between-run standard deviation | 0.6145 | 0.2669 | | Mean likelihood evaluations | 524.0 | 438.5 | | Mean elapsed time (seconds) | 75.24 | 61.48 |

The steadier search used 16.3% fewer likelihood evaluations and 18.3% less wall time. It found the better fresh-evaluated solution in four of six paired starts. Multiple starts remain important when hidden states create weakly identified parameters or local optima.

Reliability improvements around the likelihood

The fitting path also includes the following reliability improvements:

Cached likelihoods use the full parameter vector, so distinct candidate models cannot accidentally share a result.
Cache capacity is bounded more safely for long optimization runs.
Failed optimizer results are filtered before selecting the best fit.
Trait data are matched to tree tips by species name rather than assumed row position.
History identifiers work with ten or more states.
Negative trait values are shifted consistently when a model requires positive values.
Input checks and parameter bounds reject several malformed or impossible configurations earlier and more clearly.