diff --git a/R/assessment_functions.R b/R/assessment_functions.R index 4d79bb63..edb21b5b 100644 --- a/R/assessment_functions.R +++ b/R/assessment_functions.R @@ -156,7 +156,7 @@ run_assessment <- function( ctsm_ob$call.data <- NULL - + # identify which series are to be assessed in this run - defaults to all series_id <- row.names(ctsm_ob$timeSeries) @@ -276,7 +276,7 @@ assessment_engine <- function(ctsm.ob, series_id, parallel = FALSE, ...) { stations <- tibble::column_to_rownames(ctsm.ob$stations, "station_code") stations <- stations[unique(timeSeries$station_code), ] - + # set up parallel processing information @@ -305,7 +305,7 @@ assessment_engine <- function(ctsm.ob, series_id, parallel = FALSE, ...) { assessment <- pbapply::pblapply(data, ..., cl = cluster_id, FUN = function(x, ...) { # get info about the time series - + seriesID <- x$seriesID[1] seriesInfo <- sapply( timeSeries[seriesID,], diff --git a/R/import_check_functions.R b/R/import_check_functions.R index f05a40e1..b9d8a1a6 100644 --- a/R/import_check_functions.R +++ b/R/import_check_functions.R @@ -408,7 +408,7 @@ ctsm.check.species_group.biota <- function(data, info) { # check species_group appropriate for each determinand - id <- ctsm_is_contaminant(data$pargroup, exclude = "O-PAH") | data$group %in% "Auxiliary" + id <- ctsm_is_contaminant(data$pargroup, exclude = "O-PAH") | data$pargroup %in% "B-BIO" if (any(id)) data[id,] <- within(data[id,], { ok <- TRUE diff --git a/R/import_functions.R b/R/import_functions.R index 3dac46cf..f3d64d0b 100644 --- a/R/import_functions.R +++ b/R/import_functions.R @@ -2457,7 +2457,7 @@ tidy_contaminants <- function(data, info) { #' #' Cleans the data and turns it into time series structures ready for assessment #' -#' @param ctsm.obj the CTSM object, as returned from `tidy_data` +#' @param ctsm.obj the harsat object, as returned from `tidy_data` #' @param determinands the determinands to use, by default derived by #' calling `ctsm_get_determinands`, which takes values from #' the determinand reference table @@ -2842,7 +2842,7 @@ create_timeseries <- function( # merge auxiliary data with determinand data - data <- merge_auxiliary(data, info) + data <- merge_auxiliary(data, info, determinands) # impute %femalepop when missing and sex = 1 - write out remaining @@ -2852,7 +2852,7 @@ create_timeseries <- function( data <- ctsm.imposex.check.femalepop(data, info) } - + # convert data to basis of assessment data <- convert_to_target_basis(data, info, get_basis) @@ -2971,7 +2971,7 @@ create_timeseries <- function( "contact HARSAT development team") } - + # drop groups of data at stations with no data in recent years cat(" Dropping groups of compounds / stations with no data between", @@ -3437,8 +3437,16 @@ ctsm.imposex.check.femalepop <- function(data, info) { } -# utility function to get all determinand names from control structure - +#' extract determinands from control structure +#' +#' gets the names of all determinands involved in argument `determinands.control` +#' +#' @param control list of control structures passed into `create_timeseries` by +#' argument `determinand.control` +#' @param .names logical determining whether the names of the list are included +#' (`TRUE` default) or not (`FALSE`) +#' +#' @returns charcater string of determinands get_control_dets <- function(control, .names = TRUE) { if (is.null(control)) return(NULL) @@ -4216,24 +4224,41 @@ check_subseries <- function(data, info) { } -merge_auxiliary <- function(data, info) { - - # import_functions.R - # merge auxiliary variables with data +#' Merge auxiliary variables with data +#' +#' @param data a data frame containing the contaminant data in long format, both +#' the contaminants to be assessed and their auxiliary variables +#' @param info a harsat info object +#' @param determinands a character string given the identifiers of the +#' determinands that are to be assessed +#' +#' @returns a data frame containing the contaminant data in wide format, with +#' the auxiliary variables pivoted to match the determinands they are linked to +#' +#' @details +#' +#' `info$determinand` identifies which auxiliary variables are linked to each +#' determinand; `info$auxiliary` also allows the user (currently limited) some +#' control over how the auxiliary variables are linked to the determinand data +#' +#' Some variables can both be determinands to be assessed and auxiliary +#' variables (for example, CORG or DRYWT%) +#' +merge_auxiliary <- function(data, info, determinands) { control <- info$auxiliary - # identify auxiliary variables and split data set accordingly - auxiliary_var <- ctsm_get_auxiliary(data$determinand, info) + auxiliary_var <- ctsm_get_auxiliary(determinands, info) - id <- data$determinand %in% auxiliary_var - - auxiliary <- data[id, ] - data <- data[!id, ] + aux_id <- data$determinand %in% auxiliary_var + det_id <- data$determinand %in% determinands + auxiliary <- data[aux_id, ] + data <- data[det_id, ] + # ensure all auxiliary variables are present in output, by creating a # factor with levels given by auxiliary_var, and then splitting by this factor @@ -4319,12 +4344,12 @@ merge_auxiliary <- function(data, info) { } - + # finally merge data data <- merge(data, auxiliary_data, all.x = TRUE) } - + data <- droplevels(data) } @@ -4392,7 +4417,7 @@ convert_to_target_basis <- function(data, info, get_basis) { } - + # convert measurement data # print_warning gives the number of failures, which is the same for all id # so only print first time round @@ -4411,7 +4436,7 @@ convert_to_target_basis <- function(data, info, get_basis) { drywt_censoring = data[["DRYWT%.censoring"]], lipidwt = data[["LIPIDWT%"]], lipidwt_censoring = data[["LIPIDWT%.censoring"]], - exclude = data$group %in% c("Imposex", "Metabolites", "Effects") + exclude = data$group %in% c("Imposex", "Metabolites", "Effects") | data$determinand %in% c("LNMEA",'AGMEA') ), SIMPLIFY = FALSE ) @@ -4892,12 +4917,15 @@ normalise_sediment_HELCOM <- function(data, station_dictionary, info, control) { # normalises sediment concentrations - # method supplied by control + # method supplied by control - note that the value element for metals is + # hardwired by the code below (5 for AL; 52 for LI) and cannot be changed + # by the user ctsm_normalise_default <- list( metals = list(method = "pivot", normaliser = "AL", extra = NULL), copper = list(method = "hybrid", normaliser = "CORG", value = 5), - organics = list(method = "simple", normaliser = "CORG", value = 5), + organics = list(method = "simple", normaliser = "CORG", value = 5), + normalisers = list(method = "none"), exclude = NULL ) @@ -4920,8 +4948,8 @@ normalise_sediment_HELCOM <- function(data, station_dictionary, info, control) { data$normaliser <- NA_character_ data$normaliser_value <- NA_real_ data$normaliser_unit <- NA_character_ - - + + # exclude any data that do not need to be normalised # can do this globally with method = "none", but useful e.g. in the OSPAR # assessment where sediments in the Iberian Sea and Gulf of Cadiz are not @@ -4952,29 +4980,41 @@ normalise_sediment_HELCOM <- function(data, station_dictionary, info, control) { # make ad-hoc change to deal with LOIGN # must undo at the end of the code - data <- dplyr::mutate( - data, - .tmp = CORG, - .tmp.censoring = CORG.censoring, - .tmp.uncertainty = CORG.uncertainty, - CORG = dplyr::if_else(is.na(.tmp), 0.35 * LOIGN, CORG), - CORG.censoring = dplyr::if_else( - is.na(.tmp), - as.character(LOIGN.censoring), - as.character(CORG.censoring) - ), - CORG.censoring = factor(CORG.censoring), - CORG.uncertainty = dplyr::if_else(is.na(.tmp), 0.35 * LOIGN.uncertainty, CORG.uncertainty) - ) - + adjust_loign <- !is.null(data$CORG) & !is.null(data$LOIGN) + if (adjust_loign) { + data <- dplyr::mutate( + data, + .tmp = CORG, + .tmp.censoring = CORG.censoring, + .tmp.uncertainty = CORG.uncertainty, + CORG = dplyr::if_else( + is.na(.tmp), + 0.35 * LOIGN, + CORG + ), + CORG.censoring = dplyr::if_else( + is.na(.tmp), + as.character(LOIGN.censoring), + as.character(CORG.censoring) + ), + CORG.censoring = factor(CORG.censoring), + CORG.uncertainty = dplyr::if_else( + is.na(.tmp), + 0.35 * LOIGN.uncertainty, + CORG.uncertainty + ) + ) + } + # split into metals (CD, PB), copper and organics and then normalise each # with AL, CORG (LOIGN) and CORG (LOIGN) respectively# groupID <- dplyr::case_when( - data$determinand %in% c("CD", "PB") ~ "metals", - data$determinand %in% "CU" ~ "copper", - TRUE ~ "organics" + data$determinand %in% c("CD", "PB") ~ "metals", + data$determinand %in% "CU" ~ "copper", + data$determinand %in% c("CORG", "AL") ~ "normalisers", + TRUE ~ "organics" ) groupID <- factor(groupID) @@ -5214,16 +5254,17 @@ normalise_sediment_HELCOM <- function(data, station_dictionary, info, control) { data <- unsplit(data, groupID) - - data <- dplyr::mutate( - data, - CORG = .tmp, - CORG.censoring = .tmp.censoring, - CORG.uncertainty = .tmp.uncertainty, - .tmp = NULL, - .tmp.censoring = NULL, - .tmp.uncertainty = NULL - ) + if (adjust_loign) { + data <- dplyr::mutate( + data, + CORG = .tmp, + CORG.censoring = .tmp.censoring, + CORG.uncertainty = .tmp.uncertainty, + .tmp = NULL, + .tmp.censoring = NULL, + .tmp.uncertainty = NULL + ) + } if (any(exclude_id)) { data <- dplyr::bind_rows(data, excluded_data) @@ -5306,7 +5347,7 @@ normalise_biota_HELCOM <- function(data, station_dictionary, info, control) { groupID <- dplyr::if_else( data$species_group %in% "Fish" & - !(data$group %in% c("Metals", "Organofluorines", "Metabolites")), + !(data$group %in% c("Metals", "Organofluorines", "Metabolites", "Biological")), "lipid", "other" ) diff --git a/R/information_functions.R b/R/information_functions.R index 5ce64491..a02a68cc 100644 --- a/R/information_functions.R +++ b/R/information_functions.R @@ -594,27 +594,36 @@ ctsm_read_determinand <- function( }) - # check no auxiliary variables are going to be assessed - this will be - # allowed in later releases + ## check no auxiliary variables are going to be assessed - this will be + ## allowed in later releases + + # check that no determinand to be assessed has itself as its auxiliary data[paste0(compartment, "_assess")] <- lapply(compartment, function(id) { - group_id <- paste0(id, "_group") assess_id <- paste0(id, "_assess") - - not_ok <- data[[group_id]] %in% "Auxiliary" & data[[assess_id]] + aux_id <- paste0(id, "_auxiliary") + det_id <- "determinand" - if(any(not_ok)) { - det_id <- data$determinand[not_ok] + not_ok_det <- mapply( + function(det, aux) det %in% strsplit(aux, "~", fixed = TRUE)[[1]], + data[[det_id]], + data[[aux_id]], + USE.NAMES = FALSE + ) + + if(any(not_ok_det)) { + det_id <- data$determinand[not_ok_det] det_id <- sort(det_id) message( - "The following auxiliary variables have assess = TRUE which is ", - "currently not allowed.\nThese values of assess will be set to FALSE.\n", + "The following variables are set to be assessed and as auxiliary to ", + "themselves, which is not allowed.\nThese values of assess will be set ", + "to FALSE.\n", "Compartment: ", id, "\n", "Variables: ", paste(det_id, collapse = ", ") ) - - data[not_ok, assess_id] <- FALSE + + data[not_ok_det, assess_id] <- FALSE } data[[assess_id]] @@ -1105,7 +1114,7 @@ get_AC$biota <- function(data, AC, rt, export_all = FALSE) { data$datatype <- ctsm_get_datatype(data$determinand, rt) - if (!all(data$datatype %in% c("contaminant", "effect"))) { + if (!all(data$datatype %in% c("contaminant", "effect",'auxiliary'))) { stop("unrecognised datatype") } @@ -2325,7 +2334,7 @@ get_basis_default <- function(data, info) { # the exceptions are biological effects measurements where it is assumed the # data are submitted on the correct basis (or where basis isn't relevant) - + basis_id <- switch( info$compartment, biota = "W", @@ -2334,7 +2343,7 @@ get_basis_default <- function(data, info) { ) new_basis <- dplyr::if_else( - data$group %in% c("Imposex", "Metabolites", "Effects"), + data$group %in% c("Imposex", "Metabolites", "Effects") | data$determinand %in% c('LNMEA','AGMEA'), NA_character_, basis_id ) @@ -2377,6 +2386,12 @@ get_basis_most_common <- function(data, info) { return(x) } + if (unique(x$determinand) %in% c('LNMEA','AGMEA')) { + x$new_basis <- rep(NA_character_, nrow(x)) + x <- x[c(".order", "new_basis")] + return(x) + } + # check that have full basis information if (any(is.na(x$basis))) { @@ -2463,7 +2478,8 @@ get_basis_biota_OSPAR <- function(data, info) { out, .lw = .data$group %in% lw_group & !(.data$determinand %in% c("MCCP", "SCCP")), new_basis = dplyr::case_when( - .data$group %in% c("Imposex", "Effects", "Metabolites") ~ NA_character_, + .data$group %in% c("Imposex", "Effects", "Metabolites") ~ NA_character_, + .data$determinand %in% c('LNMEA','AGMEA') ~ NA_character_, .data$species_group %in% c("Bivalve", "Gastropod") ~ "D", .data$species_group %in% c("Fish", "Crustacean") & .lw & diff --git a/example_HELCOM.r b/example_HELCOM.r index 4daa2e46..c372641c 100644 --- a/example_HELCOM.r +++ b/example_HELCOM.r @@ -204,11 +204,29 @@ sediment_timeseries <- create_timeseries( normalise.control = list( metals = list(method = "pivot", normaliser = "AL", value = 5), copper = list(method = "hybrid", normaliser = "CORG", value = 5), - organics = list(method = "simple", normaliser = "CORG", value = 5) + organics = list(method = "simple", normaliser = "CORG", value = 5), + normalisers = list(method = "none") ) ) +# Note that this can be written more simply as the following, because the +# specification of normalise.control is just the default specification for +# HELCOM + +# sediment_timeseries <- create_timeseries( +# sediment_data, +# determinands.control = list( +# SBDE6 = list( +# det = c("BDE28", "BDE47", "BDE99", "BD100", "BD153", "BD154"), +# action = "sum" +# ), +# HBCD = list(det = c("HBCDA", "HBCDB", "HBCDG"), action = "sum") +# ), +# normalise = normalise_sediment_HELCOM +# ) + + # Now run the assessment. Again there is only one threshold, the EQS. This only # takes about a minute to run on my laptop. @@ -233,10 +251,11 @@ check_assessment(sediment_assessment) write_summary_table( sediment_assessment, determinandGroups = webGroups <- list( - levels = c("Metals", "Organotins", "PAH_parent", "PBDEs", "Organobromines"), + levels = c("Metals", "Organotins", "PAH_parent", "PBDEs", "Organobromines", + "Organic_constituents"), labels = c( "Metals", "Organotins", "Polycyclic aromatic hydrocarbons", - "Organobromines", "Organobromines" + "Organobromines", "Organobromines" , "Organic constituents" ) ), symbology = "default", @@ -261,7 +280,7 @@ write_summary_table( biota_data <- read_data( compartment = "biota", - purpose = "HELCOM", + purpose = "custom",#"HELCOM", contaminants = "biota.txt", stations = "stations.txt", data_dir = file.path("data", "example_HELCOM"), @@ -325,8 +344,7 @@ biota_timeseries <- create_timeseries( biota_assessment <- run_assessment( biota_timeseries, - AC = c("BAC", "EAC", "EQS", "MPC"), - parallel = TRUE + AC = c("BAC", "EAC", "EQS", "MPC") ) check_assessment(biota_assessment) @@ -355,12 +373,12 @@ write_summary_table( determinandGroups = list( levels = c( "Metals", "PAH_parent", "Metabolites", "PBDEs", "Organobromines", - "Organofluorines", "Chlorobiphenyls", "Dioxins" + "Organofluorines", "Chlorobiphenyls", "Dioxins", "Biological" ), labels = c( "Metals", "PAH compounds and metabolites", "PAH compounds and metabolites", "Organobromines", "Organobromines", "Organofluorines", - "PCBs and dioxins", "PCBs and dioxins" + "PCBs and dioxins", "PCBs and dioxins", "Biological" ) ), output_dir = file.path("output", "example_HELCOM") diff --git a/information/HELCOM_2023/determinand.csv b/information/HELCOM_2023/determinand.csv index 0b52c1a8..da226387 100644 --- a/information/HELCOM_2023/determinand.csv +++ b/information/HELCOM_2023/determinand.csv @@ -1,9 +1,9 @@ determinand,common_name,pargroup,biota_group,sediment_group,water_group,biota_assess,sediment_assess,water_assess,biota_unit,sediment_unit,water_unit,biota_auxiliary,sediment_auxiliary,water_auxiliary,biota_sd_constant,biota_sd_variable,sediment_sd_constant,sediment_sd_variable,water_sd_constant,water_sd_variable,distribution,good_status -AL,Aluminium,I-MET,Metals,Auxiliary,Metals,FALSE,FALSE,FALSE,ug/kg,%,ug/l,LNMEA~LIPIDWT%~DRYWT%,CORG~DRYWT%~LI,,,,0.000667,0.125644,,,lognormal,low +AL,Aluminium,I-MET,Metals,Metals,Metals,FALSE,FALSE,FALSE,ug/kg,%,ug/l,LNMEA~LIPIDWT%~DRYWT%,CORG~DRYWT%~LI,,,,0.000667,0.125644,,,lognormal,low CD,Cadmium,I-MET,Metals,Metals,Metals,TRUE,TRUE,TRUE,ug/kg,mg/kg,ug/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%~LI,,0.526867,0.07022,0.005,0.109988,0.002333,0.121905,lognormal,low CU,Copper,I-MET,Metals,Metals,Metals,FALSE,TRUE,FALSE,ug/kg,mg/kg,ug/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%~LI,,2,0.076316,0.158333,0.090478,0.025167,0.072003,lognormal,low HG,Mercury,I-MET,Metals,Metals,Metals,TRUE,FALSE,FALSE,ug/kg,mg/kg,ug/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%~LI,,0.541,0.087966,0.002167,0.100382,,,lognormal,low -LI,Lithium,I-MET,Metals,Auxiliary,Metals,FALSE,FALSE,FALSE,ug/kg,mg/kg,ug/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%,,,,0.033333,0.148842,,,lognormal,low +LI,Lithium,I-MET,Metals,Metals,Metals,FALSE,FALSE,FALSE,ug/kg,mg/kg,ug/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%,,,,0.033333,0.148842,,,lognormal,low PB,Lead,I-MET,Metals,Metals,Metals,TRUE,TRUE,TRUE,ug/kg,mg/kg,ug/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%~LI,,2.333333,0.10368,0.283333,0.090517,0.0115,0.08951,lognormal,low ANT,Anthracene,O-PAH,PAH_parent,PAH_parent,PAH_parent,FALSE,TRUE,FALSE,ug/kg,ug/kg,ng/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%~LOIGN,,0.066667,0.22,0.333333,0.170818,0.14445,0.1085,lognormal,low BAP,Benzo[a]pyrene,O-PAH,PAH_parent,PAH_parent,PAH_parent,TRUE,FALSE,FALSE,ug/kg,ug/kg,ng/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%~LOIGN,,0.053333,0.211934,0.333333,0.15,0.066667,0.14,lognormal,low @@ -81,14 +81,14 @@ TEQDFP,WHO-TEQ (DFP),OC-DX,Dioxins,Dioxins,,TRUE,FALSE,FALSE,ug/kg,ug/kg,,LNMEA~ PFOS,,O-FL,Organofluorines,Organofluorines,Organofluorines,TRUE,FALSE,TRUE,ug/kg,ug/kg,ng/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%~LOIGN,,0.066667,0.199772,,,0.003333,0.172725,lognormal,low N-PFOS,,O-FL,Organofluorines,Organofluorines,Organofluorines,FALSE,FALSE,FALSE,ug/kg,ug/kg,ng/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%~LOIGN,,,,,,,,lognormal,low BR-PFOS,,O-FL,Organofluorines,Organofluorines,Organofluorines,FALSE,FALSE,FALSE,ug/kg,ug/kg,ng/l,LNMEA~LIPIDWT%~DRYWT%,AL~CORG~DRYWT%~LOIGN,,,,,,,,lognormal,low -CORG,Organic Carbon,O-MAJ,,Auxiliary,,FALSE,FALSE,FALSE,,%,,,AL~DRYWT%~LOIGN,,,,0.01,0.1,,,, -LOIGN,Loss on ignition,O-MAJ,,Auxiliary,,FALSE,FALSE,FALSE,,%,,,AL~CORG~DRYWT%,,,,0,0.1,,,, -LNMEA,mean length,B-BIO,Auxiliary,,,FALSE,FALSE,FALSE,cm,,,,,,,,,,,,, -DRYWT%,dry weight,B-BIO,Auxiliary,Auxiliary,,FALSE,FALSE,FALSE,%,%,,,,,,,,,,,, -EXLIP%,extractable lipid,B-BIO,Auxiliary,,,FALSE,FALSE,FALSE,%,,,,,,,,,,,,, -LIPIDWT%,lipid weight,B-BIO,Auxiliary,,,FALSE,FALSE,FALSE,%,,,,,,,,,,,,, -FATWT%,fat weight,B-BIO,Auxiliary,,,FALSE,FALSE,FALSE,%,,,,,,,,,,,,, -%FEMALEPOP,proportion of females in population,B-END,Auxiliary,,,FALSE,FALSE,FALSE,%,,,,,,,,,,,,, +CORG,Organic Carbon,O-MAJ,,Organic_constituents,,FALSE,FALSE,FALSE,,%,,,AL~DRYWT%~LOIGN,,,,0.01,0.1,,,lognormal,low +LOIGN,Loss on ignition,O-MAJ,,Organic_constituents,,FALSE,FALSE,FALSE,,%,,,AL~CORG~DRYWT%,,,,0,0.1,,,lognormal,low +LNMEA,mean length,B-BIO,Biological,,,TRUE,FALSE,FALSE,cm,,,,,,,,,,,,lognormal,low +DRYWT%,dry weight,B-BIO,Biological,Biological,,FALSE,FALSE,FALSE,%,%,,,,,,,,,,,lognormal,low +EXLIP%,extractable lipid,B-BIO,Biological,,,FALSE,FALSE,FALSE,%,,,,,,,,,,,,, +LIPIDWT%,lipid weight,B-BIO,Biological,,,FALSE,FALSE,FALSE,%,,,,,,,,,,,,lognormal,low +FATWT%,fat weight,B-BIO,Biological,,,FALSE,FALSE,FALSE,%,,,,,,,,,,,,, +%FEMALEPOP,proportion of females in population,B-END,Imposex,,,FALSE,FALSE,FALSE,%,,,,,,,,,,,,, PYR1OH,1-hydroxy pyrene,B-MBA,Metabolites,,,TRUE,FALSE,FALSE,ng/ml,,,LNMEA~LIPIDWT%~DRYWT%,,,0.123333,0.15,,,,,lognormal,low VDS,Vas Deferens Sequence,B-END,Imposex,,,TRUE,FALSE,FALSE,st,,,%FEMALEPOP~LNMEA,,,,,,,,,multinomial,low VDSI,Vas Deferens Sequence Index,B-END,Imposex,,,FALSE,FALSE,FALSE,idx,,,%FEMALEPOP~LNMEA,,,,,,,,,quasibinomial,low diff --git a/inst/information/matrix.csv b/inst/information/matrix.csv index e650bb1f..9c648953 100644 --- a/inst/information/matrix.csv +++ b/inst/information/matrix.csv @@ -26,6 +26,7 @@ SB,soft body SEDTOT,unsieved SED20,< 20 micron fraction SED63,< 63 micron fraction +SH,shell SP,spleen TM,tail muscle UR,urine diff --git a/inst/markdown/report_assessment.Rmd b/inst/markdown/report_assessment.Rmd index e04eb34f..63cd0b00 100644 --- a/inst/markdown/report_assessment.Rmd +++ b/inst/markdown/report_assessment.Rmd @@ -717,7 +717,7 @@ ok <- info$compartment %in% c("biota", "sediment") No auxiliary variables currently plotted. ``` -```{r auxiliary_data, eval = ok, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 7} +```{r, eval = ok, include = FALSE} auxiliary_id <- switch( info$compartment, biota = c("concentration", "LNMEA", "DRYWT%", "LIPIDWT%"), @@ -732,6 +732,16 @@ auxiliary_id <- switch( } ) +auxiliary_missing <- sapply(auxiliary_id, function(x) is.null(data[[x]])) +ok <- ok & !any(auxiliary_missing) +``` + +```{asis, eval = !ok} +
+Not all required auxiliary variables are present. This is an area under development. Please contact the harsat development team. +``` + +```{r auxiliary_data, eval = ok, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 7} plot_auxiliary(data, assessment, info, assessment_object$info, auxiliary = auxiliary_id, xykey.cex = 1.2) ``` diff --git a/man/create_timeseries.Rd b/man/create_timeseries.Rd index 9feb104c..643a3083 100644 --- a/man/create_timeseries.Rd +++ b/man/create_timeseries.Rd @@ -17,7 +17,7 @@ create_timeseries( ) } \arguments{ -\item{ctsm.obj}{the CTSM object, as returned from \code{tidy_data}} +\item{ctsm.obj}{the harsat object, as returned from \code{tidy_data}} \item{determinands}{the determinands to use, by default derived by calling \code{ctsm_get_determinands}, which takes values from diff --git a/man/get_control_dets.Rd b/man/get_control_dets.Rd new file mode 100644 index 00000000..d80a5d21 --- /dev/null +++ b/man/get_control_dets.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/import_functions.R +\name{get_control_dets} +\alias{get_control_dets} +\title{extract determinands from control structure} +\usage{ +get_control_dets(control, .names = TRUE) +} +\arguments{ +\item{control}{list of control structures passed into \code{create_timeseries} by +argument \code{determinand.control}} + +\item{.names}{logical determining whether the names of the list are included +(\code{TRUE} default) or not (\code{FALSE})} +} +\value{ +charcater string of determinands +} +\description{ +gets the names of all determinands involved in argument \code{determinands.control} +} diff --git a/man/merge_auxiliary.Rd b/man/merge_auxiliary.Rd new file mode 100644 index 00000000..cf10092c --- /dev/null +++ b/man/merge_auxiliary.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/import_functions.R +\name{merge_auxiliary} +\alias{merge_auxiliary} +\title{Merge auxiliary variables with data} +\usage{ +merge_auxiliary(data, info, determinands) +} +\arguments{ +\item{data}{a data frame containing the contaminant data in long format, both +the contaminants to be assessed and their auxiliary variables} + +\item{info}{a harsat info object} + +\item{determinands}{a character string given the identifiers of the +determinands that are to be assessed} +} +\value{ +a data frame containing the contaminant data in wide format, with +the auxiliary variables pivoted to match the determinands they are linked to +} +\description{ +Merge auxiliary variables with data +} +\details{ +\code{info$determinand} identifies which auxiliary variables are linked to each +determinand; \code{info$auxiliary} also allows the user (currently limited) some +control over how the auxiliary variables are linked to the determinand data + +Some variables can both be determinands to be assessed and auxiliary +variables (for example, CORG or DRYWT\%) +} diff --git a/vignettes/example_HELCOM.Rmd.orig b/vignettes/example_HELCOM.Rmd.orig index 334d586a..8cb8977c 100644 --- a/vignettes/example_HELCOM.Rmd.orig +++ b/vignettes/example_HELCOM.Rmd.orig @@ -9,13 +9,13 @@ vignette: > # Introduction -This vignette shows how to do an assessment (mostly) following the approach -taken in HELCOM HOLAS3. +This vignette shows how to do an assessment following the approach +taken in HELCOM HOLAS3 (with a few exceptions). -The data were extracted from the ICES data base using the XHAT facilities on -the [ICES webservice](https://dome.ices.dk/api/swagger/index.html). The data -were extracted on 28 August 2023 and were filtered using is_helcom_area = TRUE -anad maxYear = 2020. The data were subsequently reduced in size to make them +The data were extracted from the ICES database on 28 August 2023 using the XHAT +facilities on the [ICES webservice](https://dome.ices.dk/api/swagger/index.html). +The data were filtered using is_helcom_area = TRUE +and maxYear = 2020. The data were subsequently reduced in size to make them more manageable for this example. **The data have not been scrutinised by data assessors so their use, and the results below, must be treated as illustrative only; in particular, they should not be used for any formal reporting**. @@ -34,22 +34,22 @@ explained in a subsequent vignette (not written yet). Second, the method for dealing with 'initial' data, unique to HELCOM, has not been implemented (this is not yet available in harsat). -We'll begin with contaminants in water which are the simplest data to assess. -We'll then move on to sediment and biota which have more features to consider. - -But before we get to any of that, we need to set up the environment. -First, load the `harsat` library (let's assume that you have already installed -it, as covered in the [Getting Started](/harsat.html) guide). +Before we start the assessments, we need to set up the environment. +First, we load the `harsat` library (if you do not have HARSAT installed go to +the [Getting Started](/harsat.html) guide and install it first). ```{r} library(harsat) ``` -We set up our main directory to find -the data files. If you use your own data files, you will need to point to -a directory containing a copy. We usually do this by putting our data files -in a directory `data`, and information files in a directory `information`, -but you can use any directory for these. +We then add the path to the main directory (the location of HARSAT, which +include the example folder) to 'working.directory', which allows us to point to +the data files used in this example (the path on Rob's computer is used as an +example). In assessments where you use your own data files, you will need to +point to the relevant directory. We usually put data files in a directory `data`, +and information files in a directory `information`, but you can use any directory +names you want. You will see the reference to 'working.directory', `data` and +`information` show up in the functions below. ```{r, echo=FALSE,results='asis'} library(rprojroot) @@ -59,16 +59,23 @@ cat("working.directory <- '", working.directory, "'\n", sep=""); cat("```\n") ``` +We will begin by assessing contaminants in water, which are the simplest to +assess. We will then do assessments of contaminants in sediment and biota, +which have more features to consider. +For sediment and biota only deviations from the water example is explained so +look through this example as well even if you are mainly interested in sediment +or biota. + # Water assessment -First, we use `read_data()` to read in the contaminant data, the station +We use the `read_data()` function to load the contaminant data, the station dictionary, and two important reference tables: the determinand and threshold reference tables. There are several things to note about the function call: * `purpose = "HELCOM"` means that the assessment configuration is set to mirror that of the HOLAS3 assessment. You can change the assessment -configuration in many ways using the control argument, but that is not -explained here. +configuration in many ways using the control argument (not +explained here). * `data_dir` identifies the directory storing the contaminant data (`water.txt`) and station dictionary (`stations.txt`). Using `file.path` prevents any difficulties using forward or backward slashes when writing file paths. @@ -79,14 +86,19 @@ explained here. * you don't have to specify the extraction date, but it can help to keep track of things. -As well as reading in the contaminant data and station dictionary, the -function matches each record in the contaminant data to a station in the -station dictionary. The process for doing this is quite complicated (and it -can take a few minutes to run) and we don't go into details here. +When the `read_data()` function is run, a number of comments are produced +indicating which data is being used, any warnings about the data, information on +station matching and max_year. + +The "Matching data with station dictionary" comment indicate that in addition +to loading the contaminant data and station dictionary, the +`read_data()` function also matches each record in the contaminant data to a +station in the station dictionary (and indicate when this is not possible). The +process is complicated and can take a few minutes to run. -Finally, a message is printed saying that the argument `max_year` is set to -2020. (This is as it should be since we set `maxYear` to be 2020 in the data -extraction.) But an important consequence is that a contaminant time series +The final comment indicates the `max_year` used for any observation in +the run. In this case it is 2020 as we set `maxYear` to be 2020 in the ICES +data extraction. An important implication is that a contaminant time series will only be assessed if it has some data in the period 2015 to 2020 (i.e. in the last six monitoring years). @@ -102,28 +114,32 @@ water_data <- read_data( ) ``` -We next simplify the data so that they are in the correct format for running -the assessment. This also involves deleting some data that do not meet the -conditions for the assessment. +The `tidy_data()` function simplifies the loaded data so that they are in the +correct format for running the assessment. This involves deleting data that +do not meet the conditions for the assessment as indicated by the printed +"Dropping .." comments. -Notice that a message appears about 'oddities'. Both tidy_data and -create_timeseries (the next function call) do a lot of checking of the data -and strange values are written to the oddities folder for you to have a look -at. This is in the hope that, if there are errors, they will get corrected and -resubmitted to the ICES database. It turns out there are no strange value at -this stage, but there are in the step that follows. +Notice that a message appears about 'oddities'. Both `tidy_data()` and +`create_timeseries()` (the next function call) do a lot of checking of the data +and strange values are written to excel files in the oddities folder. This allows +you to exmine the data problems that are found. For ICES extractions, identified +errors can be corrected by resubmission to the ICES database. If you are +subsequntly assessing your own data you can check and, if needed, correct error +in those yourself before moving on. It turns out there are no strange value in +our HELCOM test data at this stage, but there are in the step that follows. ```{r helcom-tidy-water} water_data <- tidy_data(water_data) ``` -We now do some more data cleaning and then group the data into time series. +The `create_timeseries()` function do some additional data cleaning and +group the data into time series. Each time series consists of measurements of a single determinand at a single monitoring station. Measurements in filtered and unfiltered samples are split into different time series. -The `determinands.control` argument in important to know about since it -allows related determinands to be processed in various ways. Here, it is just +The `determinands.control` argument +allows related determinands to be processed in various ways. Here, it is used for PFOS and its linear and branched components N-PFOS and BR-PFOS. PFOS measurements can be sumbitted as PFOS (the sum of the two components) or as N-PFOS and BR-PFOs (the individual components). The `determinands.control` @@ -149,14 +165,12 @@ get_timeseries(water_timeseries) |> head(10) ``` -At last it is time to run the assessment. You need to specify which thresholds -to use, otherwise the code will not use any of them. For water there is -only the EQS, but you still need to specify it. Look at the threshold -reference table to see what is available if you are unsure. Note that AC -stands for Assessment Criteria which is what thresholds are often called. The +The `run_assessment()` function runs the assessment. You need to specify which +thresholds (AC - Assessment Criteria) to use. Look at the +threshold reference table to see which types of tresholds are available if you +are unsure. For water there is only the EQS available. The parallel argument tells the code to use parallel processing. This usually -speeds things up considerably. The assessment took about 1.5 minutes to run on -my laptop. +speeds things up considerably. ```{r helcom-water-assessment} water_assessment <- run_assessment( @@ -167,24 +181,29 @@ water_assessment <- run_assessment( ``` -We now need to check whether there were any convergence issues. Lack of -convergence only rarely occurs, but when it does it is typically because there -are errors in the data (e.g. reported in incorrect units) or because there are -outliers, so the best thing to do is first check your data. However, -convergence can also be difficult if there are a lot of less-than measurements -in the time series. Describing how -to tweak the control arguments to get convergence is beyond the scope of this -vignette (need to write another vignette to discuss this). Fortunately, there -were no convergence issues here. +The `check_assessment()` function checks for any convergence issues for the time +series analysis. Lack of convergence only rarely occurs, but when it does it is +typically because there are errors in the data (e.g. reported in incorrect +units). Therefore, the first thing to do if time series have not converged is to +check your data. However, convergence problems can also relate to the presence +of outliers or a lot of less-than measurements in the time +series. Describing how to tweak the control arguments to get convergence is +beyond the scope of this vignette but an example is given in the biota +assessment below. Fortunately, there were no convergence issues here. ```{r helcom-water-check-assessment} check_assessment(water_assessment) ``` -It is time to look at the results! We can plot the data for each time series -along with the fitted model (see vignette for external data) or -print a table giving summary information about the assessment of each time -series. But first, we may need to create an output directory. +It is time to look at the results! We can print a table giving summary +information about the assessment of each time series (see below) and we +can print plots for each time series along with the fitted model (see vignette +for external data). + +We may need to create an output directory for these files if this does not +already exist. Here we create a sub-directory "example_HELCOM" under the +"output" folder to specify that this is a HELCOM assessment. You can create a +directory with any name you want for your assessment. ```{r helcom-summary-directory} summary.dir <- file.path(working.directory, "output", "example_HELCOM") @@ -194,11 +213,10 @@ if (!dir.exists(summary.dir)) { ``` Now that we have an output directory, we can write the summary table to -that directory. Of course, you can choose any directory to write to -here, so long as you are sure it exists. +that directory. -The code below prints out a CSV file giving summary information about -the assessment of each time series. This includes: +The `write_summary_tablet()` function prints a CSV file giving summary +information about the assessment of each time series. This includes: - meta-data such as the monitoring location and number of years of data for each time series @@ -230,7 +248,7 @@ The status is summarised as: Many aspects of the symbology can be controlled using the `symbology_control` argument. Users can also supply their own function to customise the symbology -to suit their needs. +to suit their needs (not explained here). ```{r helcom-water-summary} @@ -249,9 +267,9 @@ write_summary_table( # Sediment assessment -The sediment assessment is very similar, but has a few extra features related -to normalisation (to account for differences in grain size) which are described -below +The sediment assessment is very similar to the water assessment, but includes a +few extra features related to normalisation (to account for differences in grain +size) which are described below ```{r helcom-sediment-data} sediment_data <- read_data( @@ -268,24 +286,26 @@ sediment_data <- tidy_data(sediment_data) ``` The sediment data are grouped into time series which consist of the -measurements of single determinand in a single matrix (the fraction the sample -has been sieved to; e.g. `SED63` or `SEDTOT`) at a single monitoring station. +measurements of a single determinand in a single matrix (given as the fraction +the sample has been sieved to; e.g. `SED63` or `SEDTOT`) at a single monitoring +station. The `create_timeseries()` call for sediment differs from that for water in two ways. First, `determinands.control` identifies two groups of determinands that need to be summed. Second, the arguments `normalise` and `normalise.control` specify how the normalisation for grain size should be carried out. There are default functions for normalisation that will work in many cases. However, the -process for HELCOM is more complicated (because unlike other metals, copper is +process for HELCOM is more complicated (because, unlike other metals, copper is normalised to organic carbon rather than aluminium) so a customised function `normalise_sediment_HELCOM()` is provided. The argument `normalise.control` specifies that metals (apart from copper) will be normalised to 5% aluminium and copper and organics will be normalised to 5% organic carbon. The normalise -functions need a bit of work, so expect them to change. +functions need a bit of work, so they might change. All contaminant time series in sediment are assessed on a dry weight basis. The few measurements submitted on a wet weight basis are converted to a dry weight -basis using `DRYWT%` supporting information (also submitted with the data). +basis using `DRYWT%` supporting information (which needs to be submitted in +the data file). ```{r helcom-read-sediment} @@ -307,8 +327,8 @@ sediment_timeseries <- create_timeseries( ) ``` -Now run the assessment. Again there is only one threshold, the EQS. This only -takes about a minute to run on my laptop. +Now run the assessment. Similar to water there is only one threshold (AC), the +EQS. ```{r helcom-sediment-assessment} sediment_assessment <- run_assessment( @@ -318,18 +338,19 @@ sediment_assessment <- run_assessment( ) ``` -Everything has converged. +Using `check_assessment()` we see that everything has converged. ```{r helcom-sediment-convergence} check_assessment(sediment_assessment) ``` -Finally, we can plot individual time series assessments (see vignette for -external data) or print out the summary table. The `determinandGroups` argument +We print the table giving summary information for the assessment of each +time series. To print plots for each time series along with the fitted model see +vignette for external data. The `determinandGroups` argument allows the determinand groupings (as defined in the determinand reference table) to be renamed in the summary table. Below both the PBDEs and Organobromines -groups would appear as Organobromines in the summary table. +groups will appear as Organobromines in the summary table. ```{r helcom-sediment-summary} write_summary_table( @@ -353,13 +374,14 @@ write_summary_table( # Biota assessment -The main difference in the biota assessment is the inclusion of effects data. -This example has some PAH metabolite time series, but all imposex data have +The main difference between biota assessments and the other two assessments +is that they can include biological effects data in addition to contaminants. +Our example includes PAH metabolite time series, but imposex data have been excluded to keep things relatively simple. Imposex assessments have an -additional modelling stage and this will be described in another vignette (not +additional modelling stage which will be described in another vignette (not yet available). -The first two stages are just as before +The first two stages mirror that of water and sediment: ```{r helcom-biota-data} biota_data <- read_data( @@ -375,46 +397,44 @@ biota_data <- read_data( biota_data <- tidy_data(biota_data) ``` -The construction of the time series has a few more features. However, first we -need to provide the individual TEQs to allow the construction of the WHO TEQ -for dioxins, furans and planar PCBS (labelled TEQDFP). These are the values -for the human health QS. This stage won't be necessary in later releases. +The construction of the biota time series has a few more features than for water +and sediment. + +The biota data are divided into time series. Each time series consist of data +from a determinand in a single species and matrix (tissue type; +e.g. 'EH', 'LI' or 'SB') and from a single monitoring station. PAH +metabolite data are further grouped by analysis method ('method_analysis'). -The biota data are grouped into time series which consist of the -measurements of a single determinand in a single matrix (tissue type; e.g. -'EH', `LI' or `SB`) in a single species at a single monitoring station. PAH -metabolite data are further grouped by `method_analysis`. +The `determinands.control` argument contains more actions here than in the +water and sediment examples (sum, replace, bespoke). -The `determinands.control` argument does rather more here. There are five summed -variables: PFOS, SBDE6, HBCD, SCB6 and TEQDFP. There is also one variable -CB138+163 which needs to be relabeled as (replaced by) CB138. For the purposes -of the assessment, the contribution of CB163 is regarded as small. Similarly -CB138+163 is taken to be a good proxy for CB138. Note that the replacements must -be done before the six PCBs are summed to give SCB6 in order for them to be -included in the sum. +There are five variables (PFOS, SBDE6, HBCD, SCB6 and TEQDFP) that are being +summed. There is further one variable (CB138+163) which is being relabeled +(action = replace) as CB138 (for the purposes of the assessment, the +contribution of CB163 is regarded as small and CB138+163 is taken to be a good +proxy for CB138). Note that the replacements must be done before the six PCBs +are summed to give SCB6. The 'bespoke' action triggers a customised functions +that does more complicated data manipulations. Here it deals with the three +different ways in which lipid weight measurements can be submitted. The calculation of TEQDFP is more complex than that of the other summed variables. TEQDFP is the label used for the World Health Organisation Toxic Equivalent sum -for dioxins, furans and planar polychlorinated biphenyls. This is calculated -using the Toxic Equivalency Factors (TEFs) stored in `info_TEQ$HOLAS3`, with the -TEFS identified as `weights`. Note that the weights used in the example below -are those used in the HELCOM HOLAS 3 assessment, but have been superseeded; see -`help("info_TEQ")` for more details. - -There is also one 'bespoke' action in `determinands.control`. This triggers a -customised functions that does more complicated things. Here it deals with the -three different ways in which lipid weight measurements can be submitted. +for dioxins, furans and planar polychlorinated biphenyls. It is calculated +using individual Toxic Equivalency Factors (TEFs) for each substance (stored in +`info_TEQ$HOLAS3`), with the TEFs identified as `weights`. Note that the weights +used in the example below are those used in the HELCOM HOLAS 3 assessment, but +have been superseeded; see `help("info_TEQ")` for more details. -Finally, `normalise_biota_HELCOM()` is a customised function that determines which -measurements are normalised to 5% lipid in a HELCOM assessment. Again, the -normalisation functions are under active development and might well change -before the next release. +The `normalise_biota_HELCOM()` is a customised function that determines which +measurements are normalised to 5% lipid in the HELCOM assessment. The +normalisation functions are under active development and might change +in future releases. One thing that is not obvious from the function call is the choice of basis. By default, contaminant time series in biota are assessed on a wet weight basis, -with measurements submitted on a dry or lipid weight basis transformed to a wet -weight basis using supporting `DRYWT%` and `LIPIDWT%` information. Look at the -OSPAR or external data examples to see how the choice of basis can be changed. +where measurements submitted on a dry or lipid weight basis are transformed to a +wet weight basis using supporting `DRYWT%` and `LIPIDWT%` information. Look at +the OSPAR or external data examples to see how the choice of basis can be changed. ```{r helcom-biota-timeseries} biota_timeseries <- create_timeseries( @@ -446,7 +466,8 @@ biota_timeseries <- create_timeseries( ) ``` -The asssessment took about 4.3 minutes on my laptop +Now run the assessment. Where water and sediment had only an EQS threshold, +biota uses four different types of threshold for the assessment. ```{r helcom-biota-assessment} biota_assessment <- run_assessment( @@ -456,11 +477,12 @@ biota_assessment <- run_assessment( ) ``` -One time series has not converged. (The parameter estimates are fine, but the -standard errors are implausibly tight - if you look at the data, you will see why +When we run `check_assessment` we see that one time series has not converged. +(The parameter estimates are fine, but the standard errors are implausibly tight +- if you look at the data, you will see why the routines struggle with this time series.) The code below tweaks the arguments of the numerical differencing routine that calculates the standard -errors. Dealing with non-converged timeseries is a topic for a future vignette. +errors. ```{r helcom-biota-check-assessment} @@ -475,9 +497,9 @@ biota_assessment <- update_assessment( check_assessment(biota_assessment) ``` -The assessment uses four types of thresholds (BAC, EAC, EQS and MPC), but +The biota assessment uses four types of thresholds (BAC, EAC, EQS and MPC), but only one of them is applied to each timeseries, and they are all used to -delineate between good and poor status. The summary table can therefore be +delineate between good and poor status. The summary table can therefore be simplified by grouping the thresholds together using the `threshold_groups` argument and labelling them as an EQS (equivalent). The symbology then knows that, whichever threshold is applied, the status will be green if