R/run_sims_function.R

In kathalexknuts/NeutralTheorySimulations: Simulations for Neutral Theory Analysis

#' Running a Simulation for Neutral Theory analysis
#'
#' This function combines the functions mom2 and sim_communities to run a full simulation (generate data and estimate).
#' @param k This parameter sets the seed (also helpful for running this function in an apply loop)
#' @param ps this parameter adjusts the amount of variation between samples, with a range between 0 and 1. This value is the probability that a given OTU's relative abundance with be permuted in a sample. Values of 0 result in homologous taxa abundance distributions across all samples. Default in 1.
#' @param r this parameter adjusts the amount of variation between samples. Values can range from 0 to 1. r = 0 will result in a uniformly distributed relative abundance distribution within a sample. Default is 1.
#' @param n the number of samples
#' @param m true migration parameter
#' @param Nt0 total number of reads in each meta community samples. Note that if the argument for Nmeta_sd != 0, then this value will be the mean number of reads across meta community samples.
#' @param nOTU number of OTUs
#' @param Nlocal total number of reads in each local community samples. Note that if the argument for Nlocal_sd != 0, then this value will be the mean number of reads across meta community samples.
#' @param Nmeta_sd this parameter adjusts the amount of variation in the number of reads in meta community samples. Specifically, this value is the standard deviation used in a normal distribution with mean Nmeta from which samples sizes are generated
#' @param Nlocal_sd this parameter adjusts the amount of variation in the number of reads in local community samples. Specifically, this value is the standard deviation used in a normal distribution with mean Nlocal from which samples sizes are generated from.
#' @param sourc to use source data in estimation of migration parameter
#' @param pool_sim to pool meta/source community data during data simulation
#' @param pool_est to pool meta/source community data during estimation
#' @return a matrix of 4 rows: row 1 = migration parameter estimate, row 2 = mean bray curtis dissimilarity for meta community, row 3 = the estimated average shannon diversity for each sample. row 4 =  average shannon diversity after the multiple samples have been generated given the inputed value for ps and OTUs with no reads across any samples are removed
#' @export


run_sims <- function(k = 1, ps = 1, r = 1, n, m, Nt0, nOTU, Nlocal,  Nmeta_sd, Nlocal_sd, sourc = FALSE, pool_sim = FALSE, pool_est = FALSE){

  m_estimate <- matrix(0, nrow=1, ncol=1)
  bc_vals_temp <- matrix(0, nrow=1, ncol = 1)
  bc_vals_temp2 <- matrix(0, nrow=1, ncol = 1)
  if(sourc == TRUE){
    mom_cur <-get("mom2", envir = .GlobalEnv)
    set.seed(k)
    sim <- sim_communities(n=n, m=m, Nmeta=Nt0, r = r, nOTU = nOTU, ps = ps, Nlocal = Nlocal, Nmeta_sd= Nmeta_sd, Nlocal_sd = Nlocal_sd, pool_sim = pool_sim)
    est <- mom_cur(sim[[1]],sim[[2]], pool_est)
    m_estimate[1,1] <- est[[1]]
    bc_vals_temp[1, 1] <- est[[2]]
    bc_vals_temp2[1, 1] <- est[[3]]
    H <- sim[[3]]
    H_new <- sim[[4]]
  }
  
  else{
    mom_cur <-get("mom_nosource", envir = .GlobalEnv)
    set.seed(k)
    sim <- sim_communities(n=n, m=m, Nmeta=Nt0, r = r, nOTU = nOTU, ps = ps, Nlocal = Nlocal, Nmeta_sd= Nmeta_sd, Nlocal_sd = Nlocal_sd, pool_sim = pool_sim)
    est <- mom_cur(sim[[1]],sim[[2]], pool_est)
    m_estimate[1,1] <- est[[1]]
    bc_vals_temp[1, 1] <- est[[2]]
    bc_vals_temp2[1, 1] <- est[[3]]
    H <- sim[[3]]
    H_new <- sim[[4]]
  }


  return(rbind(rbind(rbind(m_estimate, rbind(bc_vals_temp, bc_vals_temp2)), H), H_new))
}