Price Equation for Diversity-Function

README.md

The Price Equation for Partitioning Diversity Effects on Ecosystem Function

This function takes a data.frame corresponding to the site-by-species "functioning" matrix (where cells contain the values of the ecosystem function), and returns the five additive components of the Price equation.

EXAMPLE

# Example 1: all species contribute equally to functioning and all occur at the baseline site
# Here, RICH_L should be negative and equal the total number of unshared species at each site
# COMP and CDE terms should be zero

ex1 <- t(data.frame(
  comm1 = c(1, 0, 1, 0, 0, 0),
  comm2 = c(0, 1, 1, 1, 0, 0),
  comm3 = c(1, 1, 1, 1, 1, 1)
))

price(ex1, relativize = F)

# Example 2: all species contribute equally but there are some that occur only at the comparison sites
# Here, RICH_G should equal the number of unshared species at each site (1)
# COMP and CDE terms should be zero

ex2 <- t(data.frame(
  comm1 = c(1, 0, 1, 0, 0, 1),
  comm2 = c(0, 1, 1, 1, 0, 0),
  comm3 = c(1, 0, 1, 1, 1, 0)
))

price(ex2, relativize = F)

# Example 3: species unique to the baseline site contribute disproportionately to functioning
# Here, the COMP_L term should be negative
# COMP_G and CDE should be zero

ex3 <- t(data.frame(
  comm1 = c(1, 0, 1, 0, 0, 0),
  comm2 = c(0, 1, 1, 1, 0, 0),
  comm3 = c(1, 0, 1, 1, 2, 2)
))

price(ex3, relativize = F)

# Example 4: species unique to the comparison sites contribute disproportionately to functioning
# Here, the COMP_G term should be positive
# COMP_L and CDE should be zero

ex4 <- t(data.frame(
  comm1 = c(0, 0, 1, 0, 2, 0),
  comm2 = c(0, 2, 1, 0, 0, 0),
  comm3 = c(1, 0, 1, 1, 0, 1)
))

price(ex4, relativize = F)

# Example 5: species shared among all sites contribute more to functioning at the baseline site
# Here, the CDE should be negative

ex5 <- t(data.frame(
  comm1 = c(1, 0, 1, 1, 0, 0),
  comm2 = c(1, 1, 1, 1, 0, 0),
  comm3 = c(2, 1, 2, 2, 1, 1)
))

price(ex5, relativize = F)

# Example 5: species shared among all sites contribute more to functioning at the comparison sites
# Here, the CDE should be positive

ex6 <- t(data.frame(
  comm1 = c(2, 0, 2, 1, 0, 0),
  comm2 = c(2, 1, 0, 0, 0, 0),
  comm3 = c(1, 1, 1, 1, 1, 1)
))

price(ex6, relativize = F)

price.R

#' `price` an R function for translating the extended Price equation to ecosystem functioning

#' Author: Jon Lefcheck
#' Last updated: 21 August 2019

#' @param mat a site-by-species "functioning" matrix (cells are values of ecosystem function)
#' @param base a vector specifying the baseline community: "best" for the highest performing community;
#' all" to compute all comparisons; or the row name(s) or number(s) of the baseline community
#' @param standardize whether results are scaled by the maximum level of functioning to the scale (-1, 1)
#' @param avg whether components should be averaged across the baseline site(s)
#' @param avglvl the top percentage of sites in terms of functioning that should be used as the baseline
#'
#' @return a data.frame 

price <- function(mat, base = "best", standardize = TRUE, avg = FALSE, avglvl = 0.1) {
  
  # Replace NAs with zeros
  mat[is.na(mat)] <- 0
  
  if(!all(apply(mat, 2, is.numeric))) stop("Remove non-numeric columns from matrix")
  
  # Remove species that do not contribute to functioning at all
  mat <- mat[, colSums(mat) != 0]
  
  # Assign rownames if not already present
  if(is.null(rownames(mat))) rownames(mat) <- 1:nrow(mat)
  
  # Identify baseline sites for comparison
  baseline <- getBaseline(mat, base = base, avg = avg, avglvl = avglvl) 
  
  # Loop over baselines
  dat <- do.call(rbind, lapply(baseline, function(i) {
    
    # Remove other baseline sets from the data
    if(length(baseline) > 1) {
      
      remove <- baseline[!baseline %in% i]
      
      m <- mat[-remove, ]
      
    } else m <- mat
    
    # Redefine baseline row
    i <- which(rownames(m) == rownames(mat)[i])
    
    # Get total number of species from baseline site
    s <- sum(m[i, ] > 0, na.rm = T)
    
    # Calculate average contribution to functioning at the baseline site
    zbar <- sum(m[i, ], na.rm = T) / s
    
    # Loop over comparisons
    dat <- do.call(rbind, lapply(which(1:nrow(m) != i), function(j) {
      
      # Get total number of species from the comparison site
      sprime <- sum(m[j, ] > 0, na.rm = T)
      
      # Get vector of species in common with both sites
      shared <- apply(m[c(i, j), ], 2, function(x) ifelse(any(x == 0), 0, 1))
      
      # Sum number of shared species in both sites
      sc <- sum(shared, na.rm = T)
      
      # Calculate mean level of functioning at comparison site
      zprimebar <- sum(m[j, ], na.rm = T) / sprime
      
      # Calculate mean level of functioning among shared species at baseline site
      zcbar <- sum(as.numeric(shared * m[i, ]), na.rm = T) /
        sum(shared * m[i, ] > 0, na.rm = T)
      
      if(is.na(zcbar)) zcbar <- 0
      
      # Calculate mean level of functioning among shared species at comparison site
      zcprimebar <- sum(as.numeric(shared * m[j, ]), na.rm = T) /
        sum(shared * m[j, ] > 0, na.rm = T)
      
      if(is.na(zcprimebar)) zcprimebar <- 0
      
      # Calculate Price equation components
      dat <- data.frame(
        baseline.site = rownames(m)[i],
        comparison.site = rownames(m)[j],
        baseline.S = s,
        comparison.S = sprime,
        shared.S = sc,
        delta_FUNC = (sprime * zprimebar) - (s * zbar),
        RICH_LOSS = (sc - s) * zbar,
        RICH_GAIN = (sprime - sc) * zprimebar,
        COMP_LOSS = sc * (zcbar - zbar),
        COMP_GAIN = -sc * (zcprimebar - zprimebar),
        CDE = sc * (zcprimebar - zcbar)
      )
      
      return(dat)
      
    } ) )
    
    return(dat)
    
  } ) )
  
  # Remove rows where site == baseline.site
  dat <- dat[as.character(dat$comparison.site) != as.character(dat$baseline.site), ]
  
  # If avg == TRUE, average across sites
  if(avg == TRUE) { 
    
    dat <- aggregate(. ~ comparison.site, "mean", data = dat[, -(1)], na.rm = TRUE)
    
    dat <- cbind(baseline.site = paste("Average of", length(baseline), "sites"), dat) 
    
  }

  # Relativize output so units are bounds (-1, 1) for each component
  if(standardize == TRUE) dat[, 7:11] <- dat[, 7:11] / rowSums(abs(dat[, 7:11]), na.rm = TRUE)
  
  # Calculate aggregate gains, losses, and across both gains and losses
  dat$TOTAL_LOSSES <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T) 
  dat$TOTAL_GAINS <- rowSums(dat[, c("RICH_GAIN", "COMP_GAIN")], na.rm = T) 
  dat$TOTAL_DIV <- rowSums(dat[, c("RICH_LOSS", "RICH_GAIN", "COMP_LOSS", "COMP_GAIN")], na.rm = T) 
  
  return(dat)
  
}

#' `getBaseline` 

#' @param mat a site-by-species "functioning" matrix (cells are values of ecosystem function),
#' @param base the row name(s) or number(s) of the baseline community; defaults to highest productivity community
#' @param avg = TRUE, whether result should be averaged across some number of baseline sites 
#' @param avglvl the top X% of sites that should be used as the baseline and then results averaged, default is top 10%
#' @return integer corresponding to the row number of the baseline community from the site-by-species functioning matrix

getBaseline <- function(mat, base = "best", avg = FALSE, avglvl = 0.10) {
  
  if(avg == TRUE & base == "best") baseline <- which(rowSums(mat, na.rm = T) >= max(rowSums(mat, na.rm = T)) * (1 - avglvl)) else
    
    if(base == "best") baseline <- which.max(rowSums(mat, na.rm = T)) else
      
      if(base == "all") baseline <- 1:nrow(mat) else
        
        if(all(base != "best") & !is.numeric(base)) baseline <- which(rownames(mat) %in% base) else
          
          baseline <- base
        
        return(baseline)
        
} 

#' `price.sim` a function to simulate communities and apply the Price equation

#' @param nspecies number of species to be used in the simulation
#' @param ncomms number of simulated communities to be generated
#' @param bscaling a vector of the % by which biomass should be scaled between the baseline and comparison sites
#' @param sloss a vector of the number of species that should be dropped from the simulated communities
#' @param directional whether the species lost should be removed directionally (vs. randomly)
#' @param top.or.bottom if directional = TRUE, whether species should be removed from the top or bottom percentage of biomass
#' @param relativize whether the components should be relativized (default is TRUE)
#' @param ... additional arguments to the function `price`

price.sim <- function(nspecies, ncomms, bscaling = 0, sloss = 0, directional = FALSE, top.or.bottom = "top", relativize = FALSE, ...) {
  
  if(sloss < 0) sloss <- abs(sloss)
  
  # Create simulated community-by-species matrix
  mat <- matrix(runif(ncomms*nspecies, 0, 1), ncomms, nspecies) 
  
  # Get combinations of sloss and bscaling
  combos <- expand.grid(sloss, bscaling)
  
  # For each row, apply scaling found in combos, compute Price equation, and return results
  do.call(rbind, lapply(1:nrow(mat), function(i) {
    
    baseline.site <- mat[i, , drop = F]
    
    do.call(rbind, lapply(1:nrow(combos), function(j) {
      
      comparison.site <- baseline.site
      
      if(directional == FALSE) {
        
        # Randomly remove species
        comparison.site[1, sample(1:ncol(mat), combos[j, 1])] <- 0
        
      } else if(directional == TRUE) {
        
        if(top.or.bottom == "top") 
          
          remove <- rev(order(comparison.site))[1:combos[j, 1]]
        
        if(top.or.bottom == "bottom") 
          
          remove <- order(comparison.site)[1:combos[j, 1]]
        
        # Directionally remove species from the top or bottom
        comparison.site[1, remove] <- 0 
        
      }
      
      # Scale biomass
      comparison.site <- comparison.site * (1 - combos[j, 2])
      
      # Compute Price components
      data.frame(
        sloss = combos[j, 1],
        bscaling = combos[j, 2],
        comparison.site = i,
        price(rbind(baseline.site, comparison.site), relativize = relativize)
      )
      
    } ) )
    
  } ) )
  
}