R/FreqPlot.R

In LewisAJones/LBGSim: Latitudinal biodiversity gradient simulator

#' Generates occurrence frequency distribution plot and dataframe
#' @param data Data.frame (colnames:species, x, y). A dataframe containing species ID and occurrences (x, y).
#' @param breaks Numeric. Number of breaks desired in plot.
#' @param reps Logical. Whether the simulation was run for more than one replication. 95 percent confidence intervals are generated by default.
#' @param write Logical. Write outputs to disk. Default is FALSE.
#' @param name Character. Name of folder to save outputs in if write = TRUE. Default is "simulation"
#' @importFrom dplyr group_by summarise bind_rows
#' @importFrom DescTools MeanCI
#' @export
FreqPlot <- function(data, breaks = 10, reps = FALSE, write = FALSE, name = "simulation"){
  if(reps == TRUE){
    nreps <- unique(data$rep)
    dat <- lapply(nreps, function(i){ #check is right set up
      subdat <- subset(data, rep == i)
      dat <- data.frame(subdat %>% dplyr::group_by(species) %>% dplyr::summarise(count=n()))
      dat})
    names(dat) <- nreps
    dat <- dplyr::bind_rows(dat)
    max <- max(dat$count, na.rm = TRUE)
    brk <- seq(from = 0, to = (max+breaks), by = breaks)
    tmp <- lapply(nreps, function(i){ #check is right set up
      subdat <- subset(data, rep == i)
      tmp <- subdat %>% dplyr::group_by(species) %>% dplyr::summarise(count=n())
      tmp <- hist(tmp$count, plot = FALSE, breaks = brk)
      tmp <- tmp$counts
      tmp
      })
    tmp <- do.call(rbind, tmp)
    rep.means <- lapply(1:ncol(tmp), function (x){
      suppressWarnings(means <- DescTools::MeanCI(tmp[,x], conf.level=0.95, na.rm = TRUE))
      means})
    tmp <- do.call(rbind, rep.means) 
    colnames(tmp) <- c("mean", "CI.Lower", "CI.Upper")
    bin <- seq(0.5, nrow(tmp), 1)
    tmp <- cbind.data.frame(tmp, bin)
    if(write == TRUE){
      suppressWarnings(dir.create("./Figures/"))
      png(paste("./Figures/", name, "_frequency_plot.png", sep = ""), width = 2400, height = 2000, res = 300)
      barplot(tmp$mean, space = 0, main = "Occurrence frequency distribution", axes = TRUE, col = "grey80", ylab = expression(bold(Frequency)), ylim = c(0, ceiling(max(tmp$CI.Upper, na.rm = TRUE)*1.15)), xlab = expression(bold(Occurrences~(italic("n")))))
      axis(side = 1, at = seq(0, length(tmp$mean), 1), labels = brk[(seq(1, length(tmp$mean)+1, 1))])
      suppressWarnings(segments(x0 = tmp$bin, y0 = tmp$CI.Lower, x1 = tmp$bin, y1 = tmp$CI.Upper))
      suppressWarnings(arrows(x0 = tmp$bin, y0 = tmp$CI.Lower, x1 = tmp$bin, y1 = tmp$CI.Upper, angle = 90, code = 3, length = 0.05))
      dev.off()
      }
    barplot(tmp$mean, space = 0, main = "Occurrence frequency distribution", axes = TRUE, col = "grey80", ylab = expression(bold(Frequency)), ylim = c(0, ceiling(max(tmp$CI.Upper, na.rm = TRUE)*1.15)), xlab = expression(bold(Occurrences~(italic("n")))))
    axis(side = 1, at = seq(0, length(tmp$mean), 1), labels = brk[(seq(1, length(tmp$mean)+1, 1))])
    suppressWarnings(segments(x0 = tmp$bin, y0 = tmp$CI.Lower, x1 = tmp$bin, y1 = tmp$CI.Upper))
    suppressWarnings(arrows(x0 = tmp$bin, y0 = tmp$CI.Lower, x1 = tmp$bin, y1 = tmp$CI.Upper, angle = 90, code = 3, length = 0.05))
    tmp$bin <- seq(breaks, (breaks*(nrow(tmp))+1), breaks)
  }
  else{tmp <- data %>% group_by(species) %>% summarise(count=n())
  max <- max(tmp$count, na.rm = TRUE)
  brk <- seq(from = 0, to = (max+breaks), by = breaks)
  if(write == TRUE){
    suppressWarnings(dir.create("./Figures/"))
    png(paste("./Figures/", name, "_frequency_plot.png", sep = ""), width = 2400, height = 2000, res = 300)
    hist(tmp$count, main = "Occurrence frequency distribution", breaks = brk, col = "grey80", freq = TRUE, ylab = expression(bold(Frequency)), xlab = expression(bold(Occurrences~(italic("n")))))
    dev.off()
  }
  p <- hist(tmp$count, main = "Occurrence frequency distribution", breaks = brk, col = "grey80", freq = TRUE, ylab = expression(bold(Frequency)), xlab = expression(bold(Occurrences~(italic("n")))))
  counts <- p$counts
  int <- seq(breaks, (breaks*length(counts))+1, breaks)
  tmp <- cbind.data.frame(counts, int)
  colnames(tmp) <- c("counts", "bin")
  }
  if(write == TRUE){
    write.csv(tmp, paste("./Figures/", name, "_frequency.csv", sep = ""))
  }
  return(tmp)
}