Compare measures of phylogenetic signal on simulated trees

phylogeneticsignal_sim.R

#### Simulations
install.packages(c("ape","reshape2","ggplot2"))
require(ape); require(reshape2); require(ggplot2)
source("http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/phylosig/v0.3/phylosig.R")
source("http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/fastBM/v0.4/fastBM.R")

# Simulation function physig_sim
# input: x = number of species in tree
# output: a vector length two with (K, lamba) 
physig_sim <- function(x) {
tree <- rcoal(x)
traits <- fastBM(tree)
physig_k <- phylosig(tree, traits, method = "K")
physig_l <- phylosig(tree, traits, method = "lambda")$lambda
sigs <- c(physig_k, physig_l)
return(sigs)
}

# Run simulation
spnumbs <- c(5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 65, 80, 95)
rands <- 1000
siglist <- alply( spnumbs, 1, function(x) t(replicate(rands, physig_sim(x))), .progress="text")
sigdf <- as.data.frame(do.call(rbind, siglist))
names(sigdf) <- c("K", "L")
sigdf$numsp_ <- rep(spnumbs, each=rands)
sigdf_m <- melt(sigdf, id = 3)

# Plot results
plotdf <- ddply(sigdf_m, .(numsp_, variable), summarise, 
mean = mean(value),
se = sd(value)/sqrt(length(value))
)
limits <- aes(ymax = mean + se, ymin = mean - se)
dodge <- position_dodge(width=0.9)
ggplot(plotdf, aes(x = numsp_, y = mean, shape = variable)) +
geom_point(position=dodge) +
geom_errorbar(limits, position=dodge, width=0.25) +
geom_smooth()
ggsave("physig_sim.jpeg")