Fit two models to two separate segments to an evolutionary sequence (time-series).

Wrapper function to find maximum likelihood solutions to two models to an evolutionary sequence.

Usage

fit.mode.shift(
  y,
  model1 = c("Stasis", "URW", "GRW", "OU"),
  model2 = c("Stasis", "URW", "GRW", "OU"),
  fit.all = FALSE,
  minb = 7,
  shift.point = NULL,
  pool = TRUE,
  silent = FALSE,
  hess = FALSE
)

Arguments

y

an univariate evoTS object.

model1

the model fitted to the first segment. Options are Stasis, URW, GRW, OU.

model2

the model fitted to the second segment. Options are Stasis, URW, GRW, OU.

fit.all

logical indicating whether to fit all pairwise combinations of the four models to the evolutionary sequence (time-series).

minb

the minimum number of samples within a segment to consider

shift.point

The sample that split the time-series into two segments. The samples are passed to the argument as a vector. Default is NULL, which means all possible shift points will be assessed constrained by how minb is defined.

pool

logical indicating whether to pool variances across samples

silent

if TRUE, less information is printed to the screen as the model is fit

hess

logical, indicating whether to calculate standard errors from the Hessian matrix.

Value

the function returns a list of all investigated models and their highest log-likelihood (and their corresponding AICc and AICc weight).

logL: the log-likelihood of the optimal solution
AICc: AIC with a correction for small sample sizes
parameters: parameter estimates
modelName: abbreviated model name
method: Joint consideration of all samples
K: number of parameters in the model
n: the number of observations/samples
all.logl: log-likelihoods for all tested partitions of the series into segments. Will return a single value if shift points have been given
GG: matrix of indices of initial samples of each tested segment configuration; each column of GG corresponds to the elements of all.logl

In addition, if fit.all=TRUE the function also returns a list of all investigated models and their highest log-likelihood (and their corresponding AICc and AICc weight).

Note

The models have been implemented to be compatible with the joint parameterization routine in the package paleoTS. The optimization is therefore fit using the actual sample values, with the autocorrelation among samples accounted for in the log-likelihood function. The joint distribution of sample means is multivariate normal, with means and variance-covariances determined by evolutionary parameters and sampling errors.

References

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

Hunt, G., Bell, M. A. & Travis, M. P. Evolution towards a new adaptive optimum: Phenotypic evolution in a fossil stickleback lineage. Evolution 62:700–710 (2008)

Author

Kjetil Lysne Voje

Examples


##Generate a paleoTS object.
x <- paleoTS::sim.GRW(30)

## Fit a mode-shift model without defining a shift point (the example may take > 5 seconds to run)
fit.mode.shift(x, model1="URW", model2="Stasis")
#> [1] "Searching all possible shift points in the evolutionary sequence"
#> Total # hypotheses:  17 
#> 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  
#> $logL
#> [1] -9.276821
#> 
#> $AICc
#> [1] 31.05364
#> 
#> $parameters
#>         anc       vstep       theta       omega      shift1 
#> -0.20433411  0.01341561 -0.25746281  0.15904385 24.00000000 
#> 
#> $modelName
#> [1] "URW-Stasis"
#> 
#> $method
#> [1] "Joint"
#> 
#> $se
#> NULL
#> 
#> $K
#> [1] 5
#> 
#> $n
#> [1] 30
#> 
#> $all.logl
#>  [1] -17.411176 -18.469248 -18.295720 -18.121415 -17.651367 -17.506559
#>  [7] -16.972842 -16.409491 -15.630322 -14.733457 -13.687043 -12.353073
#> [13] -12.440003 -11.017427 -10.941428  -9.353559  -9.276821
#> 
#> $GG
#> shift1 
#>     24 
#> 
#> attr(,"class")
#> [1] "paleoTSfit"