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) | |
| } |