R/pool_build.R

Defines functions pool_build

Documented in pool_build

#' Build the regional species pool
#' @description Regional species pool 
#' @name pool_build
#' @param seed the seed for sampling. Default = 1.
#' @param save_data whether to save species characteristic in a RData file. Logical.
#' @param data_name if `save_data` == TRUE, specify the name of the RData file. String.
#' @return  
#' @examples
#' # Load packages and example parameters
#' library(tidyverse)
#' data(example_parameter)
#' 
#' # Build regional species pool 
#' pool_build(save_data = FALSE)
#' 
#' @export

# Function that defines species abilities: m, u, w, and stoi
pool_build <- function (
  seed = 1, # Setting seed
  save_data = TRUE,
  data_name = "species_pool"
) {
  # Seed ----
  set.seed(seed)
  
  # Death rate; m ----
  ## A numeric vector with length of S
  m <- rep(mi, S_pool) 
  
  # Uptaking rate; u ----
  ## A matrix with P rows and S_pool columns
  
  if (function_group == 0) {
    # Randomly drawn from uniform distribution, where u sum is 1
    f <- rep(1, S_pool)
    u <- matrix(runif((P - 1) * S_pool), nrow = P-1, ncol = S_pool) %>%
      rbind(0, 1) %>% apply(2, sort)
    u <- matrix(u[-1, ] - u[-(P + 1), ], P, S_pool)
    
  } else if (is.numeric(function_group)) {
    # Funtional groups that have preference in uptaking resources
    # See "Emergent property SI" for how to conduct phylogenetic sampling
    suppressMessages(library(MCMCpack)) # Load this package for function that draws from dirichlet distribution
    
    # Empty matrix of uptaking rates
    u <- matrix(NA, nrow = P, ncol = S_pool)
    
    # Groups for species
    f <- rep(1:function_group, each = floor(S_pool / function_group))[1:S_pool]
    
    # Group: concentration parameter for resource in a group
    a_family <- rep(NA, S_pool)
    for (i in 1:S_pool) {
      a_preferred <- rnorm(1, group_mu, group_sigma) # Preferred resource
      a_preferred <- min(a_preferred, 1) # Set the maximum to 1
      a_other <- runif(P - 1) # Other resources
      
      u[f[i], i] <- a_preferred
      u[-f[i], i] <- (1 - a_preferred) * a_other / sum(a_other)
      a_family[i] <- a_preferred
    }
    
    # Species: total resource capacity 
    t_species <- rnorm(S_pool, 1, 0.01) 
    
    # Species: proportion of fixed budget of uptake resources
    K_family <- 100 # High K: species are very similar
    for (i in 1: S_pool) u[,i] <- rdirichlet(1, K_family * u[,i])
    
    # Uptaking rate
    u <- t_species * u * 0.9
  }
  
  # Probability of being a generalist; d_general ----
  u <- u * matrix(runif(S_pool * P) <= d_general, nrow = P, ncol = S_pool)
  
  # Conversion efficiency; w ----
  w <- matrix(wi, P, S_pool)
  
  # Stoichiometic matrix; stoi ----
  ## One available resource, no cross_feeding
  if (P == 1) {
    stoi <- array(0, dim = c(P, P, S_pool))
    for (i in 1:S_pool) diag(stoi[,,i]) <- -1
    
  ## Multiple resources, cross-feeding, different metabolism
  } else if (cross_feeding == TRUE & same_stoi == FALSE) {
    # Cross-feeding
    stoi <- array(NA, dim = c(P, P, S_pool))
    for (k in 1:S_pool) {
      for (i in 1:P) {
        w_i <- w[,k] / w[i,k] # denoted as ki in the note
        mu <- runif(1)
        d <- runif(P)
        stoi[i,,k] = (mu * d) / c(d %*% w_i)
      }
    }
    for (i in 1:S_pool) diag(stoi[,,i]) <- -1
  
  ## Multiple resources, cross-feeding, same metabolism
  } else if (cross_feeding == TRUE & same_stoi == TRUE) {
    
    # Cross-feeding
    stoi <- array(NA, dim = c(P, P, S_pool))
    for (i in 1:P) {
      w_i <- w[,1] / w[i,1] # denoted as ki in the note
      mu <- runif(1)
      d <- runif(P)
      stoi[i,,1] = (mu * d) / c(d %*% w_i)
    }
    diag(stoi[,,1]) <- -1
    stoi[,,-1] <- stoi[,,1]

  }
  
  # Species-abundance distribution; trade-off ----
  if (trade_off == TRUE) {
    abundance_species <- 
      rlnorm(n = S_pool, meanlog = 1, sdlog = 2) %>%
      round(2) %>% sort()
    
    a_preferred <- rep(NA, S_pool)
    for (i in 1:S_pool) a_preferred[i] <- u[f[i], i]
    abundance_species[rank(a_preferred)] <- abundance_species
    
  } else if (trade_off == FALSE) {
    if(species_abundance == "log-normal") {
      abundance_species <- rlnorm(n = S_pool, meanlog = 1, sdlog = 2) %>% round(2)
    } else if (species_abundance == "uniform") {
      abundance_species <- rep(1, S_pool)
    }
  }
  abundance_species <- abundance_species / sum(abundance_species)
  
  # Save the parameter in a RData ----
  if (save_data == TRUE) {
    save(m, f, u, w, stoi, abundance_species, 
         file = paste0("data/", data_name))
  }
  
  # Return the output to global environment
  m <<- m
  f <<- f
  u <<- u
  w <<- w
  stoi <<- stoi
  abundance_species <<- abundance_species
  
  # Regional pool with abundance and species names
  regional_pool <<- setNames(abundance_species, paste0("X", sprintf("%05d", 1:S_pool)))
}
Chang-Yu-Chang/MigrationCommunity documentation built on April 22, 2019, 7:36 p.m.