Title:	Analyze Paleontological Time-Series
Version:	0.6.2
Description:	Facilitates analysis of paleontological sequences of trait values. Functions are provided to fit, using maximum likelihood, simple evolutionary models (including unbiased random walks, directional evolution,stasis, Ornstein-Uhlenbeck, covariate-tracking) and complex models (punctuation, mode shifts).
Depends:	R (≥ 3.4.0)
License:	GPL-3
Encoding:	UTF-8
LazyData:	true
Imports:	mnormt, foreach, parallel, doParallel
RoxygenNote:	7.3.1
Suggests:	knitr, rmarkdown
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2024-09-10 16:33:24 UTC; hunte
Author:	Gene Hunt [aut, cre, cph], John Fricks [ctb]
Maintainer:	Gene Hunt <hunte@si.edu>
Repository:	CRAN
Date/Publication:	2024-09-10 17:20:02 UTC

paleoTS: Analyze Paleontological Time-Series

Description

Facilitates analysis of paleontological sequences of trait values. Functions are provided to fit, using maximum likelihood, simple evolutionary models (including unbiased random walks, directional evolution,stasis, Ornstein-Uhlenbeck, covariate-tracking) and complex models (punctuation, mode shifts).

Author(s)

Maintainer: Gene Hunt hunte@si.edu (ORCID) [copyright holder]

Other contributors:

• John Fricks jfricks@asu.edu [contributor]

Compute Expected Squared Divergence (ESD) for Evolutionary Models

Description

Computes for a specified model and duration of time the expected squared divergence (ESD), which is a useful measure of the magnitude or rate of change across different models.

Usage

ESD(
  y,
  dt,
  model = c("GRW", "URW", "Stasis", "allThree"),
  method = c("Joint", "AD"),
  pool = TRUE,
  ...
)

Arguments

Details

Hunt (2012) argued that rate metrics make sense only in the context of specific models of evolution. It is thus difficult to meaningfully compare rates across sequences generated by different evolutionary processes. ESD values can be used for a specified model and duration as a comparable measure of the amount of evolutionary change that is expected. Acceptable values for the model argument can be "GRW" for the general random walk (directional change), "URW" for the unbiased random walk, and "Stasis." In addition, one can also specify "allThree", in which case all these models will be fit and the resulting ESD will be the weighted average of them, using model support (Akaike weights) for the weighting (see Hunt [2012], p. 370)

Value

the ESD value

References

Hunt, G. 2012. Measuring rates of phenotypic evolution and the inseparability of tempo and mode. Paleobiology 38:351–373.

Examples

  x<- sim.GRW(ns=20)
  esd.urw<- ESD(x, dt=10, model="URW")
  esd.all<- ESD(x, dt=10, model="allThree")

Compute Information Criteria

Description

Compute Information Criteria

Usage

IC(logL, K, n = NULL, method = c("AICc", "AIC", "BIC"))

Arguments

Value

the value of the specified information criterion

Note

This function is used internally by the model-fitting functions. It will not generally be called directly by the user.

Examples

ic1 <- IC(logL = 0, K = 2, method = "AIC") # plain AIC
ic2 <- IC(logL = 0, K = 2, n = 10, method = "AICc")
ic3 <- IC(logL = 0, K = 2, n = 1000, method = "AICc") # converges to AIC with increasing n
print(rbind(ic1, ic2, ic3))

Time-varying Kalman filter calculations

Description

Time-varying Kalman filter calculations

Usage

Kfiltertv(num, y, Atv, mu0, Sigma0, Phitv, Ups, Gam, Qtv, Rtv, input)

Arguments

Details

For the dimensions of the argument arrays, n is the length of the time-series, q is the dimension of the observation variable(s), p is the dimension of the state variable(s), and r isthe dimension of the input variable(s).

This function is based on the Kfilter function of the astsa package, modified modified to allow for time-varying terms for the Kalman filter. This modification facilitates fitting a broader array of models and handling non-uniform temporal spacing of samples. See the documentation for that function, and the reference below for additional information.

Value

A list of the following elements:

• xp one-step-ahead prediction of the state

• Pp mean square prediction error

• xf filter value of the state

• Pf mean square filter error

• like log-likelihood

• innov innovation series

• sig innovation covariances

• Kn last value of the gain, needed for smoothing

Note

This function is used in the internal SSM log-likelihood functions for the models. The user will not need to use this they create their own model-fitting functions.

Author(s)

John Fricks (jfricks@asu.edu)

References

Shumway, R. H., and D. S. Stoffer. 2017. Time Series Analysis and its Applications (4th Ed.) Springer International.

Examples

y <- sim.GRW(ms = 0, vs = 1, vp = 0)
n <- length(y)
kf <- Kfiltertv(n ,y = y$mm, Atv = array(1, dim = c(1,1,n)), mu0 = y$mm[1],
                Sigma0 = y$vv[1]/y$nn[1], Phitv = array(1, dim = c(1,1,n)),
                Ups = NULL, Gam = NULL, Qtv = array(1, dim = c(1,1,n)),
                Rtv = array(0, dim = c(1,1,n)), input = NULL)

Log-rate, Log-interval (LRI) method of Gingerich

Description

Gingerich (1993) introduced a method that plots on log-log scale, the rate and interval for each pair of samples in an evolutionary sequence. On this plot, the slope is interpreted as an indicator of evolutionary mode (-1 for stasis, 0.5 for random walk, 0 for directional), and the intercept is interpreted as a measure of the rate of evolution over one generation.

Usage

LRI(y, gen.per.t = 1e+06, draw = TRUE)

Arguments

Details

Following Gingerich (1993), a robust line is fit through the points by minimizing the sum of absolute deviations. If generations are one year long and time is measured in Myr, gen.per.t= 1e6.

Value

A named vector with three elements: Intercept, slope, and GenerationalRate

Note

This method was important in early attempts to disentangle evolutionary tempo and mode. Likelihood-based methods have a more sound statistical basis, and in particular the estimation of 'Generational Rates' using LRI is compromised by sampling error; see Hunt (2012) and the example below.

References

Gingerich, P.D. 1993. Quantification and comparison of evolutionary rates. American Journal of Science 293-A:453–478.

Hunt, G. 2012. Measuring rates of phenotypic evolution and the inseparability of tempo and mode. Paleobiology 38:351–373.

See Also

lynchD

Examples

set.seed(1)
xFast <- sim.GRW(ns = 20, ms = 0.5, vs = 0.2)  # fast evolution
xSlow <- sim.Stasis(ns = 20, omega = 0)        # strict stasis (zero rates)
lri.Fast <- LRI(xFast, draw = FALSE)
lri.Slow <- LRI(xSlow, draw = FALSE)
print(lri.Fast[3], 4)
print(lri.Slow[3], 4)  # LRI thinks strict stasis rates are MUCH faster!

Compute Akaike weights from AIC scores

Description

Compute Akaike weights from AIC scores

Usage

akaike.wts(aa)

Arguments

Details

This function converts a vector of AIC or AICc scores to Akaike weights, a vector that summarizes proportional model support.

Value

a vector of Akaike weights

See Also

compareModels

Examples

akaike.wts(c(10, 14, 20))

Make a Paleontological Time-series object

Description

Combines information into an object of class paleoTS

Usage

as.paleoTS(
  mm,
  vv,
  nn,
  tt,
  MM = NULL,
  genpars = NULL,
  label = NULL,
  start.age = NULL,
  oldest = c("first", "last"),
  reset.time = TRUE
)

Arguments

Details

This function combines data into a paleoTS object. For empirical data it may be more convenient to use read.paleoTS.

If sample ages decrease through the sequence, as if given in millions of years ago, tt will automatically be converted to time elapsed from the beginning of the sequence as long as reset.time = TRUE.

Value

a paleoTS object

Note

All model-fitting functions estimate the contribution of sampling noise to the observed differences between samples. They do this assuming that the trait is represented by sample means, which have sampling variances equal to variance divided by sample size, vv/nn. If one is interested in analyzing statistics other than the sample mean (medians, quantiles, or other statistics), use the the following procedure: set the statistic in question as the mm values, replace vv with a vector of the squared standard errors for each estimate (generated by other means, for example bootstrapping), and set all values of nn to one. When fitting such time-series, be sure to set the argument pool = FALSE.

See Also

read.paleoTS

Examples

x <- as.paleoTS(mm = rnorm(20), vv = rep(1, 20), nn = rep(25, 20), tt=0:19)
plot(x) # easier to use sim.Stasis()

Create a paleoTSfit object

Description

Create a paleoTSfit object

Usage

as.paleoTSfit(
  logL,
  parameters,
  modelName,
  method,
  K,
  n,
  se,
  convergence = NULL,
  logLFunction = NULL,
  ...
)

Arguments

Value

a paleoTSfit object

Note

All model-fitting functions use this function to package the resulting data-model fits. Users will not need to call this function.

Examples

#  fake example; users won't need to use this unless they make their own model-fitting functions
w <- as.paleoTSfit(logL = 10, parameters = 2, modelName = "StrictStasis",
                   method = "Joint", K = 1, n = 25, se = NULL)

Bootstrap test to see if a complex model is significantly better than a simple model.

Description

Bootstrap test to see if a complex model is significantly better than a simple model.

Usage

bootSimpleComplex(
  y,
  simpleFit,
  complexFit,
  nboot = 99,
  minb = 7,
  ret.full.distribution = FALSE,
  parallel = FALSE,
  ...
)

Arguments

Details

Simulations suggest that AICc can be overly liberal with complex models with mode shifts or punctuations (Hunt et al., 2015). This function implements an alternative of parametric boostrapping to compare the fit of a simple model with a complex model. It proceeds in five steps:

1. Compute the observed gain in support from the simple to complex model as the likelihood ratio, LR_{obs} = -2(logL_{simple} - logL_{complex})

2. Simulate trait evolution under the specified simple model nboot times

3. Fit to each simulated sequence the specified simple and complex models

4. Measure the gain in support from simple to complex as the bootstrap likelihood ratio for each simulated sequence

5. Compute the P-value as the percentile of the bootstrap distribution corresponding to the observed LR.

Argument simpleFit should be a paleoTS object returned by the function fitSimple or similar functions (e.g., opt.joint.GRW, opt.GRW, etc.). Argument complexFit must be a paleoTS object returned by fitGpunc or fitModeShift.

Calculations can be speeded up by setting parallel = TRUE, which uses package doParallel to run the bootstrap replicates in parallel, using one fewer than the number of detected cores.

Value

A list of the observed likelihood ratio statistic, LRobs, the P-value of the test, and the number of bootstrap replicates. If ret.full.distribution = TRUE, the null distribution of likelihood ratios generated by parametric bootstrapping is also returned.

References

Hunt, G., M. J. Hopkins and S. Lidgard. 2015. Simple versus complex models of trait evolution and stasis as a response to environmental change. PNAS 112(16): 4885-4890.

See Also

sim.Stasis.RW, fitModeShift

Examples

## Not run: 
x <- sim.Stasis.RW(ns = c(15, 15), omega = 0.5, ms = 1, order = "Stasis-RW")
ws <- fitSimple(x)
wc <- fitModeShift(x, order = "Stasis-RW", rw.model = "GRW")
bootSimpleComplex(x, ws, wc, nboot = 50, minb = 7)  # nboot too low for real analysis!

## End(Not run)

Time-series of the length of lower first molar for the Cantius lineage

Description

Time-series of the length of lower first molar for the Cantius lineage

Usage

cantius_L

Format

a paleoTS object with the data

Source

Clyde, W. C. and P. D. Gingerich (1994). Rates of evolution in the dentition of early Eocene Cantius: comparison of size and shape. Paleobiology 20(4): 506-522.

Compute and (optionally) plot residuals from SSM model fit

Description

Compute and (optionally) plot residuals from SSM model fit

Usage

checkSSMresiduals(
  y,
  w,
  show.plot = TRUE,
  resid.type = c("standardized", "unstandardized")
)

Arguments

Details

It is recommended that resid.type be set to the default, "standardized", which will scale residuals by their expected standard deviation

Value

a vector of residuals, returned invisibly

Examples

y <- sim.GRW(ns = 50, ms = 0.2)
w <- fitSimple(y, model = "URW", method = "SSM")  # wrong model
checkSSMresiduals(y, w, show.plot = TRUE) # positive residuals show model mis-fit

Compare model fits for a paleontological time-series

Description

Takes output from model-fitting functions and compiles model-fit information (log-likelihood, AICc, etc.) into a convenient table

Usage

compareModels(..., silent = FALSE, sort = FALSE)

Arguments

Value

if silent = FALSE, the table is printed and nothing is returned. If silent = TRUE, printing is suppressed and a list of two objects is returned: the table of model fits, modelFits, and a list of parameter estimates, pl.

Examples

x <- sim.GRW(ns = 40, ms = 0.5, vs = 0.1)
m1 <- fitSimple(x, model = "GRW")  # the true model
m2 <- fitSimple(x, model = "URW")
plot(x, modelFit = m1)
compareModels(m1, m2)

Time-series of dorsal spine data from a fossil stickleback lineage

Description

Time-series of dorsal spine data from a fossil stickleback lineage

Usage

dorsal.spines

Format

a paleoTS object of the mean number of dorsal spines (log-transformed)

Source

Bell, M.A., M.P. Travis and D.M. Blouw 2006. Inferring natural selection in a fossil threespine stickleback. Paleobiology 32:562-577.

Hunt, G., M. A. Bell and M. P. Travis (2008). Evolution toward a new adaptive optimum: phenotypic evolution in a fossil stickleback lineage. Evolution 62(3): 700-710.

Fit a model of trait evolution with a protracted punctuation.

Description

This function fits a model of punctuated change that is is protracted enough that it is captured by multiple transitional populations. Trait evolution starts in stasis, shifts to a general random walk, and then shifts back into stasis.

Usage

fit.sgs(
  y,
  minb = 7,
  oshare = TRUE,
  pool = TRUE,
  silent = FALSE,
  hess = FALSE,
  meth = "L-BFGS-B",
  model = "GRW"
)

Arguments

Value

a paleoTSfit object

See Also

fitGpunc

Examples

x <- sim.sgs(ns = c(10, 10, 10))  # default values OK
w <- fit.sgs(x, minb = 8)  # increase minb so example takes less time; not recommended!
plot(x)
abline(v = c(16, 31), lwd = 3)  # actual shifts
abline(v = c(w$parameters[6:7]), lwd = 2, lty = 3, col = "red")  # inferred shifts

Fit a set of standard evolutionary models

Description

Fit a set of standard evolutionary models

Usage

fit3models(y, silent = FALSE, method = c("Joint", "AD", "SSM"), ...)

fit4models(y, silent = FALSE, method = c("Joint", "AD", "SSM"), ...)

Arguments

Details

Function fit3models fits the general (biased) random walk (GRW), unbiased random walk (URW), and Stasis models. In addition to these three, fit4models also fits the model of Strict Stasis.

Value

if silent = FALSE, a table of model fit statistics, also printed to the screen. if silent = TRUE, a list of the model fit statistics and model parameter values.

Functions

• fit4models(): add model of "Strict Stasis" to the three models

See Also

fitSimple

Examples

x <- sim.GRW(ns = 50, ms = 0.2)
fit4models(x)

Fit large set of models to a time-series

Description

This function fits nine models to a time-series following Hunt et al. (2015). These include the simple models fit by fit4models along with mode shift and punctuation models.

Usage

fit9models(y, silent = FALSE, method = c("Joint", "AD"), ...)

Arguments

Value

if silent = FALSE, a table of model fit statistics, also printed to the screen. if silent = TRUE, a list of the model fit statistics and model parameter values.

References

Hunt, G., M. J. Hopkins and S. Lidgard. 2015. Simple versus complex models of trait evolution and stasis as a response to environmental change. PNAS 112(16): 4885-4890.

Examples

## Not run: 
x <- sim.Stasis.RW(ns = c(15, 15), omega = 0.5, ms = 1, order = "Stasis-RW")
plot(x)
fit9models(x)

## End(Not run)

Fit trait evolution model with punctuations estimated from the data

Description

Fit trait evolution model with punctuations estimated from the data

Usage

fitGpunc(
  y,
  ng = 2,
  minb = 7,
  pool = TRUE,
  oshare = TRUE,
  method = c("Joint", "AD"),
  silent = FALSE,
  hess = FALSE,
  parallel = FALSE,
  ...
)

Arguments

Details

This function tests all possible shift points for punctuations, subject to the constraint that the number of populations in each segment is always >= minb. The shiftpoint yielding the highest log-likelihood is returned as the solution, along with the log-likelihoods (all.logl) of all tested shift points (GG).

The function uses opt.punc (if method = "AD") or opt.joint.punc (if method = "Joint") to do the fitting.

Value

a paleoTSfit object with the results of the model-fitting.

Note

Calculations can be speeded up by setting parallel = TRUE, which uses package doParallel to run the bootstrap replicates in parallel, using one fewer than the number of detected cores.

See Also

fit9models, sim.punc

Examples

x <- sim.punc(ns = c(15, 15), theta = c(0,3), omega = c(0.1, 0.1))
w.punc <- fitGpunc(x, oshare = TRUE)
plot(x, modelFit = w.punc)

Fit model in which the mode of trait evolution shifts once

Description

Trait evolution is modeled as a shift from a random walk (general or unbiased) to stasis, or vice versa.

Usage

fitModeShift(
  y,
  minb = 7,
  pool = TRUE,
  order = c("Stasis-RW", "RW-Stasis"),
  rw.model = c("URW", "GRW"),
  method = c("Joint", "AD"),
  silent = FALSE,
  hess = FALSE,
  ...
)

Arguments

Value

a paleoTSfit object

See Also

sim.Stasis.RW

Examples

x <- sim.Stasis.RW(ns = c(15, 15), omega = 0.5, ms = 1, order = "Stasis-RW")
plot(x)
w <- fitModeShift(x, order = "Stasis-RW", rw.model = "GRW")
abline(v = x$tt[15], lwd = 3)  # actual shift point
abline(v = x$tt[w$par["shift1"]], lty = 3, lwd = 2, col = "red") # inferred shift

Fit the same simple model across multiple time-series

Description

Fit the same simple model across multiple time-series

Usage

fitMult(
  yl,
  model = c("GRW", "URW", "Stasis", "covTrack"),
  method = c("Joint", "AD"),
  pool = TRUE,
  zl = NULL,
  hess = FALSE
)

Arguments

Details

This function fits a model with shared parameters across multiple trait time-series. The most likely application would be to model a common evolutionary dynamic across different sequences, perhaps representing time-series of the same trait and lineage from different localities or time intervals.

Four simple models are currently implemented:

• GRW: parameters mstep and vstep of the general random walk are shared across sequences.

• URW: parameter vstep of the unbiased random walk is shared across sequences.

• Stasis: parameter omega of stasis is shared across sequences.

• covTrack: parameters b0, b1, and evar of the covariate-tracking model are shared across sequences.

Under the joint parameterization, method = "Joint", an additional parameter, anc is fit, representing the ancestral (starting) trait value. This parameter is estimated separately in each sequence so it is not assumed that they all start at the same trait value.

Value

a paleoTSfit object with the results of the model-fitting

Note

The models are described in the help for fitSimple and the functions linked from there.

See Also

fitSimple

Examples

x1 <- sim.GRW(ms = 1, vs = 0.2)
x2 <- sim.GRW(ms = 1, vs = 0.2)
fitMult(list(x1, x2), model = "GRW")

Fit simple models of trait evolution

Description

Fit simple models of trait evolution

Usage

fitSimple(
  y,
  model = c("GRW", "URW", "Stasis", "StrictStasis", "OU", "ACDC", "covTrack"),
  method = c("Joint", "AD", "SSM"),
  pool = TRUE,
  z = NULL,
  hess = FALSE
)

Arguments

Details

This is a convenience function that calls the specific individual functions for each model and parameterization, such as opt.GRW and opt.joint.GRW. The models that this function can fit are:

• GRW: General Random Walk. Under this model, evolutionary changes, or "steps" are drawn from a distribution with a mean of mstep and variance of vstep. mstep determines directionality and vstep determines volatility (Hunt, 2006).

• URW: Unbiased Random Walk. Same as GRW with mstep = 0, and thus evolution is non-directional. For a URW, vstep is the rate parameter.

• Stasis: This parameterization follows Sheets & Mitchell (2001), with a constant mean theta and variance omega (equivalent to white noise).

• Strict Stasis: Same as Stasis with omega = 0, indicating no real evolutionary differences; all observed variation is sampling error (Hunt et al. 2015).

• OU: Ornstein-Uhlenbeck model (Hunt et al. 2008). This model is that of a population ascending a nearby peak in the adaptive landscape. The optimal trait value is theta, alpha indicates the strength of attraction to that peak (= strength of stabilizing selection around theta), vstep measures the random walk component (from genetic drift) and anc is the trait value at the start of the sequence.

• ACDC: Accelerating or decelerating evolution model (Blomberg et al. 2003). This model is that of a population undergoing a random walk with a step variance that increases or decreases over time. The initial step variance is vstep0, and the parameter r controls its rate of increase (if positive) or decrease (if negative) over time. When r < 0, the is equivalent to the "Early burst" model of Harmon et al.

• covTrack: Covariate-tracking (Hunt et al. 2010). The trait tracks a covariate with slope b1, consistent with an adaptive response. evar is the residual variance, and, under method = "Joint", b0 is the intercept of the relationship between trait and covariate. model.

Value

a paleoTSfit object with the model fitting results

Note

For the covariate-tracking model, z should be a vector of length n when method = "Joint" and n - 1 when method = "AD", where n is the number of samples in y.

Method = "Joint" is a full likelihood approach, considering each time-series as a joint sample from a multivariate normal distribution. Method = "AD" is a REML approach that uses the differences between successive samples. They perform similarly, but the Joint approach does better under some circumstances (Hunt, 2008).

References

Hunt, G. 2006. Fitting and comparing models of phyletic evolution: random walks and beyond. Paleobiology 32(4): 578-601.
Hunt, G. 2008. Evolutionary patterns within fossil lineages: model-based assessment of modes, rates, punctuations and process. p. 117-131 In From Evolution to Geobiology: Research Questions Driving Paleontology at the Start of a New Century. Bambach, R. and P. Kelley (Eds).
Hunt, G., M. A. Bell and M. P. Travis. 2008. Evolution toward a new adaptive optimum: phenotypic evolution in a fossil stickleback lineage. Evolution 62(3): 700-710.
Sheets, H. D., and C. Mitchell. 2010. Why the null matters: statistical tests, random walks and evolution. Genetica 112– 113:105–125.
Blomberg, S. P., T. Garland, and A. R. Ives. 2003. Testing for phylogenetic signal in comparative data: behavioural traits are more labile. Evolution 57(4):717-745.
Harmon, L. J. et al. 2010. Early bursts of body size and shape evolution are rare in comparative data. Evolution 64(8):2385-2396.

See Also

opt.GRW, opt.joint.GRW, opt.joint.OU, opt.covTrack

Examples

y <- sim.Stasis(ns = 20, omega = 2)
w1 <- fitSimple(y, model = "GRW")
w2 <- fitSimple(y, model = "URW")
w3 <- fitSimple(y, model = "Stasis")
compareModels(w1, w2, w3)

Approximate log-transformation of time-series data

Description

Approximate log-transformation of time-series data

Usage

ln.paleoTS(y)

Arguments

Details

For a random variable x, its approximate mean on a natural log scale is the log of its untransformed mean. The approximate variance on a log scale is equal to the squared coefficient of variation.

Value

the converted paleoTS object

Note

This transformation only makes sense for variables with dimension and a true zero point, such as lengths and areas.

References

Hunt, G. 2006. Fitting and comparing models of phyletic evolution: random walks and beyond. Paleobiology 32:578-601.
Lewontin, R. 1966. On the measurement of relative variability. Systematic Zoology 15:141-142.

Examples

x <- sim.Stasis(ns = 10, theta = 20, omega = 1)
plot(x)
xl <- ln.paleoTS(x)
plot(xl)

Compute Lynch's Delta rate metric

Description

This function computes D, the rate metric proposed by Lynch (1990). This metric derives from the random walk model, with D = Vstep/(2Vp), where Vstep is the step variance of the unbiased random walk, and Vp is the within sample variance, pooled among samples. Under mutation - drift equilibrium, D is expected to range approximately between 5e-5 and 5e-3.

Usage

lynchD(y, gen.per.t = 1e+06, pool = TRUE, method = c("Joint", "AD"), ...)

Arguments

Value

References

Lynch (1990). The rate of morphological evolution in mammals from the standpoint of the neutral expectation. The American Naturalist 136:727-741. Hunt, G. 2012. Fitting and comparing models of phyletic evolution: random walks and beyond. Paleobiology 38:351-373.

Examples

y <- sim.GRW(ns = 20, ms = 0, vs = 1e-4, tt=seq(0, 1e6, length.out=20))  # per-year simulation
lynchD(y, gen.per.t = 1)

Analytical ML estimator for random walk and stasis models

Description

Analytical ML estimator for random walk and stasis models

Usage

mle.GRW(y)

mle.URW(y)

mle.Stasis(y)

Arguments

Value

a vector of mstep and vstep for mle.GRW, vstep for mle.URW, and theta and omega for mle.Stasis

Functions

• mle.URW(): ML parameter estimates for URW model

• mle.Stasis(): ML parameter estimates for Stasis model

Note

These analytical solutions assume even spacing of samples and equal sampling variance in each, which will usually be violated in real data. They are used here mostly to generate initial parameter estimates for numerical optimization; they not likely to be called directly by the user.

See Also

fitSimple

Examples

y <- sim.GRW(ms = 1, vs = 1)
w <- mle.GRW(y)
print(w)

Fit evolutionary model using "AD" parameterization

Description

Fit evolutionary model using "AD" parameterization

Usage

opt.GRW(
  y,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE
)

opt.URW(
  y,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE
)

opt.Stasis(
  y,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE
)

opt.StrictStasis(y, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

Arguments

Details

These functions use differences between consecutive populations in the time series in order to remove temporal autocorrelation. This is referred to as the "Ancestor-Descendant" or "AD" parameterization by Hunt [2008], and it is a REML approach (like phylogenetic independent contrasts). A full ML approach, called "Joint" was found to have generally better performance (Hunt, 2008) and generally should be used instead.

Value

a paleoTSfit object with the model fitting results

Functions

• opt.URW(): fit the URW model by the AD parameterization

• opt.Stasis(): fit the Stasis model by the AD parameterization

• opt.StrictStasis(): fit the Strict Stasis model by the AD parameterization

Note

It is easier to use the convenience function fitSimple.

References

Hunt, G. 2006. Fitting and comparing models of phyletic evolution: random walks and beyond. Paleobiology 32(4): 578-601.

See Also

fitSimple, opt.joint.GRW

Examples

x <- sim.GRW(ns = 20, ms = 1)  # strong trend
plot(x)
w.grw <- opt.GRW(x)
w.urw <- opt.URW(x)
compareModels(w.grw, w.urw)

Fit random walk model with shift(s) in generating parameters

Description

Fit random walk model with shift(s) in generating parameters

Usage

opt.GRW.shift(y, ng = 2, minb = 7, model = 1, pool = TRUE, silent = FALSE)

Arguments

Details

Fits a model in which a sequence is divided into two or more segments and trait evolution proceeds as a general random walk, with each segment (potentially) getting its own generating parameters (mstep, vstep).

This function tests for shifts after each population, subject to the constraint that the number of populations in each segment is always >= minb. The shiftpoint yielding the highest log-likelihood is returned as the solution, along with the log-likelihoods (all.logl of all tested shift points (GG).

Different variants of the model can be specified by the model argument:

• model = 1: mstep is separate across segments; vstep is shared

• model = 2: mstep is shared across segments; vstep is separate

• model = 3: mstep is set to zero (unbiased random walk); vstep is separate across segments

• model = 4: mstep and vstep are both separate across segments

Value

a paleoTSfit object

See Also

sim.GRW.shift

Examples

x <- sim.GRW.shift(ns = c(15,15), ms = c(0, 1), vs = c(0.1,0.1))
w.sep <- opt.GRW.shift(x, ng = 2, model = 4)
w.sameVs <- opt.GRW.shift(x, ng = 2, model = 1)
compareModels(w.sep, w.sameVs)
plot(x)
abline(v = x$tt[16], lwd = 3)  # actual shift point
abline(v = x$tt[w.sameVs$par["shift1"]], lty = 3, col = "red", lwd = 2) # inferred shift point

Fit a model in which a trait tracks a covariate

Description

Fit a model in which a trait tracks a covariate

Usage

opt.covTrack(
  y,
  z,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE
)

opt.joint.covTrack(
  y,
  z,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE
)

Arguments

Details

In this model, changes in a trait are linearly related to changes in a covariate with a slope of b and residual variance evar: dx = b * dz + eps, where eps ~ N(0, evar). This model was described, and applied to an example in which body size changes tracked changes in temperature, by Hunt et al. (2010).

For the AD version (opt.covTrack), a trait sequence of length ns, the covariate, z, can be of length ns - 1, interpreted as the vector of changes, dx. If z is of length ns, differences are taken and these are used as the dx's, with a warning issued.

The Joint version (opt.joint.covTrack), z should be of length ns and there is an additional parameter that is the intercept of the linear relationship between trait and covariate. See warning below about using the Joint version.

Value

a paleoTSfit object with the results of the model fitting

Functions

• opt.joint.covTrack(): fits the covTrack model using the joint parameterization

Warning

The Joint parameterization of this model can be fooled by temporal autocorrelation and, especially, trends in the trait and the covariate. The latter is tested for, but the AD parameterization is generally safer for this model.

References

Hunt, G, S. Wicaksono, J. E. Brown, and K. G. Macleod. 2010. Climate-driven body size trends in the ostracod fauna of the deep Indian Ocean. Palaeontology 53(6): 1255-1268.

See Also

fitSimple

Examples

set.seed(13)
z <- c(1, 2, 2, 4, 0, 8, 2, 3, 1, 9, 4, 3)
x <- sim.covTrack(ns = 12, z = z, b = 0.5, evar = 0.2)
w.urw <- opt.URW(x)
w.cov <- opt.covTrack(x, z)
compareModels(w.urw, w.cov)

Fit evolutionary models using the "Joint" parameterization

Description

Fit evolutionary models using the "Joint" parameterization

Usage

opt.joint.GRW(
  y,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE
)

opt.joint.URW(
  y,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE
)

opt.joint.Stasis(
  y,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE
)

opt.joint.StrictStasis(y, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

Arguments

Details

These functions use the joint distribution of population means to fit models using a full maximum-likelihood approach. This approach was found to have somewhat better performance than the "AD" approach, especially for noisy trends (Hunt, 2008).

Value

a paleoTSfit object with the model fitting results

Functions

• opt.joint.URW(): fit the URW model by the Joint parameterization

• opt.joint.Stasis(): fit the Stasis model by the Joint parameterization

• opt.joint.StrictStasis(): fit the Strict Stasis model by the Joint parameterization

Note

It is easier to use the convenience function fitSimple.

References

Hunt, G., M. J. Hopkins and S. Lidgard. 2015. Simple versus complex models of trait evolution and stasis as a response to environmental change. Proc. Natl. Acad. Sci. USA 112(16): 4885-4890.

See Also

fitSimple, opt.GRW

Examples

x <- sim.GRW(ns = 20, ms = 1)  # strong trend
plot(x)
w.grw <- opt.joint.GRW(x)
w.urw <- opt.joint.URW(x)
compareModels(w.grw, w.urw)

Fit Ornstein-Uhlenbeck model using the "Joint" parameterization

Description

Fit Ornstein-Uhlenbeck model using the "Joint" parameterization

Usage

opt.joint.OU(
  y,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE
)

Arguments

Details

This function fits an Ornstein-Uhlenbeck (OU) model to time-series data. The OU model has four generating parameters: an ancestral trait value (anc), an optimum value (theta), the strength of attraction to the optimum (alpha), and a parameter that reflects the tendency of traits to diffuse (vstep). In a microevolutionary context, these parameters can be related to natural selection and genetic drift; see Hunt et al. (2008).

Value

a paleoTSfit object with the model fitting results

Note

It is easier to use the convenience function fitSimple. Note also that preliminary work found that the "AD" parameterization did not perform as well for the OU model and thus it is not implemented here.

References

Hunt, G., M. A. Bell and M. P. Travis. 2008. Evolution toward a new adaptive optimum: phenotypic evolution in a fossil stickleback lineage. Evolution 62(3): 700-710.

See Also

fitSimple, opt.joint.GRW

Examples

x <- sim.OU(vs = 0.5)  # most defaults OK
w <- opt.joint.OU(x)
plot(x, modelFit = w)

Fit a model of trait evolution with specified punctuation(s)

Description

Fit a model of trait evolution with specified punctuation(s)

Usage

opt.punc(
  y,
  gg,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE,
  oshare
)

opt.joint.punc(
  y,
  gg,
  pool = TRUE,
  cl = list(fnscale = -1),
  meth = "L-BFGS-B",
  hess = FALSE,
  oshare
)

Arguments

Details

The sequence is divided into segments, which are separated by punctuations. Means for each segment are given by the vector theta with variances given by the vector omega (or a single value if oshare = TRUE). This function calls optim to numerically fit this model to a time-series, y.

Value

a paleoTSfit object with the results of the model fitting

Functions

• opt.joint.punc(): fits the punctuation model using the joint parameterization

Note

These functions would be used in the uncommon situation in which there is a prior hypothesis as to where the punctuation(s) take place. Normally users will instead use the function fitGpunc, which uses these functions to fit a range of possible timings for the punctuations.

See Also

fitGpunc

Examples

x <- sim.punc(ns = c(15, 15), theta = c(0,3), omega = c(0.1, 0.1))
w.sta  <- fitSimple(x, model = "Stasis", method = "Joint")
w.punc <- opt.joint.punc(x, gg = rep(1:2, each = 15), oshare = TRUE)
compareModels(w.sta, w.punc)

Fit evolutionary models using state-space models (SSM)

Description

Fit evolutionary models using state-space models (SSM)

Usage

opt.ssm.GRW(y, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

opt.ssm.URW(y, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

opt.ssm.Stasis(y, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

opt.ssm.StrictStasis(y, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

opt.ssm.OU(y, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

opt.ssm.ACDC(y, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

opt.ssm.covOU(y, z, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

opt.ssm.URWshift(y, gg, pool = TRUE, cl = list(fnscale = -1), hess = FALSE)

opt.ssm.covOU_vshift(
  y,
  z,
  gg,
  pool = TRUE,
  cl = list(fnscale = -1),
  hess = FALSE
)

Arguments

Details

These functions use a state space model formulation to compute likelihoods and fit models.

Functions to fit the OU covariate tracking models (covOU, covOU_vshift) require a covariate argument, z. At present, only the OU covariate tracking with a shift in the step variance (covOU_vshift) requires the grouping vector argument (gg).

Value

a paleoTSfit object with the model fitting results

Note

For GRW, URW, Stasis, StrictStasis, ACDC and OU models, it will likely be easier to use the convenience function fitSimple with argument method = "SSM".
The grouping vector, gg, is a vector of length equal to the number of samples. It has one element for each sample and takes integer value from 1 to the number of sample groups separated by shiftpoints. See the example below.

See Also

fitSimple, opt.joint.GRW

Examples

y <- sim.GRW(ns = 30, vs = 2)
w1 <- opt.ssm.URW(y)
gg <- rep(1:2, each = 15) # shift occurs immediately after sample 15
w2 <- opt.ssm.URWshift(y, gg = gg) # test model in which the step variance shifts
compareModels(w1, w2)

Plot a paleoTS object

Description

Plot a paleoTS object

Usage

## S3 method for class 'paleoTS'
plot(
  x,
  nse = 1,
  pool = FALSE,
  add = FALSE,
  modelFit = NULL,
  pch = 21,
  pt.bg = "white",
  lwd = 1.5,
  ylim = NULL,
  xlab = NULL,
  ylab = NULL,
  add.label = TRUE,
  ...
)

Arguments

Value

none.

Examples

x <- sim.GRW(ns = 30)
w <- fitSimple(x, model = "GRW", method = "Joint")
plot(x, modelFit = w)
plot(x, xlab = "Time [Myr]", ylab = "Body length [mm]", pch = 22, pt.bg = "blue")

Compute a pooled variance

Description

Computes a pooled variance from samples in a paleontological time-series

Usage

pool.var(y, nn = NULL, minN = NULL, ret.paleoTS = FALSE)

Arguments

Details

A pooled variance of a set of populations is the weighted average of the individual variances of the populations, with the weight for each population equal to its sample size minus one.

For many kinds of traits, variation levels tend to be similar among closely related populations. When this is true and sample sizes are low, much of the observed differences in variance among samples will be due to the high noise of estimated the variances. Replacing the observed variances of all populations (or only those with nn < minN) with the estimated pooled variance can reduce this noise.

Value

if ret.paleoTS = TRUE a paleoTS object with all (or some) variances replaced with the pooled variance; otherwise the pooled variance

Examples

data(cantius_L)
cant_all <- pool.var(cantius_L, ret.paleoTS = TRUE)   # replace all variances with pooled variance
cant_n5  <- pool.var(cantius_L, minN = 5, ret.paleoTS = TRUE)  # replace only pops with n < 5

Print a paleoTSfit object

Description

Print a paleoTSfit object

Usage

## S3 method for class 'paleoTSfit'
print(x, ...)

Arguments

Value

None; this function is used only to print

Examples

y <- sim.punc(theta = c(0, 2), omega = c(0.1, 0.1))
wg <- fitGpunc(y)
print(wg)

Read a text-file with data from a paleontological time-series

Description

Read a text-file with data from a paleontological time-series

Usage

read.paleoTS(file = NULL, oldest = "first", reset.time = TRUE, ...)

Arguments

Details

This function reads a text file with a specified format and converts it into a paleoTS object. It will often be the easiest way for users to import their own data. The text file should have four columns without headers, in this order: sample size, sample means, sample variances, sample ages.

Value

a paleoTS object

See Also

as.paleoTS

Simulate random walk or directional time-series for trait evolution

Description

Simulate random walk or directional time-series for trait evolution

Usage

sim.GRW(ns = 20, ms = 0, vs = 0.1, vp = 1, nn = rep(20, ns), tt = 0:(ns - 1))

Arguments

Details

The general random walk model considers time in discrete steps. At each time step, an evolutionary change is drawn at random from a distribution of possible evolutionary "steps." It turns out that the long-term dynamics of an evolving lineage depend only on the mean and variance of this step distribution. The former, mstep, determined the directionality in a sequence and the latter, vstep, determines its volatility.

Value

a paleoTS object

Note

This function simulates an unbiased random walk if ms is equal to zero and a general (or biased) random walk otherwise.

See Also

sim.Stasis, sim.OU, as.paleoTS

Examples

x.grw <- sim.GRW(ms = 0.5)
x.urw <- sim.GRW(ms = 0)
plot(x.grw, ylim = range(c(x.grw$mm, x.urw$mm)))
plot(x.urw, add = TRUE, col = "blue")
legend(x = "topleft", c("GRW", "URW"), col = c("black", "blue"), lty = 1)

Simulate (general) random walk with shift(s) in generating parameters

Description

Simulate (general) random walk with shift(s) in generating parameters

Usage

sim.GRW.shift(
  ns = c(10, 10),
  ms = c(0, 1),
  vs = c(0.5, 0.5),
  nn = rep(30, sum(ns)),
  tt = 0:(sum(ns) - 1),
  vp = 1
)

Arguments

Details

Simulates under a model in which a sequence is divided into two or more segments. Trait evolution proceeds as a general random walk, with each segment getting its own generating parameters (mstep, vstep).

Value

a paleoTS object with the simulated time-series

See Also

sim.GRW, sim.sgs, opt.GRW.shift

Examples

x <- sim.GRW.shift(ns = c(10,10,10), ms = c(0, 1, 0), vs = c(0.1,0.1,0.1))
plot(x)
abline(v = c(9.5, 19.5), lty = 3, lwd = 2, col = "blue")  # shows where dynamics shift
text (c(5, 15, 25), c(2,2,2), paste("segement", 1:3, sep =" "), col = "blue")

Simulate an Ornstein-Uhlenbeck time-series

Description

Simulate an Ornstein-Uhlenbeck time-series

Usage

sim.OU(
  ns = 20,
  anc = 0,
  theta = 10,
  alpha = 0.3,
  vstep = 0.1,
  vp = 1,
  nn = rep(20, ns),
  tt = 0:(ns - 1)
)

Arguments

Details

This function simulates an Ornstein-Uhlenbeck (OU) process. In microevolutionary terms, this models a population ascending a nearby peak in the adaptive landscape. The optimal trait value is theta, alpha indicates the strength of attraction to that peak (= strength of stabilizing selection around theta), vstep measures the random walk component (from genetic drift) and anc is the trait value at the start of the sequence.

Value

a paleoTS object

References

Hunt, G., M. A. Bell and M. P. Travis. 2008. Evolution toward a new adaptive optimum: phenotypic evolution in a fossil stickleback lineage. Evolution 62(3): 700-710.

See Also

opt.joint.OU

Examples

x1 <- sim.OU(alpha = 0.8)  # strong alpha
x2 <- sim.OU(alpha = 0.1)  # weaker alpha
plot(x1)
plot(x2, add = TRUE, col = "blue")

Simulate Stasis time-series for trait evolution

Description

Simulate Stasis time-series for trait evolution

Usage

sim.Stasis(
  ns = 20,
  theta = 0,
  omega = 0,
  vp = 1,
  nn = rep(20, ns),
  tt = 0:(ns - 1)
)

Arguments

Value

a paleoTS object

See Also

sim.GRW, sim.OU, as.paleoTS

Examples

x <- sim.Stasis(omega = 0.5, vp = 0.1)
w.sta <- fitSimple(x, model = "Stasis")
w.ss  <- fitSimple(x, model = "StrictStasis")
compareModels(w.sta, w.ss)

Simulate trait evolution with a mode shift

Description

Trait evolution is modeled as a shift from a random walk (general or unbiased) to stasis, or vice versa.

Usage

sim.Stasis.RW(
  ns = c(20, 20),
  order = c("Stasis-RW", "RW-Stasis"),
  anc = 0,
  omega = 1,
  ms = 0,
  vs = 1,
  vp = 1,
  nn = 30,
  tt = NULL
)

Arguments

Details

The anc argument is the starting trait value, and if the first segment is stasis, this is also the value of the stasis mean. When the first segment is a random walk, the stasis mean in the second segment is equal to the true trait mean at the end of the initial random walk.

Value

a paleoTSfit object

See Also

fitModeShift

Examples

x1 <- sim.Stasis.RW(omega = 0.1, ms = 5, order = "Stasis-RW")
x2 <- sim.Stasis.RW(omega = 0.1, ms = 5, order = "RW-Stasis")
plot(x1)
plot(x2, add = TRUE, col = "blue")
abline(v = 19, lty=3)

Simulate trait evolution that tracks a covariate

Description

Simulate trait evolution that tracks a covariate

Usage

sim.covTrack(
  ns = 20,
  b = 1,
  evar = 0.1,
  z,
  nn = rep(20, times = ns),
  tt = 0:(ns - 1),
  vp = 1
)

Arguments

Details

In this model, changes in a trait are linearly related to changes in a covariate with a slope of b and residual variance evar: dx = b * dz + eps, where eps ~ N(0, evar). This model was described, and applied to an example in which body size changes tracked changes in temperature, by Hunt et al. (2010).

Value

a paleoTS object

Note

For a trait sequence of length ns, the covariate, z, can be of length ns - 1,in which case it is interpreted as the vector of changes, dz. If z is of length ns, differences are taken and these are used as the dz's.

References

Hunt, G, S. Wicaksono, J. E. Brown, and K. G. Macleod. 2010. Climate-driven body size trends in the ostracod fauna of the deep Indian Ocean. Palaeontology 53(6): 1255-1268.

Examples

set.seed(13)
z <- c(1, 2, 2, 4, 0, 8, 2, 3, 1, 9, 4, 3)
x <- sim.covTrack(ns = 12, z = z, b = 0.5, evar = 0.2)
plot(x, ylim = c(-1, 10))
lines(x$tt, z, col = "blue")

Simulate a punctuated time-series

Description

Simulates punctuated trait evolution with punctuations that are rapid relative to the spacing of samples. In practice, the time-series is divided into two or more segments, each of which has its own mean and variance.

Usage

sim.punc(
  ns = c(10, 10),
  theta = c(0, 1),
  omega = rep(0, length(theta)),
  nn = rep(30, sum(ns)),
  tt = 0:(sum(ns) - 1),
  vp = 1
)

Arguments

Details

Segments are separated by punctuations. Population means in the ith segment are drawn randomly from a normal distribution with a mean equal to ith element of theta and variance equal to the ith element of omega. The magnitudes of punctuations are determined by the differences in adjacent theta values.

Value

a paleoTS object with the simulated time-series.

See Also

fitGpunc

Examples

x <- sim.punc(ns = c(15, 15), theta = c(0,3), omega = c(0.1, 0.1))
plot(x)

Simulate protracted punctuation

Description

This function simulates a punctuated change that is is protracted enough that it is captured by multiple transitional populations. Trait evolution starts in stasis, shifts to a general random walk, and then shifts back into stasis.

Usage

sim.sgs(
  ns = c(20, 20, 20),
  theta = 0,
  omega = 1,
  ms = 1,
  vs = 0.1,
  nn = rep(30, sum(ns)),
  tt = 0:(sum(ns) - 1),
  vp = 1
)

Arguments

Details

Trait evolution proceeds in three segments: Stasis, General random walk, stasis (sgs). The initial stasis segment has a mean of theta and variance omega before shifting in the second segment to a general random walk with parameters ms and vs. Finally, the third segment is a return to stasis, centered around the trait value of the last population of the random walk.

Value

a paleoTS object

Examples

x <- sim.sgs()  # default values OK
plot(x)

Convert time-series to standard deviation units

Description

Convert time-series to standard deviation units

Usage

std.paleoTS(y, center = c("mean", "start"))

Arguments

Details

The standardization expresses each sample mean as the deviation from the overall mean, divided by the pooled within-sample standard deviation. Sample variances are also divided by the pooled sample variance.

Essentially, this converts paleontological time-series data into standard deviation units, similar to the computation of evolutionary rates in haldanes. This operation does not change the relative fit of models, but it does facilitate the comparison of parameter estimates across time-series of traits measured in different units.

If argument center = "start" the time-series is translated such that the trait mean of the first sample is zero.

Value

the converted paleoTS object

Examples

x <- sim.Stasis(ns = 8, theta = 1, omega = 4, vp = 2)
xs <- std.paleoTS(x, center = "start")
plot(x, ylim = range(c(x$mm, xs$mm)))
plot(xs, col = "red", add = TRUE)
legend(x = "topright", c("unstandardized", "standardized"), lty=1, col=c("black", "red"), cex=0.7)

Subsample a paleontological time-series

Description

Subsampling is done according to the supplied logical vector or, if none is supplied, as a proportion of samples, randomly chosen.

Usage

sub.paleoTS(y, ok = NULL, k = 0.1, reset.time = TRUE)

Arguments

Value

the sub-sampled paleoTS object

Examples

x <- sim.GRW(ns=20)
plot(x)
xs1 <- sub.paleoTS(x, k = 0.5)
plot(xs1, add = TRUE, col="green")
keep <- rep(c(TRUE, FALSE), 10)
xs2 <- sub.paleoTS(x, ok = keep)
plot(xs2, add = TRUE, col = "red")

Test for heterogeneity of variances among samples in a time-series

Description

Test for heterogeneity of variances among samples in a time-series

Usage

test.var.het(y, method = "Bartlett")

Arguments

Value

a list with the test statistic, degrees of freedom, and p-value

Note

Most often, this function will be used to assess if it is reasonable to pool variances across samples using pool.var. A significant result means that the null hypothesis of equal variances across samples is rejected. Even in this case, however, it may still be preferable to pool variances, at least for some populations, if sample sizes are quite low.

References

Sokal, R. R and F. J. Rohlf. 1995. Biometry 3rd Ed.

See Also

pool.var

Examples

data(cantius_L)
test.var.het(cantius_L)  # significant, but still may want to pool variances