#' 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)
#'
#' # 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)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.