knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(gstat)
library(sp)

Set simulation parameters

#n for the n x n landscape grid
grid_size <- 100

#number of species
n_sp <- 50

#number of environmental variables
n_env <- 3

Create gridded environmental variables

# unconditional gaussian simulation of environmental variables using gstat
grid <- expand.grid(1:grid_size, 1:grid_size)
colnames(grid) <- c("x", "y")

# define the gstat object (spatial model)
spat_mod <- gstat(formula=z~1, locations=~x+y, dummy=T, beta=1, model=vgm(psill=0.025,model="Exp",range=5), nmax=20, data = grid)

# simulate specified number of environmental variables
sims <- predict(spat_mod, newdata=grid, nsim=n_env)

# visualize
gridded(sims) = ~x+y
spplot(sims)

#convert to data frame of cells and simulated variables
env_data <- as.data.frame(sims)

Create dataframe of environmental suitability for each species

#function for getting a random range within an existing range
ran_range <- function(range){
  lower <- runif(1, range[1], range[2])

  if(lower + .5 > range[2]){
    upper <- range[2]
  } else{
    upper <- runif(1, lower + .5, range[2])}

  return(data.frame(type = c("min", "max"), range = c(lower, upper)))
}

### get ranges for each species as two min and max columns ###
#get ranges for environmental variables 
env_ranges <- map(env_data %>% select(-x, -y), ~range(.x))

#dataframe of environmental decision rules
species_env <- map_dfr(1:n_sp, 
                       ~map_dfr(env_ranges, ~ran_range(.x), .id = "env_var") %>% 
                         mutate(species = paste0("spec", .x))
                       )

Create matrix species exclusion rules

#random assignment of cooccurrence, 0 = can occur, 1 = excludes
co <- sample.int(2, n_sp^2, TRUE) - 1
dim(co) <- c(n_sp, n_sp)

#make symmetric
ind <- lower.tri(co) 
co[ind] <- t(co)[ind]

#format matrix as dataframe with row and column names
species_names <- paste0("spec", 1:n_sp)
co <- as.data.frame(co)
colnames(co) <- species_names
co <- bind_cols(co, species = species_names)

Identify potential species occupancy based on environmental requirements (environmental filtering)

get_spec_occ <- function(species_name, data = env_data, ranges = species_env){
  #use case_when to check if a cell's env value is in the species' range, assign occurence or absence accordingly
  data %>% mutate(!!species_name := case_when(sim1 > filter(ranges, species == species_name, env_var == "sim1", type == "min")$range &
                                        sim1 < filter(ranges, species == species_name, env_var == "sim1", type == "max")$range &
                                        sim2 > filter(ranges, species == species_name, env_var == "sim2", type == "min")$range &
                                        sim2 < filter(ranges, species == species_name, env_var == "sim2", type == "max")$range &
                                        sim3 > filter(ranges, species == species_name, env_var == "sim3", type == "min")$range &
                                        sim3 < filter(ranges, species == species_name, env_var == "sim3", type == "max")$range ~ 1,
                                      TRUE ~ 0)) %>%
    select(species_name) #return just the occurence column 
}

#map accross species to get an occurence column for each
occ <- map_dfc(paste0("spec", 1:n_sp), get_spec_occ)

Exclude species based on other species present (species filtering)

exclude_species <-
  function(species_name,
           exclusion_mat = co,
           occ_mat = occ) {

    #get list of species that can exclude specified species
    ex_species <- exclusion_mat %>%
      select(species, species_name) %>%
      filter(!!sym(species_name) == 1) %>%
      pull(species)

      #if there are exclusion species, remove occurence in cells where they already occur
      if (length(ex_species) > 0) {
        #if at least one species that can exclude occurs, assign as not present
        occ_mat %>% mutate(!!species_name := case_when(sum(!!!syms(ex_species)) >= 1 ~ 0, 
                                                   TRUE ~ !!sym(species_name))) %>% #otherwise accept existing value
          select(species_name)
      } else{
        return(select(occ, !!sym(species_name)))
      }

  }

data <- map_dfc(paste0("spec", 1:n_sp), exclude_species) %>% 
  bind_cols(env_data, .)

Programatically generate case example

# labels <- c("low", "mid", "high")
# breaks <- c(2, 4)
# 
# cases <- lapply(seq_along(breaks), function(i) {
#   substitute(carb <= breaks[i] ~ labels[i], list(i = i))
# })
# cases <- append(cases, substitute(TRUE ~ labels[length(labels)]))

Features to add realistic environmental variable ranges programmatically generate environmental filtering case_when sparsity of exclusion table randomize species exclusion algorithm



karinorman/communityDimensions documentation built on Jan. 22, 2020, 6:59 p.m.