R/sim_setup.R

In plthompson/mcomsimr: Simulates Metacommunity Dynamics in R

#' Generate landscape
#'
#' Generates a landscape for metacommunity simulations
#'
#' @param patches number of patches to include in landscape
#' @param xy optional dataframe with x and y columns for patch coordinates
#' @param plot option to show plot of landscape
#'
#' @return landscape with x and y coordinates
#'
#' @author Patrick L. Thompson, \email{patrick.thompson@@zoology.ubc.ca}
#'
#' @examples
#' landscape_generate()
#'
#' @export
#'
landscape_generate <- function(patches = 100, xy, plot = TRUE) {
  if (missing(xy)){

    if(patches > 10000) stop("Maximum number of patches is 10000.")
    positions_linear = sample(0:9999, patches)
    landscape = data.frame(x = floor(positions_linear / 100)+1, y = (positions_linear %% 100) + 1)

    clusters <- hclust(dist(landscape),method = "ward.D2")

    landscape <- landscape[clusters$order, ]
    rownames(landscape) <- 1:patches

  } else {
    landscape <- xy
  }
  if (plot == TRUE){
    plot(landscape, pch = 19)
  }
  return (landscape)
}

#' Generate Dispersal Matrix
#'
#' Generates dispersal matrix for metacommunity simulations
#'
#' @param landscape landscape generated by landscape_generate()
#' @param torus whether to model the landscape as a torus
#' @param disp_mat optional matrix with each column specifying the probability that an individual disperses to each other patch (row)
#' @param kernel_exp the exponential rate at which dispersal decreases as a function of the distance between patches
#' @param plot option to show plot of environmental variation
#'
#' @return matrix with dispersal probabilities
#'
#' @author Patrick L. Thompson, \email{patrick.thompson@@zoology.ubc.ca}
#'
#' @examples
#' dispersal_matrix(landscape_generate())
#'
#' @import ggplot2
#' @import dplyr
#' @importFrom magrittr %>%
#' @importFrom som.nn dist.torus
#'
#' @export
#'
dispersal_matrix <- function(landscape, torus = TRUE, disp_mat, kernel_exp = 0.1, plot = TRUE){
  if (missing(disp_mat)){
    if(torus == TRUE){
      dist_mat <- as.matrix(som.nn::dist.torus(coors = landscape))
    } else {
      dist_mat <- as.matrix(dist(landscape))
    }

    disp_mat <- exp(-kernel_exp * dist_mat)
    diag(disp_mat) <- 0
    disp_mat <- apply(disp_mat, 1, function(x) x / sum(x))
  } else {
    disp_mat <- disp_mat
    rownames(disp_mat) <- 1:nrow(disp_mat)
    colnames(disp_mat) <- 1:ncol(disp_mat)
    if (is.matrix(disp_mat) == FALSE) stop ("disp_mat is not a matrix")
    if (nrow(disp_mat) != nrow(landscape) | ncol(disp_mat) != nrow(landscape)) stop ("disp_mat does not have a row and column for each patch in landscape")
  }

  if (sum(colSums(disp_mat) > 1.001) > 0) warning ("dispersal from a patch to all others exceeds 100%. Make sure the rowSums(disp_mat) <= 1")
  if (sum(colSums(disp_mat) < 0.999) > 0) warning ("dispersal from a patch to all others is less than 100%. Some dispersing individuals will be lost from the metacommunity")

  if (plot == TRUE){
    g <- as.data.frame(disp_mat) %>%
      dplyr::mutate(to.patch = rownames(disp_mat)) %>%
      tidyr::gather(key = from.patch, value = dispersal, -to.patch) %>%
      dplyr::mutate(from.patch = as.numeric(as.character(from.patch)),
                    to.patch = as.numeric(as.character(to.patch))) %>%
      ggplot2::ggplot(ggplot2::aes(x = from.patch, y = to.patch, fill = dispersal))+
      ggplot2::geom_tile()+
      scale_fill_viridis_c()

    print(g)
  }

  return (disp_mat)
}

#' Generate Environment
#'
#' Generates density independent environmental conditions for metacommunity simulations
#'
#' @param landscape landscape generated by landscape_generate()
#' @param env.df optional dataframe with environmental conditions with columns: env1, patch, time
#' @param env1Scale scale of temporal environmental autocorrelation between -2 (anticorrelated) and 2 (correlated), default is 2
#' @param timesteps number of timesteps to simulate
#' @param plot option to show plot of environmental variation
#'
#' @return dataframe with x and y coordinates, time, and environmental conditions
#'
#' @author Patrick L. Thompson, \email{patrick.thompson@@zoology.ubc.ca}
#'
#' @examples
#' env_generate(landscape_generate())
#'
#' @importFrom synchrony phase.partnered
#' @import ggplot2
#'
#' @export
#'
env_generate <- function(landscape, env.df, env1Scale = 2, timesteps = 1000, plot = TRUE){
  if (missing(env.df)){
    repeat {
      env.df <- data.frame()
      for(i in 1:nrow(landscape)){
        env1 = phase.partnered(n = timesteps, gamma = env1Scale, mu = 0.5, sigma = 0.25)$timeseries[,1]
        env.df <- rbind(env.df, data.frame(env1 = vegan::decostand(env1,method = "range"), patch = i, time = 1:timesteps))
      }
      env.initial <- env.df[env.df$time == 1,]
      if((max(env.initial$env1)-min(env.initial$env1)) > 0.6) {break}
    }
  } else {
    if(all.equal(names(env.df), c("env1", "patch", "time")) != TRUE) stop("env.df must be a dataframe with columns: env1, patch, time")
  }

  if(plot == TRUE){
    g<-ggplot2::ggplot(env.df, aes(x = time, y = env1, group = patch, color = factor(patch)))+
      ggplot2::geom_line()+
      scale_color_viridis_d(guide = "none")
    print(g)
  }

  return(env.df)
  print("This version differs from Thompson et al. 2020 in that it does not produce spatially autocorrelated environmental variables.")
}

#' Generate Species Env. Traits
#'
#' Generates species specific traits for density independent environmental responses
#'
#' @param species number of species to simulate
#' @param max_r intrinsic growth rate in optimal environment, can be single value or vector of length species
#' @param min_env minium environmental optima
#' @param max_env minium environmental optima
#' @param env_niche_breadth standard deviation of environmental niche breadth, can be single value or vector of length species
#' @param optima optional values of environmental optima, should be a vector of length species
#' @param optima_spacing "even" or "random" to specify how optima should be distributed
#' @param plot option to show plot of environmental variation
#'
#' @return dataframe an optima and niche breadth for each species
#'
#' @author Patrick L. Thompson, \email{patrick.thompson@@zoology.ubc.ca}
#'
#' @examples
#' env_traits(species = 10)
#'
#' @export
#'
env_traits <- function(species, max_r = 5, min_env = 0, max_env = 1, env_niche_breadth = 0.5, optima, plot = TRUE, optima_spacing = "random"){
  if (missing(optima)){
    if(optima_spacing == "even"){
      optima <- seq(from = 0,to = 1,length = species)
    }
    if(optima_spacing == "random"){
      optima <- runif(n = species, min = min_env, max = max_env)
    }
  } else {
    if(length(optima)!=species) stop("optima is not a vector of length species")
    if(class(optima)!="numeric") stop("optima is not a numeric vector")
  }
  env_traits.df <- data.frame(species = 1:species, optima = optima, env_niche_breadth = env_niche_breadth, max_r = max_r)

  if(plot == TRUE){
    matplot(sapply(X = 1:species, FUN = function(x) {
      exp(-((env_traits.df$optima[x]-seq(min_env, max_env, length = 30))/(2*env_traits.df$env_niche_breadth[x]))^2)
    })*rep(max_r,each = 30), type = "l", lty = 1, ylab = "r", xlab = "environment", ylim = c(0,max(max_r)))

  }
  return(env_traits.df)
}

#' Generate Species Interaction Matrix
#'
#' Generates density dependent matrix of per capita competition
#'
#' @param species number of species to simulate
#' @param intra intraspecific competition coefficient, single value or vector of length species
#' @param min_inter min interspecific comp. coefficient
#' @param max_inter max interspecific comp. coefficient
#' @param int_mat option to supply externally generated competition matrix
#' @param comp_scaler value to multiply all competition coefficients by
#' @param plot option to show plot of competition coefficients
#'
#' @return species interaction matrix
#'
#' @author Patrick L. Thompson, \email{patrick.thompson@@zoology.ubc.ca}
#'
#' @examples
#' env_traits(species = 10)
#'
#' @import ggplot2
#' @import dplyr
#'
#' @export
#'
species_int_mat <- function(species, intra = 1, min_inter = 0, max_inter = 1.5, int_mat, comp_scaler = 0.05, plot = TRUE){
  if (missing(int_mat)){
    int_mat <- matrix(runif(n = species*species, min = min_inter, max = max_inter), nrow = species, ncol = species)
    diag(int_mat) <- intra
    int_mat <- int_mat * comp_scaler
  } else {
    if (is.matrix(int_mat) == FALSE) stop("int_mat must be a matrix")
    if (sum(dim(int_mat) != c(species,species))>0) stop("int_mat must be a matrix with a row and column for each species")
    if (is.numeric(int_mat) == FALSE) stop("int_mat must be numeric")
  }

  if (plot == TRUE){
    colnames(int_mat)<- 1:species
    g <- as.data.frame(int_mat) %>%
      dplyr::mutate(i = 1:species) %>%
      tidyr::gather(key = j, value = competition, -i) %>%
      dplyr::mutate(i = as.numeric(as.character(i)),
                    j = as.numeric(as.character(j))) %>%
      ggplot2::ggplot(ggplot2::aes(x = i, y = j, fill = competition))+
      ggplot2::geom_tile()+
      scale_fill_viridis_c(option = "E")

    print(g)
  }
  return(int_mat)
}