R/sim_result_function.R

In kathalexknuts/NeutralTheorySimulations: Simulations for Neutral Theory Analysis

#' 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)
  }))
}