#' Running a Simulation for Neutral Theory analysis given matrix of parameters
#'
#' This function combines the functions mom2 and sim_communities and run_sim functions to run a full simulation (parameter inputs, generate data and estimate).
#' @param params This should be a matrix with parameters to use in a neutral theory simultaion. Each column should be a parameter. Each row is a different combination of parameter estimates
#' @param log This should either be TRUE or FALSE, indicating whether a log file should be generated to track progress
#' @return a data frame with each row representing a simulation. Each row include parameters used in the simulations and estiamtes (including 95% CI upper and lower bounds)
#' @export
sim_result <- function(params, log){
do.call("rbind", apply(matrix(1:nrow(params)), 1, function(x){
cur_params <- params[x,]
Nt0 <- as.numeric(cur_params[1]) #
ps <- as.numeric(cur_params[2]) #
r <- as.numeric(cur_params[3]) #
m <- as.numeric(cur_params[4]) #
Nlocal <- as.numeric(cur_params[5]) ##
nOTU <- as.numeric(cur_params[6])
n <- as.numeric(cur_params[7])
Nmeta_sd <- as.numeric(cur_params[8])
Nlocal_sd <- as.numeric(cur_params[9])
num_k <- as.numeric(cur_params[10])
sourc <- as.numeric(cur_params[11])
pool_sim <- as.numeric(cur_params[12])
pool_est <- as.numeric(cur_params[13])
names(cur_params) <- c("Nt0", "ps", "r", "m", "Nlocal","nOTU", "n", "Nmeta_sd", "Nlocal_sd", "num_k", "sourc", "pool_sim", "pool_est")
if(log == TRUE){ cat(paste(x, nrow(params), sep = "/"), file="log.txt", append=TRUE)}
sim_res_list <- sapply(1:num_k, function(k){run_sims(k, ps, r, n, m, Nt0, nOTU, Nlocal, Nmeta_sd, Nlocal_sd, sourc, pool_sim, pool_est)})
H <- mean(sim_res_list[4,])
H_new <- mean(sim_res_list[5,])
BC_local <- mean(sim_res_list[2,])
BC_meta <- mean(sim_res_list[3,])
sim_res_list <- sim_res_list[1,]
mse <- (1/n)*sum((m - sim_res_list)^2) #calc mse
cur_params <- c(m, ps, r, Nlocal, Nt0, nOTU, n, Nmeta_sd, Nlocal_sd, H, pool_sim, pool_est, sourc)
df_params <- matrix(cur_params,nrow=1, ncol = length(cur_params),byrow=TRUE)
cur_result <- data.frame(df_params, mean(sim_res_list), quantile(sim_res_list,0.975), quantile(sim_res_list,0.025), BC_local, BC_meta, H_new, mse)
names(cur_result) <- c("m", "ps", "r", "Nlocal", "Nmeta", "nOTU", "n", "Nmeta_sd", "Nlocal_sd", "H", "pool_sim", "pool_est", "sourc", "Estimate", "Upper", "Lower", "meanBC_local", "meanBC_meta", "H_new", "mse")
rownames(cur_result) <- NULL
return(cur_result)
}))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.