|
#' @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) |
|
|
|
} |