Feature(#113): add calc_condition function

DESCRIPTION

@@ -38,3 +38,4 @@ Depends:
     R (>= 4.1.0)
 LazyData: true
 Config/roxygen2/version: 8.1.0
+Roxygen: list(markdown = TRUE)

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(calc_condition)
 export(calc_stratified_mean)
 export(calc_swept_area)
 export(get_area)

R/calc_condition.R

@@ -0,0 +1,233 @@
+#' Calculate species condition
+#'
+#' Calculates species [condition](https://noaa-edab.github.io/tech-doc/condition.html?q=condit#condition) from NEFSC survey data
+#'
+#' @section Species Coverage:
+#' The choice of \code{lengthweight} parameter determines species availability. When using
+#' \code{lengthweight = "Wigley"}, species coverage is constrained to the static parameters
+#' from [Wigley et al. (2003)](#references) and does not include Blueback herring or
+#' Northern Searobin. Using \code{lengthweight = "survdat"} includes all species from [Wigley et al. (2003)](#references) plus Blueback herring and Northern Searobin.
+#'
+#' @section References:
+#' Methods are outlined in the [technical documentation](https://noaa-edab.github.io/tech-doc/condition.html?q=condit#condition) for the State of the Ecosystem report
+#'
+#' Wigley SE, McBride HM, McHugh NJ. Length-weight relationships for 74 fish species collected during NEFSC research vessel bottom trawl surveys, 1992-99. 2003; Available: \url{https://repository.library.noaa.gov/view/noaa/3346}
+#'
+#'
+#' @inheritParams get_survdat_data
+#' @param surveyData Data table. NEFSC survey data generated by [get_survdat_data()]. See [examples](#ref-examples)
+#' @param lengthweight character string. Length-weight coefficients. Default value is 'Wigley'. Options = "survdat" or "Wigley". 'lengthweight = "survdat"' uses length-weight parameters from the NEFSC survdat database generated by [survdat::get_length_weight()]. 'lengthweight = "Wigley"' uses a static set of length-weight parameters from Wigley et al. (2003).
+#' @param aggregate_by vector. Default value is 'EPU'. Spatial scale to aggregate by: "EPU", "strata", or "none".
+#' @param by_sex logical. Default value is FALSE. If TRUE, calculates condition by sex. If FALSE, calculates condition across sexes.
+#' @param length_break numeric vector. Length intervals (cm) in which to calculation condition. Default (NULL) combines all lengths. Otherwise condition in calculated within provided interval. User must specify minimum and maximum lengths. For example to calculate condition for two groups \code{[0-20]} and \code{(21-70]}, use `length_breaks = c(0,20,70)`
+#'
+#' @return A data frame. Each row contains the estimated condition by year, species, sex, stratum, and length interval.
+#'
+#'
+#' \item{Species}{Species common name}
+#' \item{SVSPP}{A standard code which represents a species caught in a trawl or dredge. See [get_species()] for details}
+#' \item{sexMF}{Code indicating sex of fish. see [get_sex_fscs()] for details}
+#' \item{MeanCond}{Mean condition for year, species, stratum, sex, length interval}
+#' \item{YEAR}{Year for which condition is estimated}
+#' \item{EPU}{The Ecological Production Unit (EPU) for which condition is calculated}
+#' \item{STRATUM}{The survey stratum for which condition is calculated}
+#' \item{nCond}{Number of fish caught for Species/Sex/EPU combination in YEAR}
+#' \item{nYears}{Number of years where fish were observed for Species/Sex/EPU combination}
+#' \item{length_group}{The length interval for which condition is calulated}
+#'
+#' @examples
+#'
+#' \dontrun{
+#' # Pull data and apply conversion corrections
+#' data <- get_survdat_data(channel, getBio = TRUE)
+#' # Calculate condition aggregated by EPU, with Wigley length-weight coefficients, and with length breaks of 0, 20, and 50
+#' calc_condition(channel = channel, surveyData=data$survdat, lengthweight = "Wigley", aggregate_by = "EPU", by_sex = FALSE, length_break = c(0, 20, 50))
+#'
+#' # Calculate condition aggregated by strata, with survdat length-weight coefficients, and without length breaks
+#' calc_condition(channel = channel, surveyData=data$survdat, lengthweight = "survdat", aggregate_by = "strata", by_sex = FALSE, length_break = NULL)
+#'}
+#'
+#' @export
+
+calc_condition <- function(
+  channel,
+  surveyData,
+  lengthweight = 'Wigley',
+  aggregate_by = 'EPU',
+  by_sex = FALSE,
+  length_break = NULL
+) {
+  # add 0 to length_break if needed
+  if (!0 %in% length_break & !is.null(length_break)) {
+    length_break <- c(0, length_break)
+  }
+
+  survey.data <- surveyData |>
+    dplyr::left_join(survdat::EPUstrata, by = c("STRATUM"))
+
+  # Change sex = NA to sex = 0
+  fall <- survey.data |>
+    dplyr::filter(SEASON == "FALL") |>
+    dplyr::mutate(
+      sex = dplyr::if_else(is.na(SEX), "0", as.character(SEX))
+    )
+
+  if (lengthweight == "survdat") {
+    # pull LW data from 'get_length_weight' and filter to fall
+    survdat_lw <- survdat::get_length_weight(channel)
+    lwfall <- survdat_lw$data |>
+      dplyr::select(
+        SVSPP,
+        CATCHSEX,
+        SVLWCOEFF_FALL,
+        SVLWEXP_FALL
+      )
+  }
+
+  if (lengthweight == "Wigley") {
+    lwfall <- survdat::Wigley_LW |>
+      dplyr::filter(SEASON == "FALL") |>
+      dplyr::rename(
+        SVSPP = LW_SVSPP,
+        SVLWCOEFF_FALL = lna,
+        SVLWEXP_FALL = b
+      ) |>
+      dplyr::select(
+        SVSPP,
+        SVLWCOEFF_FALL,
+        SVLWEXP_FALL,
+        Gender
+      ) |>
+      dplyr::mutate(
+        CATCHSEX = dplyr::case_when(
+          Gender == "Combined" | Gender == "Unsexed" ~ as.character(0),
+          Gender == "Male" ~ as.character(1),
+          Gender == "Female" ~ as.character(2),
+          TRUE ~ NA
+        )
+      )
+    lwfall$CATCHSEX <- as.integer(lwfall$CATCHSEX)
+  }
+
+  # pull species common names and SVSPP from 'get_species'
+  species <- survdat::get_species(channel)
+  species_data <- species$data |>
+    dplyr::select(SVSPP, COMNAME) |>
+    tidyr::drop_na(SVSPP, COMNAME)
+
+  #group LWfall and species_data by SVSPP
+  combined_lw_species <- lwfall |>
+    dplyr::inner_join(species_data, by = "SVSPP")
+
+  # get sex info from 'get_sex' and join with combined_lw_species
+  sex_data <- survdat::get_sex(channel)
+  sex <- sex_data$data |>
+    dplyr::rename(CATCHSEX = SEX)
+
+  # add male, female, and unsexed descriptions to CATCHSEX numbers
+  ### NOTE: CATCHSEX 3-7 apply to Northern Shrimp and Lobster Only
+  LWsexed <- combined_lw_species |>
+    dplyr::inner_join(sex, by = "CATCHSEX") |>
+    dplyr::mutate(
+      SEX_DESCRIPTION = dplyr::case_when(
+        CATCHSEX == 0 ~ "Unsexed",
+        CATCHSEX == 1 ~ "Male",
+        CATCHSEX == 2 ~ "Female",
+        TRUE ~ SEX_DESCRIPTION # Keeps original value if CATCHSEX is anything else (e.g., NA)
+      )
+    )
+
+  # combine LW/sex/species data with Fall survey data
+  new_data <- LWsexed |>
+    dplyr::inner_join(fall, by = c("SVSPP", "CATCHSEX"))
+
+  # filters out values without losing rows with NAs:
+  mergewt <- dplyr::filter(new_data, is.na(INDWT) | INDWT < 900)
+  mergewtno0 <- dplyr::filter(mergewt, is.na(INDWT) | INDWT > 0.004)
+  mergelenno0 <- dplyr::filter(
+    mergewtno0,
+    is.na(LENGTH) | LENGTH > 0
+  )
+  mergelen <- dplyr::filter(mergelenno0, !is.na(LENGTH))
+  mergeindwt <- dplyr::filter(mergelen, !is.na(INDWT))
+  mergeLW <- dplyr::filter(mergeindwt, !is.na(SVLWCOEFF_FALL))
+
+  ### Calculate species condition ###
+
+  condcalc <- dplyr::mutate(
+    mergeLW,
+    predwt = (exp(SVLWCOEFF_FALL)) * LENGTH^SVLWEXP_FALL,
+    RelCond = INDWT / predwt
+  ) |>
+    dplyr::filter(is.na(RelCond) | RelCond < 300) |>
+    dplyr::group_by(SVSPP, SEX) |>
+    dplyr::mutate(
+      mean = mean(RelCond),
+      sd = stats::sd(RelCond)
+    ) |>
+    dplyr::ungroup()
+
+  condcalc <- condcalc |>
+    dplyr::filter(is.na(sex) | sex != 4) |>
+    dplyr::mutate(sexMF = sex)
+
+  species.codes <- survdat::get_species(channel = channel)$data |>
+    dplyr::rename(Species = COMNAME) |>
+    dplyr::select(Species, SVSPP) |>
+    tidyr::drop_na()
+
+  cond.epu <- dplyr::left_join(condcalc, species.codes, by = c("SVSPP"))
+
+  # Summarize annually -- parameterized groupings
+  grouping_vars <- c("Species", "YEAR")
+
+  if (aggregate_by == "strata") {
+    grouping_vars <- c(grouping_vars, "STRATUM")
+  }
+
+  if (aggregate_by == "EPU") {
+    grouping_vars <- c(grouping_vars, "EPU")
+  }
+
+  if (by_sex) {
+    grouping_vars <- c(grouping_vars, "sexMF")
+  }
+
+  if (!is.null(length_break)) {
+    cond.epu <- cond.epu |>
+      dplyr::mutate(
+        length_group = cut(
+          LENGTH,
+          breaks = length_break,
+          include.lowest = TRUE
+        )
+      )
+    grouping_vars <- c(grouping_vars, "length_group")
+  }
+
+  grouped_condition <- cond.epu |>
+    dplyr::group_by(!!!rlang::syms(grouping_vars))
+
+  if (aggregate_by == "none") {
+    grouped_condition <- grouped_condition
+  }
+
+  condition <- grouped_condition |>
+    dplyr::summarize(
+      MeanCond = mean(RelCond),
+      nCond = dplyr::n()
+    ) |>
+    dplyr::ungroup() |>
+    dplyr::filter(nCond >= 3) |>
+    # select columns
+    dplyr::select(dplyr::all_of(c(grouping_vars, "MeanCond", "nCond"))) |>
+    # group again, without YEAR
+    dplyr::group_by(
+      !!!rlang::syms(grouping_vars[-which(grouping_vars == "YEAR")])
+    ) |>
+    # filter to only species with 20+ years of data
+    dplyr::mutate(nYears = dplyr::n()) |>
+    dplyr::filter(nYears >= 20)
+
+  return(condition)
+}

R/data-EPUstrata.R

@@ -0,0 +1,13 @@
+#' Key of survey strata and corresponding EPUs
+#'
+#' Strata to EPU mappings are defined in the [technical documentation](https://noaa-edab.github.io/tech-doc/epu.html) for the State of the Ecosystem report
+#'
+#' @docType data
+#' @name EPUstrata
+#' @keywords datasets
+#' @format A tibble with 2 variables:
+#' \describe{
+#'   \item{\code{STRATUM}}{integer. A numerical code which represents a survey stratum in the NEFSC Bottom Trawl Survey. See [get_strata()] for details}
+#'   \item{\code{EPU}}{charater. An abbreviated name for the Ecological Production Unit (EPU) that corresponds to the STRATUM column}
+#' }
+"EPUstrata"

R/data-Wigley_LW.R

@@ -0,0 +1,34 @@
+#' Length-weight parameters from Wigley et al. (2003)
+#'
+#' Length-weight relationships for 74 fish species collected during NEFSC
+#' research vessel bottom trawl surveys, 1992-99 as published in [Wigley et al. (2003)](#references)
+#'
+#' @section Species Exclusions:
+#' Note that the dataset does not include length-weight parameters for Blueback herring or Northern searobin.
+#'
+#' @references
+#'
+#' Wigley SE, McBride HM, McHugh NJ. Length-weight relationships for 74 fish species collected during NEFSC research vessel bottom trawl surveys, 1992-99. 2003; Available: \url{https://repository.library.noaa.gov/view/noaa/3346}
+#'
+#'
+#' @docType data
+#' @name Wigley_LW
+#' @keywords datasets
+#' @format A tibble with 14 variables:
+#' \describe{
+#'   \item{\code{SpeciesName}}{character. Species common name}
+#'   \item{\code{LW_SVSPP}}{integer. A standard code which represents a species caught in a trawl or dredge. Refer to the SVDBS.SVSPECIES_LIST}
+#'   \item{\code{Season}}{character. Season sampled (Winter, Autumn, Spring or a combination)}
+#'   \item{\code{Gender}}{character. Gender of the fish sampled (Female, Male, or Combined)}
+#'   \item{\code{N}}{integer. Number of fish sampled}
+#'   \item{\code{ln_a}}{character. Length-weight parameter estimate}
+#'   \item{\code{SEa}}{double. Standard error of ln a}
+#'   \item{\code{b}}{double. Length-weight parameter estimate}
+#'   \item{\code{SEb}}{double. Standard error of b}
+#'   \item{\code{SE_estimate}}{double. Standard error of the weight estimate}
+#'   \item{\code{r2}}{double. Regression correlation coefficient}
+#'   \item{\code{lna1}}{double. Length-weight parameter estimate}
+#'   \item{\code{lna}}{double. Inverse of ln_a and lna1}
+#'   \item{\code{SEASON}}{character. Parsed season data into Winter, Spring, or Fall}
+#' }
+"Wigley_LW"

R/get_survdat_clam_data.r

@@ -12,6 +12,14 @@
 #'
 #' @return A list containing a Data frame (data.table) (n x 21) and a list of SQL queries used to pull the data, the date of the pull, and the call expression
 #'
+#' @family survdat
+#'
+#'@examples
+#'\dontrun{
+#' channel <- dbutils::connect_to_database("serverName","userName")
+#' get_survdat_clam_data(channel)
+#' }
+#'
 #'@export
 
 #-------------------------------------------------------------------------------

data-raw/create_EPUstrata.R

@@ -0,0 +1,23 @@
+#' Create survdat::EPUstrata data object
+
+create_EPUstrata <- function() {
+  # Define vectors containing strata values for each EPU
+  MAB <- c(1010:1080, 1100:1120, 1600:1750, 3010:3450, 3470, 3500, 3510)
+  GB <- c(1090, 1130:1210, 1230, 1250, 3460, 3480, 3490, 3520:3550)
+  GOM <- c(1220, 1240, 1260:1290, 1360:1400, 3560:3830)
+  SS <- c(1300:1352, 3840:3990)
+
+  EPU_df <- data.frame(
+    STRATUM = c(MAB, GB, GOM, SS),
+    EPU = rep(
+      c("MAB", "GB", "GOM", "SS"),
+      times = c(length(MAB), length(GB), length(GOM), length(SS))
+    )
+  )
+
+  EPUstrata <- tibble::as_tibble(EPU_df)
+
+  usethis::use_data(EPUstrata, overwrite = TRUE)
+
+  return(EPUstrata)
+}

data-raw/create_Wigley_LW.R

@@ -0,0 +1,14 @@
+#' Create survdat::Wigley_LW data object
+#' Dataset read in as csv ("tech_memo_parameters_table_format.csv") in original condition calculation here (https://github.com/Laurels1/Condition/blob/master/R/RelConditionEPU.R
+
+create_Wigley_LW <- function() {
+  # Read in the CSV file
+  Wigley_LW <- read.csv(
+    here::here('data-raw/data/Wigley_LW.csv'),
+    stringsAsFactors = FALSE
+  )
+
+  usethis::use_data(Wigley_LW, overwrite = TRUE)
+
+  return(Wigley_LW)
+}