Skip to content

Instantly share code, notes, and snippets.

@jslefche
Last active January 26, 2023 15:45
Show Gist options
  • Save jslefche/543f4282b8e72cb9da2c451075ee0128 to your computer and use it in GitHub Desktop.
Save jslefche/543f4282b8e72cb9da2c451075ee0128 to your computer and use it in GitHub Desktop.
Decomposition of community thermal index into individual species contributions

A function that takes individual species thermal indices (STIs) and their abundances, and partitions their contributions to the change in the community thermal index through at least two time points.

example <- data.frame(
  site = "site1",
  date = c(rep(as.Date("2001-01-01"), 5), rep(as.Date("2002-01-01"), 5)),
  species = paste0("species", LETTERS[1:5]),
  STI = c(21.2, 23.4, 19.2, 26.1, 22.0),
  abundance = c(0, 0, 10, 8, 3, 1, 3, 5, 0, 0)
  )

cti_decomp(example)  
#' Decomposing temporal change in CTI into individual species contributions
#'
#' @param x a community matrix in long-form (observations as rows) with the following columns:
#' "site" = the site(s) at which the survey was conducted
#' "date" = the date(s) of the survey
#' "species" = the species observed at each site and date
#' "abundance" = the abundance of the species at each site and date
#' "STI" = the species' thermal index
#' @param scale whether the output should be scaled by the absolute change in CTI. Default is FALSE
#' @param log whether abundances should be log+1-transformed before proceeding. Default is FALSE
#' @param span whether the decomposition should be performed for all timepoints ("all"),
#' consecutive time points ("consecutive") or for just the two endpoints ("endpoints").
#' Default is "all"
#' @.progressBar whether a progress bar is printed. Default is TRUE
#'
#' @return data.frame containing abundance ("abund") of each species in each time period,
#' start and end dates, average STI across the assemblage at each site,
#' proportion ("p") of each species in each time period, and each species' contribution
#' to CTI ("species_contribution")
#'
#' @requires tidyverse
#'
#' @author Jon Lefcheck <lefcheckj@@si.edu>
#'
#' @example
#'
#' example <- data.frame(
#' site = "site1",
#' date = c(rep(as.Date("2001-01-01"), 5), rep(as.Date("2002-01-01"), 5)),
#' species = paste0("species", LETTERS[1:5]),
#' STI = c(21.2, 23.4, 19.2, 26.1, 22.0),
#' abundance = c(0, 0, 10, 8, 3, 1, 3, 5, 0, 0))
#'
#' cti_decomp(example)
#'
cti_decomp <- function(x, log = FALSE, scale = FALSE, span = "consecutive", .progressBar = TRUE) {
# Set progress bar
if(.progressBar == TRUE)
pb <- txtProgressBar(min = 0, max = length(unique(x$site)), style = 3, width = 50)
# Loop over site codes and compute date pairs and look at change in STI per species
out <- do.call(rbind, lapply(unique(x$site), function(i) {
# Subset site
x <- dplyr::filter(x, site == i)
if(log == TRUE) x$abundance <- log10(x$abundance + 1)
if(nrow(x) == 0 | length(unique(x$date)) == 1) data.frame() else {
# Generate date pairs
dates <- unique(x$date)[order(unique(x$date))]
if(span == "all") {
pairs <- as.data.frame(t(combn(as.character(dates), 2)))
colnames(pairs) <- c("start", "end")
pairs$start <- as.Date(pairs$start)
pairs$end <- as.Date(pairs$end)
} else if(span == "consecutive") {
pairs <- data.frame(start = dates[-length(dates)], end = dates[-1])
} else if(span == "endpoints") {
pairs <- data.frame(start = min(x$date), end = max(x$date))
} else stop("Span must be 'all' or 'endpoints.'")
# Kick back empty data.frame if there is only 1 survey date at the site
if(nrow(pairs) == 0 | all(pairs[1] == pairs[2])) data.frame() else {
# Loop over date pairs
out <- do.call(rbind, lapply(1:nrow(pairs), function(j) {
y <- dplyr::filter(x, date %in% c(pairs[j, 1:2]))
# Order dates
y <- y[order(y$date), ]
# Cast dates as rows
y <- tidyr::pivot_wider(y, names_from = "date", values_from = "abundance", values_fill = 0)
if(any(duplicated(y$species))) stop("The same species has different STIs in the dataset.")
# Add dates as rows
y$start_date <- colnames(y)[ncol(y)-1]
y$end_date <- colnames(y)[ncol(y)-1]
# Rename columns for abundance at time t1 and t2
colnames(y)[(ncol(y)-3):(ncol(y)-2)] <- c("abund_t1", "abund_t2")
# Compute mean STI
y$mean_STI <- mean(y$STI, na.rm = T)
# Compute proportional abundances
y$p1 <- y$abund_t1 / sum(y$abund_t1)
y$p2 <- y$abund_t2 / sum(y$abund_t2)
# Compute species contributions to delta CtI based on difference in proportion weighted by deviance in STI from community mean
y$species_contribution <- (y$p2 - y$p1) * (y$STI - y$mean_STI)
# Scale output by total change in CTI to sum to 1
if(scale == TRUE) y$species_contribution <- y$species_contribution / sum(y$species_contribution)
# Return output
return(y)
} ) )
# Update progress bar
if(.progressBar == TRUE) setTxtProgressBar(pb, which(i == unique(x$site)))
# Return data
return(out)
}
}
} ) )
return(out)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment