R/sim_result_function.R

In kathalexknuts/ntsims:

#' 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])
    Nmeta_sd <- as.numeric(cur_params[7])
    Nlocal_sd <- as.numeric(cur_params[8])
    num_k <- as.numeric(cur_params[9])
    sourc <- as.numeric(cur_params[10])

    names(cur_params) <- c("Nt0", "ps", "r", "m", "Nlocal","nOTU", "Nmeta_sd", "Nlocal_sd", "num_k", "sourc")
    
    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, m, Nt0, nOTU, Nlocal, Nmeta_sd, Nlocal_sd, sourc)})
    H <- mean(sim_res_list[2,])
    H_new <- 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, Nmeta_sd, Nlocal_sd, H, 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), H_new, mse)
    names(cur_result) <- c("m", "ps", "r", "Nlocal", "Nmeta", "nOTU", "Nmeta_sd", "Nlocal_sd", "H", "sourc", "Estimate", "Upper", "Lower", "H_new", "MSE")
    rownames(cur_result) <- NULL
    
    return(cur_result)
  }))
}