simulate_div {diverge}

Simulate trait divergence

Description

Simulate replicate sets of lineage pair trait divergence under one of ten evolutionary models.

Usage

simulate_div(model, pars, ages, me_prop=NULL, GRAD = NULL, cats=NULL, 
  breakpoint = NULL, Nsets = 1)

Arguments

model

Character string defining one of 11 models of trait divergence (options: "BM_null", "BM_linear", "OU_null", "OU_linear", "DA_null", "DA_linear", "DA_cat", "DA_wt", "DA_bp", "DA_wt_linear", "DA_bp_linear", "DA_OU", "DA_OU_linear", "DA_OU_cat", "DA_BM", "DA_BM_linear", "DA_BM_cat", "OU_BM", "OU_BM_linear", "OU_BM_cat"). See model_select for model descriptions.

pars

Vector of model parameters to use in sumulation. Parameters must be in the correct order in the vector for each model. The parameter orders are
BM_null: sig2,
BM_linear: sig2_slope, sig2_int
OU_null: alpha, sig2
OU_linear: alpha_int, alpha_slope, sig2
DA_null: alpha, sig2, psi
DA_linear: alpha, sig2, psi_slope, psi_int
DA_cat: alpha, sig2, psi1, psi2, psi3(optional)
DA_wt: alpha, sig2, psi1, psi2, wt
DA_bp: alpha, sig2, psi1, psi2
DA_wt_linear: alpha, sig2, psi1_slope, psi1_int, psi2_slope, psi2_int, wt
DA_bp_linear: alpha, sig2, psi1_slope, psi1_int, psi2_slope, psi2_int

ages

Vector of 'ages' (i.e. times since divergence) of lineage pairs for which to simulate trait divergence.

me_prop

Induced measurement error as a proportion of the expected variance. Can be a single value applied to all simulated species or a vector of values of length equal to length(ages).

GRAD

Vector containing the gradient position of each pair to be simulated. This is the value of a continuous variable such as latitude or body size across which parameters will vary. Required for all models with 'linear' suffix.

cats

Vector containing the category code (0, 1, or 2) for each pair (see model_select help page for details).

breakpoint

Vector of breakpoint times (times after 0 at which a discrete shift occurs in psi) for each lineage pair to be simulated. Required for all "bp" models. Must align with "ages" vector such that ages[i] and breakpoint[i] correspond to the same lineage pair. When using bp models, datasets can and should be simulated with a mix of pairs that have and have not experienced a shift in psi (see model_select for details on bp models). Set bp=0 for pairs that experience no shift.

Nsets

Number of datasets to simulate.

Details

Simulates replicate sets of trait divegence under one of ten evolutionary models given a user-defined set of lineage pair ages. Since all evolutionary models are stochastic process models, replicate datasets simulated under the same model parameters are not identical. Users can simulate the addition of measurement error using the ms_prop argument, which denotes measurement error as a proportion of the expected variance (e.g. ms_prop = 0.1 corresponds to a measurement error in each species equal to 10 percent of the expected variance in the trait in that species, given the parameters).

Value

Returns a vector (if Nsets=1) representing divergences for one dataset or a list of vectors (if Nsets>1) for replicate datasets of lineage-pair trait divergences.

Author(s)

Sean A. S. Anderson

Examples

## Simulate datasets of trait divergence for 150 lineage pairs 
# under different evolutionary models

# Define vector of 150 lineage pair ages
ages = rep(c(0.5, 1, 1.5, 2, 3, 8), 25)

## simulate a dataset of trait divergence under a single-process Brownian Motion model
sig2 = 0.2
sis_div_BM = simulate_div(model="BM_null", ages=ages, pars=c(sig2))
summary(sis_div_BM)
sis_div_BM

## simulate 1000 divergence datasets under the null model of divergent adaptation
sig2 = 0.2
alpha = 0.8
psi = 0.8
sis_div_DA = simulate_div(model="DA_null", ages=ages, pars=c(alpha, sig2, psi), N=1000)
length(sis_div_DA)
sis_div_DA[[1]]

## Simulate 1000 divergence datasets under the breakpoint model
# note: a breakpoint vector must be created
# pairs that experience no epoch shift (no shift in the psi parameter) are assigned bp=0
# for all two-epoch pairs, the breakpoint time must be lower than the age of the pair
# here we make half the dataset into two-epoch pairs, half into one-epoch pairs
# we set the breakpoint time equal to half the age of each two-epoch pair
bp = c(ages[1:75]/2, rep(0, 75))

# run the simulation
sig2 = 0.2
alpha = 0.8
psi1 = 0.3
psi2 = 0.9
sis_div_bp = simulate_div(model="DA_bp", ages=ages, breakpoint=bp, 
  pars=c(alpha, sig2, psi1, psi2), N=1000)
length(sis_div_bp)
sis_div_bp[[1]]