Simulate communities for BEF

README.md

simulateComm: simulate communities with predefined properties

Simulate community-by-species matrix of "functioning" (e.g., biomass) where the total community richness has a predefined correlation with total community biomass (i.e., the row sums).

# Generate community data
ncomms <- 100

nspecies <- 20

rho <- 0.5

mat <- simulateComm(rho, nspecies, ncomms, propzero = 0.8)

# Get correlation
cor(rowSums(mat), rowSums(mat > 0))

simulateComm.R

#' Generate species-by-site matrix with prespecified degree of across-site correlation 
#' between community richness and total community biomass
#' 
#' @param rho = the bivariate correlation between the rowSums and the rowSums > 0
#' @param nspecies = the number of species in the simulated data (ncol)
#' @param ncomms = the number of communities in the simulated data (nrow)
#' @param propzero = the proportion of the matrix that is zeros (default is 0.5)
#' 
simulateComm <- function(rho, nspecies, ncomms, propzero = 0.5) {
  
  # Create vector of 1s and replace with the specified proprtion of 0s
  mat <- rep(1, ncomms * nspecies)
  
  remove <- sample(sum(mat), sum(mat) * propzero)
  
  mat[remove] <- 0
  
  # Convert vector to matrix
  mat <- matrix(mat, ncomms, nspecies)
  
  # For each row select random number of species (columns) to populate
  richness <- rowSums(mat > 0)
  
  # Get vector of biomasses with predefined correlation
  theta <- acos(rho)
  
  x <- runif(nrow(mat))
  
  X <- cbind(richness, x)
  
  Xctr <- scale(X)    
  
  Id <- diag(nrow(mat))
  
  Q <- qr.Q(qr(Xctr[, 1, drop = FALSE]))
  
  P <- tcrossprod(Q)
  
  x2o <- (Id - P) %*% Xctr[, 2]  
  
  Xc2 <- cbind(Xctr[, 1], x2o)
  
  Y <- Xc2 %*% diag(1 / sqrt(colSums(Xc2^2)))
  
  biomass <- Y[, 2] + (1 / tan(theta)) * Y[, 1]
  
  # Add minimum to biomass to set lower bound at zero
  biomass <- (biomass + abs(min(biomass)) + 0.1) * 1000
  
  # Randomly apportion biomass to species
  for(j in 1:nrow(mat)) {
    
    # Get columns to populate
    bcols <- richness[j]
    
    # Get biomass sum
    bsum <- biomass[j]
    
    # Divide biomass into random bins corresponding to the number of columns
    bvec <- bsum * (sample(1:bcols) / sum(sample(1:bcols)))
    
    # Assign to community matrix
    mat[j, mat[j,] == 1] <- bvec
    
  }
  
  # print(cor(rowSums(mat), apply(mat, 1, function(x) sum(x > 0)), use = "complete.obs"))
  
  return(mat)
  
}