The Price Equation: Creating the Optimal Society

priceEquation.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 values are the values of the ecosystem function that are to be compared), and returns the five additive components of the Price equation.

EXAMPLE

# Create example communities
func = t(data.frame(
  baseline = c(0, 10, 4, 0, 1),
  random.gain = c(0, 10, 4, 5, 1),
  random.loss = c(0, 10, 0, 0, 1),
  compositional.gain = c(20, 10, 4, 0, 1),
  compositional.loss = c(0, 0, 4, 0, 1),
  abundance = c(0, 8, 20, 0, 50)
))

# In the above examples, the value should load exclusively on only one component of the Price equation
price(func, base = "baseline")

                site baseline.site deltaS RICH_L    RICH_G COMP_L COMP_G ABUN    T_Tprime RICH_L_COMP_L RICH_G_COMP_G RICH_G_L_COMP_G_L
1        random.gain      baseline     -1      0 0.5714286    0.0      0    0  0.07936508           0.0          0.25              0.25
2        random.loss      baseline      1     -1 0.0000000    0.2      0    0 -0.06349206          -0.4          0.00             -0.20
3 compositional.gain      baseline     -1      0 1.0000000    0.0      1    0  0.31746032           0.0          1.00              1.00
4 compositional.loss      baseline      1     -1 0.0000000   -1.0      0    0 -0.15873016          -1.0          0.00             -0.50
5          abundance      baseline      0      0 0.0000000    0.0      0    1  1.00000000           0.0          0.00              0.00

priceEquation.R

#' @param func = a site-by-species "functioning" matrix (cells are values of ecosystem function),
#' @param base = the row name or number of the baseline community; defaults to highest productivity community
#' @param standardize = TRUE, whether results are scaled by the maximum to units (-1, 1), 
#' @param avg = TRUE, whether 
#' @param avglvl = The top X% of sites that should be used as the baseline and then results averaged, default is top 90%
#' @return data.frame of Price components & summed effects

price = function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl = 0.90) {
  
  # Replace NAs with zeros
  func[is.na(func)] = 0
  
  # Remove species that do not contribute to functioning at all
  func = func[, colSums(func) != 0]
  
  # Identify baseline sites for comparison
  baseline = get.baseline(func, 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(avg == TRUE & length(baseline) > 1) {
        
        remove = baseline[!baseline %in% i]
        
        func.temp = func[-remove, ]
        
      } else func.temp = func 
      
      # Get total number of species from baseline site
      s = sum(func.temp[i, ] > 0, na.rm = T)
      
      # Calculate mean level of functioning at baseline site
      zbar = sum(func.temp[i, ], na.rm = T) / sum(func.temp[i, ] > 0, na.rm = T)
      
      # Determine new baseline
      new.baseline = get.baseline(func.temp, base = base, avg = avg, avglvl = avglvl)
      
      # Next, loop over comparisons
      dat = do.call(rbind, lapply((1:nrow(func.temp))[-new.baseline], function(j) {
        
        # Get total number of species from the comparison site
        sprime = sum(func.temp[j, ] > 0, na.rm = T)
        
        # Get vector of species in common with both sites
        shared = apply(func.temp[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(func.temp[j, ], na.rm = T) / sum(func.temp[j, ] > 0, na.rm = T)
        
        # Calculate mean level of functioning among shared species at baseline site
        zcbar = sum(as.numeric(shared * func.temp[i, ]), na.rm = T) /
          sum(shared * func.temp[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 * func.temp[j, ]), na.rm = T) /
          sum(shared * func.temp[j, ] > 0, na.rm = T)
        
        if(is.na(zcprimebar)) zcprimebar = 0
        
        # Calculate Price equation components
        datf = data.frame(
          site = rownames(func.temp)[j],
          baseline.site = rownames(func.temp)[i],
          deltaS = s - sprime,
          RICH_L = (sc - s) * zbar,
          RICH_G = (sprime - sc) * zprimebar,
          COMP_L = sc * (zcbar - zbar),
          COMP_G = -sc * (zcprimebar - zprimebar),
          ABUN = sc * (zcprimebar - zcbar),
          T_Tprime = (sprime * zprimebar) - (s * zbar)
        )
        
        return(datf)
        
      } ) )
      
      # Calculate aggregate gains, losses, and across both gains and losses
      dat$RICH_L_COMP_L =rowSums(dat[, c("RICH_L", "COMP_L")], na.rm = T) 
      dat$RICH_G_COMP_G = rowSums(dat[, c("RICH_G", "COMP_G")], na.rm = T) 
      dat$RICH_G_L_COMP_G_L = rowSums(dat[, c("RICH_L", "RICH_G", "COMP_L", "COMP_G")], na.rm = T) 
      
      return(dat)
      
    } ) )
    
    # Remove rows where site == baseline.site
    dat = dat[as.character(dat$site) != as.character(dat$baseline.site), ]
    
    # If avg == TRUE, average across sites
    if(avg == TRUE) { 
      
      dat = aggregate(. ~ site + baseline.site, "mean", data = dat)
      
      if(length(baseline) > 1) dat$baseline.site = NA
      
    }
    
    # Relativize output so units are bounds (-1, 1) for each component
    if(standardize == TRUE) dat[, -(1:3)] = apply(dat[, -(1:3)], 2, function(x) x / max(abs(x), na.rm = T))
    
    return(dat)
    
}