R/LatCRCurve.R

In LewisAJones/LBGSim: Latitudinal biodiversity gradient simulator

#' Latitudinal classical rarefaction curve
#'
#' @param data Dataframe (colnames: rep, species, x, y). A dataframe containing species ID, occurrences and (optional) replications.
#' @param bandsize Integer. Size of latitudinal bins/bands desired for analyses. Defaults to 15. 
#' @param knots Integer. Defaults to 100.
#' @param iterations Integer. Defaults to 100, number of iterations to sample each rep for.
#' @param write Logical. Write dataframe of latitudinal diversity.
#' @param name Character. File name to write. 
#' @return A dataframe containing data for rarefaction curve. 
#' @importFrom matrixStats rowMeans2 rowSds rowSums2
#' @importFrom dplyr bind_rows

#' @export
LatCRCurve <- function(data, bandsize = 15, knots = 100, iterations = 100, write = FALSE, name = "simulation"){
  nbins <- 180/bandsize
  bin <- as.data.frame(seq(1, nbins, 1))
  max <- as.data.frame(rev(seq((-90+bandsize), 90, bandsize))) 
  min <- as.data.frame(rev(seq(-90, (90-bandsize), bandsize)))
  mid <- as.data.frame(rev(seq((-90+(bandsize/2)), (90-(bandsize/2)), bandsize)))
  hemisphere <- c("N", "S")
  hemisphere <- as.data.frame(rep(hemisphere, each = nbins/2))
  latbin <- cbind(bin, max, min, mid, hemisphere)
  colnames(latbin) <- c("bin", "max", "min", "mid", "hemisphere")
  binned <- data.frame()
  maxDiv <- 0
  for(i in 1:length(latbin$bin)){
      tmp <- subset(data, y <= latbin[i, "max"] & y >= latbin[i, "min"])
      if(nrow(tmp) == 0){next}
      if(length(unique(tmp$species)) > maxDiv){
        maxDiv <- as.numeric(length(unique(tmp$species)))
      }}
    k <- ceiling(maxDiv/knots)
    k <- seq(from = 1, to = maxDiv*1.25, by = k)
    tmpdat <- lapply(k, function (j){
      iter <- lapply(1:iterations, function (x){
        species_per_latbin <- vector("numeric") #count occurrences within each latitudinal bin
        for (i in 1:nrow(latbin)) {
          out <- subset(data, y <= latbin[i, "max"] & y >= latbin[i, "min"])
          if(j > nrow(out)){
            species_per_latbin[i] <- NA
          }
          else{
            out <- out[sample(nrow(out), size = j, replace = FALSE),]
            species_per_latbin[i] <- length(unique(out$species))
          }
        }
        species_per_latbin
      })
      names(iter) <- 1:iterations
      iter <- do.call(cbind.data.frame, iter)
      iter
    })
    names(tmpdat) <- k
    tmpdat <- dplyr::bind_rows(tmpdat)
    tmpdat <- data.matrix(tmpdat)
    mean <- matrixStats::rowMeans2(tmpdat, na.rm = TRUE)
    sd <- matrixStats::rowSds(tmpdat, na.rm = TRUE)
    n <- matrixStats::rowSums2(!is.na(tmpdat))
    error <- qnorm(0.975)*sd/sqrt(n)
    CI.Lower <- mean - error
    CI.Upper <- mean + error
    knots <- rep(k, each = nrow(bin))
    tmplat <- latbin[rep(seq_len(nrow(latbin)), length(k)), ]
    master <- cbind.data.frame(tmplat, knots, mean, sd, n, error, CI.Lower, CI.Upper)
    rownames(master) <- c()
    master[is.na(master)] <-NA
    tmpdat <- master
    col <- viridisLite::viridis(nrow(bin))
    tmpdat$col <- rep(col, length(k))
    tmpdat <- subset(tmpdat, !mean == 0)
    plot(mean~knots, tmpdat, type = "l", col = NA, main = "Rarefaction curve", ylab = "Species richness", xlab = "Number of occurrences", pch = 20, ylim = c(0, maxDiv))
    for(x in 1:nrow(bin)){
      tmp <- subset(tmpdat, bin == x)
      #points(mean~knots, tmp, pch = 20)
      #arrows(tmp$knots, tmp$CI.Lower, tmp$knots, tmp$CI.Upper, length=0.05, angle=90, code=3)
      suppressWarnings(polygon(c(tmp$knots,rev(tmp$knots)),c(tmp$CI.Lower,rev(tmp$CI.Upper)), col = adjustcolor(tmp$col, alpha.f=0.5), border = FALSE))
      lines(mean~knots, tmp, pch = 20, col = tmp$col)
    }
    legend("topleft", title = "Latitudinal midpoint", legend=unique(tmpdat$mid),
           fill=unique(tmpdat$col), cex=0.8, bg = NA,
           box.lty=0)
    if(write == TRUE){
      suppressWarnings(dir.create("./CR/"))
      suppressWarnings(dir.create("./CR/Figures/"))
      png(paste("./CR/Figures/", name, "_CR_curve_LBG_plot.png", sep = ""), width = 2400, height = 2000, res = 300)
      plot(mean~knots, tmpdat, type = "l", col = NA, main = "Rarefaction curve", ylab = "Species richness", xlab = "Number of occurrences", pch = 20, ylim = c(0, maxDiv))
      for(x in 1:nrow(bin)){
        tmp <- subset(tmpdat, bin == x)
        #points(mean~knots, tmp, pch = 20)
        #arrows(tmp$knots, tmp$CI.Lower, tmp$knots, tmp$CI.Upper, length=0.05, angle=90, code=3)
        suppressWarnings(polygon(c(tmp$knots,rev(tmp$knots)),c(tmp$CI.Lower,rev(tmp$CI.Upper)), col = adjustcolor(tmp$col, alpha.f=0.5), border = FALSE))
        lines(mean~knots, tmp, pch = 20, col = tmp$col)
      }
      legend("topleft", inset = c(0.02, 0.02), title = "Latitudinal midpoint", legend=unique(tmpdat$mid),
             fill=unique(tmpdat$col), cex=0.8, bg = NA,
             box.lty=0)
      dev.off()
      write.csv(master, paste("./CR/", name, "_rarefaction_curve.csv", sep = ""))
    }
    
  return(master)
}