R/utils.R

In ebirdst: Access and Analyze eBird Status and Trends Data Products

#' Calculate MCC and F1 score
#'
#' Given binary observed and predicted response, estimate Matthews correlation
#' coefficient (MCC) and the F1 score.
#'
#' @param observed logical or 0/1; the observed binary response.
#' @param predicted logical or 0/1; the predicted binary response. This will typically
#'   need to be generated by applying a threshold to the continuous predicted
#'   response.
#'
#' @return A list with two elements: `mcc` and `f1`.
#' @export
#' @examples
#' obs <- c(rep(1L, 1000L), rep(0L, 10000L))
#' pred <- c(rbeta(300L, 12, 2), rbeta(700L, 3, 4), rbeta(10000L, 2, 3))
#' calculate_mcc_f1(obs > 0, pred > 0.5)
calculate_mcc_f1 <- function(observed, predicted) {
  if (!requireNamespace("PresenceAbsence", quietly = TRUE)) {
    stop("Package 'PresenceAbsence' must be installed to use this function.")
  }
  stopifnot(
    is.logical(observed) || all(observed %in% c(0L, 1L)),
    is.logical(predicted) || all(predicted %in% c(0L, 1L)),
    length(observed) == length(predicted)
  )

  # confusion matrix
  obs_pred <- data.frame(
    blank = "x",
    obs = as.numeric(observed),
    pred = as.numeric(predicted)
  )
  cmx <- PresenceAbsence::cmx(obs_pred, na.rm = TRUE)
  tp <- cmx[1, 1]
  fp <- cmx[1, 2]
  fn <- cmx[2, 1]
  tn <- cmx[2, 2]

  # f1 score
  p_pat <- if ((tp + fp) == 0) 0 else tp / (tp + fp)
  r_pat <- if ((tp + fn) == 0) 0 else tp / (tp + fn)
  f1 <- if ((p_pat + r_pat) == 0) 0 else 2 * ((p_pat * r_pat) / (p_pat + r_pat))

  # mcc
  mcc_denom <- sqrt((tp + fp) * (tp + fn) * (tn + fp) * (tn + fn))
  mcc <- if (mcc_denom == 0) 0 else ((tp * tn) - (fp * fn)) / mcc_denom
  return(list(f1 = f1, mcc = mcc))
}


#' Get the Status and Trends week that a date falls into
#'
#' @param dates a vector of dates.
#' @param version One of `2021` for the date scheme used for the 2021 and prior
#' data releases or `2022` for the date scheme used in the 2022 and subsequent
#' releases.  Default is `2022`.
#' @return An integer vector of weeks numbers from 1-52.
#' @export
#' @examples
#' d <- as.Date(c("2016-04-08", "2018-12-31", "2014-01-01", "2018-09-04"))
#' date_to_st_week(d)
date_to_st_week <- function(dates, version = 2022) {
  stopifnot(version %in% c(2021, 2022))

  if (version == 2021) {
    # as used in ebirdst 2021 data release
    dv <- seq(from = 0, to = 1, length.out = 52 + 1)
    days <- (as.POSIXlt(dates)$yday + 0.5) / 366
    check_d <- function(x) {
      which(x >= dv[-length(dv)] & x < dv[-1])
    }
    return(vapply(days, check_d, FUN.VALUE = integer(length = 1)))
  }

  # 2022 date scheme
  # Each week is defined by it's central Julian Date which are: seq(4, 366, 7)
  # Dates are converted to Julian dates and assigned to nearest week.
  breaks <- c(-Inf, seq(7.5, 357.5, 7), Inf) # breaks are midway between weeks
  jd <- as.POSIXlt(dates)$yday + 1 # Julian Day, yday starts with 0 for Jan 1
  return(findInterval(jd, breaks))
}


#' Get eBird species code for a set of species
#'
#' Give a vector of species codes, common names, and/or scientific names, return
#' a vector of 6-letter eBird species codes. This function will only look up
#' codes for species for which eBird Status and Trends results exist.
#'
#' @param x character; vector of species codes, common names, and/or scientific
#'   names.
#'
#' @return A character vector of eBird species codes.
#' @export
#'
#' @examples
#' get_species(c("Black-capped Chickadee", "Poecile gambeli", "carchi"))
get_species <- function(x) {
  stopifnot(is.character(x))
  r <- ebirdst::ebirdst_runs
  r <- r[!stringr::str_detect(r$species_code, "example"), ]
  x <- tolower(trimws(x))

  # species code
  code <- match(x, tolower(r$species_code))
  # scientific name
  sci <- match(x, tolower(r$scientific_name))
  # common names
  com <- match(x, tolower(r$common_name))
  # combine
  codes <- r$species_code[dplyr::coalesce(code, sci, com)]
  # adjust for example dataset
  codes[x == "yebsap-example"] <- "yebsap-example"
  return(codes)
}


# internal ----

is_integer <- function(x) {
  return(isTRUE(is.integer(x) || (is.numeric(x) && all(x == as.integer(x)))))
}

is_count <- function(x) {
  return(isTRUE(is_integer(x) && length(x) == 1 && x >= 0))
}

is_flag <- function(x) {
  return(isTRUE(is.logical(x) && length(x) == 1 && !is.na(x)))
}