R/metaweb_build.R

In alaindanet/SizeTrophicInteractions: What the Package Does (One Line, Title Case)

#' Metaweb building
#'
#'
#' @param data a dataframe containing species and size variable 
#' @param diet_shift data.frame containing species and the predation window.
#' @inheritParams compute_prey_size 
#' @inheritParams compute_classes 
#' @param fish_resource_method character Either "overlap", "willem" or "midpoint"
#'
#' @details class_method = "quantile" use the `quantile` function to create size
#'  class in nb_class . When class_method = "percentile", the range of the size distribution
#'  is divided by nb_class. The difference is that `quantile` start the first
#'  class at the minimum value of the size distribution while percentile begins
#'  at 0. It therefore gives slightly different results.
#'  When the predation window definition is defined as overlap or midpoint, there
#'  will be possibility of overlap between size class of a predator and a prey.
#'  It means that for trophic link established between a prey and a predator, the
#'  smallest predator of a given size class is smaller than the biggest prey
#'  fish. When pred_win_method is set to "no_overlap", the upper limit of the predation window of
#'  a predator size class is defined by the smallest predator of the size class,
#'  so there is no overlap between predator and prey. With midpoint option, the
#'  existence of a trophic link is defined by the middle of the predator size
#'  class.
#'
#' @return matrix
#'
#' @export
build_metaweb <- function(data, species, size, pred_win, beta_min, beta_max,
  fish_diet_shift, low_bound, upper_bound, fish,
  resource_diet_shift, class_method = "quantile",
  nb_class = 9, pred_win_method = "midpoint", na.rm = FALSE,
  replace_min_by_one = FALSE, fish_resource_method = "overlap") {

  #Capture var:
  species <- rlang::enquo(species)
  size <- rlang::enquo(size)
  low_bound <- rlang::enquo(low_bound)
  upper_bound <- rlang::enquo(upper_bound)
  fish <- rlang::enquo(fish)
  beta_min <- rlang::enquo(beta_min)
  beta_max <- rlang::enquo(beta_max)

  # Check
  data <- sanatize_metaweb(
    data = data,
    species = !!species,
    fish_diet_shift = fish_diet_shift,
    resource_diet_shift = resource_diet_shift,
    nb_class = nb_class
  )
  # TODO: Check concordance of column names between datasets
  ## eg. species, (life_stage), low_bound, upper_bound, fish


  # Build the size classes
  size_class <- compute_classes(data, group_var = !!species, var = !!size,
    class_method = class_method, na.rm = na.rm, replace_min_by_one =
      replace_min_by_one)


  ## Compute the th_prey size min & max + mid-point
  th_prey_size <- compute_prey_size(size_class, pred_win, !!species,
    !!beta_min, !!beta_max, pred_win_method = pred_win_method)

  ## Get piscivory index for each size class  
  piscivory_index <- compute_piscivory(size_class, fish_diet_shift,
    species = !!species, low_bound = !!low_bound, upper_bound = !!upper_bound,
    fish = !!fish)

  ## Fill the matrix
  int_matrices <- compute_links(size_class, th_prey_size, piscivory_index,
    fish_diet_shift, resource_diet_shift, !!species, !!fish, !!low_bound, !!upper_bound,
    fish_resource_method = fish_resource_method, pred_win_method = pred_win_method
  )
  fish_fish_int <- int_matrices$fish_fish_int
  fish_resource_int <- int_matrices$fish_resource_int

  # Resource-Resource matrix
  resource_list <- int_matrices$resource_list
  resource_resource_int <- t(resource_diet_shift[, resource_list])
  # Resource-Fish matrix
  resource_fish_int <- matrix(0, ncol = ncol(resource_resource_int),
    nrow = nrow(fish_fish_int)) #0 because no resources are eating fishes
  colnames(resource_resource_int) <- colnames(resource_fish_int) <- resource_list

  # Merge the matrix
  metaweb <- rbind(
    cbind(fish_fish_int, resource_fish_int),
    cbind(fish_resource_int, resource_resource_int)
    )

  list(
  metaweb = metaweb,
  species = int_matrices$species_list,
  resource = resource_list,
  nb_class = nb_class,
  size_class = size_class,
  piscivory_index = piscivory_index,
  th_prey_size = th_prey_size
  )

}

#' Compute species links
#'
#' Compute feeding links between species based on their size 
#'
#' @inheritParams compute_piscivory
#' @param th_prey_size data.frame generated by compute_prey_size
#' @param piscivory_index data.frame generated by compute_piscivory
#' @param resource_diet_shift data.frame generated by compute_prey_size
#'
#' @details
#'
#' @return list containing the species-species and the
#' species-resource interaction matrix
#' @export
compute_links <- function(size_class, th_prey_size, piscivory_index,
  fish_diet_shift, resource_diet_shift, species, fish, low_bound, upper_bound,
  fish_resource_method = "overlap", pred_win_method = "midpoint") {

  species <- rlang::enquo(species)
  low_bound <- rlang::enquo(low_bound)
  upper_bound <- rlang::enquo(upper_bound)
  fish <- rlang::enquo(fish)

  #Check
  stopifnot(fish_resource_method %in% c("overlap", "midpoint", "willem"))

  # Build the empty matrices
  ## Fish-Fish
  species_list <- unique(dplyr::select(size_class, !!species)) %>% unlist
  nb_species <- length(species_list)
  nb_class <- length(unique(size_class$class_id))
  fish_fish_int <- matrix(rep(0, (nb_species * nb_class) ^ 2), ncol = nb_species * nb_class)
  fish_class_names <- rep(species_list, each = nb_class) %>%
    paste(., seq(1, nb_class), sep = "_")
  rownames(fish_fish_int) <- colnames(fish_fish_int)  <- fish_class_names
  ## Fish-Resource
  resource_list <- dplyr::select(resource_diet_shift, !!species) %>% unlist
  names(resource_list) <- NULL
  nb_resource <- length(resource_list)
  fish_resource_int <- matrix(rep(0, nb_resource * nb_class * nb_species), ncol = nb_class * nb_species)
  rownames(fish_resource_int) <- resource_list
  colnames(fish_resource_int) <- fish_class_names
  
  # Set method of resource assignement 
    if (fish_resource_method == "midpoint") {
      condition <- rlang::quo(
    #size class midpoint in stage interval
    midpoint >= !!low_bound & midpoint <= !!upper_bound
      )
    } else if (fish_resource_method == "willem") {
      condition <- rlang::quo(
    # min stage in size class: ]min_pred;max_pred]
    (min_pred < !!low_bound & !!low_bound < max_pred) |#<-diff with overlap
      # max stage in size class:]min_pred;max_pred]
    (min_pred < !!upper_bound & !!upper_bound <= max_pred) |
    #stage class strictly around size class:]min_pred;max_pred]
    ## stage class
    (!!low_bound < min_pred & max_pred <= !!upper_bound) |#<-diff with overlap
    #size class strictly in stage class:]low_bound;upper_bound]
    (!!low_bound > min_pred & max_pred >= !!upper_bound)
      )
    } else {
     condition <- rlang::quo(
    # min stage in size class: ]min_pred;max_pred]
    (min_pred < !!low_bound & !!low_bound <= max_pred) |
      # max stage in size class:]min_pred;max_pred]
    (min_pred < !!upper_bound & !!upper_bound <= max_pred) |
    #stage class strictly around size class:]min_pred;max_pred]
    ## stage class
    (!!low_bound <= min_pred & max_pred < !!upper_bound) |
    #size class strictly in stage class:]low_bound;upper_bound]
    (!!low_bound > min_pred & max_pred >= !!upper_bound)
      )
    }

  # merge dataset 
  trophic_data <- dplyr::left_join(size_class, th_prey_size,
    by = c(rlang::quo_name(species), "class_id")) %>%
    dplyr::left_join(., piscivory_index, by = c(rlang::quo_name(species), "class_id")) %>%
    tidyr::unite(sp_class, !!species, class_id)
    
  pred_classes <- colnames(fish_fish_int)
  trophic_data <- trophic_data[order(match(unlist(trophic_data[, "sp_class"]), rownames(fish_fish_int))), ]
  prey_data <- dplyr::select(trophic_data, sp_class, lower, upper)

  for (j in 1:ncol(fish_fish_int)) {#Predators
    pred_data <- dplyr::filter(trophic_data, sp_class == pred_classes[j]) %>%
      unlist

    if (pred_data["pisc_index"] != 0) {

      min_prey   <- as.numeric(pred_data["min_prey"])
      max_prey   <- as.numeric(pred_data["max_prey"])
      pisc_index <- as.integer(pred_data["pisc_index"])

      ## Check if the middle of the prey size class is in the range of prey of the predator
      if (pred_win_method == "midpoint") {
      troph_link <- dplyr::mutate(prey_data,
    prey_mid = (lower + upper) / 2,
    troph_link = (min_prey <= prey_mid & max_prey  >= prey_mid) * pisc_index)
      } else {
    stop("pred_win_method is different from midpoint. Other methods are not yet implemented")
      }

      fish_int_values <- troph_link %>%
    dplyr::select(troph_link) %>% unlist
      fish_fish_int[, j] <- fish_int_values
    }

    #Fish-Resource matrix
    ## What the predator size class is eating: 
    ### Get predator size:
    min_pred <- as.numeric(pred_data["lower"])
    max_pred <- as.numeric(pred_data["upper"])
    ### Filter diet:
    #### get species name
    pred_diet <- fish_diet_shift %>%
      dplyr::filter(!!species ==
        stringr::str_remove(pred_classes[j], "_\\d+")
      ) %>%
      dplyr::select(-!!fish)

    if (fish_resource_method == "midpoint") {
      midpoint <- (min_pred + max_pred) / 2
    }
    pred_diet_class_match <- pred_diet %>%
      dplyr::filter(!!condition)
    if (nrow(pred_diet_class_match) == 0) {
      message(paste("Species/class", pred_classes[j], "had no matches with life stages."))

      resource_int_values <- rep(0, length(resource_list))
    } else {
      resource_int_values <- pred_diet_class_match %>%
    dplyr::select(!! resource_list) %>%
    tidyr::gather(resource, troph_index) %>%
    dplyr::group_by(resource) %>%
    dplyr::summarise(troph_index = (sum(troph_index) > 0) * 1) %>%
    tidyr::spread(resource, troph_index) %>%
    unlist
      good_order <- match(names(resource_int_values), resource_list)
      resource_int_values <- resource_int_values[good_order]
    }
    fish_resource_int[, j] <- resource_int_values
  }

  list(
    fish_fish_int = fish_fish_int,
    fish_resource_int = fish_resource_int,
    species_list = species_list,
    resource_list = resource_list
  )
}

#' Compute piscivory link 
#'
#' Determine which size class is piscivorous
#'
#' @param size_class a data.frame generated by compute_classes
#' @param fish_diet_shift a data.frame containing species, life stage, lower and
#' upper bound size for each stage
#' @param species variable name
#' @param low_bound variable containing lower limit of the stage 
#' @param upper_bound variable containing upper limit of the stage. Not used.
#' @param fish name of logical piscivory variable. 0 = no piscivory; 1 = piscivory
#'  
#' @details If given species class is considered as piscivor if the upper bound
#' of the class is greater than the mininum size to reach piscivory.
#'
#' @return a data.frame
#' @export
compute_piscivory <- function (size_class, fish_diet_shift, species, low_bound, upper_bound, fish) {

  species <- rlang::enquo(species)
  lower_stage_bound <- rlang::enquo(low_bound)
  upper_stage_bound <- rlang::enquo(upper_bound) #not use until now
  fish <- rlang::enquo(fish)

  # Check variables
  type_variable <- fish_diet_shift %>%
    #dplyr::select(!!lower_stage_bound) %>%#, !!upper_stage_bound
    dplyr::summarise(type = typeof(!!lower_stage_bound)) %>%
    dplyr::select(type) %>% unlist(.)
  if (any(!(type_variable %in% c("double", "integer")))) {
    stop("In fish_diet_shift, the low_bound variable is not numeric")
  }
  

  get_piscivory <- function(fish, lower) {
    if (sum(fish) > 0) { # Is there at least one picivory life stage?
      # Get min size at piscivory
      min(lower[which(fish != 0)])
    } else {
      NA
    }
  }
  ## Get minimum size for each species where there is piscivory:
  min_pisc <- fish_diet_shift %>%
    dplyr::group_by(!!species) %>%
    dplyr::summarise(min_pisc = get_piscivory(!!fish, !!lower_stage_bound)) %>%
    dplyr::select(!!species, min_pisc)

  ## Piscivory if upper bound of the size class is sup to the size min at which
  ## the species begins to be piscivorous (if it does)
  dplyr::left_join(size_class, min_pisc, by = rlang::quo_name(species)) %>%
    dplyr::mutate(pisc_index = ifelse(is.na(min_pisc), 0, (upper > min_pisc) * 1)) %>%
    dplyr::select(-min_pisc, - lower, - upper)
}

#' Compute prey size
#'
#' Compute prey size for class size species can eat. It is theoretical, indeep,
#' the information about real diet of the fish is not taken in account here. 
#'
#' @param class_size a data.frame generated by compute_classes   
#' @param pred_win data.frame containing species and the predation window
#'  parameters (\alpha and \beta).
#' @param pred_win_method character midpoint, overlap or no_overlap. Default to midpoint.
#' @param species blabla
#' @param beta_min column name containing beta_min parameters
#' @param beta_max column name contains beta_max parameters
#'
#' @details For now, the function compute the mininum prey size as follow:
#' beta_min * min_predator_size. But the general formula is in hte form:
#' alpha_min + beta_min * min_predator_size. But here, we do not have data about
#' alpha. We will add it if necessary.
#' 
#' @return matrix
#'
#' @export
compute_prey_size <- function(class_size, pred_win, species, beta_min, beta_max, pred_win_method = "midpoint") {

  #Get var:
  species  <- rlang::enquo(species)
  beta_min <- rlang::enquo(beta_min)
  beta_max <- rlang::enquo(beta_max)
  
  if (pred_win_method != "midpoint") {
    message("Other methods are not yet implemented")
  }

  # check species concondance
  class_sp <- dplyr::select(class_size, !!species)[[1]] %>% unique()
  pred_sp <- dplyr::select(pred_win, !!species)[[1]] %>% unique()
  missing_sp <- which(!class_sp %in% pred_sp)
  if (length(missing_sp) > 0) {
    msg <- paste0(
      "The following species ", 
      class_sp[missing_sp],
      " have no predation windows defined.", collapse = "\n")
    stop(msg)
  }

  dplyr::left_join(class_size, pred_win, by = rlang::quo_name(species)) %>%
    dplyr::mutate(
      min_prey = !!beta_min * ( (lower + upper) / 2),
      max_prey = !!beta_max * ( (lower + upper) / 2)
      ) %>%
  dplyr::select(!!species, class_id, min_prey, max_prey)
}

#' Compute classes
#'
#' @param size data.frame containing species and length.
#' @param group_var variable to group the data. E.g. species.
#' @param var variable containing the values to class. 
#' @param na.rm logical. Should the NAs be removed ? 
#' @param replace_min_by_one logical. In the first class, replace low bound by
#' 1.    
#' @inheritParams split_in_classes 
#' @details The replace_min_by_one option is to reflect the choice made by
#' Bonnafé et al. to replace the low bound of the first class (= lowest value)
#' by one. It leads to slighty different results (env. 0.6% of the links).
#'
#' @return blabla
#' @export
compute_classes <- function(size, group_var, var, class_method = "quantile",
  nb_class = 9, na.rm = FALSE, replace_min_by_one = FALSE) {

  stopifnot(is.numeric(nb_class))

  #Capture variables:
  group_var <- rlang::enquo(group_var)
  var <- rlang::enquo(var)

  if (na.rm) {
    size %<>% na.omit
    message("NAs has been removed with na.omit")
  } else {
    if (length(which(is.na(size))) > 0) {
      stop("There are NAs in your dataset. Please set na.rm = TRUE")
    }
  }

  size %<>%
    dplyr::group_by(!!group_var) %>%
    dplyr::select(!!group_var, !!var) # ensure that there is only the variable
  #of interest
  nested_size <- size %>%
    tidyr::nest()

  ## Check that there are at least 2 different values
  check_nb_data <- size %>%
    dplyr::summarise(nb = length(unique(!!var))) %>%
    dplyr::mutate(good = ifelse(nb >= 2, TRUE, FALSE)) %>%
    dplyr::select(-nb)

  ## Warn which species species had been removed
  if (!all(check_nb_data$good)) {
  species_not_good <- dplyr::filter(check_nb_data, good == FALSE) %>%
    dplyr::select(!!group_var) %>% unlist(.)
    

  nested_size %<>% dplyr::filter(!(!!group_var %in% species_not_good))

  if (nrow(nested_size) == 0) {
    stop("None of the species got more of two unique values. Check your dataset.")
  }

  warning("The following species had less than two unique size values, so we got rid of them:", species_not_good)

  }


  ## Compute classes
  nested_size %<>%
    dplyr::mutate(
      classes = purrr::map(data, split_in_classes, nb_class = nb_class,
    class_method = class_method)
      ) %>%
    dplyr::select(-data)

  output <- nested_size %>%
    tidyr::unnest(classes) %>%
    dplyr::ungroup()

  if (replace_min_by_one) {
    output[output$class_id == 1, "lower"] <- 1
  }

  return(output)
}

#' Split in classes
#'
#' @param to_class a numeric vector
#' @param class_method character percentile or quantile. Default to percentile. 
#' @param nb_class integer number of size class to create. Default to 9. 
#' @param round_limits logical Does the class limits should be rounded ?.
#' @export
split_in_classes <- function(to_class, class_method = "quantile",
  nb_class = 9, round_limits = TRUE) {

  stopifnot(class_method %in% c("percentile", "quantile"))

  if (class_method == "quantile") {
    if (is.list(to_class)) {
      to_class %<>% unlist
      stopifnot(is.numeric(to_class))
    }
    classified <- quantile(to_class, probs = seq(0, 1, by = 1 / nb_class), names = FALSE)
  } else {
    classified <- seq(min(to_class), max(to_class),
      by = (max(to_class) - min(to_class)) / nb_class)
  }

  if (round_limits) {
    classified %<>% round
  }
    # Determine the lower and upper limit of each size class:
    tibble::tibble(
      class_id = seq(1, nb_class),
      lower = classified[-length(classified)],
      upper = classified[-1]
    )
  }

#' Sanatize fish length for metaweb
#'
#' @param data blabla
#' @param species blabla
#' @param fish_diet_shift blabla
#' @param nb_class blabla
#'
#' @export
sanatize_metaweb <- function(
  data = NULL,
  species = NULL,
  fish_diet_shift = NULL,
  resource_diet_shift = NULL,
  nb_class = NULL) {

  #Capture var:
  species <- rlang::enquo(species)
  sp_chr <- rlang::quo_name(species)

  data_sp <- unique(data[[sp_chr]])
  onto_sp <- unique(fish_diet_shift[[sp_chr]])

  # Stop if resource species are not found in fish diet data
  mask_res_missing <-
    !resource_diet_shift[[sp_chr]] %in% colnames(fish_diet_shift)
  if (any(mask_res_missing)) {
    res_missing <- resource_diet_shift[[sp_chr]][mask_res_missing]
    msg <- c(
      "The following resources are not found in fish_diet_shift: ",
      paste(res_missing, "\n", sep = ", "),
      "They should be present to build the metaweb !")
    stop(msg)
  }

  ## Del species for which life traits are not filled
  sp_missing <- data_sp[!data_sp %in% onto_sp]
  if (length(sp_missing) != 0) {
    msg <- c(
      "For the following species, ",
      "the life stage are not found in fish_diet_shift are absent: ",
      paste(sp_missing, "\n", sep = ", "),
    "They have been removed from data")
    message(msg)

    data <- data[! data[[sp_chr]] %in% sp_missing, ]
  }

  # Del species for which the number of individuals is inferior to the number of
  # classes:
  nb_ind_sp <- data %>%
    dplyr::group_by(!!species) %>%
    dplyr::summarise(nind = dplyr::n())

  species_chr <- rlang::quo_name(species)
  mask <- nb_ind_sp$nind <= nb_class
  sp_to_low <- nb_ind_sp[mask, ][[species_chr]]
  if (length(sp_to_low) != 0) {
    msg <- c(
      "For the following species, ",
      "# obs is not superior to nb_class of the metaweb: ",
      paste(sp_to_low, "\n", sep = ", "),
    "They have been removed from data")
    message(msg)


    data <- data[! data[[sp_chr]] %in% sp_to_low, ]
  }
  data
}