Skip to content

Instantly share code, notes, and snippets.

@jslefche

jslefche/README.md

Last active August 21, 2019 17:57
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save jslefche/66d461cdf6ffbbc81c31 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save jslefche/66d461cdf6ffbbc81c31 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. Jon Lefcheck revised this gist Aug 21, 2019. 1 changed file with 5 additions and 9 deletions.
    14 changes: 5 additions & 9 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,18 +1,18 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 06 July 2018
    #' 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 relativize whether results are scaled by the maximum level of functioning to the scale (-1, 1)
    #' @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", relativize = TRUE, avg = FALSE, avglvl = 0.1) {
    price <- function(mat, base = "best", standardize = TRUE, avg = FALSE, avglvl = 0.1) {

    # Replace NAs with zeros
    mat[is.na(mat)] <- 0
    @@ -110,13 +110,9 @@ price <- function(mat, base = "best", relativize = TRUE, avg = FALSE, avglvl = 0
    dat <- cbind(baseline.site = paste("Average of", length(baseline), "sites"), dat)

    }

    # Relativize output so units are bounds (-1, 1) for each component
    if(relativize == TRUE) {

    dat[, (7:11)] <- t(apply(dat[, (7:11)], 1, function(x) x / max(abs(x), na.rm = TRUE)))

    }
    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)
  2. Jon Lefcheck revised this gist Jul 6, 2018. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -56,7 +56,7 @@ price <- function(mat, base = "best", relativize = TRUE, avg = FALSE, avglvl = 0
    sprime <- sum(m[j, ] > 0, na.rm = T)

    # Get vector of species in common with both sites
    shared <- colSums(m[c(i, j), ] > 0) / 2
    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)
  3. Jon Lefcheck revised this gist Jul 6, 2018. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -38,7 +38,7 @@ price <- function(mat, base = "best", relativize = TRUE, avg = FALSE, avglvl = 0

    m <- mat[-remove, ]

    } m <- mat
    } else m <- mat

    # Redefine baseline row
    i <- which(rownames(m) == rownames(mat)[i])
  4. Jon Lefcheck revised this gist Jul 6, 2018. 1 changed file with 28 additions and 16 deletions.
    44 changes: 28 additions & 16 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -31,43 +31,55 @@ price <- function(mat, base = "best", relativize = TRUE, avg = FALSE, avglvl = 0
    # 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, ]

    } m <- mat

    # Redefine baseline row
    i <- which(rownames(m) == rownames(mat)[i])

    # Get total number of species from baseline site
    s <- sum(mat[i, ] > 0, na.rm = T)
    s <- sum(m[i, ] > 0, na.rm = T)

    # Calculate average contribution to functioning at the baseline site
    zbar <- sum(mat[i, ], na.rm = T) / s
    zbar <- sum(m[i, ], na.rm = T) / s

    # Loop over comparisons
    dat <- do.call(rbind, lapply(which(1:nrow(mat) != i), function(j) {

    dat <- do.call(rbind, lapply(which(1:nrow(m) != i), function(j) {
    # Get total number of species from the comparison site
    sprime <- sum(mat[j, ] > 0, na.rm = T)
    sprime <- sum(m[j, ] > 0, na.rm = T)

    # Get vector of species in common with both sites
    shared <- colSums(mat[c(i, j), ] > 0) / 2
    shared <- colSums(m[c(i, j), ] > 0) / 2

    # Sum number of shared species in both sites
    sc <- sum(shared, na.rm = T)

    # Calculate mean level of functioning at comparison site
    zprimebar <- sum(mat[j, ], na.rm = T) / sprime
    zprimebar <- sum(m[j, ], na.rm = T) / sprime

    # Calculate mean level of functioning among shared species at baseline site
    zcbar <- sum(as.numeric(shared * mat[i, ]), na.rm = T) /
    sum(shared * mat[j, ] > 0, na.rm = T)
    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 * mat[j, ]), na.rm = T) /
    sum(shared * mat[j, ] > 0, na.rm = T)
    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(mat)[i],
    comparison.site = rownames(mat)[j],
    baseline.site = rownames(m)[i],
    comparison.site = rownames(m)[j],
    baseline.S = s,
    comparison.S = sprime,
    shared.S = sc,
    @@ -134,9 +146,9 @@ getBaseline <- function(mat, base = "best", avg = FALSE, avglvl = 0.10) {
    if(all(base != "best") & !is.numeric(base)) baseline <- which(rownames(mat) %in% base) else

    baseline <- base

    return(baseline)

    return(baseline)
    }

    #' `price.sim` a function to simulate communities and apply the Price equation
  5. Jon Lefcheck revised this gist Jul 6, 2018. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -102,7 +102,7 @@ price <- function(mat, base = "best", relativize = TRUE, avg = FALSE, avglvl = 0
    # Relativize output so units are bounds (-1, 1) for each component
    if(relativize == TRUE) {

    dat[, (6:10)] <- t(apply(dat[, (6:10)], 1, function(x) x / max(abs(x), na.rm = TRUE)))
    dat[, (7:11)] <- t(apply(dat[, (7:11)], 1, function(x) x / max(abs(x), na.rm = TRUE)))

    }

  6. Jon Lefcheck revised this gist Jul 6, 2018. 1 changed file with 18 additions and 16 deletions.
    34 changes: 18 additions & 16 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -3,15 +3,16 @@
    #' Author: Jon Lefcheck
    #' Last updated: 06 July 2018

    #' @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, or "best" for the highest performing community, or "all"
    #' to compute all comparisons; default is "all"
    #' @param relativize = TRUE, whether results are scaled by the maximum to be on the scale (-1, 1)
    #' @param avg = FALSE, whether result should be averaged across some number of baseline sites
    #' @param avglvl = 0.1, the top X% of sites in terms of functinoing that should be used as the baseline, and then results averaged
    #' @return data.frame of Price components & summed effects
    #' @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 relativize 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 = "all", relativize = TRUE, avg = FALSE, avglvl = 0.1) {
    price <- function(mat, base = "best", relativize = TRUE, avg = FALSE, avglvl = 0.1) {

    # Replace NAs with zeros
    mat[is.na(mat)] <- 0
    @@ -36,14 +37,14 @@ price <- function(mat, base = "all", relativize = TRUE, avg = FALSE, avglvl = 0.
    # Calculate average contribution to functioning at the baseline site
    zbar <- sum(mat[i, ], na.rm = T) / s

    # Next, loop over comparisons
    dat <- do.call(rbind, lapply(which(1:nrow(j) != i), function(j) {
    # Loop over comparisons
    dat <- do.call(rbind, lapply(which(1:nrow(mat) != i), function(j) {

    # Get total number of species from the comparison site
    sprime <- sum(mat[j, ] > 0, na.rm = T)

    # Get vector of species in common with both sites
    shared <- apply(mat[c(i, j), ], 2, function(x) ifelse(any(x == 0), 0, 1))
    shared <- colSums(mat[c(i, j), ] > 0) / 2

    # Sum number of shared species in both sites
    sc <- sum(shared, na.rm = T)
    @@ -53,7 +54,7 @@ price <- function(mat, base = "all", relativize = TRUE, avg = FALSE, avglvl = 0.

    # Calculate mean level of functioning among shared species at baseline site
    zcbar <- sum(as.numeric(shared * mat[i, ]), na.rm = T) /
    sum(shared * mat[i, ] > 0, na.rm = T)
    sum(shared * mat[j, ] > 0, na.rm = T)

    if(is.na(zcbar)) zcbar <- 0

    @@ -67,8 +68,9 @@ price <- function(mat, base = "all", relativize = TRUE, avg = FALSE, avglvl = 0.
    dat <- data.frame(
    baseline.site = rownames(mat)[i],
    comparison.site = rownames(mat)[j],
    delta_S = sprime - s,
    shared_S = sc,
    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,
    @@ -107,7 +109,7 @@ price <- function(mat, base = "all", relativize = TRUE, avg = FALSE, avglvl = 0.
    # 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_DIVERSITY <- rowSums(dat[, c("RICH_LOSS", "RICH_GAIN", "COMP_LOSS", "COMP_GAIN")], na.rm = T)
    dat$TOTAL_DIV <- rowSums(dat[, c("RICH_LOSS", "RICH_GAIN", "COMP_LOSS", "COMP_GAIN")], na.rm = T)

    return(dat)

  7. Jon Lefcheck revised this gist Jul 5, 2018. No changes.
  8. Jon Lefcheck revised this gist Jul 5, 2018. 1 changed file with 3 additions and 3 deletions.
    6 changes: 3 additions & 3 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -37,7 +37,7 @@ price <- function(mat, base = "all", relativize = TRUE, avg = FALSE, avglvl = 0.
    zbar <- sum(mat[i, ], na.rm = T) / s

    # Next, loop over comparisons
    dat <- do.call(rbind, lapply(baseline[baseline != i], function(j) {
    dat <- do.call(rbind, lapply(which(1:nrow(j) != i), function(j) {

    # Get total number of species from the comparison site
    sprime <- sum(mat[j, ] > 0, na.rm = T)
    @@ -91,7 +91,7 @@ price <- function(mat, base = "all", relativize = TRUE, avg = FALSE, avglvl = 0.
    # If avg == TRUE, average across sites
    if(avg == TRUE) {

    dat <- aggregate(. ~ comparison.site, "mean", data = dat[, -(1)])
    dat <- aggregate(. ~ comparison.site, "mean", data = dat[, -(1)], na.rm = TRUE)

    dat <- cbind(baseline.site = paste("Average of", length(baseline), "sites"), dat)

    @@ -100,7 +100,7 @@ price <- function(mat, base = "all", relativize = TRUE, avg = FALSE, avglvl = 0.
    # Relativize output so units are bounds (-1, 1) for each component
    if(relativize == TRUE) {

    dat[, (6:10)] <- t(apply(dat[, (6:10)], 1, function(x) x / max(abs(x))))
    dat[, (6:10)] <- t(apply(dat[, (6:10)], 1, function(x) x / max(abs(x), na.rm = TRUE)))

    }

  9. Jon Lefcheck revised this gist Jul 5, 2018. 1 changed file with 3 additions and 3 deletions.
    6 changes: 3 additions & 3 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -105,9 +105,9 @@ price <- function(mat, base = "all", relativize = TRUE, avg = FALSE, avglvl = 0.
    }

    # Calculate aggregate gains, losses, and across both gains and losses
    dat$TOTAL_DIV_LOSS <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_DIV_GAIN <- 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)
    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_DIVERSITY <- rowSums(dat[, c("RICH_LOSS", "RICH_GAIN", "COMP_LOSS", "COMP_GAIN")], na.rm = T)

    return(dat)

  10. Jon Lefcheck revised this gist Jul 5, 2018. 1 changed file with 6 additions and 6 deletions.
    12 changes: 6 additions & 6 deletions README.md
    View file
    Original file line number Diff line number Diff line change
    @@ -14,7 +14,7 @@ ex1 <- t(data.frame(
    comm3 = c(1, 1, 1, 1, 1, 1)
    ))
    price(ex1, standardize = F)
    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)
    @@ -26,7 +26,7 @@ ex2 <- t(data.frame(
    comm3 = c(1, 0, 1, 1, 1, 0)
    ))
    price(ex2, standardize = F)
    price(ex2, relativize = F)
    # Example 3: species unique to the baseline site contribute disproportionately to functioning
    # Here, the COMP_L term should be negative
    @@ -38,7 +38,7 @@ ex3 <- t(data.frame(
    comm3 = c(1, 0, 1, 1, 2, 2)
    ))
    price(ex3, standardize = F)
    price(ex3, relativize = F)
    # Example 4: species unique to the comparison sites contribute disproportionately to functioning
    # Here, the COMP_G term should be positive
    @@ -50,7 +50,7 @@ ex4 <- t(data.frame(
    comm3 = c(1, 0, 1, 1, 0, 1)
    ))
    price(ex4, standardize = F)
    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
    @@ -61,7 +61,7 @@ ex5 <- t(data.frame(
    comm3 = c(2, 1, 2, 2, 1, 1)
    ))
    price(ex5, standardize = F)
    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
    @@ -72,5 +72,5 @@ ex6 <- t(data.frame(
    comm3 = c(1, 1, 1, 1, 1, 1)
    ))
    price(ex6, standardize = F)
    price(ex6, relativize = F)
    ```
  11. Jon Lefcheck revised this gist Jul 5, 2018. 1 changed file with 31 additions and 40 deletions.
    71 changes: 31 additions & 40 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,83 +1,72 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 18 June 2018
    #' Last updated: 06 July 2018

    #' @param func 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 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, or "best" for the highest performing community, or "all"
    #' to compute all comparisons; default is "all"
    #' @param relativize = TRUE, whether results are scaled by the maximum to be on the scale (-1, 1)
    #' @param avg = FALSE, whether result should be averaged across some number of baseline sites
    #' @param avglvl = 0.1, the top X% of sites in terms of functinoing that should be used as the baseline, and then results averaged
    #' @return data.frame of Price components & summed effects

    price <- function(func, base = "best", relativize = TRUE, avg = FALSE, avglvl = 0.1) {
    price <- function(mat, base = "all", relativize = TRUE, avg = FALSE, avglvl = 0.1) {

    # Replace NAs with zeros
    func[is.na(func)] <- 0
    mat[is.na(mat)] <- 0

    if(!all(apply(func, 2, is.numeric))) stop("Remove non-numeric columns from matrix")
    if(!all(apply(mat, 2, is.numeric))) stop("Remove non-numeric columns from matrix")

    # Remove species that do not contribute to functioning at all
    func <- func[, colSums(func) != 0]
    mat <- mat[, colSums(mat) != 0]

    # Assign rownames if not already present
    if(is.null(rownames(func))) rownames(func) <- 1:nrow(func)
    if(is.null(rownames(mat))) rownames(mat) <- 1:nrow(mat)

    # Identify baseline sites for comparison
    baseline <- getBaseline(func, base = base, avg = avg, avglvl = avglvl)
    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]

    func.temp <- func[-remove, ]

    } else func.temp <- func

    # Redefine baseline row
    brow <- which(rownames(func.temp) == rownames(func)[i])

    # Get total number of species from baseline site
    s <- sum(func.temp[brow, ] > 0, na.rm = T)
    s <- sum(mat[i, ] > 0, na.rm = T)

    # Calculate average contribution to functioning at the baseline site
    zbar <- sum(func.temp[brow, ], na.rm = T) / s
    zbar <- sum(mat[i, ], na.rm = T) / s

    # Next, loop over comparisons
    dat <- do.call(rbind, lapply((1:nrow(func.temp))[-brow], function(j) {
    dat <- do.call(rbind, lapply(baseline[baseline != i], function(j) {

    # Get total number of species from the comparison site
    sprime <- sum(func.temp[j, ] > 0, na.rm = T)
    sprime <- sum(mat[j, ] > 0, na.rm = T)

    # Get vector of species in common with both sites
    shared <- apply(func.temp[c(brow, j), ], 2, function(x) ifelse(any(x == 0), 0, 1))
    shared <- apply(mat[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) / sprime
    zprimebar <- sum(mat[j, ], na.rm = T) / sprime

    # Calculate mean level of functioning among shared species at baseline site
    zcbar <- sum(as.numeric(shared * func.temp[brow, ]), na.rm = T) /
    sum(shared * func.temp[brow, ] > 0, na.rm = T)
    zcbar <- sum(as.numeric(shared * mat[i, ]), na.rm = T) /
    sum(shared * mat[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)
    zcprimebar <- sum(as.numeric(shared * mat[j, ]), na.rm = T) /
    sum(shared * mat[j, ] > 0, na.rm = T)

    if(is.na(zcprimebar)) zcprimebar <- 0

    # Calculate Price equation components
    dat <- data.frame(
    baseline.site = rownames(func.temp)[brow],
    comparison.site = rownames(func.temp)[j],
    baseline.site = rownames(mat)[i],
    comparison.site = rownames(mat)[j],
    delta_S = sprime - s,
    shared_S = sc,
    delta_FUNC = (sprime * zprimebar) - (s * zbar),
    @@ -126,21 +115,23 @@ price <- function(func, base = "best", relativize = TRUE, avg = FALSE, avglvl =

    #' `getBaseline`

    #' @param func a site-by-species "functioning" matrix (cells are values of ecosystem function),
    #' @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(func, base = "best", avg = FALSE, avglvl = 0.10) {
    getBaseline <- function(mat, base = "best", avg = FALSE, avglvl = 0.10) {

    if(avg == TRUE & all(base == "best")) baseline <- which(rowSums(func, na.rm = T) >= max(rowSums(func, na.rm = T)) * (1 - avglvl)) else
    if(avg == TRUE & base == "best") baseline <- which(rowSums(mat, na.rm = T) >= max(rowSums(mat, na.rm = T)) * (1 - avglvl)) else

    if(all(base == "best")) baseline <- which.max(rowSums(func, na.rm = T)) else
    if(base == "best") baseline <- which.max(rowSums(mat, na.rm = T)) else

    if(all(base != "best") & !is.numeric(base)) baseline <- which(rownames(func) %in% base) else
    if(base == "all") baseline <- 1:nrow(mat) else

    baseline <- base
    if(all(base != "best") & !is.numeric(base)) baseline <- which(rownames(mat) %in% base) else

    baseline <- base

    return(baseline)

    @@ -211,4 +202,4 @@ price.sim <- function(nspecies, ncomms, bscaling = 0, sloss = 0, directional = F

    } ) )

    }
    }
  12. Jon Lefcheck revised this gist Jun 18, 2018. 1 changed file with 5 additions and 3 deletions.
    8 changes: 5 additions & 3 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,7 +1,7 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 05 February 2018
    #' Last updated: 18 June 2018

    #' @param func 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
    @@ -116,8 +116,8 @@ price <- function(func, base = "best", relativize = TRUE, avg = FALSE, avglvl =
    }

    # Calculate aggregate gains, losses, and across both gains and losses
    dat$TOTAL_LOSS <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_GAIN <- rowSums(dat[, c("RICH_GAIN", "COMP_GAIN")], na.rm = T)
    dat$TOTAL_DIV_LOSS <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_DIV_GAIN <- 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)
    @@ -159,6 +159,8 @@ getBaseline <- function(func, base = "best", avg = FALSE, avglvl = 0.10) {

    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)

  13. Jon Lefcheck revised this gist Feb 5, 2018. 1 changed file with 2 additions and 2 deletions.
    4 changes: 2 additions & 2 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,7 +1,7 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 02 February 2018
    #' Last updated: 05 February 2018

    #' @param func 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
    @@ -111,7 +111,7 @@ price <- function(func, base = "best", relativize = TRUE, avg = FALSE, avglvl =
    # Relativize output so units are bounds (-1, 1) for each component
    if(relativize == TRUE) {

    dat[, (5:9)] <- t(apply(dat[, (5:9)], 1, function(x) x / max(abs(x))))
    dat[, (6:10)] <- t(apply(dat[, (6:10)], 1, function(x) x / max(abs(x))))

    }

  14. Jon Lefcheck revised this gist Feb 2, 2018. 1 changed file with 3 additions and 2 deletions.
    5 changes: 3 additions & 2 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,7 +1,7 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 27 November 2017
    #' Last updated: 02 February 2018

    #' @param func 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
    @@ -79,6 +79,7 @@ price <- function(func, base = "best", relativize = TRUE, avg = FALSE, avglvl =
    baseline.site = rownames(func.temp)[brow],
    comparison.site = rownames(func.temp)[j],
    delta_S = sprime - s,
    shared_S = sc,
    delta_FUNC = (sprime * zprimebar) - (s * zbar),
    RICH_LOSS = (sc - s) * zbar,
    RICH_GAIN = (sprime - sc) * zprimebar,
    @@ -208,4 +209,4 @@ price.sim <- function(nspecies, ncomms, bscaling = 0, sloss = 0, directional = F

    } ) )

    }
    }
  15. Jon Lefcheck revised this gist Nov 27, 2017. 1 changed file with 80 additions and 17 deletions.
    97 changes: 80 additions & 17 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,16 +1,16 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 18 October 2017
    #' Last updated: 27 November 2017

    #' @param func 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 standardize = TRUE, whether results are scaled by the maximum to be on the scale (-1, 1)
    #' @param relativize = TRUE, whether results are scaled by the maximum to be on the scale (-1, 1)
    #' @param avg = FALSE, whether result should be averaged across some number of baseline sites
    #' @param avglvl = 0.1, the top X% of sites in terms of functinoing that should be used as the baseline, and then results averaged
    #' @return data.frame of Price components & summed effects

    price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl = 0.1) {
    price <- function(func, base = "best", relativize = TRUE, avg = FALSE, avglvl = 0.1) {

    # Replace NAs with zeros
    func[is.na(func)] <- 0
    @@ -95,11 +95,6 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =

    } ) )

    # Calculate aggregate gains, losses, and across both gains and losses
    dat$TOTAL_LOSS <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_GAIN <- 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)

    # Remove rows where site == baseline.site
    dat <- dat[as.character(dat$comparison.site) != as.character(dat$baseline.site), ]

    @@ -113,14 +108,17 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    }

    # Relativize output so units are bounds (-1, 1) for each component
    if(standardize == TRUE) {
    if(relativize == TRUE) {

    dat[, -(1:3)] <- apply(dat[, -(1:3)], 2, function(x) x / max(abs(x), na.rm = T))

    dat[is.na(dat)] <- 0
    dat[, (5:9)] <- t(apply(dat[, (5:9)], 1, function(x) x / max(abs(x))))

    }

    # Calculate aggregate gains, losses, and across both gains and losses
    dat$TOTAL_LOSS <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_GAIN <- 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)

    }
    @@ -136,13 +134,78 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    getBaseline <- function(func, base = "best", avg = FALSE, avglvl = 0.10) {

    if(avg == TRUE & all(base == "best")) baseline <- which(rowSums(func, na.rm = T) >= max(rowSums(func, na.rm = T)) * (1 - avglvl)) else

    if(all(base == "best")) baseline <- which.max(rowSums(func, na.rm = T)) else

    if(all(base != "best") & !is.numeric(base)) baseline <- which(rownames(func) %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, ...) {

    # 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) {

    if(all(base == "best")) baseline <- which.max(rowSums(func, na.rm = T)) else
    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(all(base != "best") & !is.numeric(base)) baseline <- which(rownames(func) %in% base) else
    if(top.or.bottom == "top")

    baseline <- base
    remove <- rev(order(comparison.site))[1:combos[j, 1]]

    return(baseline)
    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)
    )

    } ) )

    } ) )

    }
  16. Jon Lefcheck revised this gist Oct 18, 2017. 1 changed file with 10 additions and 10 deletions.
    20 changes: 10 additions & 10 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,16 +1,16 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 12 October 2017
    #' Last updated: 18 October 2017

    #' @param func 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 standardize = TRUE, whether results are scaled by the maximum to be on the scale (-1, 1)
    #' @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 90%
    #' @param avg = FALSE, whether result should be averaged across some number of baseline sites
    #' @param avglvl = 0.1, the top X% of sites in terms of functinoing that should be used as the baseline, and then results averaged
    #' @return data.frame of Price components & summed effects

    price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl = 0.90) {
    price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl = 0.1) {

    # Replace NAs with zeros
    func[is.na(func)] <- 0
    @@ -130,19 +130,19 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    #' @param func 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 90%
    #' @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(func, base = "best", avg = FALSE, avglvl = 0.90) {
    getBaseline <- function(func, base = "best", avg = FALSE, avglvl = 0.10) {

    if(avg == TRUE & all(base == "best"))

    baseline <- which(rowSums(func, na.rm = T) >= max(rowSums(func, na.rm = T)) * avglvl) else
    if(avg == TRUE & all(base == "best")) baseline <- which(rowSums(func, na.rm = T) >= max(rowSums(func, na.rm = T)) * (1 - avglvl)) else

    if(all(base == "best")) baseline <- which.max(rowSums(func, na.rm = T)) else

    if(all(base != "best") & !is.numeric(base)) baseline <- which(rownames(func) %in% base) else

    if(all(base != "best") & is.numeric(base)) baseline <- base
    baseline <- base

    return(baseline)

    }
  17. Jon Lefcheck revised this gist Oct 12, 2017. 1 changed file with 4 additions and 4 deletions.
    8 changes: 4 additions & 4 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,7 +1,7 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 0p October 2017
    #' Last updated: 12 October 2017

    #' @param func 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
    @@ -78,7 +78,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    dat <- data.frame(
    baseline.site = rownames(func.temp)[brow],
    comparison.site = rownames(func.temp)[j],
    delta_S = s - sprime,
    delta_S = sprime - s,
    delta_FUNC = (sprime * zprimebar) - (s * zbar),
    RICH_LOSS = (sc - s) * zbar,
    RICH_GAIN = (sprime - sc) * zprimebar,
    @@ -111,7 +111,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    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) {

    @@ -120,7 +120,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    dat[is.na(dat)] <- 0

    }

    return(dat)

    }
  18. Jon Lefcheck revised this gist Oct 9, 2017. No changes.
  19. Jon Lefcheck revised this gist Oct 9, 2017. 1 changed file with 67 additions and 17 deletions.
    84 changes: 67 additions & 17 deletions README.md
    View file
    Original file line number Diff line number Diff line change
    @@ -1,26 +1,76 @@
    ## 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.
    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
    ```
    # 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)
    # 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, standardize = 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, standardize = 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, standardize = 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, standardize = 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)
    ))
    # In the above examples, the value should load exclusively on only one component of the Price equation
    price(func, base = "baseline")
    price(ex5, standardize = 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)
    ))
    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
    price(ex6, standardize = F)
    ```
  20. Jon Lefcheck revised this gist Oct 9, 2017. 1 changed file with 3 additions and 3 deletions.
    6 changes: 3 additions & 3 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -44,8 +44,8 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    # Get total number of species from baseline site
    s <- sum(func.temp[brow, ] > 0, na.rm = T)

    # Calculate mean level of functioning at baseline site
    zbar <- sum(func.temp[brow, ], na.rm = T) / sum(func.temp[brow, ] > 0, na.rm = T)
    # Calculate average contribution to functioning at the baseline site
    zbar <- sum(func.temp[brow, ], na.rm = T) / s

    # Next, loop over comparisons
    dat <- do.call(rbind, lapply((1:nrow(func.temp))[-brow], function(j) {
    @@ -60,7 +60,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    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)
    zprimebar <- sum(func.temp[j, ], na.rm = T) / sprime

    # Calculate mean level of functioning among shared species at baseline site
    zcbar <- sum(as.numeric(shared * func.temp[brow, ]), na.rm = T) /
  21. Jon Lefcheck revised this gist Oct 9, 2017. 1 changed file with 3 additions and 3 deletions.
    6 changes: 3 additions & 3 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,7 +1,7 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 06 October 2017
    #' Last updated: 0p October 2017

    #' @param func 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
    @@ -106,9 +106,9 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    # If avg == TRUE, average across sites
    if(avg == TRUE) {

    dat <- cbind(baseline.site = dat[, 1], aggregate(. ~ comparison.site, "mean", data = dat[, -(1)]))
    dat <- aggregate(. ~ comparison.site, "mean", data = dat[, -(1)])

    if(length(baseline) > 1) dat$baseline.site = paste("Average of", length(baseline), "sites")
    dat <- cbind(baseline.site = paste("Average of", length(baseline), "sites"), dat)

    }

  22. Jon Lefcheck revised this gist Oct 6, 2017. 1 changed file with 4 additions and 4 deletions.
    8 changes: 4 additions & 4 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -30,7 +30,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    dat <- do.call(rbind, lapply(baseline, function(i) {

    # Remove other baseline sets from the data
    if(avg == TRUE | length(baseline) > 1) {
    if(length(baseline) > 1) {

    remove <- baseline[!baseline %in% i]

    @@ -76,8 +76,8 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =

    # Calculate Price equation components
    dat <- data.frame(
    site = rownames(func.temp)[j],
    baseline.site = rownames(func.temp)[brow],
    comparison.site = rownames(func.temp)[j],
    delta_S = s - sprime,
    delta_FUNC = (sprime * zprimebar) - (s * zbar),
    RICH_LOSS = (sc - s) * zbar,
    @@ -101,12 +101,12 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    dat$TOTAL_DIV <- rowSums(dat[, c("RICH_LOSS", "RICH_GAIN", "COMP_LOSS", "COMP_GAIN")], na.rm = T)

    # Remove rows where site == baseline.site
    dat <- dat[as.character(dat$site) != as.character(dat$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(. ~ site, "mean", data = dat)
    dat <- cbind(baseline.site = dat[, 1], aggregate(. ~ comparison.site, "mean", data = dat[, -(1)]))

    if(length(baseline) > 1) dat$baseline.site = paste("Average of", length(baseline), "sites")

  23. Jon Lefcheck revised this gist Oct 6, 2017. 1 changed file with 5 additions and 5 deletions.
    10 changes: 5 additions & 5 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -30,7 +30,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    dat <- do.call(rbind, lapply(baseline, function(i) {

    # Remove other baseline sets from the data
    if(avg == TRUE & length(baseline) > 1) {
    if(avg == TRUE | length(baseline) > 1) {

    remove <- baseline[!baseline %in% i]

    @@ -99,7 +99,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    dat$TOTAL_LOSS <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_GAIN <- 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)

    # Remove rows where site == baseline.site
    dat <- dat[as.character(dat$site) != as.character(dat$baseline.site), ]

    @@ -108,10 +108,10 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =

    dat <- aggregate(. ~ site, "mean", data = dat)

    if(length(baseline) > 1) dat$baseline.site = ">1"
    if(length(baseline) > 1) dat$baseline.site = paste("Average of", length(baseline), "sites")

    }

    # Relativize output so units are bounds (-1, 1) for each component
    if(standardize == TRUE) {

    @@ -120,7 +120,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    dat[is.na(dat)] <- 0

    }

    return(dat)

    }
  24. Jon Lefcheck revised this gist Oct 6, 2017. 1 changed file with 3 additions and 3 deletions.
    6 changes: 3 additions & 3 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -104,7 +104,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    dat <- dat[as.character(dat$site) != as.character(dat$baseline.site), ]

    # If avg == TRUE, average across sites
    if(all(base != "base") | length(base) > 1 | avg == TRUE) {
    if(avg == TRUE) {

    dat <- aggregate(. ~ site, "mean", data = dat)

    @@ -117,7 +117,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =

    dat[, -(1:3)] <- apply(dat[, -(1:3)], 2, function(x) x / max(abs(x), na.rm = T))

    dat[is.na(dat)] = 0
    dat[is.na(dat)] <- 0

    }

    @@ -145,4 +145,4 @@ getBaseline <- function(func, base = "best", avg = FALSE, avglvl = 0.90) {

    if(all(base != "best") & is.numeric(base)) baseline <- base

    }
    }
  25. Jon Lefcheck revised this gist Oct 6, 2017. 1 changed file with 33 additions and 31 deletions.
    64 changes: 33 additions & 31 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,10 +1,10 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 04 October 2017
    #' Last updated: 06 October 2017

    #' @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 base the row name(s) or number(s) of the baseline community; defaults to highest productivity community
    #' @param standardize = TRUE, whether results are scaled by the maximum to be on the scale (-1, 1)
    #' @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 90%
    @@ -24,7 +24,7 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    if(is.null(rownames(func))) rownames(func) <- 1:nrow(func)

    # Identify baseline sites for comparison
    baseline = get.baseline(func, base = base, avg = avg, avglvl = avglvl)
    baseline <- getBaseline(func, base = base, avg = avg, avglvl = avglvl)

    # Loop over baselines
    dat <- do.call(rbind, lapply(baseline, function(i) {
    @@ -36,25 +36,25 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =

    func.temp <- func[-remove, ]

    } else func.temp = func
    } else func.temp <- func

    # Redefine baseline row
    brow <- which(rownames(func.temp) == rownames(func)[i])

    # Get total number of species from baseline site
    s <- sum(func.temp[i, ] > 0, na.rm = T)
    s <- sum(func.temp[brow, ] > 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)
    zbar <- sum(func.temp[brow, ], na.rm = T) / sum(func.temp[brow, ] > 0, na.rm = T)

    # Next, loop over comparisons
    dat <- do.call(rbind, lapply((1:nrow(func.temp))[-new.baseline], function(j) {
    dat <- do.call(rbind, lapply((1:nrow(func.temp))[-brow], 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))
    shared <- apply(func.temp[c(brow, j), ], 2, function(x) ifelse(any(x == 0), 0, 1))

    # Sum number of shared species in both sites
    sc <- sum(shared, na.rm = T)
    @@ -63,8 +63,8 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    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)
    zcbar <- sum(as.numeric(shared * func.temp[brow, ]), na.rm = T) /
    sum(shared * func.temp[brow, ] > 0, na.rm = T)

    if(is.na(zcbar)) zcbar <- 0

    @@ -75,9 +75,9 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    if(is.na(zcprimebar)) zcprimebar <- 0

    # Calculate Price equation components
    datf <- data.frame(
    dat <- data.frame(
    site = rownames(func.temp)[j],
    baseline.site = rownames(func.temp)[i],
    baseline.site = rownames(func.temp)[brow],
    delta_S = s - sprime,
    delta_FUNC = (sprime * zprimebar) - (s * zbar),
    RICH_LOSS = (sc - s) * zbar,
    @@ -87,28 +87,28 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    CDE = sc * (zcprimebar - zcbar)
    )

    return(datf)
    return(dat)

    } ) )

    # Calculate aggregate gains, losses, and across both gains and losses
    dat$TOTAL_LOSS <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_GAIN <- 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)

    } ) )

    # Calculate aggregate gains, losses, and across both gains and losses
    dat$TOTAL_LOSS <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_GAIN <- 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)

    # 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) {
    if(all(base != "base") | length(base) > 1 | avg == TRUE) {

    dat <- aggregate(. ~ site + baseline.site, "mean", data = dat)
    dat <- aggregate(. ~ site, "mean", data = dat)

    if(length(baseline) > 1) dat$baseline.site = NA
    if(length(baseline) > 1) dat$baseline.site = ">1"

    }

    @@ -125,22 +125,24 @@ price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =

    }

    #' `getBaseline`

    #' @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 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 90%
    #' @return integer corresponding to the row number of the baseline community from the site-by-species functioning matrix

    get.baseline <- function(func, base = "best", avg = FALSE, avglvl = 0.90) {
    getBaseline <- function(func, base = "best", avg = FALSE, avglvl = 0.90) {

    if(avg == TRUE & base == "best")
    if(avg == TRUE & all(base == "best"))

    baseline <- which(rowSums(func, na.rm = T) >= max(rowSums(func, na.rm = T)) * avglvl) else

    if(base == "best") baseline <- which.max(rowSums(func, na.rm = T)) else
    if(all(base == "best")) baseline <- which.max(rowSums(func, na.rm = T)) else

    if(base != "best" & !is.numeric(base)) baseline <- which(rownames(func) == base) else
    if(all(base != "best") & !is.numeric(base)) baseline <- which(rownames(func) %in% base) else

    if(base != "best" & is.numeric(base)) baseline <- base
    if(all(base != "best") & is.numeric(base)) baseline <- base

    }
    }
  26. Jon Lefcheck revised this gist Oct 4, 2017. 1 changed file with 55 additions and 55 deletions.
    110 changes: 55 additions & 55 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,101 +1,81 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 07 October 2016
    #' Last updated: 04 October 2017

    #' @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 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 90%
    #' @return integer corresponding to the row number of the baseline community from the site-by-species functioning matrix

    get.baseline = function(func, base = "best", avg = FALSE, avglvl = 0.90) {

    if(avg == TRUE & base == "best")

    baseline = which(rowSums(func, na.rm = T) >= max(rowSums(func, na.rm = T)) * avglvl) else

    if(base == "best") baseline = which.max(rowSums(func, na.rm = T)) else

    if(base != "best" & !is.numeric(base)) baseline = which(rownames(func) == base) else

    if(base != "best" & is.numeric(base)) baseline = base

    }

    #' @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%
    #' @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 be on the scale (-1, 1)
    #' @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 90%
    #' @return data.frame of Price components & summed effects

    price = function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl = 0.90) {
    price <- function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl = 0.90) {

    # Replace NAs with zeros
    func[is.na(func)] = 0

    func[is.na(func)] <- 0
    if(!all(apply(func, 2, is.numeric))) stop("Remove non-numeric columns from matrix")

    # Remove species that do not contribute to functioning at all
    func = func[, colSums(func) != 0]
    func <- func[, colSums(func) != 0]

    # Assign rownames if not already present
    if(is.null(rownames(func))) rownames(func) = 1:nrow(func)
    if(is.null(rownames(func))) rownames(func) <- 1:nrow(func)

    # 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) {
    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]
    remove <- baseline[!baseline %in% i]

    func.temp = func[-remove, ]
    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)
    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)
    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)
    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) {
    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)
    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))
    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)
    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)
    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) /
    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
    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) /
    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
    if(is.na(zcprimebar)) zcprimebar <- 0

    # Calculate Price equation components
    datf = data.frame(
    datf <- data.frame(
    site = rownames(func.temp)[j],
    baseline.site = rownames(func.temp)[i],
    delta_S = s - sprime,
    @@ -112,21 +92,21 @@ price = function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    } ) )

    # Calculate aggregate gains, losses, and across both gains and losses
    dat$TOTAL_LOSS =rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_GAIN = 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)
    dat$TOTAL_LOSS <- rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_GAIN <- 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)

    } ) )

    # Remove rows where site == baseline.site
    dat = dat[as.character(dat$site) != as.character(dat$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)
    dat <- aggregate(. ~ site + baseline.site, "mean", data = dat)

    if(length(baseline) > 1) dat$baseline.site = NA

    @@ -135,12 +115,32 @@ price = function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    # 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))
    dat[, -(1:3)] <- apply(dat[, -(1:3)], 2, function(x) x / max(abs(x), na.rm = T))

    dat[is.na(dat)] = 0

    }

    return(dat)

    }
    }

    #' @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 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 90%
    #' @return integer corresponding to the row number of the baseline community from the site-by-species functioning matrix

    get.baseline <- function(func, base = "best", avg = FALSE, avglvl = 0.90) {

    if(avg == TRUE & base == "best")

    baseline <- which(rowSums(func, na.rm = T) >= max(rowSums(func, na.rm = T)) * avglvl) else

    if(base == "best") baseline <- which.max(rowSums(func, na.rm = T)) else

    if(base != "best" & !is.numeric(base)) baseline <- which(rownames(func) == base) else

    if(base != "best" & is.numeric(base)) baseline <- base

    }
  27. Jon Lefcheck revised this gist Oct 7, 2016. 1 changed file with 15 additions and 13 deletions.
    28 changes: 15 additions & 13 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,11 +1,11 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 04 October 2016
    #' Last updated: 07 October 2016

    #' @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 avg = TRUE, whether
    #' @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 90%
    #' @return integer corresponding to the row number of the baseline community from the site-by-species functioning matrix

    @@ -26,14 +26,16 @@ get.baseline = function(func, base = "best", avg = FALSE, avglvl = 0.90) {
    #' @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 result should be averaged across some number of baseline sites
    #' @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

    if(!all(apply(func, 2, is.numeric))) stop("Remove non-numeric columns from matrix")

    # Remove species that do not contribute to functioning at all
    func = func[, colSums(func) != 0]
    @@ -96,23 +98,23 @@ price = function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    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),
    CDE = sc * (zcprimebar - zcbar),
    T_Tprime = (sprime * zprimebar) - (s * zbar)
    delta_S = s - sprime,
    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(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)
    dat$TOTAL_LOSS =rowSums(dat[, c("RICH_LOSS", "COMP_LOSS")], na.rm = T)
    dat$TOTAL_GAIN = 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)

  28. Jon Lefcheck revised this gist Oct 5, 2016. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -26,7 +26,7 @@ get.baseline = function(func, base = "best", avg = FALSE, avglvl = 0.90) {
    #' @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 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 90%
    #' @return data.frame of Price components & summed effects

  29. Jon Lefcheck revised this gist Oct 5, 2016. 1 changed file with 86 additions and 77 deletions.
    163 changes: 86 additions & 77 deletions price.R
    View file
    Original file line number Diff line number Diff line change
    @@ -1,7 +1,7 @@
    #' `price` an R function for translating the extended Price equation to ecosystem functioning

    #' Author: Jon Lefcheck
    #' Last updated: 06 Apr 2016
    #' Last updated: 04 October 2016

    #' @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
    @@ -38,98 +38,107 @@ price = function(func, base = "best", standardize = TRUE, avg = FALSE, avglvl =
    # Remove species that do not contribute to functioning at all
    func = func[, colSums(func) != 0]

    # Assign rownames if not already present
    if(is.null(rownames(func))) rownames(func) = 1:nrow(func)

    # 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) {
    # 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 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
    remove = baseline[!baseline %in% i]

    # Get total number of species from baseline site
    s = sum(func.temp[i, ] > 0, na.rm = T)
    func.temp = func[-remove, ]

    # Calculate mean level of functioning at baseline site
    zbar = sum(func.temp[i, ], na.rm = T) / sum(func.temp[i, ] > 0, na.rm = T)
    } 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) {

    # Determine new baseline
    new.baseline = get.baseline(func.temp, base = base, avg = avg, avglvl = avglvl)
    # Get total number of species from the comparison site
    sprime = sum(func.temp[j, ] > 0, na.rm = T)

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

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

    # 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)
    # Sum number of shared species in both sites
    sc = sum(shared, na.rm = T)

    return(dat)
    # Calculate mean level of functioning at comparison site
    zprimebar = sum(func.temp[j, ], na.rm = T) / sum(func.temp[j, ] > 0, na.rm = T)

    } ) )

    # 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) {
    # 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

    dat = aggregate(. ~ site + baseline.site, "mean", data = dat)
    # 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(length(baseline) > 1) dat$baseline.site = NA
    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),
    CDE = sc * (zcprimebar - zcbar),
    T_Tprime = (sprime * zprimebar) - (s * zbar)
    )

    return(datf)

    } ) )

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

    dat[is.na(dat)] = 0

    }

    return(dat)

    }
  30. Jon Lefcheck revised this gist Apr 6, 2016. 1 changed file with 7 additions and 0 deletions.
    7 changes: 7 additions & 0 deletions README.md
    View file
    Original file line number Diff line number Diff line change
    @@ -16,4 +16,11 @@ func = t(data.frame(
    # 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
    ```