R/SawlogLiklihood.R

In ryanmismith/inventoryfunctions: Forestry Inventory Functions

#' Likelihood of Having A Sawlog
#'
#' For explanation of form/risk analysis see Castle. et al 2017 . This function
#' is meant to supplement the Form.Risk merchandising function.
#'
#'
#'@details
#'Will work for AB, BT, PB, QA, RM, RO, SM, WA, WO, YB, EH, OC, RS, WS, WP
#'###
#'These values can be used in a multitude of ways, but if you are looking for binary TRUE/FALSE results
#'for if a tree is likely to produce a sawlog, the optimal cutoff point for determining whether or not a tree
#'has a sawlog can be found using the coords function
#'in the pRoc package. Cutpoints that have been used with exceptionally high AUC values (greater than .9)
#'in Northern Maine.
#'
#'@param SPP The species identification using FVS codes: ex 'RO' = Red Oak
#'@param DBH Diameter at breast height in cm
#'@param HT Height in m
#'@param Form Form classes 1-8 or 'GF', 'AF', 'PF' values are accepted.
#'@param Risk Risk class may be entered using 1-4 values or 'HR' or 'LR'.
#'@family Merchandising Functions
#'@return
#'Returns a predictive value to be used in AUC analysis to determine the likelihood that a tree contains a sawlog.
#'
#'@author Ryan Smith
#'
#'@examples
#'Sawlog.Likelihood("RO", 42, 20, 1, 1)
#'@export


Sawlog.Likelihood <- function(SPP, DBH, HT, Form, Risk){
  if(SPP %in% c("AB", "BA", "BT", "PB", "QA", "RM", "RO", "SM", "WA", "WO", "YB", "OH")){
  # Intercept and DBH Coefficients
  I <- -77.78576
  blogDBH <- 24.51994
  bHT <- 3.69166

  # Species Coefficients
  if (SPP == "AB") { # American Beech
    b4 <- -0.58341
  }
  if (SPP == "BA") { # Black Ash
    b4 <- -0.53830
  }
  if (SPP == "BT") { # Bigtooth
    b4 <- -0.55497
  }
  if (SPP == "OH") { # Other Hardwood
    b4 <- -1.00221
  }
  if (SPP == "PB") { # Paper Birch
    b4 <- -0.52068
  }
  if (SPP == "QA") { # Quaking Aspen
    b4 <- -0.54344
  }
  if (SPP == "RM") { # Red Maple
    b4 <- -0.55394
  }
  if (SPP == "RO") { # Red Oak
    b4 <- -0.49415
  }
  if (SPP == "SM") { # Sugar Maple
    b4 <- -0.53914
  }
  if (SPP == "WA") { # White Ash
    b4 <- -0.52310
  }
  if (SPP == "WO") { # White Oak
    b4 <- -0.53612
  }
  if (SPP == "YB") { # Yellow Birch
    b4 <- -0.53755
  }


  # Convert Form and Risk
  if (Form == 1) {
    b3 <- 0.89539
  }
  if (Form %in% c(2, 3, 5, 6, 7, 8)) {
    b3 <- 0
  }
  if (Form == 4) {
    b3 <- -1.76630
  }
  if (Risk == 1) {
    b2 <- 0
  }
  if (Risk == 2) {
    b2 <-  -1.04725
  }
  if (Risk == 3 | Risk == 4) {
    b2 <- -3.28995
  }


  Liklihood <- (exp(I + (blogDBH*log(DBH)) + b2 + b3 + b4*DBH + bHT*log(HT))) /
    (1 + exp(I + (blogDBH*log(DBH)) + b2 + b3 + b4*DBH + bHT*log(HT)))

  return(Liklihood)
  } else if(SPP %in% c("RS", "WS", "BF", "BS", "OC", "WC", "WP", "OS")){
    I <- -42.95
    blogDBH <- 18.03
    bHT <-  2.770

    if(Risk == 1){
      b2 = 0
    }
    if(Risk == 2){
      b2 = -2.077
    }
    if(Risk == 3){
      b2 = -4.575
    }
    if(Risk == 4){
      b2 = -4.575
    }


    if(SPP == "BF"){
      b3 = 0
      b4 = -0.7318
    }
    if(SPP == "BS"){
      b3 = -0.6871
      b4 = -0.7504
    }
    if(SPP == "EH"){
      b3 = -5.667
      b4 = -0.4432
    }
    if(SPP == "OC" | SPP == "WC"){
      b3 = -12.07
      b4 = -0.46986
    }
    if(SPP == "OS"){
      b3 = -43.01
      b4 = 0
    }
    if(SPP == "RS" | SPP == "WS"){
      b3 = -1.469
      b4 = -0.6444
    }
    if(SPP == "WP"){
      b3 = -10.02
      b4 = -0.3996
    }

    Liklihood <- (exp(I + (blogDBH*log(DBH)) + b2 + b3 + b4*DBH + bHT*log(HT))) /
      (1 + exp(I + (blogDBH*log(DBH)) + b2 + b3 + b4*DBH + bHT*log(HT)))

    return(Liklihood)
  } else {
    Liklihood <- 0
    return(Liklihood)
  }
}