################## Simulate trees using function simbal
install.packages(c("plyr","ggplot2","ape","apTreeshape","bipartite","reshape2"))
library(plyr); library(ape); library(apTreeshape); library(bipartite); library(ggplot2); library(reshape2)
source("fastBM.r"); source("nodeHeights.r"); source("pbtree.r")
### plant trees
foo <- function(tips) {
trees_colless_plants <- simbal(t=tips, "colless", 1000, cutlow = -0.8, cuthigh = 0.8)
trees_colless_plants_bal <- compact(lapply(trees_colless_plants, function(x) x$colless_bal)) # get the balanced trees
trees_colless_plants_unbal <- compact(lapply(trees_colless_plants, function(x) x$colless_unbal)) # get the unbalanced trees
### animal trees
trees_colless_anim <- simbal(t=tips, "colless", 1000, cutlow = -0.8, cuthigh = 0.8)
trees_colless_anim_bal <- compact(lapply(trees_colless_anim, function(x) x$colless_bal)) # get the balanced trees
trees_colless_anim_unbal <- compact(lapply(trees_colless_anim, function(x) x$colless_unbal)) # get the unbalanced trees
################## Simulate traits on each tree
# Plants
t1_col_plants_bal <- llply(trees_colless_plants_bal, fastBM, a = 10, bounds=c(0,100))
t1_col_plants_unbal <- llply(trees_colless_plants_unbal, fastBM, a = 10, bounds=c(0,100))
t2_col_plants_bal <- llply(trees_colless_plants_bal, fastBM, a = 10, bounds=c(0,100))
t2_col_plants_unbal <- llply(trees_colless_plants_unbal, fastBM, a = 10, bounds=c(0,100))
abd_col_plants_bal <- round(rlnorm(length(trees_colless_plants_bal), meanlog=5, 1), 0)
abd_col_plants_unbal <- round(rlnorm(length(trees_colless_plants_unbal), meanlog=5, 1), 0)
# Animals
t1_col_anim_bal <- llply(trees_colless_anim_bal, fastBM, a = 10, bounds=c(0,100))
t1_col_anim_unbal <- llply(trees_colless_anim_unbal, fastBM, a = 10, bounds=c(0,100))
t2_col_anim_bal <- llply(trees_colless_anim_bal, fastBM, a = 10, bounds=c(0,100))
t2_col_anim_unbal <- llply(trees_colless_anim_unbal, fastBM, a = 10, bounds=c(0,100))
abd_col_anim_bal <- round(rlnorm(length(trees_colless_anim_bal), meanlog=5, 1), 0)
abd_col_anim_unbal <- round(rlnorm(length(trees_colless_anim_unbal), meanlog=5, 1), 0)
################## Get plant-animal phylogenetic tree pairs
# Balanced trees
all_bal <- list(trees_colless_plants_bal, trees_colless_anim_bal)
smaller_bal <- which.min(sapply(all_bal, length))
larger_bal <- which.max(sapply(all_bal, length))
tree_pairs_bal <- list(all_bal[[smaller_bal]], all_bal[[larger_bal]][1:length(all_bal[[smaller_bal]])])
# Unbalanced trees
all_unbal <- list(trees_colless_plants_unbal, trees_colless_anim_unbal)
smaller_unbal <- which.min(sapply(all_unbal, length))
larger_unbal <- which.max(sapply(all_unbal, length))
tree_pairs_unbal <- list(all_unbal[[smaller_unbal]], all_unbal[[larger_unbal]][1:length(all_unbal[[smaller_unbal]])])
# Pairs of traits from balanced trees
all_t1_bal <- list(t1_col_plants_bal[1:length(tree_pairs_bal[[1]])],
t1_col_anim_bal[1:length(tree_pairs_bal[[1]])])
all_t2_bal <- list(t2_col_plants_bal[1:length(tree_pairs_bal[[1]])],
t2_col_anim_bal[1:length(tree_pairs_bal[[1]])])
all_abd_bal <- list(abd_col_plants_bal[1:length(tree_pairs_bal[[1]])],
abd_col_anim_bal[1:length(tree_pairs_bal[[1]])])
# Pairs of traits form unbalanced trees
all_t1_unbal <- list(t1_col_plants_unbal[1:length(tree_pairs_unbal[[1]])],
t1_col_anim_unbal[1:length(tree_pairs_unbal[[1]])])
all_t2_unbal <- list(t2_col_plants_unbal[1:length(tree_pairs_unbal[[1]])],
t2_col_anim_unbal[1:length(tree_pairs_unbal[[1]])])
all_abd_unbal <- list(abd_col_plants_unbal[1:length(tree_pairs_unbal[[1]])],
abd_col_anim_unbal[1:length(tree_pairs_unbal[[1]])])
################## Simulate networks on each tree
# Simulate completely random networks, regardless of traits, etc.
mats_rand_bal <- sim_rand_nets(tree_pairs_bal)
mats_rand_unbal <- sim_rand_nets(tree_pairs_unbal)
# Simulate networks, with interactions propoprtional to trait matching
# ratio
mats_traits1_bal_ratio <- sim_traits_nets(all_t1_bal, method = "r", value = 1.5)
mats_traits1_unbal_ratio <- sim_traits_nets(all_t1_unbal, method = "r", value = 1.5)
mats_traits2_bal_ratio <- sim_traits_nets(all_t2_bal, method = "r", value = 1.5)
mats_traits2_unbal_ratio <- sim_traits_nets(all_t2_unbal, method = "r", value = 1.5)
mats_traitsabd_bal_ratio <- sim_traits_nets(all_abd_bal, method = "r", value = 1.5)
mats_traitsabd_unbal_ratio <- sim_traits_nets(all_abd_unbal, method = "r", value = 1.5)
# complementarity
mats_traits1_bal_comp <- sim_traits_nets(all_t1_bal, method = "c", value = 2)
mats_traits1_unbal_comp <- sim_traits_nets(all_t1_unbal, method = "c", value = 2)
mats_traits2_bal_comp <- sim_traits_nets(all_t2_bal, method = "c", value = 2)
mats_traits2_unbal_comp <- sim_traits_nets(all_t2_unbal, method = "c", value = 2)
mats_traitsabd_bal_comp <- sim_traits_nets(all_abd_bal, method = "c", value = 2)
mats_traitsabd_unbal_comp <- sim_traits_nets(all_abd_unbal, method = "c", value = 2)
# barrier (no need to give value parameter)
mats_traits1_bal_barr <- sim_traits_nets(all_t1_bal, method = 'b')
mats_traits1_unbal_barr <- sim_traits_nets(all_t1_unbal, method = 'b')
mats_traits2_bal_barr <- sim_traits_nets(all_t2_bal, method = 'b')
mats_traits2_unbal_barr <- sim_traits_nets(all_t2_unbal, method = 'b')
mats_traitsabd_bal_barr <- sim_traits_nets(all_abd_bal, method = 'b')
mats_traitsabd_unbal_barr <- sim_traits_nets(all_abd_unbal, method = 'b')
################## Calculate network metrics on matrices
df_rand <- getnetmets(mats_rand_bal, mats_rand_unbal) # random networks
df_traits1_ratio <- getnetmets(mats_traits1_bal_ratio, mats_traits1_unbal_ratio) # ratio traits
df_traits1_comp <- getnetmets(mats_traits1_bal_comp, mats_traits1_unbal_comp) # complementarity traits
df_traits1_barr <- getnetmets(mats_traits1_bal_barr, mats_traits1_unbal_barr) # barrier traits
df_traits2_ratio <- getnetmets(mats_traits2_bal_ratio, mats_traits2_unbal_ratio) # ratio traits
df_traits2_comp <- getnetmets(mats_traits2_bal_comp, mats_traits2_unbal_comp) # complementarity traits
df_traits2_barr <- getnetmets(mats_traits2_bal_barr, mats_traits2_unbal_barr) # barrier traits
alldat <- rbind(df_rand, df_traits1_ratio, df_traits2_ratio,
df_traits1_comp, df_traits2_comp,
df_traits1_barr, df_traits2_barr)
alldat$model <- c( rep("Random",nrow(df_rand)),
rep("Ratio 1",nrow(df_traits1_ratio)),
rep("Ratio 2",nrow(df_traits2_ratio)),
rep("Complementarity 1",nrow(df_traits1_comp)),
rep("Complementarity 2",nrow(df_traits2_comp)),
rep("Barrier 1",nrow(df_traits1_barr)),
rep("Barrier 2",nrow(df_traits2_barr))
)
alldat$numsp <- rep(tips, nrow(alldat))
alldat
}
out <- suppressMessages(llply(c(20, 40, 60, 80, 100, 140, 160), foo, .progress="text"))
outdf <- ldply(out)
outdf_melt <- melt(outdf, id.vars=c(1,3,4))
outdf_melt_ <- ddply(outdf_melt, .(type, model, numsp, variable), summarise,
mean = mean(value, na.rm=T),
se = sd(value, na.rm=T)/sqrt(na.omit(length(value)))
)
##### Visualize stuff
ggplot(outdf_melt_, aes(numsp, mean, colour = type)) +
geom_point(size=4, alpha=0.5) +
theme_bw() +
facet_wrap(~ model)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.