# simulation.R
rm(list=ls())
require(wrightscape)
require(snowfall)
require(ggplot2)
# store the unique id of this script version
require(socialR)
gitaddr <- gitcommit("simulation.R")
id <- gitlog()$shortID
print(id)
data(labrids)
# rename the regimes less technically
levels(pharyngeal) = c("wrasses", "parrotfish")
regime <- pharyngeal
# Create a dataset by simulation, where the parrotfish all have lower alpha
a1_spec <- list(alpha = "indep", sigma = "global", theta = "global")
a1 <- multiTypeOU(data = dat["close"], tree = tree, regimes = regime,
model_spec = a1_spec, control = list(maxit=5000))
# expected "startree" standard deviations: wrasses 1.25, parrotfish 5
a1$alpha[1] <- 10
a1$alpha[2] <- 1e-10
a1$sigma <- c(5, 5)
a1$theta <- c(0,0)
dat[["constraint release"]] <-simulate(a1)[[1]]
# Check out the variance in the relative groups -- it should be larger in parrotfish
# in an extreme example, but need not be...
testcase <- dat[["constraint release"]]
testcase[regime =="wrasses" & !is.na(testcase) & testcase != 0] -> lowvar
testcase[regime !="wrasses" & !is.na(testcase) & testcase != 0] -> highvar
print(c(var(lowvar), var(highvar)))
## We can repeat the whole thing with a model based on differnt sigmas, to make sure
s1_spec <- list(alpha = "global", sigma = "indep", theta = "global")
s1 <- multiTypeOU(data = dat["close"], tree = tree, regimes = pharyngeal,
model_spec = s1_spec, control = list(maxit=5000))
# Order of entries in sigma is the order regime names are given by levels (alphabetical)
names(s1$sigma) <- levels(regime)
s1$sigma[1] <- sqrt(2*5*1.25)
s1$sigma[2] <- sqrt(2*5*5)
s1$alpha <- c(5, 5) # We can keep those parameters estimated from data or update them
a1$theta <- c(0,0)
dat[["faster evolution"]] <-simulate(s1)[[1]]
testcase <- dat[["faster evolution"]]
testcase[regime == "wrasses" & !is.na(testcase) & testcase != 0 ] -> lowvar
testcase[regime != "wrasses" & !is.na(testcase) & testcase != 0 ] -> highvar
print(c(var(lowvar), var(highvar)))
## Now we have a trait where change in alpha is responsible,
## and one in which sigma change is responsible.
## Can we correctly identify each??
traits <- c("constraint release", "faster evolution")
sfInit(par=T, 10) # for debugging locally
sfLibrary(wrightscape)
sfLibrary(pmc)
fits <- lapply(traits, function(trait){
multi <- function(modelspec)
m <- multiTypeOU(data = dat[[trait]], tree = tree, regimes = pharyngeal,
model_spec = modelspec,
control = list(maxit=5000)
)
s1 <- multi(list(alpha = "global", sigma = "indep", theta = "global"))
a1 <- multi(list(alpha = "indep", sigma = "global", theta = "global"))
sfExportAll()
mc <- montecarlotest(s1, a1, nboot=100)
})
save(list=ls(), file=sprintf("sim_modelchoice_%s.Rdat", id))
print(id)
#upload(sprintf("%s_*.png", id), gitaddr=gitaddr, tag="phylogenetics", save=FALSE)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.