A cautionary note on the use of Ornstein Uhlenbeck models in macroevolutionary studies

1Phylogenetic comparative methods (PCMs) are powerful tools for identifying patterns in the evolution of species traits, and for potentially inferring the evolutionary processes that underlie them (e.g. Freckleton, 2009; Nunn, 2011; O'Meara, 2012; Pennell & Harmon, 2013). These approaches have been used, for example, to infer potential rates of species responses to climate change (Quintero & Wiens, 2013), test the role of ecological niche as a driver of morphological evolution (Pienaar et al., 2013) and test for constraints in adaptive radiations (Blackburn et al., 2013).

2The majority of PCMs use an explicit evolutionary model to characterize trait evolution (Freckleton et al. 2011). Most model‐based methods for characterizing trait evolution are based on the Brownian constant variance model (for exceptions see Price, 1997; Harvey & Rambaut, 2000; Freckleton & Harvey, 2006). The Brownian model, first applied in a phylogenetic context by Cavalli‐Sforza & Edwards (1967) and to across‐species data by Felsenstein (1973), is a simple model of trait evolution in which trait variance accrues as a linear function of time, and makes the prediction that traits of closely related species are more similar than those of distantly related ones. The Brownian model has been modified in various ways to account for a suite of ecological and evolutionary processes (e.g. Grafen, 1989; Hansen, 1997; Pagel, 1997, 1999). Most of these involve a transformation of the tree and thereby fitting a model with one or more extra parameters. These modified Brownian models tend to fit better and often have links to process‐based interpretations.

3One of the most commonly used Brownian‐like models is the Ornstein Uhlenbeck (OU) model. The OU model was introduced to population genetics by Lande (1976) to model stabilizing selection in which the trait is drawn towards a fitness optimum on an adaptive landscape. The process operating in comparative data is analogous to but distinct from stabilizing selection. The phylogenetic OU model is a modification of the Brownian model with an additional parameter α that measures the strength of return towards a theoretical optimum (Hansen, 1997) that is shared across a clade or subset of species. Although widely used, the properties of the OU model, and other direct extensions of the Brownian model, are poorly understood leading to the potential for inappropriate use and misinterpretation of results.

4In this paper we present an introduction to the OU model, its general properties and some issues with its use in ecology, evolution and palaeontology. We use simulations to demonstrate the inherent bias in estimating the core parameter of the OU model, α, that describes the strength of pull towards a central value (typically referred to as the trait or selective optima). We discuss the intricacies of interpreting OU models biologically, and provide advice for appropriate use of OU models in phylogenetic comparative analyses. We also show that very small amounts of intraspecific trait variation (including measurement error) can profoundly affect the performance of models. These findings will be applicable to other models of evolution, but we focus on the OU model because of its widespread use and because of the ambiguity in the link between pattern and process when interpreting estimates of the α parameter. We are not the first to describe some of these problems (e.g. Ives & Garland, 2010; Boettiger, Coop & Ralph, 2012; Hansen & Bartoszek, 2012; Ho & Ané, 2013, 2014). However, widespread use of the model is clear evidence that many are unaware of the potential problems. We use a simulation approach to summarize the problems and to generate practical recommendations of how to deal with them.

5The popularity of the OU model has grown extensively in recent years (Fig. 1); even just between 2012 and 2014 over 2500 ecology, evolution and palaeontology papers containing the phrase ‘Ornstein Uhlenbeck’ were published (Google Scholar search 15 March 2015; see Supporting Information). This may partly be because these models are now easy to apply via packages in R (e.g. ouch, GEIGER and OUwie; Butler & King, 2004; Harmon et al., 2008; Beaulieu & O'Meara, 2012). Additionally, although the OU model is pattern‐based, it has several attractive biological interpretations. For example, fit to an OU model is used as evidence for processes such as phylogenetic niche conservatism, convergent evolution and stabilizing selection (e.g. Wiens et al., 2010; Christin et al., 2013; Ingram & Mahler, 2013).

6It is important to note, however, that although the OU model is frequently described and interpreted as a model of ‘stabilizing selection’, this is inaccurate and misleading. As formulated by Hansen (1997), a trait has a primary optimum that is the mean of individual species optima for that trait. Under this formulation, α can be considered as the strength of the pull towards a central trait value (the primary optimum; Hansen, 2012). However, this is not an estimate of stabilizing selection in the population genetics sense, where it is a measure of selection within a population towards a fitness optimum on an adaptive landscape (Lande, 1976). This is a qualitatively different process to trait evolution among species which is more akin to a trait tracking movement of the adaptive optima itself.

7The OU model is most commonly used to model the evolution of a single continuous character (Table 1; Supporting Information). Usually several models of evolution (Brownian motion, OU, Early Burst, etc.; Harmon et al., 2010; Cooper & Purvis, 2010; Cardillo, 2015; Slater, 2015) are fit to the same continuous character and model selection is then used to determine which model best fits the data. OU models can be fit with one optimal trait value, or multiple different optima (Butler & King, 2004; Beaulieu et al., 2012). The latter represents evolution under multiple selective regimes, and may be more biologically realistic for many datasets. OU models with various numbers of optima are often included in the pool of evolutionary models being compared (e.g. Christin et al., 2013; Table 1). OU models are also commonly used to model phylogenetically structured residual error in evolutionary correlations (Revell, 2010; often colloquially referred to as ‘controlling for phylogeny’; Table 1; Supporting Information).

8As an extension to modelling single traits, phylogenetic generalized least squares (PGLS) models incorporate information about the relationships among species into the error term of a generalized least squares model. This error term generally consists of a variance–covariance matrix of the phylogeny, but various transformations are used (e.g. Pagel's λ; Pagel, 1997) to improve the fit of the model to the data. The α parameter from an OU model can also be used to transform the tree and the variance–covariance matrix. This is rarely interpreted as corresponding to any kind of process; instead it improves the fit of the PGLS models (e.g. Blankers, Adams & Wiens, 2012). However, it is not clear that it was originally intended that the OU model would merely be used in such a context. Finally, the OU model is also used to reconstruct ancestral states (Martins, 1999) and to detect clade‐wide convergent evolution (Ingram & Mahler, 2013; Uyeda & Harmon, 2014).

9Most papers use the OU model to model the evolution of a continuous character (Table 1; Supporting Information), so we focus on this use of the OU model in our simulations below. The principles here also apply to a range of other macroevolutionary models that can be fit to continuous data and compared using model testing procedures (e.g. κ, λ, δ, ACDC, Early Burst; Pagel, 1997, 1999; Blomberg, Garland & Ives, 2003; Harmon et al., 2008). Note that we focus on OU models with a single stationary optimum trait value because these are more commonly used (Table 1; Supporting Information) and easier to simulate. For discussions on the performance of multiple optima OU models we refer the reader to Beaulieu & O'Meara (2012).

10According to the Brownian model (Cavalli‐Sforza & Edwards, 1967; Felsenstein, 1973), a trait X evolves at random at a rate σ:where W(t) is drawn at random from a normal distribution with mean 0 and variance σ2. The model assumes that there is no overall drift in the direction of evolution (hence the expectation of W(t) is zero) and that the rate of evolution is constant. Because the direction of change in trait values at each step is random, Brownian motion is often described as a ‘random walk’. The model assumes the correlation structure among trait values is proportional to the extent of shared ancestry for pairs of species. This means that close relatives will be more similar in their trait values than more distant relatives. It also means that variance in the trait will increase (linearly) in proportion to time. The model has two parameters, the Brownian rate parameter, σ2, and the state of the root at time zero, X(0).

11The OU model (Hansen, 1997; Butler & King, 2004) is a random walk in which trait values revert back towards some ‘optimal’ value, μ (also called θ), with an attraction strength proportional to the parameter α. The model has the following form:

12Note that this model has two parameters in addition to those of the Brownian model: α and μ. α is the strength of evolutionary force that returns traits back towards the long‐term mean, μ, if they evolve away from it. α is sometimes referred to as the ‘rubber band’ parameter because of the way it pulls traits back towards μ. The parameter μ is a long‐term mean, and it is assumed that species traits evolve around this value. For more details see the Appendix.

13In many implementations of the OU model (e.g. GEIGER; Harmon et al., 2008), μ is the same as the state of the root at time zero, X(0). This is referred to as a ‘single stationary peak’ or SSP model (sensu Harmon et al., 2010). Some implementations also allow users to estimate X(0) (e.g. OUwie; Beaulieu & O'Meara, 2012). Where X(0) ≠ μ this is referred to as a single peak OU process. However, estimating X(0) can sometimes lead to nonsensical values of μ, i.e. far outside the range of values for the trait, so results should be interpreted with caution. X(0) can also be defined a priori, but this is only appropriate when fossil data, or other independent evidence, allow confident estimates of the root state.

14When α is close to zero, evolution is approximately Brownian (but note that in the special case of a single peak OU model with X(0) ≠ μ, when μ is zero, evolution approximates Brownian motion with a trend; Hansen, 1997; Benson et al., 2014), then as α gets larger the non‐Brownian behaviour of the model starts to become apparent. Eventually, when α is really large, all imprint of history is lost and the trait evolution is essentially a rapid burst at the present. Note that α scales with tree height (i.e. the maximum distance from the root of the tree to the tips); taller trees will have lower α values, all else being equal, because there is more time for traits to return to the optimum value, and thus the strength of the pull towards the optimum, α, can be smaller. Thus, α values need to be interpreted relative to tree height. Generally the simplest solution is to rescale tree heights to 1 (e.g. Ives & Garland, 2010; see simulations below). Ives & Garland (2010) suggest interpreting log (α) after rescaling trees to a height of 1, rather than raw α. They equate –log (α) = 4 as a very low, almost Brownian, value and –log (α) = −4 as a very high value. Others (e.g. Hansen & Bartoszek, 2012; Slater, 2015) prefer to use the phylogenetic half‐life: t12 (see ‘Recommendations for interpreting α’ below).

15To explore some issues with the OU model in more detail, we ran a number of simulations designed to mirror the use of OU models in the literature.

16We simulated phylogenies with 25, 50, 100, 150, 200, 500, or 1000 tips under pure birth, constant‐rate Birth–Death (extinction fractions of 0.25, 0.5 and 0.75) or temporally varying speciation rate (speciation rate modelled as time from the root raised to the power 0.2, 0.5, 2 and 5) models. We simulated 1000 phylogenies for each combination of tips and models resulting in 56 000 simulated phylogenies in total. Trees were simulated using the R package TESS (Hohna, 2013). We then simulated the evolution of a single trait under a Brownian motion model on each phylogeny using the R package MOTMOT (Thomas & Freckleton, 2011). All our simulated trees and data are available on GitHub: https://github.com/nhcooper123/OhYou.

17To determine whether α is biased under conditions where it should be zero, and whether Likelihood ratio tests are appropriate for use with the OU model, we estimated α for each of our simulated phylogenies and data, and compared the fit of a Brownian model with that of an OU model using a Likelihood ratio test with 1 degree of freedom with the transformPhylo.ML function in MOTMOT (Thomas & Freckleton, 2011; https://github.com/ghthomas/motmot). This mirrors the common situation where researchers fit Brownian and OU models and then use Likelihood ratio tests to select the ‘best’ model. We then estimated the rejection rate of the null (Brownian) model for each set of simulations. We refer to this measure as the Type I error rate for the OU model.

18We find that Type I error rates are unacceptably high when tree size is small (Tables 2, 3), i.e. the OU model is often favoured, even though the Brownian model was used to generate the data. For some tree shapes, particularly where speciation rates accelerate towards the present, Type I error remains > 0.05 even for trees with 1000 tips (Table 3). This shows that, in general, analyses based on small datasets are prone to biases that decrease only slowly as the size of the dataset increases. Unfortunately, OU models are often fitted to phylogenies with fewer than 100 taxa (mean = 166.97 ± 43.86, median = 58; Fig. 3; see Supporting Information). We provide recommendations related to these findings below.

19For small trees the confidence limits on parameter estimates are broad (Fig. 2). However, the median is typically low irrespective of tree size (see Tables 2, 3) and even when the OU model is favoured, scrutiny of the model parameters suggest that the favoured OU model is biologically indistinguishable from Brownian motion. Indeed, examination of the effect of α on the expected covariances among taxa (see Fig. 4) confirms that covariances are barely distinguishable when α < 1 for unit height trees. This suggests that, notwithstanding the elevated Type I errors of OU models, their interpretation requires examination of parameters. We return to the issue of parameter interpretation below, specifically with reference to the increasingly commonly used phylogenetic half‐life.

20Measurement error can strongly affect the results of comparative analyses (Silvestro et al., 2015), and appears to influence whether OU is the favoured model across a range of datasets (see Pennell et al., 2015). Therefore, we also investigated whether adding error to our simulated data influenced estimates of α.

21We used the same procedure as above to simulate trait data under a Brownian motion model with known error. Specifically, we simulated trees under a Yule model with 25, 50, 100, 150, 200, 500 or 1000 tips and added branch length of (1) 1%, (2) 5% or (3) 10% of the tree height to the tips of the simulated trees. We then simulated data under a Brownian model on each tree. All our simulated trees and data are available on GitHub: https://github.com/nhcooper123/OhYou.

22We compared the fit of a Brownian model with that of an OU model using Likelihood ratio tests as described above, using the original trees without the addition of extra branch length to the tip edges. Our expectation is that the OU model should fit the data better than the Brownian model because the α parameter should account for much of the error. In this case, rejection of the Brownian model does not represent Type I error because it would be quite correct to reject the Brownian model. However, the reason for the better fit would be entirely unrelated to any macroevolutionary process.

23Table 4 shows the proportion of data sets in which the OU model is favoured over the Brownian model for data simulated under Brownian motion with error. The expectation is that the OU model should fit better because the branch length transformation partially captures the non‐Brownian component (the error). There are two points worth noting. First, the frequency with which the OU model is favoured increases with tree size. With as little as 5% error, the OU model becomes extremely difficult to reject, even for trees with just 100 species. This is important for the interpretation of the OU model; we cannot conclude anything about the evolutionary process from a single optimum OU model unless error is adequately accounted for. Second, for moderate amounts of error (5–10%), estimates of α are consistently > 1. Large values of α are similarly difficult to interpret because they are indicative that the signal of the past has been overwritten.

24Although it is possible to create and implement new models for comparative data that encompass a range of processes, we have to be aware that such models are statistically complex and may behave in unexpected ways. Transformations of the variance–covariance matrix in the Brownian model (e.g. λ; Pagel, 1997) are an attractive and computationally simple way to modify the basic model to include evolutionary processes. But, as first pointed out by Grafen (1989), the statistical consequences of these modifications can include biases and problems with interpretation. The results of our simulations illustrate that such problems can occur under conditions that closely match the size and type of datasets that are commonly used (Fig. 3, Supporting Information, Table S1).

25In the case of the OU model, there are several limitations worth highlighting: Type I error rates are high when sample size is low. The results of the simulations indicate that, in general, analyses based on small datasets are prone to biases that decrease only slowly as the size of the dataset increases. Likelihood ratio tests are untrustworthy. All Likelihood ratio tests assume that the Likelihood ratio statistic is asymptotically (i.e. as sample sizes become large) χ2‐distributed. In the OU model α is bounded (i.e. it cannot be any smaller than zero) and has a non‐linear effect on the expected variances, so the assumptions of the Likelihood ratio test are not likely to be upheld for small samples. Our simulations indicate that Likelihood ratio tests should not be relied upon for analyses with small sample sizes, and that for robust inference and testing, alternatives, such as simulation and Monte Carlo Markov chains (MCMC), should be considered [see Recommendations and potential solutions (1) and (2) below]. Note that these issues will also apply when using the Akaike information criterion (AIC) to compare models, as this is essentially the same as performing a Likelihood ratio test. AIC also presents other difficulties when used for complex models such as those used in phylogenetic comparative analyses, so we recommend avoiding them for these reasons. Measurement error increases Type I error rates. Our results show that a simple Brownian process can be mistaken for an OU process when a small amount of error is added to the data. The effects of measurement error become more severe with increasing tree size. These limitations mean that when evidence for the OU model is found, the results should be interpreted with caution, particularly where there is likely to be intraspecific variation or measurement error in the data.

26Type I error rates are high when sample size is low. The results of the simulations indicate that, in general, analyses based on small datasets are prone to biases that decrease only slowly as the size of the dataset increases.

27Likelihood ratio tests are untrustworthy. All Likelihood ratio tests assume that the Likelihood ratio statistic is asymptotically (i.e. as sample sizes become large) χ2‐distributed. In the OU model α is bounded (i.e. it cannot be any smaller than zero) and has a non‐linear effect on the expected variances, so the assumptions of the Likelihood ratio test are not likely to be upheld for small samples. Our simulations indicate that Likelihood ratio tests should not be relied upon for analyses with small sample sizes, and that for robust inference and testing, alternatives, such as simulation and Monte Carlo Markov chains (MCMC), should be considered [see Recommendations and potential solutions (1) and (2) below]. Note that these issues will also apply when using the Akaike information criterion (AIC) to compare models, as this is essentially the same as performing a Likelihood ratio test. AIC also presents other difficulties when used for complex models such as those used in phylogenetic comparative analyses, so we recommend avoiding them for these reasons.

28Measurement error increases Type I error rates. Our results show that a simple Brownian process can be mistaken for an OU process when a small amount of error is added to the data. The effects of measurement error become more severe with increasing tree size. These limitations mean that when evidence for the OU model is found, the results should be interpreted with caution, particularly where there is likely to be intraspecific variation or measurement error in the data.

29We have highlighted several issues with the OU model above. These can be broadly divided into two classes: (1) fitting the model and (2) interpreting model parameters. These issues can be at least partially addressed with additional or alternative analytical approaches. Below we make recommendations for best practice in fitting the model and approaches to interpret model parameters. We highlight areas of future research that may be important in alleviating outstanding issues.

30Simulate under the null model. OU models should not be applied to small trees. Our simulations indicate that trees with > 200 tips are necessary to obtain acceptable Type I error rates. However, we are cautious about recommending a minimum tree size because the performance of the OU model may vary among datasets for reasons other than tree size. We note also that large trees are particularly susceptible to issues arising from unaccounted measurement error in the data. Instead, we suggest using simulations to assess the model fit. Simulating data under the Brownian model will generate null distributions (e.g. Boettiger et al., 2012) and allow straightforward assessment of model fit, for example by generating appropriate critical values for Likelihood ratio tests. Boettiger et al. (2012) discuss this at length and other papers use similar approaches (e.g. Martins & Garland, 1991; Freckleton, Harvey & Pagel, 2002). Consider Bayesian approaches. An alternative but less explored approach is to use Bayesian methods. Bayesian model fitting has not been widely available for OU model fitting until recently but is now possible in R packages including diversitree (FitzJohn, 2012) and GEIGER (the Single Stationary Peak model in fitContinuousMCMC; Harmon et al., 2008) and stand‐alone software including BayesTraits (Pagel & Meade, 2013).

31Simulate under the null model. OU models should not be applied to small trees. Our simulations indicate that trees with > 200 tips are necessary to obtain acceptable Type I error rates. However, we are cautious about recommending a minimum tree size because the performance of the OU model may vary among datasets for reasons other than tree size. We note also that large trees are particularly susceptible to issues arising from unaccounted measurement error in the data. Instead, we suggest using simulations to assess the model fit. Simulating data under the Brownian model will generate null distributions (e.g. Boettiger et al., 2012) and allow straightforward assessment of model fit, for example by generating appropriate critical values for Likelihood ratio tests. Boettiger et al. (2012) discuss this at length and other papers use similar approaches (e.g. Martins & Garland, 1991; Freckleton, Harvey & Pagel, 2002).

32Consider Bayesian approaches. An alternative but less explored approach is to use Bayesian methods. Bayesian model fitting has not been widely available for OU model fitting until recently but is now possible in R packages including diversitree (FitzJohn, 2012) and GEIGER (the Single Stationary Peak model in fitContinuousMCMC; Harmon et al., 2008) and stand‐alone software including BayesTraits (Pagel & Meade, 2013).

33We explored Bayesian methods as a possible remedy to the limitations of fitting OU models in a Likelihood framework. We repeated our simulations in a Bayesian framework implemented in BayesTraits (Pagel & Meade, 2013). We determined the fit to an OU model using Bayes factors estimated from a stepping stone sampling procedure (Xie et al., 2010). The marginal likelihoods of the models were calculated using a stepping stone sampler in which 50 stones were drawn from a beta distribution with α = 0.4 and β = 1. Each stone was sampled for 20 000 iterations (with the first 5000 iterations discarded). We use stepping stone sampling as it has been shown to estimate the marginal likelihood better than the harmonic mean (Baele et al., 2012). We treated Bayes factors > 2 as evidence favouring the OU model (Kass & Raftery, 1995). We ran the MCMCs for 1 × 106 iterations, disregarding the first 1 × 104 as burn‐in. Following burn‐in the chains were sampled every 1000 iterations to ensure independence of each consecutive sample. Multiple independent chains were run for each analysis to ensure convergence was reached. An important challenge for Bayesian approaches is selection of appropriate priors. We used three alternative sets of priors on α: (1) an exponential distribution with mean = 1, (2) an exponential distribution with mean = 10 and (3) a uniform distribution bounded at 0 and 20. For all analyses we used a uniform −100 to 100 prior for μ and uniform 0–100 prior for σ2.

34Table 5 shows the results from the exponential prior with mean = 10. This is a broad, liberal prior, but our results are similar regardless of priors. The Bayesian approach results in highly conservative rejection rates regardless of tree shape in the absence of measurement error. While this is encouraging from the perspective of falsely rejecting the Brownian null model, it is also indicative of potentially low statistical power, although testing would be required to confirm this. The Bayesian approach also appears to more readily handle low levels of measurement error (Table 6). With 1% measurement error the Bayesian approach retains acceptable rejection rates (< 0.05) for trees of up to 150 tips. However, more error or larger trees result in the frequent rejection of the Brownian model. As noted above, this is not an issue of Type I error and it is entirely correct that the Brownian model is rejected. In such cases emphasis should shift to how the α parameter is interpreted. Bayesian analysis is a promising approach for OU model fitting. However, further testing of the Bayesian approach to simulated data sets under a wide range of values of α is necessary to fully characterize performance.

35Consider plausible alternative hypotheses. Because the OU model was proposed with evolutionary processes in mind, alternative and arguably more parsimonious explanations for favouring OU models and interpreting non‐zero α are often overlooked. The effects of measurement error in particular suggest that α must always be interpreted with caution. Many issues of misinterpretation can be solved by carefully inspecting the α parameter when an OU model is favoured. Often, when Likelihood ratio tests suggest the OU model should be favoured over the Brownian model, estimates of α are actually very small and biologically indistinguishable from Brownian (e.g. examples in Harmon et al., 2010). In these circumstances, it is likely that measurement error, intraspecific variation or phylogenetic uncertainty are generating noise that is more effectively modelled by the extra parameters in the OU model than by the Brownian model alone. Thus, this does not reflect any kind of OU process underlying the data. The similarity between Brownian and OU models with small α is demonstrated in Fig. 4. At the other extreme, as values of α become larger, the effects of changing α on model predictions are increasingly small and large values of α are indistinguishable from white‐noise.

36Consider plausible alternative hypotheses. Because the OU model was proposed with evolutionary processes in mind, alternative and arguably more parsimonious explanations for favouring OU models and interpreting non‐zero α are often overlooked. The effects of measurement error in particular suggest that α must always be interpreted with caution. Many issues of misinterpretation can be solved by carefully inspecting the α parameter when an OU model is favoured. Often, when Likelihood ratio tests suggest the OU model should be favoured over the Brownian model, estimates of α are actually very small and biologically indistinguishable from Brownian (e.g. examples in Harmon et al., 2010). In these circumstances, it is likely that measurement error, intraspecific variation or phylogenetic uncertainty are generating noise that is more effectively modelled by the extra parameters in the OU model than by the Brownian model alone. Thus, this does not reflect any kind of OU process underlying the data. The similarity between Brownian and OU models with small α is demonstrated in Fig. 4. At the other extreme, as values of α become larger, the effects of changing α on model predictions are increasingly small and large values of α are indistinguishable from white‐noise.

37One good strategy for data exploration would be to simulate data under Brownian and the favoured OU model to generate distributions of parameters under known values. These can then be compared with results for your dataset (see Slater, 2014; Slater & Pennell, 2014; for a related approach). This is important because we have shown that the shape of a phylogeny has consequences for parameter biases and hypothesis tests. Any given tree will therefore generate unique parameter estimates. Generating data under the favoured OU model will allow an assessment of whether it is possible to retrieve known values, or whether there is evidence of bias.

38Calculate the phylogenetic half‐life. The α parameter ranges from zero to infinity (although in practice an upper bound is often set; for example, in GEIGER this is α = 150; Harmon et al., 2008), thus recognizing ‘small’ or ‘large’ values may not be intuitive. α can sometimes be interpreted more easily by using it to estimate the ‘phylogenetic half‐life’ (t12) of a trait, i.e. the time it takes for a species entering a new niche to evolve halfway toward its new expected optimum (Hansen, 1997), as follows:

39Calculate the phylogenetic half‐life. The α parameter ranges from zero to infinity (although in practice an upper bound is often set; for example, in GEIGER this is α = 150; Harmon et al., 2008), thus recognizing ‘small’ or ‘large’ values may not be intuitive. α can sometimes be interpreted more easily by using it to estimate the ‘phylogenetic half‐life’ (t12) of a trait, i.e. the time it takes for a species entering a new niche to evolve halfway toward its new expected optimum (Hansen, 1997), as follows:

40If t12 is short relative to the branch lengths of the phylogeny, evolution towards the optimum trait value is fast, residual phylogenetic correlations are weak and there is little influence of the past on trait values (Hansen, 1997). t12 equal to the height of the phylogeny is a moderate value (Hansen & Bartoszek, 2012). We would not advise interpreting t12 as literally being ‘the time it takes for a species entering a new niche to evolve halfway toward its new expected optimum’ (Hansen, 1997). However, if t12 is extremely large relative to tree height, it suggests that if an OU process is acting, it is extremely weak (e.g. for a clade that is 50 Myr old, a t12 of 100 Myr suggests a species will not approach the optimum within the temporal range of the clade), and thus should not be interpreted as evidence of any kind of process. As a further note of caution, it is important to recognize that biases in the estimation of α would lead to similar biases in t12. If possible, include ancillary data. If data on fossils are available then these could be incorporated into the analysis (Slater, Harmon & Alfaro, 2012). Indeed, if fossil taxa can be reliably placed in the phylogeny then they may improve model accuracy and a strong case can be made that fossils should be included. On the other hand, if there is substantial uncertainty in fossil placement then they should be treated cautiously. Placement of fossil taxa is implemented in GEIGER (Harmon et al., 2008) and BayesTraits (Pagel & Meade, 2013). A caution here is that the OU model for non‐ultrametric trees has to be carefully parameterized because for non‐ultrametric trees the co‐variances depend on both the shared distances between species and the distance of a node to the nearest tip (see eqn A6). This creates potential problems in parameterization and in interpretation because the variance–covariance matrix is no longer tree‐like; for example, related species can effectively become more similar to one another than to themselves – an inherently non tree‐like pattern (Slater, 2014). Some current implementations of the OU model are based on transforming the tree directly, rather than transforming the variance–covariance matrix (e.g. MOTMOT; Thomas & Freckleton, 2011). These implementations should not be used with fossil data.

41If possible, include ancillary data. If data on fossils are available then these could be incorporated into the analysis (Slater, Harmon & Alfaro, 2012). Indeed, if fossil taxa can be reliably placed in the phylogeny then they may improve model accuracy and a strong case can be made that fossils should be included. On the other hand, if there is substantial uncertainty in fossil placement then they should be treated cautiously. Placement of fossil taxa is implemented in GEIGER (Harmon et al., 2008) and BayesTraits (Pagel & Meade, 2013). A caution here is that the OU model for non‐ultrametric trees has to be carefully parameterized because for non‐ultrametric trees the co‐variances depend on both the shared distances between species and the distance of a node to the nearest tip (see eqn A6). This creates potential problems in parameterization and in interpretation because the variance–covariance matrix is no longer tree‐like; for example, related species can effectively become more similar to one another than to themselves – an inherently non tree‐like pattern (Slater, 2014). Some current implementations of the OU model are based on transforming the tree directly, rather than transforming the variance–covariance matrix (e.g. MOTMOT; Thomas & Freckleton, 2011). These implementations should not be used with fossil data.

42A recurring theme from our simulations is that interpretation of OU models is not straightforward. More focus is needed on the interpretation of α rather than simply model fit. Even when the OU model is favoured, α may be so small as to be indistinguishable from Brownian motion in any biological sense. It would clearly be very useful to have estimates of measurement error for all species traits, although the inclusion of species‐specific variances has to be done carefully (e.g. Grafen, 1989). Several approaches for accounting for error have been proposed and warrant wider implementation for OU models and other models of trait evolution (e.g. Lynch, 1991; Martins & Hansen, 1997; Ives, Midford & Garland, 2007; Hansen & Bartoszek, 2012; Rohlfs, Harrigan & Nielsen, 2014), but it is outside our scope to explore these here.

43Indeed, the problems that we report are not limited to OU models. Any model of trait evolution that attempts to account for non‐Brownian components of trait variation is susceptible to being misled by measurement error, and in some scenarios measurement error can also incorrectly favour Brownian motion over the true model, e.g. if the true model is Early Burst. The fundamental problem is that rejection of the Brownian model in favour of another model does not necessarily say anything about process. This problem can be alleviated to some extent if model comparisons are set in a firm hypothesis‐testing framework in which alternative hypotheses make clear predictions of emerging patterns that can be unambiguously associated with particular models (e.g. Cooper, Freckleton & Jetz, 2011), although this does not appear to be possible for comparisons of the single stationary peak OU model with noisy Brownian processes. We should therefore not use any statistical model without thinking carefully about the limits in terms of both data and interpretation.

44At the symposium that generated this special issue, one of us (N.C.) gave a talk on Open Science and reproducibility. We have therefore tried to make this paper as open and reproducible as possible. All simulated phylogenies, data and R code are available on GitHub: https://github.com/nhcooper123/OhYou. However, although our R code is provided, it is disorganized and thus difficult to use. We also do not provide an automated way of reproducing the data collection for our literature review. Nor do we use tools such as Travis CI (travis‐ci.org), Docker (www.docker.com) or packrat (Ushey et al., 2015) to increase the reproducibility of our analyses. Additionally, we use BayesTraits to run our Bayesian analyses. BayesTraits is free to download as a binary executable for various platforms but it is not Open Source. This means that while BayesTraits analyses are reproducible, the software has limitations with respect to long‐term development of the code. We included this epilogue to highlight that fully reproducible and Open Science is challenging, but we are trying to improve. With a little effort most people should be able to produce something vaguely reproducible, and to provide their data and code, moving us all slightly closer to truly reproducible and Open Science.