require(wrightscape)
require(snowfall)
require(ggplot2)
require(reshape2)
require(Hmisc)
data(labrids)
traits <- c("bodymass", "close", "open", "kt", "gape.y", "prot.y", "AM.y", "SH.y", "LP.y")
regimes <- two_shifts
Setup parallel environment
sfInit(par=T, 4)
R Version: R version 2.14.1 (2011-12-22)
sfLibrary(wrightscape)
Library wrightscape loaded.
sfExportAll()
Bootstrap the fits by replication (bootstraps parameter values without bootstrapping the likelihood ratios)
fits <- sfLapply(traits, function(trait){
multi <- function(modelspec, reps = 100){
m <- multiTypeOU(data = dat[trait], tree = tree, regimes = regimes,
model_spec = modelspec,
control = list(maxit=8000)
)
replicate(reps, bootstrap(m))
}
bm <- multi(list(alpha = "fixed", sigma = "global", theta = "global"))
ou <- multi(list(alpha = "global", sigma = "global", theta = "global"))
bm2 <- multi(list(alpha = "fixed", sigma = "indep", theta = "global"))
a2 <- multi(list(alpha = "indep", sigma = "global", theta = "global"))
t2 <- multi(list(alpha = "global", sigma = "global", theta = "indep"))
list(bm=bm,brownie=bm2, ou=ou, ouch=t2, alphas=a2)
})
Reformat and label data for plotting
names(fits) <- traits # each fit is a different trait (so use it for a label)
data <- melt(fits)
names(data) <- c("regimes", "param", "rep", "value", "model", "trait")
subdat <- subset(data, param %in% c("alpha")
& model %in% c("alphas")
& value < 20)
r <- cast(subdat, regimes ~ model ~ trait ~ param, smedian.hilow, conf.int=.5, na.rm=T)
upper <- sapply(c("alpha"), function(t) max(r[, , , t]))
Generate a plot of the relevant parameters
p0 <-
ggplot(subdat, aes(model, value, fill=regimes)) +
stat_summary(fun.data=median_hilow, geom="pointrange", aes(color=regimes),
position = position_dodge(width=0.90), conf.int=.5) +
facet_grid(param ~ trait, scales = "free_y") +
coord_cartesian(ylim=c(0,upper["alpha"]), wise=TRUE) +
opts(title = "alpha")
p0
save(list=ls, file = "~/public_html/data/labrids_par_bootstrap.rda")
Error: invalid first argument
Additional plots, just for variety
range <- cast(data, regimes ~ model ~ trait ~ param, smedian.hilow,
conf.int=.5, na.rm=T)
upper <- sapply(c("sigma", "alpha"), function(t) max(range[, , , t]))
p1 <- ggplot(subset(data, param=="loglik")) +
geom_boxplot(aes(model, value)) +
facet_wrap(~ trait, scales="free_y")
p1
p2 <- ggplot(subset(data, param %in% c("sigma", "alpha") & model != "bm"),
aes(model, value, fill=regimes)) +
stat_summary(fun.y=mean, geom="bar", position="dodge") +
stat_summary(fun.data=median_hilow, geom="pointrange",
aes(color=regimes),
position = position_dodge(width=0.90), conf.int=.5) +
scale_y_log() +
facet_grid(param ~ trait, scales = "free_y") +
coord_cartesian(ylim=c(0,max(upper)), wise=TRUE)
p2
p3 <- ggplot(subset(data, param %in% c("sigma") ),
aes(model, value, fill=regimes)) +
stat_summary(fun.y=mean, geom="bar", position="dodge", alpha=.5) +
stat_summary(fun.data=median_hilow, geom="pointrange",
aes(color=regimes),
position = position_dodge(width=0.90), conf.int=.5) +
facet_grid(param ~ trait, scales = "free_y") +
coord_cartesian(ylim=c(0,upper["sigma"]), wise=TRUE) +
opts(title="sigma")
p3
p4 <- ggplot(subset(data, param %in% c("alpha") ),
aes(model, value, fill=regimes)) +
stat_summary(fun.y=mean, geom="bar", position="dodge", alpha=.5) +
stat_summary(fun.data=median_hilow, geom="pointrange",
aes(color=regimes),
position = position_dodge(width=0.90), conf.int=.5) +
facet_grid(param ~ trait, scales = "free_y") +
coord_cartesian(ylim=c(0,upper["alpha"]), wise=TRUE) +
opts(title = "alpha")
p4
save(list=ls, file = "~/public_html/data/labrids_par_bootstrap.rda")
Error: invalid first argument
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.