demo/sims_in_parallel.R

In kathalexknuts/NeutralTheorySimulations: Simulations for Neutral Theory Analysis

# install.packages("gtools", dependencies=TRUE, repos='http://cran.rstudio.com/')
# install.packages("vegan", dependencies=TRUE, repos='http://cran.rstudio.com/') 
# install.packages("doParallel", dependencies=TRUE, repos='http://cran.rstudio.com/')
# install.packages("foreach", dependencies=TRUE, repos='http://cran.rstudio.com/')
# install.packages("doSNOW", dependencies=TRUE, repos='http://cran.rstudio.com/')
# install.packages("dplyr", dependencies=TRUE, repos='http://cran.rstudio.com/')
# install.packages("parallel", dependencies=TRUE, repos='http://cran.rstudio.com/') #
# install.packages("devtools", dependencies=TRUE, repos='http://cran.rstudio.com/') 
library(gtools)
library(vegan)
library(doParallel)
library(foreach)
library(doSNOW)
library(dplyr)
library(parallel)
library(devtools)

install_github("kathalexknuts/NeutralTheorySimulations")
numCores <- detectCores() - 2

###############################
Nmeta <- 1000
n <- 1 #seq(30, 300, by = 30)  #number of samples
nOTU <-30 #number of OTUs in meta community
num_k <- 50 #just how many times to loop
ps <- seq(0, 0.9, by = 0.1)#c(0.2, 0.4, 0.6)
r <- seq(0.1, 1, by = 0.1)#seq(0.2, 0.8, by = 0.2)
Nlocal <- 1000
Nmeta_sd <-0
Nlocal_sd <- 0
m <- 0.3
sourc <- c(TRUE, FALSE)
pool_sim <- c(TRUE, FALSE)
pool_est <- c(TRUE, FALSE)
################

params <- NeutralTheorySims::param_mat(Nmeta, ps, r, m, Nlocal, nOTU, n, Nmeta_sd, Nlocal_sd, num_k, sourc, pool_sim, pool_est)
sim_result_multi <- function(cur_params){
  return(sim_result(cur_params, FALSE))
}

groups <- rep(1:numCores, each = nrow(params)/numCores,length.out = nrow(params))
param_split <- split(params, as.factor(groups))



cl <- makeCluster(numCores)
registerDoParallel(cl)
registerDoSNOW(cl)



mcoptions <- list(preschedule=FALSE, set.seed=FALSE, cores = numCores)
out <- foreach(i=1:numCores, .combine=rbind,.verbose=T, .packages = c("NeutralTheorySims", "vegan", "gtools"),
               .export=c("sim_result_multi", "param_split")) %dopar% {
                 sink("./log.txt", append=TRUE)
                 sim_result_multi(param_split[[i]]) }

#out <- sim_result(params, FALSE)
write.csv(out, "~/NTSims/pooled_source/test4.csv")
stopCluster(cl)


require(ggplot2)
require(dplyr)
out <- read.csv("~/NTSims/pooled_source/test3.csv", header = TRUE)
out$H_new <- as.factor(out$H_new)
out$r <- as.factor(out$r)
out$H <- as.factor(out$H)
#out$nOTU <- as.factor(out$nOTU)
out$Nlocal <- as.factor(out$Nlocal)
out$H_new <- as.numeric(substr(as.character(out$H_new),1, 4))
out$H <- as.numeric(substr(as.character(out$H),1, 4))
out[which(out$H_new < 0.05),"H_new"] <- 0
out$Bias <- abs(as.numeric(out$Estimate - out$m))
out$m <- as.character(out$m)
out$pool_sim <- ifelse(as.character(out$pool_sim) == "1", "Pooled", "Not Pooled")
out$pool_est <- ifelse(as.character(out$pool_est) == "1", "Pooled", "Not Pooled")
out$sourc <- ifelse(as.character(out$sourc) == "1", "Source", "No Source")
out$combine <- paste(paste(paste("Source", out$sourc, sep = ": "), paste("Sim", out$pool_sim, sep = ": "), sep = "\n"), paste("Est", out$pool_est, sep = ": "), sep = "\n")

jpeg("~/NTSims/pooled_source/figure6.jpeg", units="in", width=7, height=5, res=300)
out %>% filter(r %in% c(0.1, 0.3, 0.5, 0.7, 0.9)) %>% ggplot(aes(x = meanBC, y = Estimate)) + geom_point(aes(color = combine)) + geom_line(aes(color = combine)) + facet_grid(H ~ combine) + 
  labs(title = "Migration Parameter Estimates for Simulated Data, Pooled Source", 
       subtitle = "Shannon Diversity Index in rows", y= "Estimate", x = "Mean Bray Curtis Dissimilarity")+
  geom_ribbon(aes(ymin=Lower,ymax=Upper),alpha=0.3) + geom_hline(aes(yintercept = as.numeric(m)), color = "red", linetype = 3)+
  theme(strip.text.x = element_blank())+
  theme(legend.position="bottom")+
  theme(legend.title=element_blank())
dev.off()