45 handling of relative errors

NAMESPACE

@@ -8,13 +8,15 @@ S3method(get_analytical_columns,ASTR)
 S3method(get_concentration_columns,ASTR)
 S3method(get_contextual_columns,ASTR)
 S3method(get_element_columns,ASTR)
+S3method(get_error_columns,ASTR)
 S3method(get_isotope_columns,ASTR)
 S3method(get_ratio_columns,ASTR)
 S3method(get_unit_columns,ASTR)
 S3method(print,ASTR)
 S3method(remove_units,ASTR)
 S3method(validate,ASTR)
 S3method(validate,default)
+export(abs_to_rel)
 export(albarede_juteau_1984)
 export(as_ASTR)
 export(at_to_wt)
@@ -30,13 +32,15 @@ export(get_analytical_columns)
 export(get_concentration_columns)
 export(get_contextual_columns)
 export(get_element_columns)
+export(get_error_columns)
 export(get_isotope_columns)
 export(get_ratio_columns)
 export(get_unit_columns)
 export(oxide_to_element)
 export(pb_iso_age_model)
 export(pointcloud_distribution)
 export(read_ASTR)
+export(rel_to_abs)
 export(remove_units)
 export(stacey_kramers_1975)
 export(unify_concentration_unit)

NEWS.md

@@ -2,3 +2,4 @@
 
 * Implementation of ASTR schema.
 * Support and conversions for geochemical non-SI units: *at%*, *wt%* (elements and oxides) 
+* Support and conversions for relative and absolute analytical precisions. 

R/ASTR_basic.R

@@ -174,17 +174,19 @@ as_ASTR <- function(
   df4 <- remove_unit_substrings(df3)
   # turn into tibble-derived object
   df5 <- tibble::new_tibble(df4, nrow = nrow(df4), class = "ASTR")
+  # convert relative to absolute errors
+  df6 <- rel_to_abs(df5)
   # post-reading validation
   if (validate) {
-    validation_output <- validate(df5, quiet = FALSE)
+    validation_output <- validate(df6, quiet = FALSE)
     if (nrow(validation_output) > 0) {
       warning(
         "See the full list of validation output with: ",
         "ASTR::validate(<your ASTR object>)."
       )
     }
   }
-  return(df5)
+  return(df6)
 }
 
 # helper function to rename column names

R/ASTR_colname_parser.R

@@ -282,7 +282,7 @@ is_err_percent <- function(colname) {
   grepl(err_percent(), colname, perl = TRUE)
 }
 is_err_abs <- function(colname) {
-  grepl(err_abs(), colname, perl = TRUE)
+  grepl(paste0("(", err_abs(), ")(?!%)"), colname, perl = TRUE)
 }
 is_isotope_ratio <- function(colname) {
   grepl(isotope_ratio(), colname, perl = TRUE)
@@ -351,6 +351,8 @@ err_percent <- function() {
 }
 err_2sd_percent <- function() "\\_err2SD%"
 err_sd_percent <- function() "\\_errSD%"
+err_2se_percent <- function() "\\_err2SE%"
+err_se_percent <- function() "\\_errSE%"
 
 err_abs <- function() {
   paste0(c(

R/ASTR_column_select.R

@@ -80,6 +80,16 @@ get_concentration_columns.ASTR <- function(x, ...) {
   get_cols_with_ac_class(x, c("ASTR_id", "ASTR_concentration"))
 }
 
+#' @rdname ASTR
+#' @export
+get_error_columns <- function(x, ...) {
+  UseMethod("get_error_columns")
+}
+#' @export
+get_error_columns.ASTR <- function(x, ...) {
+  get_cols_with_ac_class(x, c("ASTR_id", "ASTR_error"))
+}
+
 get_cols_with_unit <- function(x, units) {
 
   units <- sapply(units, function(unit) transform_notation(unit))

R/ASTR_conversion_error.R

@@ -0,0 +1,100 @@
+#' Convert between relative and absolute analytical uncertainties
+#'
+#' Convert relative to absolute analytical uncertainties and vice versa in ASTR
+#' objects. Work only for objects of class `ASTR`.
+#'
+#' @param df An ASTR object
+#' @return An ASTR object with converted analytical precision columns. The
+#'   unchanged input, if it does not contain columns of the respective
+#'   analytical precision type
+#'
+#'
+#' @name error_conversion
+#' @export
+#'
+#' @examples
+#' test_file <- system.file("extdata", "test_data_input_good.csv", package = "ASTR")
+#' arch <- read_ASTR(test_file, id_column = "Sample", context = 1:7)
+#'
+#' arch2 <- abs_to_rel(arch)
+#'
+#' arch3 <- rel_to_abs(arch2)
+#'
+#' # Conversion is lossless
+#' all.equal(arch$`SiO2_errSD`, arch3$`SiO2_errSD`)
+#'
+rel_to_abs <- function(df) {
+
+  # Basic checks
+  checkmate::assert_class(df, "ASTR")
+
+  # Find all error columns
+  error_cols <- colnames(df)[is_err_percent(colnames(df))]
+
+  if (length(error_cols) == 0) {
+    return(df)
+  }
+
+  df_old <- df
+
+  # Process all error columns
+  for (err_col in error_cols) {
+    base_name <- remove_suffix(err_col)
+
+    if (base_name %in% names(df)) {
+      # absolute = relative * measured_value
+      df[[err_col]] <- df[[err_col]] * df[[base_name]] / ifelse(inherits(df[[base_name]], "units"), 1, 100)
+
+      # Set units to match the concentration column
+      if (inherits(df[[base_name]], "units")) {
+        units(df[[err_col]]) <- units(df[[base_name]])
+      } else {
+        units(df[[err_col]]) <- NULL
+      }
+    }
+  }
+
+  # assign ASTR class
+  df <- preserve_ASTR_attrs(df, df_old)
+
+  # rename column names
+  colnames(df)[colnames(df) %in% error_cols] <- gsub("%", "", error_cols)
+
+  return(df)
+}
+
+#' @rdname error_conversion
+#' @export
+abs_to_rel <- function(df) {
+
+  # Basic checks
+  checkmate::assert_class(df, "ASTR")
+
+  # Find all error columns
+  error_cols <- colnames(df)[is_err_abs(colnames(df))]
+
+  if (length(error_cols) == 0) {
+    return(df)
+  }
+
+  df_old <- df
+
+  # Process all error columns
+  for (err_col in error_cols) {
+    base_name <- remove_suffix(err_col)
+
+    # relative = (absolute / measured_value)
+    df[[err_col]] <- df[[err_col]] / df[[base_name]] * ifelse(inherits(df[[base_name]], "units"), 1, 100)
+
+    # Set units to percent
+    units(df[[err_col]]) <- units::as_units("%")
+  }
+
+  # assign ASTR class
+  df <- preserve_ASTR_attrs(df, df_old)
+
+  # rename columns
+  colnames(df)[colnames(df) %in% error_cols] <- paste0(error_cols, "%")
+
+  return(df)
+}

_pkgdown.yml

@@ -23,6 +23,7 @@ reference:
     Functions for creating and handling ASTR objects.
   contents:
   - ASTR
+  - error_conversion
 - title: Data subsetting
   desc: >
     Pre-compiled lists for easier subsetting of data

tests/testthat/_snaps/ASTR_basic.md

@@ -45,21 +45,21 @@
       12      ICP-MS    0.512991  2.04       0.8000 4.5 [wtP] 0.070 [wtP]  1.56 [wtP]
       13      ICP-MS          NA -0.24      -0.0504 3.2 [wtP] 0.040 [wtP]  3.64 [wtP]
       14      ICP-MS          NA -0.42      -0.0420 3.7 [wtP] 0.060 [wtP]  7.51 [wtP]
-                MgO      Al2O3        SiO2 SiO2_errSD%       P2O5          S
-      1  0.77 [wtP] 3.92 [wtP] 31.63 [wtP]    4.30 [%] 0.15 [wtP] 2.57 [atP]
-      2  0.55 [wtP] 3.46 [wtP] 29.06 [wtP]    2.88 [%]   NA [wtP]   NA [atP]
-      3  0.33 [wtP] 5.87 [wtP] 30.50 [wtP]    5.71 [%] 0.11 [wtP] 3.68 [atP]
-      4  0.76 [wtP] 4.33 [wtP] 25.73 [wtP]    2.22 [%] 0.23 [wtP] 2.76 [atP]
-      5  0.52 [wtP] 4.17 [wtP] 43.50 [wtP]    4.65 [%] 0.40 [wtP] 0.57 [atP]
-      6  0.58 [wtP] 3.58 [wtP] 33.83 [wtP]    3.67 [%] 0.29 [wtP] 0.61 [atP]
-      7  0.64 [wtP] 5.60 [wtP] 33.81 [wtP]    3.00 [%] 0.43 [wtP] 0.77 [atP]
-      8  0.59 [wtP] 4.47 [wtP] 26.39 [wtP]    2.87 [%] 0.62 [wtP] 3.04 [atP]
-      9  0.53 [wtP] 3.73 [wtP] 21.18 [wtP]    3.44 [%] 0.24 [wtP] 3.93 [atP]
-      10 0.70 [wtP] 6.27 [wtP] 31.73 [wtP]    6.12 [%] 0.28 [wtP] 0.21 [atP]
-      11 0.46 [wtP] 2.91 [wtP] 16.91 [wtP]    3.73 [%] 0.22 [wtP] 3.57 [atP]
-      12 0.58 [wtP] 3.60 [wtP] 28.49 [wtP]    4.89 [%] 0.19 [wtP] 1.95 [atP]
-      13 0.88 [wtP] 5.24 [wtP] 38.04 [wtP]   10.12 [%] 0.30 [wtP] 1.35 [atP]
-      14 0.84 [wtP] 3.97 [wtP] 31.67 [wtP]    4.30 [%] 0.18 [wtP] 0.83 [atP]
+                MgO      Al2O3        SiO2     SiO2_errSD       P2O5          S
+      1  0.77 [wtP] 3.92 [wtP] 31.63 [wtP] 1.360090 [wtP] 0.15 [wtP] 2.57 [atP]
+      2  0.55 [wtP] 3.46 [wtP] 29.06 [wtP] 0.836928 [wtP]   NA [wtP]   NA [atP]
+      3  0.33 [wtP] 5.87 [wtP] 30.50 [wtP] 1.741550 [wtP] 0.11 [wtP] 3.68 [atP]
+      4  0.76 [wtP] 4.33 [wtP] 25.73 [wtP] 0.571206 [wtP] 0.23 [wtP] 2.76 [atP]
+      5  0.52 [wtP] 4.17 [wtP] 43.50 [wtP] 2.022750 [wtP] 0.40 [wtP] 0.57 [atP]
+      6  0.58 [wtP] 3.58 [wtP] 33.83 [wtP] 1.241561 [wtP] 0.29 [wtP] 0.61 [atP]
+      7  0.64 [wtP] 5.60 [wtP] 33.81 [wtP] 1.014300 [wtP] 0.43 [wtP] 0.77 [atP]
+      8  0.59 [wtP] 4.47 [wtP] 26.39 [wtP] 0.757393 [wtP] 0.62 [wtP] 3.04 [atP]
+      9  0.53 [wtP] 3.73 [wtP] 21.18 [wtP] 0.728592 [wtP] 0.24 [wtP] 3.93 [atP]
+      10 0.70 [wtP] 6.27 [wtP] 31.73 [wtP] 1.941876 [wtP] 0.28 [wtP] 0.21 [atP]
+      11 0.46 [wtP] 2.91 [wtP] 16.91 [wtP] 0.630743 [wtP] 0.22 [wtP] 3.57 [atP]
+      12 0.58 [wtP] 3.60 [wtP] 28.49 [wtP] 1.393161 [wtP] 0.19 [wtP] 1.95 [atP]
+      13 0.88 [wtP] 5.24 [wtP] 38.04 [wtP] 3.849648 [wtP] 0.30 [wtP] 1.35 [atP]
+      14 0.84 [wtP] 3.97 [wtP] 31.67 [wtP] 1.361810 [wtP] 0.18 [wtP] 0.83 [atP]
                 CaO       TiO2        MnO      FeOtot FeOtot_err2SD       ZnO
       1  2.11 [wtP] 0.52 [wtP] 0.20 [wtP] 43.83 [wtP]   4.120 [wtP]  5.64 [%]
       2  2.34 [wtP] 0.49 [wtP] 0.54 [wtP] 51.02 [wtP]   3.890 [wtP]  4.07 [%]

tests/testthat/test_ASTR_column_selection.R

@@ -32,8 +32,14 @@ test_that("column selection based on ASTR column types", {
         "U", "V", "Cr", "Co", "Ni", "Sr", "Se"
       )
   )
+  expect_all_true(
+    colnames(get_error_columns(test_input)) ==
+      c("ID", "d65Cu_err2SD", "SiO2_errSD", "FeOtot_err2SD", "206Pb/204Pb_err2SD",
+        "207Pb/204Pb_err2SD", "208Pb/204Pb_err2SD", "207Pb/206Pb_err2SD", "208Pb/206Pb_err2SD"
+      )
+  )
   expect_all_true(
     colnames(get_unit_columns(test_input, c("ng/g", "µg/ml", "%"))) ==
-      c("ID", "SiO2_errSD%", "ZnO", "Ag", "Sn")
+      c("ID", "ZnO", "Ag", "Sn")
   )
 })

tests/testthat/test_conversion_error.R

@@ -0,0 +1,77 @@
+test_input <- suppressWarnings(
+  read_ASTR(
+    system.file("extdata", "test_data_input_good.csv", package = "ASTR"),
+    id_column = "Sample",
+    context = c("Lab no.", "Site", "latitude", "longitude", "Type", "method_comp")
+  )
+)
+
+# --- abs_to_rel ---
+
+result <- abs_to_rel(test_input)
+
+test_that("abs_to_rel: numeric conversion is correct", {
+  expect_equal(as.numeric(result[["d65Cu_err2SD%"]]),
+               as.numeric(test_input[["d65Cu_err2SD"]] / test_input[["d65Cu"]] * 100), tolerance = 1e-6)
+})
+
+test_that("abs_to_rel: error col units are percent after conversion", {
+  expect_equal(as.character(units(result[["SiO2_errSD%"]])), "%")
+})
+
+test_that("abs_to_rel: input returned unchanged if no absolute error cols present", {
+  df_no_err <- test_input[, !is_err_abs(colnames(test_input))]
+  expect_identical(abs_to_rel(df_no_err), df_no_err)
+})
+
+test_that("abs_to_rel: errors on non-ASTR input", {
+  expect_error(abs_to_rel(data.frame(x = 1)))
+})
+
+# --- rel_to_abs ---
+
+result2 <- rel_to_abs(result)
+
+test_that("rel_to_abs: numeric conversion is correct", {
+  expect_equal(as.numeric(result2[["SiO2_errSD"]]),
+               as.numeric(result[["SiO2_errSD%"]] * result[["SiO2"]] / 100), tolerance = 1e-6)
+})
+
+test_that("rel_to_abs: column renamed to absolute error", {
+  expect_all_false(grepl("_err.{3}%", colnames(result2), perl = TRUE))
+})
+
+test_that("rel_to_abs: error col units match base col after conversion", {
+  expect_equal(units(result2[["SiO2_errSD"]]),
+               units(result2[["SiO2"]]))
+})
+
+test_that("rel_to_abs: base col and class unchanged", {
+  expect_equal(result2[["SiO2"]], result[["SiO2"]])
+  expect_s3_class(result2, "ASTR")
+})
+
+test_that("rel_to_abs: input returned unchanged if no relative error cols present", {
+  df_no_err <- test_input[, !is_err_percent(colnames(test_input))]
+  expect_identical(rel_to_abs(df_no_err), df_no_err)
+})
+
+# test not valid because column does no include NAs
+#test_that("rel_to_abs: NA in base col produces NA in error col, no crash", {
+#  na_rows <- is.na(as.numeric(test_input[["SiO2"]]))
+#  expect_true(all(is.na(as.numeric(result[["SiO2_errSD"]])[na_rows])))
+#})
+
+test_that("rel_to_abs: errors on non-ASTR input", {
+  expect_error(rel_to_abs(data.frame(x = 1)))
+})
+
+# --- conversion is reversible ---
+
+error_reversed <- rel_to_abs(abs_to_rel(test_input))
+
+test_that("error conversion is reversible and unitless values handled correctly", {
+  expect_equal(test_input[["SiO_err2SD"]], error_reversed[["SiO2_err2SD"]])
+  expect_equal(test_input[["d65Cu_err2SD"]], error_reversed[["d65Cu_err2SD"]])
+  expect_equal(test_input[["206Pb/204Pb_err2SD"]], error_reversed[["206Pb/204Pb_err2SD"]])
+})