R/RSFIA_analysis.R

In bmcnellis/RSFIA: FIA Tree Mortality Analysis

#' @title Secondary analysis functions for RSFIA
#' @description
#'
#' Additional analysis might be needed on the main mortality dataset, this
#' set of functions aids that.
#'
#' @examples RSFIA_analysis()
#' @export
RSFIA_analysis <- function() {
  invisible()
}
#' @describeIn RSFIA_analysis Bins mortality rates
#' @export
CategoricalMort <- function(mort_vec,
                            breaks = c(0.2, 0.5),
                            names = c('Low', 'Moderate', 'High')) {

  if (!is.vector(mort_vec)) stop('Mort vec needs to be 1d')
  if (!is.numeric(mort_vec)) stop('Mort vec needs to be numeric')
  if (!(length(breaks) + 1 == length(names))) names <- rep('', length(breaks) + 1)

  for (i in 1:(length(names) - 1)) {
    tag0 <- paste0(breaks[i - 1] * 100, '% <')
    tag1 <- ifelse(i == 1, '0% <', tag0)
    tag2 <- paste0('< ', breaks[i] * 100, '%')
    names[i] <- paste(tag1, names[i], tag2)
  }
  tag1 <- paste0(breaks[length(breaks)] * 100, '% <')
  tag2 <- '< 100%'
  names[length(names)] <- paste(tag1, names[length(names)], tag2)

  out_vec <- rep(names[1], length(mort_vec))
  for (i in 2:length(names)) {
    out_vec <- ifelse(mort_vec > breaks[i - 1], names[i], out_vec)
  }
  return(out_vec)
}
#' @describeIn RSFIA_analysis Generates a summary table for binned mortality
#' @export
CalcCategoricalMortByBreaks <- function(x, group, breaks) {
  bins <- 1:length(breaks)
  y <- rep(bins[1], length(x))
  for (i in 2:length(bins)) {
    y[which(x > breaks[i])] <- bins[i]
  }
  z <- table(y)
  za <- aggregate(x, by = list(group, y), FUN = length)
  return(list(z, za))
}
#' @describeIn RSFIA_analysis Generates a summary table OR summary plot of zero-mort plots
#' @export
SummarizeAndPlotZeroMortPlotsBySection <- function(agg = F) {

  # Prepare plot environment:
  df <- PrepPlotEnvir()
  ytag <- 'Number of Plots With No Observed Mortality'
  mort_df <- FIA_mortality_with_explanatory
  mortality <- mort_df$mort_rate
  df <- data.frame(mortality, mort_df$Cleland_section, mort_df$Cleland_province,
                   stringsAsFactors = F)
  df[, 2] <- KeyClelandCode(df[, 2], lvl = 'section')
  df[, 3] <- KeyClelandCode(df[, 3], lvl = 'province')
  NA_cols <- aggregate(df[, 1], by = list(df[, 2]),
                       FUN = function(x) all(is.na(x)))
  NA_cols <- NA_cols$Group.1[which(NA_cols[, 'x'])]
  df[which(df[, 2] %in% NA_cols), 1] <- 0
  x1 <- aggregate(mortality, by = list(df[, 2]), FUN = function(x) {
    ret <- ifelse('0' %in% names(table(x)), table(x)[['0']], 0)
    return(ret)
  })
  x2 <- aggregate(mortality, by = list(df[, 2]), FUN = function(x) {
    ret <- ifelse('0' %in% names(table(x)), table(x)[['0']], 0)
    ret0 <- round(ret / sum(table(x), na.rm = T), 3) * 100
    return(ret0)
  })
  out <- dplyr::full_join(x1, x2, by = 'Group.1')
  colnames(out) <- c('Cle_sect', 'n_zero', 'perc_zero')
  if (agg) {
    return(out)
  }
  for (i in 1:nrow(out)) {
    out[i, 'Cle_prov'] <- df[which(df[, 2] == out$Cle_sect[i])[1], 3]
  }
  out_plot <- ggplot() + theme_bw() +
    geom_boxplot(aes(x = out[, 4], y = out[, 3])) +
    theme(axis.text.x = element_text(angle = -45, hjust = 0)) +
    theme(plot.margin = unit(c(1, 2.5, 1, 1), 'cm')) +
    xlab('') + ylab(ytag)

  # Return:
  return(out_plot)
}
#' @describeIn RSFIA_analysis Assess relationships between 10-year weather and FIA mort
#' @export
AssessYearBinWeatherMort <- function(source, bin = '10', model = 'lm') {
  # Sanity checks:
  stopifnot(bin %in% c('2', '4', '6', '8', '10')) # available bins
  stopifnot(model %in% c('lm', 'two_stage', 'tweedie', 'tobit'))

  # Prepare data:
  mort <- FIA_mortality_with_explanatory
  mort_cols0 <- c('mort_rate', 'non_harv_fire_beet_mort_rate',
                 'Fire_rate_wrt_ntree', 'Insect_rate_wrt_ntree')
  mort_cols <- c('LAT', 'LON', 'Cleland_province', 'Cleland_section', mort_cols0)
  mort <- mort[, which(colnames(mort) %in% mort_cols)]
  colnames(mort)[which(colnames(mort) == source)] <- 'mortality'
  wth_10 <- FIA_weather_by_year_bin
  wth_10 <- wth_10[, c(1, 2, grep(colnames(wth_10), pattern = bin))]
  data <- dplyr::left_join(mort, wth_10, by = c('LAT', 'LON'))
  data$Cleland_section <- RSFIA::KeyClelandCode(data$Cleland_section, lvl = 'section')

  # Declare:
  out_df <- data.frame(unique(data$Cleland_section), stringsAsFactors = F)
  out_list <- list(out_df, out_df, out_df)
  names(out_list) <- c('p_val', 'shapiro_p_val', 'r_2')
  for (i in 1:nrow(out_df)) { # for each Cleland section...
    row_indx <- which(data$Cleland_section == out_df[i, 1])
    y <- data[row_indx, 'mortality']
    bin_y <- ifelse(y == 0, 0, 1)
    cnt <- 0
    for (j in c('mean_max_min', 'min_temp_by_Q_val', 'max_temp_by_Q_val')) { #
      cnt <- cnt + 1
      x <- data[row_indx, grep(j, colnames(data))]
      if (model == 'lm') {
        mod <- lm(y ~ x)
        out_list[[1]][, j] <- round(summary(mod)$coefficients[2, 4], 2)
        out_list[[2]][, j] <- round(shapiro.test(mod$residuals)$p.value, 2)
        out_list[[3]][, j] <- round(summary(mod)$r.squared, 2)
      }
      if (model == 'two_stage') {
        zinf <- pscl::zeroinfl(bin_y ~ x)
        browser()
      }

      browser()
    }
  }
  return(out_list)
}