Modern Sciences is a premier science journal that bridges the gap between science and its application to society.
LATEST STORIES

Using state space models to understand trait evolution in fossil lineages

Gene Hunt1*, Wilmer Martinez Rivera2, Melanie Hopkins3, John Fricks4 and Beckett Sterner5

    1. 1The Smithsonian National Museum of Natural History, 10th Street and Constitution Avenue NW, Washington, DC 20560
    2. 2Banco de la República de Colombia, Cra. 7 #14-78, Santa Fé, Bogotá, Columbia
    3. 3American Museum of Natural History, 200 Central Park West, New York, NY 10024
    4. 4School of Mathematical and Statistical Sciences, Arizona State University, 901 South Palm Walk, Tempe, AZ 85287-1804, USA
    5. 5School of Life Sciences, Arizona State University, 427 East Tyler Mall, Tempe, AZ 85287-1804, USA
    6. *Corresponding author; email: hunte{at}si.edu

    bioRxiv preprint DOI: https://doi.org/10.1101/2024.07.17.603977

    Posted: July 19, 2024, Version 1

    Copyright: This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available for use under a CC0 license

    Abstract

    Linear state space models provide a useful framework for investigating phenotypic evolution in fossil lineages in a wide variety of models including Brownian motion, Ornstein-Uhlenbeck processes, and models that incorporate potentially explanatory environmental covariates. A state space framework also provides access to residuals for the predicted and observed values at each time point as well as improved numerical stability. We illustrate the value of the state space approach by re-analyzing a classic dataset of diatom evolution in Yellowstone Lake. We find that number of spines is best explained by adaptation to changing solar insolation as an exogenous environmental covariate.

    1 Introduction

    The project of documenting and explaining historical patterns of evolutionary change has had enduring significance since George Gaylord Simpson introduced the ideas of evolutionary tempo and mode as a tool for bridging observation and theory across micro- and macro-evolutionary time scales (Simpson 1944). Fossil trait series provide a sequence of phenotypic measurements drawn from multiple organisms in the same lineage over a period of time, typically on the scale of tens to hundreds of thousands of years. In contrast to comparative methods, which rely on measurements of different species at particular moments, fossil trait series document historical patterns of evolutionary change in a single species over much longer durations than we can access in the present. This makes fossil trait series a uniquely valuable data source for addressing fundamental questions in evolutionary biology occurring at those scales, including why evolutionary rates of change appear to show a robust scaling with time (Harmon et al. 2021) and whether evolutionary divergence is driven by the gradual accumulation of small changes or short, rapid bursts of change that punctuate long periods of stasis (Hunt, Hopkins, and Lidgard 2015). In addition, paleobiologists have looked for signatures of other evolutionary and ecological processes, seeking to gain insights into processes of adaptation (Voje 2020; Kearns et al. 2021), extinction (Brombacher et al. 2017; Brombacher et al. 2018), stasis (Antell et al. 2021) and parallel evolution (Stuart, Travis, and Bell 2020).

    Since evolution is stochastic and involves the interaction of multiple processes such as selection and drift, statistical modeling is essential to reliably estimating evolutionary rates and classifying fossil trait series into biologically significant patterns or modes of change (Hunt 2012; Gingerich 2019). In addition, the sample sizes of fossil lineage data are often small relative to the potential complexity of the system, so it is especially important to make effective use of the available information. The appropriate interpretation of evolutionary rates in a system, for example, is sensitive to the underlying mode of change: for a lineage evolving according to classical Brownian motion, we can understand its rate of change as the variance of random fluctuations it undergoes between points of time, but for a lineage exhibiting a linear directional trend in addition to Brownian motion we also need to incorporate the magnitude of that linear tendency (Hunt 2012). Numerical estimation of evolutionary rates is also sensitive to model choice, and using an inadequate model for the data can generate scaling artifacts in the magnitude of rates in relation to the absolute time duration of the trait series (Hunt 2012; Gingerich 2019; Harmon et al. 2021).

    In this paper, we present some practical tools and background theory for the use of linear state space models, also called dynamic linear models, in the analysis of phenotypic time series from fossil lineages. This approach provides some key advances over previous methodological approaches. First, state space models, which apply to time series with observations at discrete time points, are grounded in continuous-time models, which are important for knowing how to handle time-varying parameters when there are uneven time steps between observations. Second, the state space framework allows for the ready evaluation of exogenous environmental covariates. Specifically, the framework provides access to residuals for the predicted and observed values at each time point, a powerful model diagnostic tool that is in standard use across other scientific disciplines. We illustrate the value of these residuals to propose a new interpretation of a classic dataset of diatom evolution in Yellowstone Lake (Edward C. Theriot et al. 2006). We also present some simulation studies emphasizing the quality of parameter estimation and model selection within this framework in comparison to previous frameworks for estimation.

    We focus on univariate evolutionary models from a single species because this foundation is essential for numerically efficient and stable parameter estimation as paleontologists start to increasingly take advantage of high-throughput specimen processing, including machine-learning methods for image segmentation and trait extraction (Porto and Voje 2020; Goswami and Clavel 2024; He et al. 2024). However, the modeling framework we develop here is readily generalized to analyzing multivariate measurements from variable numbers of individual specimens at each time point. This approach also demonstrates the broad relevance of stochastic integral models to analyzing fossil lineage datasets, building on prior work by Reitan (Reitan, Schweder, and Henderiks 2012) and is complementary to stochastic integral models for species richness (Hannisdal and Liow 2018; Reitan and Liow 2019) and comparative phylogenetic datasets (Blomberg, Rathnayake, and Moreau 2020; Turley 2020).

    2 Methods

    In this section, we briefly present the definition of a linear state space model, sometimes known as a dynamic linear model. We introduce the general structure of the model along with some adaptations and quantities that this formulation facilitates. Then, we show how observing an Ornstein-Uhlenbeck (OU) process at discrete time points can be formulated as a linear state-space model. In particular, we make clear the relationship between the continuously-defined model and the discretely-sampled observations. We also discuss the relationship of other familiar models, such as an undirected random walk (RW), random walk with a directional trend (DT), and stasis (ST) to the linear state space models in the appendix.

    2.1 Definition of a Linear State Space Model

    The linear state space model consists of two recursively defined equations: the system equation and the observation equation. We will focus on the uni-dimensional case for each equation, although multivariate extension is straightforward. The system equation is a linear stochastic equation describing the underlying dynamics of a phenomena, such as the change of the level of a trait measurement for a population at a given time. The system equation can be written asEmbedded Imagefor t = 1, …, T where W1, …,WT are a sequence of independent and identically distributed normal random variables with zero-mean and variance Embedded Image. This is an autoregressive model with coefficient Φ. The model is called autoregressive because the sequence could be viewed as having the current time observation depend on the previous observation as a covariate. The observation equation is then defined asEmbedded Imagefor t = 1, …, T where V1, …,VT is a sequence of independent and identically distributed normal random variables with zero-mean and variance Embedded Image which is also independent of the sequence W1, …,WT. This sequence of Y1, …YT will model the observed data and is a linear transformation by A of the system process, Xt, plus an additional observation error, Vt. Note that the notation here follows the book by Shumway and Stoffer, where detailed calculations associated with this model can be found (Shumway, Stoffer, and Stoffer 2000).

    In the sub-section that follows and in the appendix, we will show how a number of familiar trait evolution models fit into this framework, especially when a few features are added. This linear state space model is the framework for calculating model likelihoods using the renowned Kalman filter. The Kalman filter recursively calculates the conditional distribution of Xt given observations Y1, …,Yt; this means that if you have the Kalman filter calculated up to time t − 1 and receive a new observation Yt, then the conditional distribution of Xt given the data up to time t can be calculated without passing back through all of the previous data. A by-product of calculating the Kalman filter is that the likelihood function can be calculated with one pass through the data, so that the computational efficiency of this calculation is of order T. The likelihood function for a linear state space model can also be calculated using a multivariate normal distribution of size T, which has been the standard approach in previous paleontological research (Hunt, Bell, and Travis 2008; Voje 2020); however, that approach requires calculating the inverse of a T ×T matrix. Hence there can be a significant speed up and improved stability using the Kalman filter algorithm to calculate the relevant likelihood function, especially for time-series with many samples.

    The Kalman filter approach also makes prediction relatively straightforward. The filter calculates the conditional mean and variance of Xt given Y1, …,Yt, giving us access to a “best guess” for Xt with the data up to that time. This also allows us to predict, i.e. to find the conditional mean and variance, of the next observation Yt+1 given the observations up to time t (i.e. Y1, …Yt), which is known as the one-step ahead predictor of the data and also allows us to construct residuals for our data after having fit it to a specific state space model. Again, using the notation of Shumway and Stoffer, the standardized residuals would be defined asEmbedded Imagefor t = 1, …, T −1, where Σt+1 is the conditional variance of Yt+1 given Y1, …,Yt, which is calculated as part of the Kalman filter algorithm. Residuals in a time series context are an important tool for discerning the quality of fit of a model, in a manner similar to their use in regression analysis. It is a way to “approximate” the sequence Vt, and we therefore expect them to be approximately independent and identically distributed.

    In addition to facilitating likelihood function calculations and predictions of both the system equation model and future observations, the linear state space model facilitates regression with exogenous variables, which is especially important when considering the effects of environmental variables on trait mean dynamics. These can enter either through the system equation or the observation equations. In addition, the parameters defining the linear state space model can be time-varying. These modifications giveEmbedded ImageandEmbedded Imagewhere ut is an r × 1 column vector of exogenous variables, typically what we think of as covariates, for each t, and ϒ and Γ are the coefficients that convert changes in the input variables into changes in Xt and Yt. In this context, we can still use the Kalman filter to calculate the likelihood and to create predictions and thus residuals to evaluate model performance. Visually examining the residuals is a useful complement for formal model misspecification tests for common fossil lineage models introduced in (Voje 2018).

    2.2 Unbiased Random Walk as a Linear State Space model

    Many commonly used trait evolution models can be represented in this state space framework. Perhaps the simplest case is that of an unbiased random walk, which we develop here as an example. Under this model, the trait value at one time step (Xt) is equal to the value at the previous time step (Xt−1), plus an evolutionary perturbation (Wt). The autoregressive coefficient is unity (Φ = 1), and the variance of the Wt is usually called the step variance in paleontology. The state equation represents the true trait mean of the population at each time step. Since sample sizes are finite, however, we can never know the true means of our samples. The calculated trait means that we observe, Yt, reflects the true population mean (Xt), plus sampling noise (Vt). In most paleontological cases, our measurements of traits are noisy but unbiased, and so A = 1. For trait means, the observational variance of Vt for each sample will be the within-sample trait variance divided by the number of individuals measured in that sample.

    2.3 Ornstein-Uhlenbeck as a Linear State Space model

    The Ornstein-Uhlenbeck process is used in evolutionary studies to model adaptation (Hunt, Bell, and Travis 2008; Hansen 1997). It is often defined by its stochastic differential equation formulation.Embedded Image

    This type of differential is a short hand to define an integral equation. Such an equation is defined as a solution to sequences of difference equations:Embedded Imageas Δ goes to zero and nΔ converges to some terminal time T and with X (0) defining an initial value for the equations. (This sequence is on an evenly spaced grid, which is not necessary in general as long as the grid spacing distance shrinks to zero.) The Ornstein-Uhlenbeck process has a solution (Oksendal 2013) in the form of a stochastic integral with deterministic integrand,Embedded Image

    The definition of these SDEs as a limit of solutions to certain difference equations point to one possible discretization. Namely, we could re-arrange the difference equation to arrive atEmbedded Imagebut this is only an approximate solution that accurately represents the model definition when Δ is small. However, if we take the integral solution of Equation 3 and manipulate it appropriately, we arrive atEmbedded Image

    We can use Ito’s isometry to calculate the variance of the last term, which is Gaussian and has a zero mean. The variance is thenEmbedded Image

    We solve the integral to arrive at the variance of the last termEmbedded Image

    Note that when Δ is small, then a Taylor approximation will verify that this expression is approximately equal to σ2Δ corresponding to equation 4. So, we can express this discretization asEmbedded Imagewhere εn is a sequence of independent standard normal random variables.

    The Ornstein-Uhlenbeck model will eventually settle around zero regardless of the initial condition. We can modify this part of the model for the OU process to be centered around another constant, θ. The exact discrete version would beEmbedded Image

    This formulation allows for one way to introduce covariates with a linear combination of covariates in place of θ.

    2.4 Model estimation and selection simulations

    We used simulations to validate model fitting in the state space framework using the Kalman filter compared to prior work. We also explored parameter estimation and model selection performance in the state space framework. In particular, there are two inherent observational time scales in trait series: the time step between observed samples (Δ above) and the total duration between the first and last observations. Parameter convergence may therefore depend on how increasing the sample size alters these time scales. For example, subdividing a fixed interval of time with more observations does not lead the linear trend parameter, µ, to converge asymptotically in the DT model. In contrast, smaller time steps are valuable for estimating the Appendix 7 provides more details on parameter estimation and model selection performance.

    3 Application: Stephanodiscus yellowstonensis trait evolution

    To illustrate the linear state space framework, we re-analyzed the Stephanodiscus yellowstonensis fossil trait series published in (Voje 2020), originally created by (Edward C. Theriot et al. 2006), using the discretized models described in Section 2.1 and implemented in the state space frame-work. S. yellowstonensis is a species of diatom endemic to Yellowstone Lake, Wyoming, United States and likely descended from the pre-existing species S. niagarae, which is still extant through the region (Edward C. Theriot et al. 2006). The fossil trait series is derived from 63 samples from a sediment core collected from the lake’s central basin, and it covers approximately 14,000 years ago until the present. For each sample, Theriot et al. measured 50 individuals, occasionally fewer if this number of specimens was not available. They measured three traits on each individual: valve diameter, the number of costae per valve, and the number of spines per valve. All three traits show a relatively rapid increase in values from about 12,000 to 10,000 years ago, with a slower and fluctuating decrease thereafter. As noted by (Voje 2020), all three traits are considered to be ecologically important for diatoms; valve spines, in particular, may enhance nutrient uptake and photosynthetic rate by affecting how diatoms sink through the water column. To illustrate the model-fitting methods, we focus on just one variable, spine count (Fig 1) Before analysis, spine counts were log-transformed because we consider a proportional scale to be more appropriate for the evolution of this trait. This transformation also has the effect of removing a strong correlation between the mean and variance among samples (r = 0.70, P < 1e − 9 for untransformed spine counts, r = 0.04, P = 0.77 after log-transformation).

    Figure 1:

    Evolutionary time-series of log spine counts for Stephanodiscus yellowstonensis. Open circles and vertical bars show mean values +/-one standard error on those means. The red line is the model-predicted trajectory of the fitness optimum for the best supported model (see text).

    Theriot et al (Edward C. Theriot et al. 2006) posit that morphological changes in the S. yellow-stonensis lineage track environmental changes, and these authors synthesize the available records of regional change through the study interval. Here, we quantitatively analyze two of these records. The first is the proportion of the dominant pollen type, attributable to Pinus contortus, as reflecting floral change in the area (Fig 2, digitized from Fig. 10 in ref (Edward C. Theriot et al. 2006)). The second environmental record we analyzed is solar insolation, which peaked around 11,000 years ago and decreased to the present day. Insolation values were taken from model output Lasker et al. (2004), using the web interface vo.imce.fr. Insolation values were computed in W/m2 at the latitude and longitude of Lake Yellowstone. As S. yellowstonensis blooms in summer (Edward C. Theriot et al. 2006), we used insolation values for the month of June (Fig 2).

    Figure 2:

    Measured environmental covariates, including June solar insolation (red solid line) and proportion of pollen attributable to Pinus contortus (blue dashed; note reversed axis).

    The sampling times of the environmental records did not precisely match those of the trait time-series. We used linear interpolation to produce time-series of the environmental records, sampled at the same times as the traits. Both environmental records were mean-centered prior to analysis to facilitate model fitting.

    Our overall model-fitting strategy started with the five models considered by Voje (Voje 2020): stasis (ST), random walk (RW), directional trend (DT), Ornstein-Uhlenbeck (OU), and decelerated evolution (DE). To test Theriot et al.’s suggestion that environmental changes influenced morphological evolution, we added OU models in which the trait optimum linearly tracks the two environmental covariates described above, June solar insolation (OUinsol) and the proportional abundance of Pinus contortus pollen (OUpollen). Because model fits and residuals indicated a decrease in stochastic evolutionary change through the core (see Results), we considered additional models that allowed for a one-time decrease in the step variance. We specified this time to be 10,000 years ago, following Theriot et al.’s observation [(Edward C. Theriot et al. 2006), p. 45] that environmental conditions were much more stable after this date.

    All the above models were fit using functions in the R package paleoTS, the recent update of which (v. 0.6-1) allows for fitting models via the Kalman filter and a state space model approach. Confidence intervals on parameter estimates were generated using the dentist package to compute approximate profile confidence intervals.

    Of the five models fit by Voje, we found that DE was best supported by AICc (Table 1, models 1-5), consistent with Voje’s findings. The maximum-likelihood parameter estimate for the decay parameter of the DE model implies a roughly 7-fold decrease in the step variance over the course of the 14 kyr sequence (Embedded Image, Fig. 3). Examination of residuals indicate that decrease in the stochastic component over time is an important signal in this dataset. Models without this dynamic show residuals with elevated spread early in the sequence (Fig. 4 as an example from the RW model). In contrast, residuals from the DE model are not structured in this way, showing a pattern closer to the ideal uniform spread (Fig.5). A model with a single step-down in variance at 10 ka is also consistent with these data, fitting very slightly better than the DE model (Table 1, model 6, is about 0.3 units of AICc better than DE). Parameter estimates from this model indicate that the step variance decreases from 2.77 * 10−5 to 1.40 * 10−5, about a 2-fold drop (Fig. 3).

    Table 1:

    Model fits to spine counts of the Stephanodiscus yellowstonensis lineage. From left to right, columns give model abbreviations, log-likelihoods, number of parameters, AICc scores, ΔAICc scores, and Akaike weights. Model abbreviations: RW = random walk, DT = directional trend, OU = Ornstein-Uhlenbeck, DE = decelerating evolution, RWshift = random walk with a shift in the step variance parameter at 10 ka, OUpollen = OU model in which the trait optimum tracks the proportion of the pollen comprised of Pinus contortus, OUinsol = OU model in which the trait optimum tracks solar insolation, OUinsol−shift = OU model in which the trait optimum tracks solar insolation with a shift in the step variance at 10 ka.

    Table 2:

    Maximum-likelihood estimates (MLE) and confidence intervals (CI) for the best-fitting model: an OU process in which the position of the optimum depends linearly on the value of summer insolation and the step variance is estimated separately before  and after  10,000 years ago. b0 and b1 are the intercept and slope of the relationship between solar insolation and the trait optimum, α represents the force of attraction to that optimum, and anc is the estimated trait value at the start of the time-series.

    Figure 3:

    Modeled changes in the step variance predicted by the two models for which this parameter varies over time. The DE model has a step variance that exponentially decreases, whereas the RWshift model posits a single decrease in step variance that occurs at 10 Ka.

    Figure 4:

    Residuals from the RWmodel. Note the greater spread of residuals early in the sequence.

    Figure 5:

    Residuals from the DE model. Note they are much less structured over time compare to those from the RW model.

    Comparing the two covariate tracking models, it is more plausible that spine counts tracked June solar insolation than the pollen data (Table 1, models 8 and 7, AICc difference of 4.9). The OUinso model shares the features of the OU model, except that the fitness optimum varies with solar insolation instead of being constant over time. The large increase in support between the OU to OUinso (Table 1, ΔAICc = 6.1) is therefore a measure of the importance of solar insolation in accounting for these observations. Combining this insolation-tracking dynamic with a step decrease in stochastic variance results in a model that is best supported overall (Table 1, model 9). This model implies a dynamic where spine counts deterministically follow summer insolation, with overlaid stochastic evolution that is initially high, but then decreases later on.

    4 Discussion

    4.1 Using the State Space Modeling Framework

    We have presented a novel framework based on stochastic integrals and linear state space models for describing, simulating, and analyzing five models for univariate trait evolution in fossil lineages. We have shown how the stochastic integral approach provides a clearer conceptual basis for relating underlying parameters stated in continuous time to models incorporating discretized sampling and observational error. In particular, we showed how a property of the sampling regime, the duration between observed time points (Δ), enters into the system equations of the Ornstein-Uhlenbeck model. Looking forward, the framework is naturally generalizable to multivariate systems.

    The five base models considered here have all been implemented before for paleontological time-series, and they have been fit using maximum likelihood (hunt 2006; Voje 2020). The present approach, using state space models and the Kalman filter, offers an alternative means to compute these same model likelihoods. The two approaches will yield log-likelihoods that are the same (within a constant) and the resulting maximum likelihood parameter estimates are equivalent, within precision of the hill-climbing algorithm used to search for the best parameter values. Two practical benefits of using the state space approach are that the Kalman filter calculations (i) naturally produce residuals useful for assessing model adequacy, and (ii) do not require inverting a large T by T matrix, where T is the number of samples in the time-series. The second benefit applies mostly to rather long time-series (T > 100), as this inversion becomes slow and is prone to fail for very large matrices.

    The analysis of spine counts in Stephanodiscus yellowstonensis is a good example of a typical workflow with the state space approach. An initial set of models were considered, drawn from existing theory and prior interpretations of the system under study. Model fits, as well as examination of residuals, suggest that there are two important signals in the data captured by these models: (i) a decrease over time in the stochastic component of evolutionary change, and (ii) a correlation between diatom phenotypes and summer solar insolation. The Kalman filter calculations allowed us to quickly implement a model with both of these components, which turned out to be the best supported among those considered. The modular nature of the Kalman filter can thus facilitate model development, as it is allows users to easily combine evolutionary components into new models of trait evolution.

    We added a single step-down in variance, rather than the exponential decrease of the DE model, to the OU covariate-tracking model. The step decrease in variance is slightly favored over the DE model, but the decision was also a practical one; its incorporation into the OU with covariate tracking is more straightforward. Given the near-equivalance in model support between the DE and this discrete shift (Table 1), it is unlikely that these data could discriminate between the two different ways of modeling a reduction in the step variance over time.

    4.2 Microevolution in S. yellowstonensis

    The set of evolutionary models considered here are usually interpreted as phenomenological, not mechanistic. Phenomenological models are useful for representing, using just a few parameters, qualitatively different kinds of dynamics, such as meandering change (random walk), fluctuations around a stable mean (stasis), or directional trends. In some cases, these models can be shown to be the expected outcome of specific microevolutionary scenarios (e.g., neutral genetic drift will produce a random walk). But these models are usually interpreted descriptively, rather than as the outcome of specific microevolutionary mechanisms.

    One potential exception is the OU model. Under a set of simplified but reasonable assumptions, this model describes the expected outcome of a population evolving in the vicinity of a peak in the adaptive landscape (Lande 1976). The peak corresponds to the trait value that results in highest mean population fitness. This peak is stable in the OU model, but changes with extrinsic variables in the covariate-tracking versions implemented here (OUinsol, OUpollen). If the best-fitting of these OUinsol – were a complete description of the microevolutionary process, its parameters can be related to population genetic parameters related to the strength of natural selection (from the α parameter) and the effective population size (Ne, from σ2) as described by Hunt et al. 2008 (Hunt, Bell, and Travis 2008)[page=10].

    Under the above set of assumptions, the natural selection is inferred to be rather weak. Mean fitness decreases only 1% or less for population means three standard deviations away from the optimum (Table 3). Some caution should be exercised here because the timescale of adaptation is inferred to be rapid relative to temporal sampling resolution. Its half-life – the amount of time it would take the population to progress half-way to the optimum – is only about 280 years (Table 3), which is close to the median spacing between samples (258 yr). As a result, selection could be much stronger than what is estimated but we would not be able to detect it without finer temporal resolution.

    Table 3:

    Estimates of microevolutionary parameters calculated from the OU model in which the trait optimum follows summer insolation, with a step down in step variance at 10 ka (model 9 in Table 1). Shown are calculations assuming low (0.1) and high (0.7) plausible values of trait heritability and 10 generations per year. ω is the computed variance of the population fitness function; larger values indicate broader fitness curves and therefore weaker stabilizing selection. Fitness reduction is the resulting decrease in population mean fitness between the optimal trait value and the trait values corresponding to three population standard deviations away from the optimum. Effective population size (Ne) is computed separately before and after 10 ka.

    The other population genetic parameter that can be calculated is the effective population size, Ne. This parameter determines the magnitudes of change due to genetic drift in these models. Drift is more potent in smaller populations and thus lower Ne correspond to larger stochastic changes (= higher step variances) around the adaptive trajectory of OU models. The fit of the best model implies about a 3x increase in effective population size at 10 ka. This is consistent with a stepwise increase in the abundance of S. yellowstonensis observed at this time [(Edward C. Theriot et al. 2006)].

    Although the direction of this change is consistent with an increase in the observed absolute abundances of this lineage, the magnitudes of estimated effective population sizes, ranging from 102 to 104 (Table 3), seem rather low for these unicellular algae, which can be found living at abundances high enough to produce that many individuals in just 100 liters of water, or less. (See p. 679 of (Interlandi, Kilham, and Edward C Theriot 1999)). However, it is important to note that effective population size is generally much lower than census population size, with the discrepancy between the two increasing with fluctuations in population size, differences in fitness across individuals, inbreeding, and other factors. The two literature estimates of Ne for diatoms are for widespread marine species and are about 107, which although very high, is still orders of magnitude lower than their peak absolute abundances (Krasovec, Sanchez-Brosseau, and Piganeau 2019). It is likely that these Ne estimates are unrealistically low, but the population genetics of lake diatoms is not well enough investigated to be sure of this.

    Assuming that the Ne values computed from the best model are unrealistically low, then genetic drift would not be sufficient to account for the stochastic component of spine count evolution. Therefore other factors, in addition to June insolation, will have caused changes in the position of the adaptive optimum for spine counts. Theriot et al.’s (Edward C. Theriot et al. 2006) presentation of the paleoenvironmental record provides a detailed account of environmental variation that might contribute to these insolation-independent evolutionary change. In particular, the decrease in the stochastic evolutionary component after 10 ka may be explained by the shift to more stable conditions at this time, leading to more modest selective fluctuations diatom morphology. Voje (Voje 2020) offers an alternative explanation in which ecological opportunity is initially high, perhaps because of the phenotypic changes associated with the origin of the new species S. yellowstonensis. With high ecological opportunity, stabilizing selection may be weakened, permitting greater variation and larger evolutionary steps. This interpretation is consistent with the observation that standing variation in spine counts is initially high and decreases steadily for the first 3 or 4 kyr of the sequence (Fig. 6). These two explanations for the reduction in stochastic evolutionary change – decreasing environmental variation and decreasing ecological opportunity – are not mutually exclusive.

    Figure 6:

    Variance in spine count decreases in the initial part of the S. yellowstonensis sequence.

    4.3 Handling covariates in models of trait evolution

    Although a model in which spine counts follow solar insolation as an OU process is the best supported among those considered, we caution that aspects of this model may make it less suitable for some situations. Our implementation requires an assumption that the position of the trait optimum is constant in between the time points at which we have observations. This is reasonable when, as is the case here, the covariates show point-to-point changes that are small compared to the total range of the time-series. This assumption will be less realistic for covariates that fluctuate widely on short time scales. In addition, the microevolutionary intepretation of the OU dynamics is only tenable if the sampling resolution is fine enough to potentially capture the adaptive dynamics of the population chasing the moving fitness peak. Even with the exceptional temporal resolution of the S. yellowstonensis data, the evolutionary dynamics may be too rapid to well-constrain the microevolutionary dynamics. In addition, this modeling approach may be prone to receiving spuriously strong support when applied to traits that show clearly directional change if analyzed with covariates that are also trended. Including the simpler model of a trend, as done here, may protect against this effect as the fewer parameters of the trend model give it an AICc advantage.

    State space models are flexible enough to allow for other approaches to incorporate the effects of exogenous covariates that may be more suitable in other circumstances. For example, one could model a trait as an unbiased random walk, with an additional pulse of change that is proportional to changes in a covariate, implemented through the ϒ term of the state equation. Such a modeling approach does not attempt to capture the dynamics of a population climbing an adaptive peak and instead would be consistent with an assumption that enough time has elapsed between samples for the population to have reached the adaptive optimum. This approach would therefore would be more appropriate for trait time-series at more typical paleontological resolutions. And the use of changes in covariates as input variables, rather than covariates themselves, would render this approach less susceptible to trended sequences as described above.

    5 Conclusion

    The state space framework provides a practical approach for analyzing phenotypic evolution in fossil lineages that facilitates model development incorporating exogenous environmental variables based on easier access to residuals as a diagnostic tool for model fit. We highlighted some additional useful features, especially the ease of accessing time series residuals and enhanced numerical stability and efficiency. Our analysis suggested a novel biological interpretation of evolution in spine count for S. yellowstonensis based on stabilizing selection to changing solar insolation levels. Our focus on univariate trait models in this paper provides a foundation for expanding into more complex, multi-variate models that, for example, allow for estimation of trait covariances in a time series setting. This is essential to investigate how solar insolation that may jointly affect all three traits measured by Theriot et al.

    6. Appendix 1

    Common models as linear state space models

    6.1 Stochastic Integral Models and Their Discretizations

    In this section, we will discuss a number of well-known models of trait evolution in fossil lineages and give a corresponding continuous time equivalent. We will then show that each of these models can be exactly discretized to correspond to discretely observed data.

    Each of these standard models may be expressed as an Ito integral with a deterministic integrand. We can then look at how each can then be expressed as a linear state space model when observed at discrete time points.

    First, the definition of a Ito integral for a deterministic integrand is:Embedded Imagewhere f (t, s) is a deterministic function and W (s) is a standard Weiner process (Oksendal 2013). In other words, for each s, W (s) is a normal random variable with zero mean and variance equal to s. This process also has the property of independent increments, implying that W (t) −W (s) is independent W (v) −W (u) as long as (s, t] and (u, v] do not overlap. The fact that W (s) has zero mean for each s implies that X (t) does also. Each of the models that we examine in this manuscript is a Gaussian process and using this integral representation allows us to express the models in a unified way.

    An important property of this integral, especially when calculating variances, is the Ito isometry (Oksendal 2013).Embedded Image

    Now, we write some basic models in the form of such an integral along with their exact discretizations and approximate discretizations where appropriate.

    6.1.1 Random Walk

    For the RW model, the deterministic function is simply a constant, σ, andEmbedded Image

    Note then, that a discretely observed version of this model, assuming equal spacing in time, would beEmbedded Imagewhich we could we write as follows by using the fact of independent increments of Brownian motion:Embedded Imagewhere εn are a sequence of independent standard normal random variable and Δ is the amount of time between observations. This can be derived directly from the integral definition above with the Embedded Image corresponding to W (nΔ) −W ((n − 1)Δ). In this discrete version, Xn corresponds to X (nΔ).

    6.1.2 Directed Random Walk

    For this model, we add a deterministic linear directionality to the Random Walk model.Embedded Image

    In a similar way then, a discretely observed version of this model, again assuming equal spacing in time, would beEmbedded Imagewith similar interpretations as for the Random Walk.

    6.1.3 Decelerated Evolution

    Voje introduced a model for the evolution of a trait where the step variance of a random walk declines exponentially. In other words, a model that could be describe with the recursionEmbedded Imagewhere σV and r are positive parameters and εn are a sequence of standard normal random variables. A natural way to write this as a stochastic integral would beEmbedded Image

    An Euler approximation of this integral would then be defined by the recursionEmbedded Imageand if we identify Embedded Image we notice that this approximation corresponds to Voje’s original definition. However, we can write down an exact discretization of the stochastic integral.Embedded Image

    We can again look at the second term on the right and calculate the variance using Ito’s isometryEmbedded Image

    So, solving the integral on the right hand side of the above equation, we findEmbedded Image

    Applying a Taylor expansion with Δ small, we see thatEmbedded Imageand this corresponds to the Euler approximation above.

    6.1.4 Stasis

    The stasis model assumes that each observation is independent and identically distributed, typically with a normal distribution. So, Xn is normal with mean c and variance Embedded Image, which we could write asEmbedded Imagewhere the εn are a sequence of independent and identically distributed random variables. Effectively, there is no continuous version of this model. One way to think of this model, however, is as a discretely observed Ornstein-Uhlenbeck process with sufficient spacing between the observations relative to the parameter α. The Ornstein-Uhlenbeck process is stationary, meaning that after a sufficient time, the initial condition is trivially relevant and the relationship between observations at two time depends only on the distance between those observations. Since the Ornstein-Uhlenbeck process is Gaussian, the variance and covariances define the process. The covariance between two observations of a (zero-mean) Ornstein-Uhlenbeck process is given byEmbedded Image

    So, if the observations are sufficiently spaced this covariance is effectively zero (as long as α is not too small). So, the stasis model could be viewed as an Ornstein-Uhlenbeck process that is approximately stationary and sufficiently spaced with the relationship between the variances being Embedded Image where σ2 is the infinitesimal variance of the Ornstein-Uhlenbeck model. If we look at equation 5, with αΔ being large, then the we see that the OU model is effectively the stasis model.

    7 Appendix

    Simulations and Model Selection

    To ensure comparability of results, we re-used code from previous studies to simulate data from the five models, including the paleoTS R package maintained by Hunt and supplementary materials from (Voje 2020). We used paleoTS to simulate data for the ST, RW, DT, and OU models and code from (Voje 2020) to simulate data for the DE model. However, we found that prior simulation studies explored parameter values and data sampling regimes in a way that varied both observational and biological time scales simultaneously (Hunt 2006; Hunt 2008; Voje 2020), making their separate effects on performance difficult to disentangle. We chose simulation parameters to feature substantial levels of model uncertainty in order to illustrate dependence on sampling scales.

    For model selection, we calculated model goodness-of-fit using the corrected Akaike Information Criterion (AICc), which is modified for better performance in small sample sizes (Hunt 2008). In general, the AIC provides an unbiased estimator of a model’s expected likelihood, and picking the model with the best (i.e. lowest) AIC score will asymptotically converge on the true distribution when it is unique to a single model (Burnham and Anderson 2002). We also calculated Akaike weights for each model using the AICc scores (Wagenmakers and Farrell 2004), which approximate the probability that a model is the best out of the candidates considered. For model estimation, we used the linear Kalman filter in the state space models and compared the results of our procedures to those in PaleoTS and Voje’s code.

    7.1 Parameter estimation

    We find that estimation accuracy and uncertainty are not uniformly influenced by sampling time scales. We show simulation results for the RW and DT models (Fig. 7), OU model (Fig. 8), and DE model (Fig. 9). We explored three scenarios for modifying sample time scales: first, increasing sampling with constant total duration (shrinking Δ); second, increasing duration and increasing sample size (uniform Δ); and third, increasing duration while holding sample size constant (increasing Δ). The ST model is not shown but is effectively an independent, identically distributed (i.i.d.) process with time-indexed observations, and so estimation of Embedded Image will depend on sample size but not total duration or Δ.

    Figure 7:

    Parameter estimation for random walk and directional trend models using varying ratios of stepwise and total observational time scales. The left column shows estimation of the diffusion parameter for the RW model. The middle and right columns show estimation of the diffusion and the drift parameters, respectively, for the DT model. The true parameter values are σ2 = 25 and µ = 5. We use different values for the sample size N = {20, 40, 80, 160}, the size of the increments Δ = {1/8, 1/4, 1/2, 1}, and the terminal time T = {2.5, 5, 10, 20}. Box plots show 100 replicates

    Figure 8:

    Parameter estimation on simulated data from the Ornstein-Uhlenbeck model for different combinations of sample size, time step, and total duration. Each of the columns shows parameter estimation for α, σ, Z0, and θ, respectively. True values are α =0.50, σ2 = 20, Z0 = 40, and θ = 50, as represented with dashed red lines. We varied the sample size {N = 20, 40, 80, 160 }, the size of the increments {Δ = 1/8, 1/4, 1/2, 1}, and the terminal time T = {2.5, 5, 10, 20}. Box plots show 100 replicates

    Figure 9:

    Parameter estimation for Decelerated Evolution model on simulated data using varying combinations of sample size, time step, and total duration. The columns show parameter estimation of r, σ and Z0, respectively. Dashed red lines show the true values of r =-1, Embedded Image, and Z0 = 40. We used sample sizes N = {20, 40, 80, 160}, time increments of Δ = {1/8, 1/4, 1/2, 1}, and terminal times T = {2.5, 5, 10, 20}. Box plots show 100 replicates.

    Figure 10:

    Model selection performance using the AIC criterion on simulated data for varying ratios of stepwise and total observational time scales. In each panel, the true model is labeled on the x-axis, and the stacked histogram shows the average Akaike weight for each model. Perfect model performance would show each bar as completely filled by the corresponding true model’s shade on the legend (e.g. ST as white, RW as light gray, DT as medium gray, etc.). The true model parameters are ST: θ = 50, ω = 20; RW: Embedded Image DT: Dri f t = 5, Embedded Image OU: θ = 50, α = 20, Embedded Image and DE: r = −1, Embedded Image. The initial condition for all the models is Z0 = 40, the variance of the evolutionary step is Vp = 5, and the vector of population sample size is m = 50. We varied sample size N = {20, 40, 80, 160}, the size of the increments Δ = {1/8, 1/4, 1/2, 1}, and the terminal time T = {2.5, 5, 10, 20}.

    For the RW and DT models, the results in Fig. 7 show that increasing total time while keeping the sample size fixed doesn’t affect Embedded Image but does improve Embedded Image. Note that the sufficient statistic for Embedded Image is the value of the process at the terminal point of the time series, so intermediate values don’t matter for estimation, only the end point. Increasing duration with constant Δ improves both Embedded Image and Embedded Image, but for different reasons: Embedded Image is improving because the total time observed is increasing, but Embedded Image is improving because there are more steps observed.

    In contrast, the OU model parameters in Fig. 8 show several different types of response to time scales. Both Embedded Image and Embedded Image improve for both the scenarios with increasing sampling with constant duration or increasing duration. As Δ becomes larger in row 3, however, estimation gets worse because exp(−Δ) goes to zero and the process starts to look i.i.d., so that Embedded Image and Embedded Image are both being fit to a normal distribution with mean θ and variance σ2 /(2α). The initial value parameter, Embedded Image, is unaffected in all three scenarios because better estimation of the restoring force, α can only improve precision for the initial condition up to a point. Similar to Embedded Image ‘s behavior in the RW and DT models, the OU equilibrium value parameter, Embedded Image, does not converge under the increased sampling intensity.

    Neighboring observations are positively correlated, so adding more time points within a fixed interval gives diminishing returns for estimating the mean, but if T is increasing, the observations are spaced further apart and so are more independent.

    For the DE model, the variance decay parameter Embedded Image in Fig. 9 shows improved precision in the increasing sampling intensity and increasing duration, constant sampling scenarios but remains biased below the true value for the simulation setups we examined. The DE model shows phenomenologically the same behavior in Embedded Image and Embedded Image as the OU model but for different reasons. As the step variance of the process decays exponentially to zero with time, observing the process over a longer duration provides progressively less information.

    7.2 Model selection performance

    Figure 10 shows how model selection performance, measured in terms of the average Akaike weight of the true model, varies with sampling. Columns in the figure show increasing total duration. Rows show denser sampling as total duration shrinks. Diagonals from top-left to bottom-right show increasingly dense sampling within a fixed total duration.

    ST is almost exclusively conflated with OU. Increased data appears to slightly worsen false positives for OU when stasis is true, likely because the AICc has a bias for nested models toward the model with more parameters. RW is most frequently confused with DT and OU models. The average Akaike weight does not vary significantly with sampling time scales, again likely because of the AICc’s bias toward complex models. As expected, the evidence for DT improves significantly as total duration grows. OU shows improvement with greater duration and sampling density. The same is true for DE, which is mainly competitive with ST for small sample sizes.

    Footnotes

    References

    1. Antell, Gwen S. et al. (May 2021). “Thermal niches of planktonic foraminifera are static throughout glacial201A00C400ECinterglacial climate change”. en. In: Proceedings of the National Academy of Sciences 118.18, e2017105118. ISSN: 0027-8424, 1091-6490. DOI:10.1073/pnas.2017105118. URL: http://www.pnas.org/lookup/doi/10.1073/pnas.2017105118 (visited on 08/10/2021).
    2. Blomberg, Simone P., Suren I. Rathnayake, and Cheyenne M. Moreau (Feb. 2020). “Beyond Brownian Motion and the Ornstein-Uhlenbeck Process: Stochastic Diffusion Models for the Evolution of Quantitative Characters”. In: The American Naturalist 195.2. Publisher: The University of Chicago Press, pp. 145–165. ISSN: 0003-0147. DOI:10.1086/706339. URL: https://www.journals.uchicago.edu/doi/full/10.1086/706339 (visited on 07/10/2020).
    3. Brombacher, Anieke et al. (June 2017). “The Breakdown of Static and Evolutionary Allometries during Climatic Upheaval”. In: The American Naturalist 190.3. Publisher: The University of Chicago Press, pp. 350–362. ISSN: 0003-0147. DOI:10.1086/692570. URL: https://www.journals.uchicago.edu/doi/full/10.1086/692570 (visited on 01/05/2021).
    4. Brombacher, Anieke et al. (July 2018). “Temperature is a poor proxy for synergistic climate forcing of plankton evolution”. In: Proceedings of the Royal Society B: Biological Sciences 285.1883. Publisher: Royal Society, p. 20180665. DOI:10.1098/rspb.2018.0665. URL: https://royalsocietypublishing.org/doi/full/10.1098/rspb.2018.0665 (visited on 05/18/2020).
    5. Burnham, K P and David R Anderson (2002). Model Selection and Inference: A Practical Information-Theoretic Approach. New York: Springer-Verlag.
    6. Gingerich, Philip D. (2019). Rates of Evolution: A Quantitative Synthesis. Cambridge: Cambridge University Press. ISBN: 978-1-107-16724-7. DOI:10.1017/9781316711644. URL: https://www.cambridge.org/core/books/rates-of-evolution/D1CC64790C70A249E39D034EFABAA551 (visited on 01/04/2021).
    7. Goswami, Anjali and Julien Clavel (Jan. 2024). “Morphological evolution in a time of Phenomics”. en. In: Publisher: EcoEvoRxiv. URL: https://ecoevorxiv.org/repository/view/6507/ (visited on 02/24/2024).
    8. Hannisdal, Bjarte and Lee Hsiang Liow (July 2018). “Causality from palaeontological time series”. en. In: Palaeontology 61.4. Ed. by Andrew Smith, pp. 495–509. ISSN: 00310239. DOI:10.1111/pala.12370. URL: http://doi.wiley.com/10.1111/pala.12370 (visited on 03/09/2020).
    9. Hansen, Thomas F (1997). “Stabilizing selection and the comparative analysis of adaptation”. In: Evolution 51.5, pp. 1341–1351.
    10. Harmon, Luke J. et al. (Nov. 2021). “Causes and Consequences of Apparent Timescaling Across All Estimated Evolutionary Rates”. In: Annual Review of Ecology, Evolution, and Systematics. Publisher: Annual Reviews. ISSN: 1543-592X. DOI:10.1146/annurev-ecolsys-011921-023644. URL: https://www.annualreviews.org/doi/10.1146/annurev-ecolsys-011921-023644 (visited on 09/14/2021).
    11. He, Yichen et al. (Jan. 2024). “Challenges and opportunities in applying AI to evolutionary morphology”. en. In: Publisher: EcoEvoRxiv. URL: https://ecoevorxiv.org/repository/view/6525/ (visited on 02/24/2024).
    12. Hunt, Gene (2006). “Fitting and comparing models of phyletic evolution: random walks and beyond”. en. In: Paleobiology 32.4, pp. 578–601. ISSN: 0094-8373, 1938-5331. DOI:10.1666/05070.1. URL: https://www.cambridge.org/core/product/identifier/S0094837300018832/type/journal_article (visited on 03/09/2020).
    13. Hunt, Gene (Oct. 2008). “Evolutionary Patterns Within Fossil Lineages: Model-Based Assessment of Modes, Rates, Punctuations and Process”. en. In: The Paleontological Society Papers 14, pp. 117–131. ISSN: 1089-3326, 2399-7575. DOI:10.1017/S1089332600001649. URL: https://www.cambridge.org/core/product/identifier/S1089332600001649/type/journal_article (visited on 03/09/2020).
    14. Hunt, Gene (2012). “Measuring rates of phenotypic evolution and the inseparability of tempo and mode”. en. In: Paleobiology 38.3, pp. 351–373. ISSN: 0094-8373, 1938-5331. DOI:10.1666/11047.1. URL: https://www.cambridge.org/core/product/identifier/S0094837300000695/type/journal_article (visited on 03/09/2020).
    15. Hunt, Gene, Michael A. Bell, and Matthew P. Travis (Mar. 2008). “Evolution Toward A New Adaptive Optimum: Phenotypic Evolution in A Fossil Stickleback Lineage”. en. In: Evolution 62.3, pp. 700–710. ISSN: 00143820, 15585646. DOI:10.1111/j.1558-5646.2007.00310.x. URL: http://doi.wiley.com/10.1111/j.1558-5646.2007.00310.x (visited on 03/09/2020).
    16. Hunt, Gene, Melanie J. Hopkins, and Scott Lidgard (Apr. 2015). “Simple versus complex models of trait evolution and stasis as a response to environmental change”. en. In: Proceedings of the National Academy of Sciences 112.16, pp. 4885–4890. ISSN: 0027-8424, 1091-6490. DOI:10.1073/pnas.1403662111. URL: http://www.pnas.org/lookup/doi/10.1073/pnas.1403662111 (visited on 03/09/2020).
    17. Interlandi, Sebastian J, Susan S Kilham, and Edward C Theriot (1999). “Responses of phytoplankton to varied resource availability in large lakes of the Greater Yellowstone Ecosystem”. In: Limnology and oceanography 44.3, pp. 668–682.
    18. Kearns, Lorna E. et al. (July 2021). “Searching for Function: Reconstructing Adaptive Niche Changes Using Geochemical and Morphological Data in Planktonic Foraminifera”. en. In: Frontiers in Ecology and Evolution 9, p. 679722. ISSN: 2296-701X. DOI:10.3389/fevo.2021.679722. URL: https://www.frontiersin.org/articles/10.3389/fevo.2021.679722/full (visited on 11/16/2021).
    19. Krasovec, Marc, Sophie Sanchez-Brosseau, and Gwenael Piganeau (2019). “First estimation of the spontaneous mutation rate in diatoms”. In: Genome biology and evolution 11.7, pp. 1829–1837.
    20. Lande, Russell (1976). “Natural Selection and Random Genetic Drift in Phenotypic Evolution”. In: Evolution 30.2. Publisher: [Society for the Study of Evolution, Wiley], pp. 314–334. ISSN: 0014-3820. DOI:10.2307/2407703. URL: https://www.jstor.org/stable/2407703 (visited on 07/10/2023).
    21. Oksendal, Bernt (2013). Stochastic differential equations: an introduction with applications. Springer Science & Business Media.
    22. Porto, Arthur and Kjetil L. Voje (2020). “ML-morph: A fast, accurate and general approach for automated detection and landmarking of biological structures in images”. en. In: Methods in Ecology and Evolution 11.4. eprint: https://besjournals.onlinelibrary.wiley.com/doi/pdf/10.1111/2041-210X.13373, pp. 500–512. ISSN: 2041-210X. DOI: 10.1111/2041-210X.13373. URL: https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.13373 (visited on 02/15/2021).
    23. Reitan, Trond and Lee Hsiang Liow (Sept. 2019). “layeranalyzer: inferring correlative and causal connections from time series data in R”. In: Methods in Ecology and Evolution 10. DOI:10.1111/2041-210X.13299.
    24. Reitan, Trond, Tore Schweder, and Jorijntje Henderiks (Dec. 2012). “Phenotypic evolution studied by layered stochastic differential equations”. EN. In: Annals of Applied Statistics 6.4. Publisher: Institute of Mathematical Statistics, pp. 1531–1551. ISSN: 1932-6157, 1941-7330. DOI:10.1214/12-AOAS559. URL: https://projecteuclid.org/euclid.aoas/1356629050 (visited on 01/04/2021).
    25. Shumway, Robert H, David S Stoffer, and David S Stoffer (2000). Time series analysis and its applications. Vol. 3. Springer.
    26. Simpson, George Gaylord (1944). Tempo and mode in evolution. 15. Columbia University Press.
    27. Stuart, Yoel E., Matthew P. Travis, and Michael A. Bell (Nov. 2020). “Inferred genetic architecture underlying evolution in a fossil stickleback lineage”. en. In: Nature Ecology & Evolution 4.11. Number: 11 Publisher: Nature Publishing Group, pp. 1549–1557. ISSN: 2397-334X. DOI:10.1038/s41559-020-01287-x. URL: https://www.nature.com/articles/s41559-020-01287-x (visited on 01/26/2021).
    28. Theriot, Edward C. et al. (Jan. 2006). “Late Quaternary rapid morphological evolution of an endemic diatom in Yellowstone Lake, Wyoming”. In: Paleobiology 32.1. Publisher: The Paleontological Society, pp. 38–54. ISSN: 0094-8373, 1938-5331. DOI:10.1666/02075.1. URL: https://bioone.org/journals/paleobiology/volume-32/issue-1/02075.1/Late-Quaternary-rapid-morphological-evolution-of-an-endemic-diatom-in/10.1666/02075.1.full (visited on 02/15/2021).
    29. Turley, Lydia Marissa (2020). “Ornstein Uhlenbeck-Jump Models of Evolution”. en. Accepted: 2020-09-22T21:55:39Z. Thesis. University of Otago. URL: https://ourarchive.otago.ac.nz/handle/10523/10401 (visited on 04/29/2022).
    30. Voje, Kjetil L. (2018). “Assessing adequacy of models of phyletic evolution in the fossil record”. en. In: Methods in Ecology and Evolution 9.12. eprint: https://besjournals.onlinelibrary.wiley.com/doi/pdf/10.210X.13083, pp. 2402–2413. ISSN: 2041-210X. DOI:10.1111/2041-210X.13083. URL: https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.13083 (visited on 08/26/2020).
    31. Voje, Kjetil L. (Jan. 2020). “Testing eco-evolutionary predictions using fossil data: Phyletic evolution following ecological opportunity*”. en. In: Evolution 74.1, pp. 188–200. ISSN: 0014-3820, 1558-5646. DOI:10.1111/evo.13838. URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/evo.13838 (visited on 06/05/2020).
    32. Wagenmakers, Eric-Jan and Simon Farrell (Feb. 2004). “AIC model selection using Akaike weights”. en. In: Psychonomic Bulletin & Review 11.1, pp. 192–196. ISSN: 1069-9384, 1531-5320. DOI:10.3758/BF03206482. URL: http://link.springer.com/10.3758/BF03206482 (visited on 11/07/2021).

    Journal of Modern Sciences

    info@modernsciences.org

    About
    Modern Sciences is a media company that aims to provide the academics, researches, and science enthusiasts with the latest news and stories in the world of Applied and Pure Science.