Calculation Update

Priscilla Lau, Bjørn Kopperud, Sebastian Höhna, Jeremy M. Beaulieu, and Brian C. O’Meara

03/27/2026

There was an issue in OUMA and OUMVA models in OUwie versions from its origin through 2.17. This has now been fixed. Details below:

Timeline

Priscilla Lau reached out to Jeremy Beaulieu in Aug, 2025, at a conference about an issue she and her collaborators (Bjørn Kopperud and Sebastian Höhna) found in how OUwie calculates the variance-covariance matrix under the OU model when there are different $\alpha$ parameters, as well as an error in one of the equations in the original Beaulieu et al. (2012) paper. She followed up with an email on Sept 3, 2025; after a discussion of how best to handle this, we put a warning on the use of models with different $\alpha$ in OUwie (on Sept. 6, 2025). A further email from Priscilla and colleagues had additional details on the issues.

Brian spent some time trying to fix the code; Jeremy swooped in and actually fixed the code.

Rendered equation error in original Beaulieu et al. (2012)

In the appendix on page 2383, when it describes the covariation under an OUMVA model, the right-most calculation in the final term is rendered as $$\left(\sum_{\gamma=1}^{k(i,j)} \sigma_\gamma^2 \frac{e^{2s_{ij,\gamma}} - e^{2s_{ij,\gamma-1}}}{2\alpha_\gamma}\right),$$ when it should be $$\left(\sum_{\gamma=1}^{k(i,j)} \sigma_\gamma^2 \frac{e^{2 \alpha s_{ij,\gamma}} - e^{2 \alpha s_{ij,\gamma-1}}}{2\alpha_\gamma}\right),$$ as it is correctly written in Eq. A.26. This is unrelated to the issue below but still worth noting.

Variance-covariance matrix calculation error in OUwie

For model inference the way OUwie calculated the variance-covariance matrix under the OU model with multiple $\alpha$ parameters was indeed incorrect. In an Ornstein–Uhlenbeck process, shared variance decays through time as a function of the $\alpha$ parameter: higher values of $\alpha$ cause trait values to lose memory of their shared history more rapidly, making the optima toward which lineages are pulled more important than their ancestral states.

The key issue is how variance decay is handled. The recursive approach of Lau et al. (In review) builds variance step-by-step along each branch, correctly applying decay at each stage. In contrast, Eq. 13 on page 2371 gives a closed-form expression, but does not fully account for how decay accumulates across segments when $\alpha$ varies. As a result, it is not equivalent in the general multi-$\alpha$ case, although the two agree when $\alpha$ is constant (i.e., the BM1, BMS, OU1, OUM, and OUMV cases).

Adopting the corrections proposed by Priscilla Lau and colleagues, we have updated the code in OUwie version 3.01 (released March 6, 2026) to correctly calculate the variance-covariance matrix under the OU model with multiple $\alpha$ parameters. We thank Priscilla Lau, Bjørn Kopperud, and Sebastian Höhna for their careful work in identifying this issue and proposing a solution.

FAQ

How big of a difference does this make?

The exact impact will depend on the tree, the $\alpha$ values, and the data. In some cases it could be substantial. There is no predicted effect of whether this makes OUMVA or OUMA models seem to fit better or worse than they should otherwise.

Can I verify how this affects my analyses?

Yes, with OUwie() and similar functions we have added a revert.old argument that defaults to FALSE. If you set revert.old=TRUE, it will always use the old, incorrect method for calculating the variance-covariance matrix. You can compare results with revert.old=TRUE and revert.old=FALSE to see how much of a difference it makes for your specific analysis. We wanted to make sure that any problems are transparent and reproducible.

Was hOUwie affected?

Yes, in the same way for the multiple $\alpha$ models only, and corrected in the same way.

Can I see more details on the changes?

Yes, compare, for example, the code around line 56 of varcov.ou.R with the code around line 166 of the original calculation for trees with stochastic maps

for(regimeindex in 1:length(currentmap)){
  dt <- as.numeric(currentmap[regimeindex])
    regimenumber <- which(colnames(phy$mapped.edge) == names(currentmap)[regimeindex])
    a <- alpha[regimenumber]
    s2 <- sigma[regimenumber]
    nodevar1[i] <- nodevar1[i] + a * dt
    nodevar2[i] <- nodevar2[i] + s2 * exp(2 * Acur) * (exp(2 * a * dt) - 1) / (2 * a)
    Acur <- Acur + a * dt
}
Ato[desc] <- Acur

for (regimeindex in 1:length(currentmap)){
  regimeduration <- currentmap[regimeindex]
  newtime <- oldtime+regimeduration
  regimenumber <- which(colnames(phy$mapped.edge)==names(currentmap)[regimeindex])
    nodevar1[i] <- nodevar1[i]+alpha[regimenumber]*(newtime-oldtime)
    nodevar2[i] <- nodevar2[i]+sigma[regimenumber]*((exp(2*alpha[regimenumber]*newtime)-exp(2*alpha[regimenumber]*oldtime))/(2*alpha[regimenumber]))
    oldtime <- newtime
  newregime <- regimenumber
}

Some of the variables have changed names, but a key change is in how the time accumulates along the tree, with it being scaled by alpha in the new code.