## 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]]
Run the code above in your browser using DataLab