Skip to content

Instantly share code, notes, and snippets.

@jslefche
Last active August 21, 2019 17:57
Show Gist options
  • Save jslefche/66d461cdf6ffbbc81c31 to your computer and use it in GitHub Desktop.
Save jslefche/66d461cdf6ffbbc81c31 to your computer and use it in GitHub Desktop.
Price Equation for Diversity-Function

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` 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)
)
} ) )
} ) )
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment