New package version that can handle mass and volume separately

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: pyroresp
-Version: 0.1.1.9001
+Version: 0.1.1.9002
 Title: Analyse respirometry data captured by firesting loggers
 Authors@R: person("Hugo", "Flávio", role = c("aut", "cre"), email = "hflavio@wlu.ca", comment = c(ORCID = "0000-0002-5174-1197"))
 Description: Designed to analyse intermittent-flow respirometry data. This is a sibling package of the FishResp R package.
@@ -8,7 +8,7 @@ BugReports: https://github.com/hugomflavio/pyroresp/issues
 License: GPL-3
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Language: en-GB
 NeedsCompilation: no
 VignetteBuilder: knitr,

NAMESPACE

@@ -18,7 +18,7 @@ export(filter_r2)
 export(killen_summary)
 export(load_experiment)
 export(load_phases)
-export(load_pyro_data)
+export(load_pyro_folder)
 export(melt_resp)
 export(patch_NAs)
 export(plot_bg)

NEWS.md

@@ -2,6 +2,12 @@
 
 Find out the main highlights of each update.
 
+## pyroresp dev
+
+Enhancements:
+  * Separate animal mass and volume into different variables so the user may define those individually.
+  * Calculate background progression as a function of time rather than a function of measurement cycle: Fixes under/over estimations when bg doesn't perfectly bracket the data being analysed.
+
 ## pyroresp 0.1.1
 
 Version used for microtag respirometry paper

R/aux_functions.R

@@ -263,3 +263,73 @@ auc <- function(x, y, zero = 0) {
 
 	return(df)
 }
+
+#' needs a better place in the future.
+#' needs to be exported in the future.
+#' 
+#' @keywords internal
+#' 
+process_probe_info <- function(input, vol_unit = "ml", mass_unit = "g") {
+  required_cols <- c("id", "chamber_vol", "probe")
+  cols_missing <- !(required_cols %in% colnames(input))
+  if (any(cols_missing)) {
+    stop("The following required columns are missing ",
+       "from the input: ",
+       paste0(required_cols[cols_missing], collapse = ", "),
+       call. = FALSE)
+  }
+  if (any(is.na(input[, required_cols]))) {
+  	stop("NAs found in required probe_info columns", call. = FALSE)
+  }
+  if ("animal_vol" %in% colnames(input) &
+      !"animal_mass" %in% colnames(input)) {
+    warning("Column 'animal_mass' not found in the input.",
+            " Won't be able to calculate mass-corrected MO2.",
+            immediate. = TRUE, call. = FALSE)
+  }
+  if (!"animal_vol" %in% colnames(input) &
+      "animal_mass" %in% colnames(input)) {
+    warning("Column 'animal_vol' not found in the input.",
+            " Will assume animal density is 1.",
+            immediate. = TRUE, call. = FALSE)
+    input$animal_vol <- input$animal_mass
+  }
+  if (!"animal_vol" %in% colnames(input) &
+      !"animal_mass" %in% colnames(input)) {
+    warning("Neither 'animal_vol' nor 'animal_mass' found in the",
+            " input Won't be able to correct chamber volume",
+            " nor calculate mass-corrected MO2.",
+            immediate. = TRUE, call. = FALSE)
+  }
+  if ("animal_mass" %in% colnames(input) &
+  	any(is.na(input$animal_mass))) {
+  	stop("probe_info contains animal_mass but some data is missing.",
+  			 call. = FALSE)
+  }
+  if ("animal_vol" %in% colnames(input) &
+  	any(is.na(input$animal_vol))) {
+  	stop("probe_info contains animal_vol but some data is missing.",
+  			 call. = FALSE)
+  }
+
+  units(input$chamber_vol) <- vol_unit
+
+  check <- c("animal_vol", "animal_mass") %in% colnames(input)
+  if (check[1]) {
+    units(input$animal_vol) <- vol_unit
+    input$water_vol <- input$chamber_vol - input$animal_vol
+    input$volvol_ratio <- input$water_vol / input$animal_vol
+  } else {
+    input$water_vol <- input$chamber_vol
+  }
+
+  if (check[2]) {
+    units(input$animal_mass) <- mass_unit
+    input$volmass_ratio <- input$water_vol / conv_w_to_ml(input$animal_mass)
+  }
+
+  if (all(check)) {
+    input$animal_density <- input$animal_mass / input$animal_vol
+  }    
+  return(input)
+}

R/bg_functions.R

@@ -11,7 +11,7 @@
 #'
 #' @export
 #'
-calc_bg <- function(input, method = c('mean', 'first', 'last')){
+calc_bg <- function(input, method = c('mean', 'first', 'last')) {
 
   method <- match.arg(method)
 
@@ -44,7 +44,8 @@ calc_bg <- function(input, method = c('mean', 'first', 'last')){
     bg_lm$coefficients[1] <- 0
 
     output <- data.frame(slope = bg_lm$coefficients[2],
-                         R2 = summary(bg_lm)$adj.r.squared)
+                         R2 = summary(bg_lm)$adj.r.squared,
+                         date_time = tail(trimmed$date_time, 1))
     # done
     return(output)
   })
@@ -151,7 +152,7 @@ replace_bg <- function(input, replace, with) {
 subtract_bg <- function(input, pre, post,
                         method = c("pre", "post", "average",
                                    "linear", "parallel", "none"),
-                        ref_probe){
+                        ref_probe) {
 
   method <- match.arg(method)
 
@@ -201,7 +202,8 @@ subtract_bg <- function(input, pre, post,
       }
       link <- match(input$probe_info$ref, input$probe_info$probe)
       if (any(is.na(link))) {
-        stop("method = 'parallel' but not all values in probe_info$ref match probe names")
+        stop("method = 'parallel' but not all values in probe_info$ref",
+             " match probe names")
       }
     }
 
@@ -302,25 +304,31 @@ calc_linear_bg <- function(input, pre, post) {
     # cat(probe, "\n")
     pre_slope <- pre$bg$slope[pre$bg$probe == probe]
     post_slope <- post$bg$slope[post$bg$probe == probe]
+    pre_time <- pre$bg$date_time[pre$bg$probe == probe]
+    post_time <- post$bg$date_time[post$bg$probe == probe]
 
     slope_diff <- post_slope - pre_slope
+    time_diff <- as.numeric(difftime(post_time, pre_time, units = "s"))
+    units(time_diff) <- "s"
+    bg_progression <- slope_diff/time_diff
 
     if (as.numeric(slope_diff) == 0) {
       stop("The pre-bg and post-bg are exactly the same!",
            " Cannot calculate linear progression.", call. = FALSE)
     }
 
     # linear slope_incr.
-    slope_incr <- slope_diff / (cycles - 1)
-
-    slope_bg <- seq(from = pre_slope,
-                   to = post_slope,
-                   by = slope_incr)
+    the_slopes <- input$slopes[input$slopes$probe == probe, ]
 
     output <- data.frame(probe = probe,
-                         cycle = 1:cycles,
-                         slope_bg = slope_bg)
-
+                         date_time = the_slopes$date_time,
+                         cycle = the_slopes$cycle)
+
+    output$time_since_bg <- as.numeric(difftime(output$date_time,
+                                                pre_time,
+                                                units = "s"))
+    units(output$time_since_bg) <- "s"
+    output$slope_bg <- pre_slope + bg_progression * output$time_since_bg
     return(output)
   })
   output <- do.call(rbind, my_bg)

R/epoc_functions.R

@@ -60,11 +60,18 @@ calc_auc <- function(mmr, smr, smr_col, smr_buffer = 0.1) {
     the_probe <- the_mr$probe[1]
     the_smr <- smr$smr[smr$smr$probe == the_probe, smr_col]
 
-    the_auc <- auc(y = the_mr$mr_over_smr,
-                   x = the_mr$time_since_chase,
-                   zero = the_smr * smr_buffer)
-    colnames(the_auc)[1:2] <- c("time_delta", "mr_delta")
-    the_auc$probe <- the_probe
+    if (nrow(the_mr) > 1) {
+      the_auc <- auc(y = the_mr$mr_over_smr,
+                     x = the_mr$time_since_chase,
+                     zero = the_smr * smr_buffer)
+      colnames(the_auc)[1:2] <- c("time_delta", "mr_delta")
+      the_auc$probe <- the_probe
+    } else {
+      warning("Not enough MR points to calculate AUC for probe ",
+              the_mr$probe[1], ". Skipping.", call. = FALSE,
+              immediate. = TRUE)
+      the_auc <- NULL
+    }
     return(the_auc)
   })
 
@@ -83,6 +90,8 @@ calc_auc <- function(mmr, smr, smr_col, smr_buffer = 0.1) {
 #'
 #' @param input a resp object with an "auc" sub-object.
 #'   The output of \code{\link{calc_auc}}
+#' @param min_duration minimum time during which auc mas continue being counted
+#'   even if it reaches 0 before it. Useful for animals that can stop breathing.
 #' @param max_duration maximum time allowed for the auc to reach 0 (in hours).
 #'   if AUC does not reach 0 before this threshold, the cumulative AUC for the
 #'   whole duration is used.
@@ -91,24 +100,29 @@ calc_auc <- function(mmr, smr, smr_col, smr_buffer = 0.1) {
 #'
 #' @export
 #' 
-extract_epoc <- function(input, max_duration = Inf) {
+extract_epoc <- function(input, min_duration = 0, max_duration = Inf) {
   if (is.null(input$auc)) {
     stop("Couldn't find object 'auc' inside input.",
          " Have you run calc_auc?", call. = FALSE)
   }
   units(max_duration) <- units(input$auc$time_delta)
+  units(min_duration) <- units(input$auc$time_delta)
 
   by_probe <- split(input$auc, input$auc$probe)
 
   recipient <- lapply(by_probe, function(x) {
-    x <- x[x$time_delta <= max_duration, ]
+    x <- x[x$time_delta >= min_duration & x$time_delta <= max_duration, ]
     epoc_ended <- min(which(as.numeric(x$mr_delta) < 0), nrow(x))
     output <- data.frame(probe = x$probe[1],
                          epoc = x$cumauc[epoc_ended],
                          epoc_duration = x$time_delta[epoc_ended])
-    this_probe <- input$mr$probe == x$probe[1]
-    this_time <- input$mr$time_since_chase == x$time_delta[epoc_ended]
-    output$fold_over_smr_at_end <- input$mr$fold_over_smr[this_probe & this_time]
+    if (drop_units(output$epoc_duration) == 0) {
+      output$fold_over_smr_at_end <- NA
+    } else {
+      this_probe <- input$mr$probe == x$probe[1]
+      this_time <- input$mr$time_since_chase == x$time_delta[epoc_ended]
+      output$fold_over_smr_at_end <- input$mr$fold_over_smr[this_probe & this_time]
+    }
     return(output)
   })
   the_epoc <- do.call(rbind, recipient)

R/killen_summary.R

@@ -14,9 +14,8 @@
 #' @return a list of summary information
 #'
 killen_summary <- function(input) {
-  # 04
-  pt04 <- with(input$probe_info,
-               (volume - conv_w_to_ml(mass)) / conv_w_to_ml(mass))
+  # 04: water:mass ratio
+  pt04 <- input$probe_info$volmass_ratio
   # 15
   recipient <- lapply(input$pyro$source_data, function(i) {
     prb <- data.frame(probe = paste0(attributes(i)$device, attributes(i)$ch))
@@ -85,9 +84,9 @@ killen_summary <- function(input) {
   # bring it all together
   output <- list(
     pt01 = list(description = "mass of the animals",
-                value = input$probe_info$mass),
+                value = input$probe_info$animal_mass),
     pt02 = list(description = "volume of empty respirometer",
-                value = input$probe_info$volume),
+                value = input$probe_info$chamber_vol),
     pt04 = list(description = "ratio of net resp volume",
                 value = pt04),
     pt14 = list(description = "wait time (in number of data points)",

R/melt_resp.R

@@ -80,11 +80,8 @@ melt_resp <- function(input) {
   if (!is.null(input$probe_info)) {
     # Include the extra information
     link <- match(pre_output$probe, input$probe_info$probe)
-    if ("mass" %in% colnames(input$probe_info)) {
-      to_transfer <- c("id", "mass", "volume")
-    } else {
-      to_transfer <- c("id", "volume")
-    }
+    to_transfer <- c("id", "animal_mass", "water_vol")
+    to_transfer <- to_transfer[to_transfer %in% colnames(input$probe_info)]
     output <- cbind(pre_output,
                     input$probe_info[link, to_transfer])
   } else {

R/mmr_functions.R

@@ -10,20 +10,17 @@ extract_mmr <- function(mr){
   by_probe <- split(mr, mr$probe)
 
   recipient <- lapply(by_probe, function(the_probe) {
-      if ("mass" %in% colnames(the_probe)) {
-        index <- head(order(the_probe$mr_g, decreasing = TRUE), 1)
-        the_probe$mr <- the_probe$mr_g
-      } else {
-        index <- head(order(the_probe$mr_abs, decreasing = TRUE), 1)
-        the_probe$mr <- the_probe$mr_abs
-      }
-      output <- the_probe[index, ]
-    })
+    target <- ifelse ("animal_mass" %in% colnames(the_probe),
+                      "mr_g", "mr_abs")
+    the_probe$mmr <- the_probe[, target]
+    index <- order(the_probe$mmr, decreasing = TRUE)[1]
+    output <- the_probe[index, ]
+  })
 
   mmr <- as.data.frame(data.table::rbindlist(recipient))
 
   # Keep only needed columns
-    mmr <- mmr[, c("probe", "date_time", "cycle", "mr")]
-  
+    mmr <- mmr[, c("probe", "date_time", "cycle", "mmr")]
+
   return(mmr)
 }