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

More importantly 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 derivation presented in Beaulieu et al. (2012) is correct; however, the appendix does not explicitly describe how the exponential decay term should accumulate when $\alpha$ varies across regimes. In this situation, the decay of shared evolutionary history must accumulate along lineage segments using the $\alpha$ associated with each regime. For developers, the decay factor is computed as $\exp\left(-\sum_{\gamma} \alpha_\gamma \Delta t_\gamma\right)$, where $\Delta t_\gamma$ is the duration of regime $\gamma$ along the lineage. Earlier versions of OUwie effectively applied a single $\alpha$ when calculating this decay while traversing the tree. As a result, the variance-covariance matrix and weight vectors were calculated incorrectly for models that allow multiple $\alpha$ parameters (i.e., OUMA and OUMVA models).

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.

Are the equations in the Beaulieu et al. (2012) paper correct?

Other than the one noted above, yes. So other software that implemented the equations correctly should be fine, this is just an issue with how OUwie implemented the calculations.

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 58 of weight.mat.R with the code around line 174 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]
    nodevar.root.tot[i] <- nodevar.root.tot[i] - a * dt
    if(regimenumber == j){
        nodevar.k[i] <- nodevar.k[i] + (exp(Acur + a * dt) - exp(Acur))
    }
    Acur <- Acur + a * dt
}

for (regimeindex in 1:length(currentmap)){
    regimeduration <- currentmap[regimeindex]
    newtime <- oldtime + regimeduration
    regimenumber <- which(colnames(phy$mapped.edge)==names(currentmap)[regimeindex])
    if(regimenumber == j){
        nodevar.root.tot[i] <- -alpha[regimenumber]*(newtime-oldtime)
        nodevar.k[i] <- exp(alpha[regimenumber]*newtime)-exp(alpha[regimenumber]*oldtime)
    }else{
        nodevar.root.tot[i] <- -alpha[regimenumber]*(newtime-oldtime)
        nodevar.k[i] <- nodevar.k[i] + 0
    }
    oldtime <- newtime
}

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.