diff --git a/DESCRIPTION b/DESCRIPTION index b63ffd9..d584679 100644 --- a/DESCRIPTION +++ b/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, diff --git a/NAMESPACE b/NAMESPACE index 10fca6c..8a6e704 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index 452a0a3..cae781d 100644 --- a/NEWS.md +++ b/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 diff --git a/R/aux_functions.R b/R/aux_functions.R index 6df4c86..48a1d04 100644 --- a/R/aux_functions.R +++ b/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) +} diff --git a/R/bg_functions.R b/R/bg_functions.R index aa89960..a2997a9 100644 --- a/R/bg_functions.R +++ b/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,8 +304,13 @@ 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!", @@ -311,16 +318,17 @@ calc_linear_bg <- function(input, pre, post) { } # 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) diff --git a/R/epoc_functions.R b/R/epoc_functions.R index 36179e6..61cd7b3 100644 --- a/R/epoc_functions.R +++ b/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) diff --git a/R/killen_summary.R b/R/killen_summary.R index 98b22cb..df270f2 100644 --- a/R/killen_summary.R +++ b/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)", diff --git a/R/melt_resp.R b/R/melt_resp.R index af874c9..bca590e 100644 --- a/R/melt_resp.R +++ b/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 { diff --git a/R/mmr_functions.R b/R/mmr_functions.R index c6e93cb..0eed961 100644 --- a/R/mmr_functions.R +++ b/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) } diff --git a/R/mr_functions.R b/R/mr_functions.R index 933291f..a2c413d 100644 --- a/R/mr_functions.R +++ b/R/mr_functions.R @@ -1,33 +1,25 @@ #' Calculate metabolic rates #' -#' Calculates the absolute and mass-specific metabolic rates, -#' as well as the relative importance of background (in percentage). +#' Calculates absolute and mass-specific metabolic rates. #' -#' @param slope_data a data frame obtained by using either the -#' function \code{\link{calc_slopes}} or \code{\link{filter_r2}} -#' @inheritParams conv_w_to_ml +#' @param input a data frame with columns slope_cor (mandatory), +#' water_vol (mandatory), and animal_mass (optional). #' -#' @return A data frame with mr_abs, bg, and mr_cor. +#' @return The input with columns mr_abs and mr_g if animal_mass provided. #' #' @export #' -calc_mr <- function(slope_data, density = 1){ - if ("mass" %in% colnames(slope_data)) { - vol <- slope_data$volume - conv_w_to_ml(slope_data$mass, density) - slope_data$mr_abs <- -(slope_data$slope_cor * vol) - - # temporarily drop and reassign units. - # this is needed to avoid {units} merging oxygen and animal - # weight when O2 is measured in weight (e.g. mg). +calc_mr <- function(input){ + input$mr_abs <- -(input$slope_cor * input$water_vol) + if ("animal_mass" %in% colnames(input)) { # see: https://github.com/r-quantities/units/issues/411 - slope_data$mr_g <- drop_units(slope_data$mr_abs) / drop_units(slope_data$mass) - units(slope_data$mr_g) <- paste0(units(slope_data$mr_abs), - "/", units(slope_data$mass)) - } else { - slope_data$mr_abs <- -slope_data$slope_cor + prev_option <- units_options("simplify") + units_options(simplify = FALSE) + input$mr_g <- input$mr_abs / input$animal_mass + units_options(simplify = prev_option) } - return(slope_data) + return(input) } #' Break down MR from a single cycle @@ -77,6 +69,8 @@ roll_mr <- function(input, probe, cycle, smoothing, density = 1, r2 = 0.95) { data = this_data[(i - smoothing + 1):i, ]) output <- data.frame(id = this_id, probe = probe, + animal_mass = this_data$animal_mass[1], + water_vol = this_data$water_vol[1], cycle = cycle, phase_time = i, smoothing = smoothing, @@ -98,9 +92,7 @@ roll_mr <- function(input, probe, cycle, smoothing, density = 1, r2 = 0.95) { output$slope_cor <- output$slope - output$slope_bg # convert to metabolic rate - the_vol <- this_data$volume[1] - conv_w_to_ml(this_data$mass[1], density) - output$mr_abs <- -(output$slope_cor * the_vol) - output$mr_g <- output$mr_abs / this_data$mass[1] + output <- calc_mr(output) # look out for bad slopes r2_link <- output$r2 >= r2 @@ -111,7 +103,12 @@ roll_mr <- function(input, probe, cycle, smoothing, density = 1, r2 = 0.95) { } # find the max value trim_mr <- output[r2_link, ] - index <- which.max(trim_mr$mr_g) + + if ("mr_g" %in% colnames(trim_mr)) { + index <- which.max(trim_mr$mr_g) + } else { + index <- which.max(trim_mr$mr_abs) + } max_mr <- trim_mr[index, , drop = FALSE] if (is.null(input$rolling_mr)) { diff --git a/R/o2_calcs_functions.R b/R/o2_calcs_functions.R index 20d9efe..1cb8686 100644 --- a/R/o2_calcs_functions.R +++ b/R/o2_calcs_functions.R @@ -45,7 +45,7 @@ trim_resp <- function(input, meas_max = Inf, meas_min = 60, link <- !(target_probes %in% names(first_cycle)) if (any(link)) { stop("Length of argument 'first_cycle' is not 1 but named values", - " not supplied for all probes (missing", + " not supplied for all probes (missing ", paste0(target_probes[link], collapse = ", "), ").", call. = FALSE) } diff --git a/R/plot_functions.R b/R/plot_functions.R index 5ad05e3..e5a7792 100644 --- a/R/plot_functions.R +++ b/R/plot_functions.R @@ -1,7 +1,7 @@ #' Plot background readings #' #' @param input An experiment list with calculated background. -#' The output of \code{\link{calc_bg}}. +#' The output of \code{\link{calc_bg}}. #' @param linewidth The width of the averaged/linear bg. #' @inheritParams plot_meas #' @@ -10,58 +10,58 @@ #' @export #' plot_bg <- function(input, probes, linewidth = 1.5) { - # ggplot variables - o2_delta <- cycle <- o2_bg_delta <- phase_time <- NULL - - if(is.null(input$bg)) { - stop("Could not find a bg object in the input") - } - - if (!missing(probes)) { - probes <- check_arg_in_data(probes, input$trimmed$probe, "probes") - input$trimmed <- input$trimmed[input$trimmed$probe %in% probes, ] - input$bg <- input$bg[input$bg$probe %in% probes, ] - } - - input$trimmed$cycle <- as.factor(input$trimmed$cycle) - input$probe <- factor(input$idprobe, levels = input$probe_info$probe) - - bg_lines <- lapply(1:nrow(input$bg), function(i) { - probe <- input$bg$probe[i] - link <- input$trimmed$probe == probe - nsecs <- max(input$trimmed$phase_time[link]) - output <- data.frame(probe = probe, - phase_time = c(0, nsecs), - o2_delta = c(0, input$bg$slope[i] * nsecs)) - units(output$phase_time) <- units(input$trimmed$phase_time) - units(output$o2_delta) <- units(input$trimmed$o2_delta) - return(output) - }) - bg_lines <- do.call(rbind, bg_lines) - - p <- ggplot2::ggplot(data = input$trimmed) - p <- p + ggplot2::aes(x = phase_time, y = o2_delta) - - p <- p + ggplot2::geom_line(ggplot2::aes(group = cycle, colour = cycle)) - p <- p + ggplot2::geom_line(data = bg_lines, colour = 'red', - linewidth = linewidth) - # if (!is.null(input$sim_params)) { - # p <- plot_sim_watermark(p = p, - # x = input$trimmed$phase_time, - # y = input$trimmed$o2_delta) - # } - p <- p + ggplot2::theme_bw() - p <- p + ggplot2::labs(x = "Phase time", y = "Delta~O[2]") - p <- p + ggplot2::facet_wrap(. ~ probe, ncol = 4) - - p + # ggplot variables + o2_delta <- cycle <- o2_bg_delta <- phase_time <- NULL + + if(is.null(input$bg)) { + stop("Could not find a bg object in the input") + } + + if (!missing(probes)) { + probes <- check_arg_in_data(probes, input$trimmed$probe, "probes") + input$trimmed <- input$trimmed[input$trimmed$probe %in% probes, ] + input$bg <- input$bg[input$bg$probe %in% probes, ] + } + + input$trimmed$cycle <- as.factor(input$trimmed$cycle) + input$probe <- factor(input$idprobe, levels = input$probe_info$probe) + + bg_lines <- lapply(1:nrow(input$bg), function(i) { + probe <- input$bg$probe[i] + link <- input$trimmed$probe == probe + nsecs <- max(input$trimmed$phase_time[link]) + output <- data.frame(probe = probe, + phase_time = c(0, nsecs), + o2_delta = c(0, input$bg$slope[i] * nsecs)) + units(output$phase_time) <- units(input$trimmed$phase_time) + units(output$o2_delta) <- units(input$trimmed$o2_delta) + return(output) + }) + bg_lines <- do.call(rbind, bg_lines) + + p <- ggplot2::ggplot(data = input$trimmed) + p <- p + ggplot2::aes(x = phase_time, y = o2_delta) + + p <- p + ggplot2::geom_line(ggplot2::aes(group = cycle, colour = cycle)) + p <- p + ggplot2::geom_line(data = bg_lines, colour = 'red', + linewidth = linewidth) + # if (!is.null(input$sim_params)) { + # p <- plot_sim_watermark(p = p, + # x = input$trimmed$phase_time, + # y = input$trimmed$o2_delta) + # } + p <- p + ggplot2::theme_bw() + p <- p + ggplot2::labs(x = "Phase time", y = "Delta~O[2]") + p <- p + ggplot2::facet_wrap(. ~ probe, ncol = 4) + + p } #' Plot raw measurements of oxygen and temp #' #' @param input An experiment list with melted oxygen and temperature data. -#' The output of \code{\link{process_experiment}} or -#' \code{\link{melt_resp}}. +#' The output of \code{\link{process_experiment}} or +#' \code{\link{melt_resp}}. #' @param cycles A numeric vector of which cycles to plot #' @param probes A string of which probes to plot #' @param type Should the plot show only "lines", "points", or "both"? @@ -74,155 +74,160 @@ plot_bg <- function(input, probes, linewidth = 1.5) { #' plot_meas <- function(input, cycles, probes, type = c("lines", "points", "both"), - show_temp = FALSE, verbose = TRUE) { - # ggplot variables - phase <- date_time <- temp <- o2 <- probe <- NULL - xmin <- xmax <- ymin <- ymax <- NULL - - type <- match.arg(type) - - meas <- input$melted - meas$idprobe <- paste0(meas$id, " (", meas$probe, ")") - idprobe_levels <-paste0(input$probe_info$id, - " (", input$probe_info$probe, ")") - meas$idprobe <- factor(meas$idprobe, levels = idprobe_levels) - - # strip units to avoid conflicts - o2_unit <- units(meas$o2) - meas$o2 <- as.numeric(meas$o2) - temp_unit <- units(meas$temp) - meas$temp <- as.numeric(meas$temp) - - if (!missing(cycles)) { - cycles <- check_arg_in_data(cycles, meas$cycle, - "cycles", verbose = verbose) - meas <- meas[meas$cycle %in% cycles, ] - } - - if (!missing(probes)) { - probes <- check_arg_in_data(probes, meas$probe, - "probes", verbose = verbose) - meas <- meas[meas$probe %in% probes, ] - meas$idprobe <- droplevels(meas$idprobe) - } - - if (any(grepl("F|W", meas$phase))) { - paint_flush <- TRUE - - aux <- split(meas, meas$idprobe) - - recipient <- lapply(names(aux), function(the_idprobe) { - # instead of using split, I need to manually break these - # because F0 can appear multiple times. split would - # group them together and mess up the plot - phase_info <- rle(aux[[the_idprobe]]$phase) - row_breaks <- c(1, cumsum(phase_info$lengths)) - aux2 <- lapply(2:length(row_breaks), function(i) { - meas[row_breaks[i-1]:row_breaks[i], ] - }) - # -- - names(aux2) <- phase_info$values - aux2 <- aux2[grepl("F|W", names(aux2))] - # now extract start and end of phase - aux2 <- lapply(aux2, function(the_phase) { - data.frame(idprobe = the_idprobe, - xmin = the_phase$date_time[1], - xmax = the_phase$date_time[nrow(the_phase)], - ymin = -Inf, - ymax = Inf) - }) - return(as.data.frame(data.table::rbindlist(aux2))) - }) - rm(aux) - - flush_times <- as.data.frame(data.table::rbindlist(recipient)) - } else { - paint_flush <- FALSE - } - - - p <- ggplot2::ggplot(data = meas) - - if (paint_flush) { - p <- p + ggplot2::geom_rect(data = flush_times, - ggplot2::aes(xmin = xmin, xmax = xmax, - ymin = ymin, ymax = ymax), - fill = "black", alpha = 0.1) - } - - if (show_temp) { - aux <- range(meas$o2, na.rm = TRUE) - ox_range <- aux[2] - aux[1] - ox_mid <- mean(aux) - - aux <- range(meas$temp, na.rm = TRUE) - temp_range <- aux[2] - aux[1] - temp_mid <- mean(aux) - - mid_dif <- temp_mid - ox_mid - - if (ox_range == 0 | temp_range == 0) { - compress_factor <- 1 - } else { - compress_factor <- ox_range/temp_range - } - - data.link <- function(x, a = temp_mid, - b = ox_mid, c_factor = compress_factor) { - b + ((x - a) * c_factor) # = y - } - - axis_link <- function(y, a = temp_mid, - b = ox_mid, c_factor = compress_factor) { - (y + a * c_factor - b) / c_factor # = x - } - - if (type %in% c("lines", "both")) { - p <- p + ggplot2::geom_line(ggplot2::aes(x = date_time, - y = data.link(temp), - group = phase, - colour = "temp")) - } - if (type %in% c("points", "both")) { - p <- p + ggplot2::geom_point(ggplot2::aes(x = date_time, - y = data.link(temp), - group = phase, - colour = "temp")) - } - p <- p + ggplot2::scale_y_continuous( - sec.axis = ggplot2::sec_axis(trans = axis_link, - name = units::make_unit_label("Temperature", temp_unit)) - ) - p <- p + ggplot2::scale_colour_manual(values = c("royalblue", "red")) - } else { - p <- p + ggplot2::scale_colour_manual(values = "royalblue") - } - - if (type %in% c("lines", "both")) { - p <- p + ggplot2::geom_line(ggplot2::aes(x = date_time, - y = o2, group = phase, - colour = "o2")) - } - if (type %in% c("points", "both")) { - p <- p + ggplot2::geom_point(ggplot2::aes(x = date_time, - y = o2, group = phase, - colour = "o2")) - } - p <- p + ggplot2::theme_bw() - - p <- p + ggplot2::theme(legend.position = "none") - p <- p + ggplot2::labs(y = units::make_unit_label("O[2]", o2_unit), - x = "Time") - p <- p + ggplot2::facet_wrap(idprobe ~ ., ncol = 1) - - return(p) + show_temp = FALSE, verbose = TRUE) { + # ggplot variables + phase <- date_time <- temp <- o2 <- probe <- NULL + xmin <- xmax <- ymin <- ymax <- NULL + + type <- match.arg(type) + + meas <- input$melted + if (is.null(input$probe_info)) { + meas$idprobe <- meas$probe + idprobe_levels <- sort(unique(meas$probe)) + } else { + meas$idprobe <- paste0(meas$id, " (", meas$probe, ")") + idprobe_levels <- paste0(input$probe_info$id, + " (", input$probe_info$probe, ")") + } + meas$idprobe <- factor(meas$idprobe, levels = idprobe_levels) + + # strip units to avoid conflicts + o2_unit <- units(meas$o2) + meas$o2 <- as.numeric(meas$o2) + temp_unit <- units(meas$temp) + meas$temp <- as.numeric(meas$temp) + + if (!missing(cycles)) { + cycles <- check_arg_in_data(cycles, meas$cycle, + "cycles", verbose = verbose) + meas <- meas[meas$cycle %in% cycles, ] + } + + if (!missing(probes)) { + probes <- check_arg_in_data(probes, meas$probe, + "probes", verbose = verbose) + meas <- meas[meas$probe %in% probes, ] + meas$idprobe <- droplevels(meas$idprobe) + } + + if (any(grepl("F|W", meas$phase))) { + paint_flush <- TRUE + + aux <- split(meas, meas$idprobe) + + recipient <- lapply(names(aux), function(the_idprobe) { + # instead of using split, I need to manually break these + # because F0 can appear multiple times. split would + # group them together and mess up the plot + phase_info <- rle(aux[[the_idprobe]]$phase) + row_breaks <- c(1, cumsum(phase_info$lengths)) + aux2 <- lapply(2:length(row_breaks), function(i) { + meas[row_breaks[i-1]:row_breaks[i], ] + }) + # -- + names(aux2) <- phase_info$values + aux2 <- aux2[grepl("F|W", names(aux2))] + # now extract start and end of phase + aux2 <- lapply(aux2, function(the_phase) { + data.frame(idprobe = the_idprobe, + xmin = the_phase$date_time[1], + xmax = the_phase$date_time[nrow(the_phase)], + ymin = -Inf, + ymax = Inf) + }) + return(as.data.frame(data.table::rbindlist(aux2))) + }) + rm(aux) + + flush_times <- as.data.frame(data.table::rbindlist(recipient)) + } else { + paint_flush <- FALSE + } + + + p <- ggplot2::ggplot(data = meas) + + if (paint_flush) { + p <- p + ggplot2::geom_rect(data = flush_times, + ggplot2::aes(xmin = xmin, xmax = xmax, + ymin = ymin, ymax = ymax), + fill = "black", alpha = 0.1) + } + + if (show_temp) { + aux <- range(meas$o2, na.rm = TRUE) + ox_range <- aux[2] - aux[1] + ox_mid <- mean(aux) + + aux <- range(meas$temp, na.rm = TRUE) + temp_range <- aux[2] - aux[1] + temp_mid <- mean(aux) + + mid_dif <- temp_mid - ox_mid + + if (ox_range == 0 | temp_range == 0) { + compress_factor <- 1 + } else { + compress_factor <- ox_range/temp_range + } + + data.link <- function(x, a = temp_mid, + b = ox_mid, c_factor = compress_factor) { + b + ((x - a) * c_factor) # = y + } + + axis_link <- function(y, a = temp_mid, + b = ox_mid, c_factor = compress_factor) { + (y + a * c_factor - b) / c_factor # = x + } + + if (type %in% c("lines", "both")) { + p <- p + ggplot2::geom_line(ggplot2::aes(x = date_time, + y = data.link(temp), + group = phase, + colour = "temp")) + } + if (type %in% c("points", "both")) { + p <- p + ggplot2::geom_point(ggplot2::aes(x = date_time, + y = data.link(temp), + group = phase, + colour = "temp")) + } + p <- p + ggplot2::scale_y_continuous( + sec.axis = ggplot2::sec_axis(trans = axis_link, + name = units::make_unit_label("Temperature", temp_unit)) + ) + p <- p + ggplot2::scale_colour_manual(values = c("royalblue", "red")) + } else { + p <- p + ggplot2::scale_colour_manual(values = "royalblue") + } + + if (type %in% c("lines", "both")) { + p <- p + ggplot2::geom_line(ggplot2::aes(x = date_time, + y = o2, group = phase, + colour = "o2")) + } + if (type %in% c("points", "both")) { + p <- p + ggplot2::geom_point(ggplot2::aes(x = date_time, + y = o2, group = phase, + colour = "o2")) + } + p <- p + ggplot2::theme_bw() + + p <- p + ggplot2::theme(legend.position = "none") + p <- p + ggplot2::labs(y = units::make_unit_label("O[2]", o2_unit), + x = "Time") + p <- p + ggplot2::facet_wrap(idprobe ~ ., ncol = 1) + + return(p) } #' Plot the calculated oxygen deltas #' #' @param input An experiment list with calculated deltas. -#' The output of \code{\link{process_experiment}} or -#' \code{\link{calc_delta}}. +#' The output of \code{\link{process_experiment}} or +#' \code{\link{calc_delta}}. #' @inheritParams plot_meas #' #' @return a ggplot object @@ -230,50 +235,50 @@ plot_meas <- function(input, cycles, probes, #' @export #' plot_deltas <- function(input, cycles, probes, verbose = TRUE) { - # ggplot variables - o2_delta <- o2_cordelta <- o2_bg_delta <- date_time <- NULL - phase <- plot_this <- NULL - - trimmed <- input$trimmed - trimmed$idprobe <- paste0(trimmed$id, " (", trimmed$probe, ")") - idprobe_levels <-paste0(input$probe_info$id, - " (", input$probe_info$probe, ")") - trimmed$idprobe <- factor(trimmed$idprobe, levels = idprobe_levels) - - if (!missing(cycles)) { - cycles <- check_arg_in_data(cycles, trimmed$cycle, - "cycles", verbose = verbose) - trimmed <- trimmed[trimmed$cycle %in% cycles, ] - } - - if (!missing(probes)) { - probes <- check_arg_in_data(probes, trimmed$probe, - "probes", verbose = verbose) - trimmed <- trimmed[trimmed$probe %in% probes, ] - trimmed$idprobe <- droplevels(trimmed$idprobe) - } - - p <- ggplot2::ggplot(data = trimmed, - ggplot2::aes(x = date_time, Group = phase)) - - # if (!is.null(input$sim_params)) { - # p <- plot_sim_watermark(p = p, - # x = trimmed$date_time, - # y = trimmed$o2_delta) - # } - p <- p + ggplot2::geom_path(ggplot2::aes(y = o2_delta), col = "royalblue") - p <- p + ggplot2::theme_bw() - p <- p + ggplot2::facet_wrap(.~idprobe) - p <- p + ggplot2::labs(x = "", y = expression(Delta~O[2])) - - return(p) + # ggplot variables + o2_delta <- o2_cordelta <- o2_bg_delta <- date_time <- NULL + phase <- plot_this <- NULL + + trimmed <- input$trimmed + trimmed$idprobe <- paste0(trimmed$id, " (", trimmed$probe, ")") + idprobe_levels <-paste0(input$probe_info$id, + " (", input$probe_info$probe, ")") + trimmed$idprobe <- factor(trimmed$idprobe, levels = idprobe_levels) + + if (!missing(cycles)) { + cycles <- check_arg_in_data(cycles, trimmed$cycle, + "cycles", verbose = verbose) + trimmed <- trimmed[trimmed$cycle %in% cycles, ] + } + + if (!missing(probes)) { + probes <- check_arg_in_data(probes, trimmed$probe, + "probes", verbose = verbose) + trimmed <- trimmed[trimmed$probe %in% probes, ] + trimmed$idprobe <- droplevels(trimmed$idprobe) + } + + p <- ggplot2::ggplot(data = trimmed, + ggplot2::aes(x = date_time, Group = phase)) + + # if (!is.null(input$sim_params)) { + # p <- plot_sim_watermark(p = p, + # x = trimmed$date_time, + # y = trimmed$o2_delta) + # } + p <- p + ggplot2::geom_path(ggplot2::aes(y = o2_delta), col = "royalblue") + p <- p + ggplot2::theme_bw() + p <- p + ggplot2::facet_wrap(.~idprobe) + p <- p + ggplot2::labs(x = "", y = expression(Delta~O[2])) + + return(p) } #' Plot the calculated slopes #' #' @param input An experiment list with calculated slopes. -#' The output of \code{\link{process_experiment}} or -#' \code{\link{calc_slopes}}. +#' The output of \code{\link{process_experiment}} or +#' \code{\link{calc_slopes}}. #' @param r2 show the respective R2 for each slope? #' @inheritParams plot_mr #' @@ -282,105 +287,105 @@ plot_deltas <- function(input, cycles, probes, verbose = TRUE) { #' @export #' plot_slopes <- function(input, cycles, probes, r2 = TRUE, verbose = TRUE) { - # ggplot variables - date_time <- slope_cor <- valid <- NULL - - slopes <- input$slopes - slopes$idprobe <- paste0(slopes$id, " (", slopes$probe, ")") - idprobe_levels <-paste0(input$probe_info$id, - " (", input$probe_info$probe, ")") - slopes$idprobe <- factor(slopes$idprobe, levels = idprobe_levels) - - # remove units to make the double y-axis easy - slope_unit <- units(slopes$slope_cor) - slopes$slope_cor <- as.numeric(slopes$slope_cor) - r2_threshold <- attributes(input$good_slopes)$r2_threshold - slopes$valid <- slopes$r2 > r2_threshold - slopes$valid[is.na(slopes$valid)] <- FALSE - - if (!missing(cycles)) { - cycles <- check_arg_in_data(cycles, slopes$cycle, - "cycles", verbose = verbose) - slopes <- slopes[slopes$cycle %in% cycles, ] - } - - if (!missing(probes)) { - probes <- check_arg_in_data(probes, slopes$probe, - "probes", verbose = verbose) - slopes <- slopes[slopes$probe %in% probes, ] - slopes$idprobe <- droplevels(slopes$idprobe) - } - - p <- ggplot2::ggplot(data = slopes, ggplot2::aes(x = date_time)) - - # if (!is.null(input$sim_params)) { - # aux <- range(slopes$slope_cor, na.rm = TRUE) - # if (r2) { - # # compensate for artificial y axis raise set below - # aux[2] <- aux[1] + ((((aux[2] - aux[1]) / 2) * 1.15) * 2) - # } - # p <- plot_sim_watermark(p = p, - # x = slopes$date_time, - # y = aux) - # } - - p <- p + ggplot2::geom_line(ggplot2::aes(y = slope_cor, col = 'slope_cor')) - p <- p + ggplot2::geom_point(ggplot2::aes(y = slope_cor, col = 'slope_cor')) - - - if (r2) { - aux <- range(slopes$slope_cor, na.rm = TRUE) - # artificially raise centre so R2 tends to stay about the slopes - slope_centre <- aux[1] + (((aux[2] - aux[1]) / 2) * 1.3) - slope_range <- aux[2] - aux[1] - r2_centre <- 0.5 - r2_range <- 1 - - if (slope_range == 0 | r2_range == 0) { - compress_factor <- 1 - } else { - compress_factor <- slope_range/r2_range - } - - data.link <- function(x, a = r2_centre, b = slope_centre, - c_factor = compress_factor) { - b + ((x - a) * c_factor) # = y - } - - axis_link <- function(y, a = r2_centre, b = slope_centre, - c_factor = compress_factor) { - (y + a * c_factor - b) / c_factor # = x - } - - p <- p + ggplot2::geom_hline(yintercept = data.link(r2_threshold)) - p <- p + ggplot2::geom_line(ggplot2::aes(y = data.link(r2)), - col = "grey") - p <- p + ggplot2::geom_point(ggplot2::aes(y = data.link(r2), - col = valid)) - p <- p + ggplot2::scale_y_continuous(sec.axis = - ggplot2::sec_axis(trans = axis_link, name = "R2") - ) - if (any(!slopes$valid)) { - p <- p + ggplot2::scale_colour_manual( - values = c("Red", "Black", "Grey")) - } else { - p <- p + ggplot2::scale_colour_manual(values = c("Black", "Grey")) - } - } - - p <- p + ggplot2::theme_bw() - p <- p + ggplot2::theme(legend.position = "none") - p <- p + ggplot2::labs(y = units::make_unit_label("Slope", slope_unit), - x = "Time") - p <- p + ggplot2::facet_wrap(idprobe ~ ., ncol = 1) - return(p) + # ggplot variables + date_time <- slope_cor <- valid <- NULL + + slopes <- input$slopes + slopes$idprobe <- paste0(slopes$id, " (", slopes$probe, ")") + idprobe_levels <-paste0(input$probe_info$id, + " (", input$probe_info$probe, ")") + slopes$idprobe <- factor(slopes$idprobe, levels = idprobe_levels) + + # remove units to make the double y-axis easy + slope_unit <- units(slopes$slope_cor) + slopes$slope_cor <- as.numeric(slopes$slope_cor) + r2_threshold <- attributes(input$good_slopes)$r2_threshold + slopes$valid <- slopes$r2 > r2_threshold + slopes$valid[is.na(slopes$valid)] <- FALSE + + if (!missing(cycles)) { + cycles <- check_arg_in_data(cycles, slopes$cycle, + "cycles", verbose = verbose) + slopes <- slopes[slopes$cycle %in% cycles, ] + } + + if (!missing(probes)) { + probes <- check_arg_in_data(probes, slopes$probe, + "probes", verbose = verbose) + slopes <- slopes[slopes$probe %in% probes, ] + slopes$idprobe <- droplevels(slopes$idprobe) + } + + p <- ggplot2::ggplot(data = slopes, ggplot2::aes(x = date_time)) + + # if (!is.null(input$sim_params)) { + # aux <- range(slopes$slope_cor, na.rm = TRUE) + # if (r2) { + # # compensate for artificial y axis raise set below + # aux[2] <- aux[1] + ((((aux[2] - aux[1]) / 2) * 1.15) * 2) + # } + # p <- plot_sim_watermark(p = p, + # x = slopes$date_time, + # y = aux) + # } + + p <- p + ggplot2::geom_line(ggplot2::aes(y = slope_cor, col = 'slope_cor')) + p <- p + ggplot2::geom_point(ggplot2::aes(y = slope_cor, col = 'slope_cor')) + + + if (r2) { + aux <- range(slopes$slope_cor, na.rm = TRUE) + # artificially raise centre so R2 tends to stay about the slopes + slope_centre <- aux[1] + (((aux[2] - aux[1]) / 2) * 1.3) + slope_range <- aux[2] - aux[1] + r2_centre <- 0.5 + r2_range <- 1 + + if (slope_range == 0 | r2_range == 0) { + compress_factor <- 1 + } else { + compress_factor <- slope_range/r2_range + } + + data.link <- function(x, a = r2_centre, b = slope_centre, + c_factor = compress_factor) { + b + ((x - a) * c_factor) # = y + } + + axis_link <- function(y, a = r2_centre, b = slope_centre, + c_factor = compress_factor) { + (y + a * c_factor - b) / c_factor # = x + } + + p <- p + ggplot2::geom_hline(yintercept = data.link(r2_threshold)) + p <- p + ggplot2::geom_line(ggplot2::aes(y = data.link(r2)), + col = "grey") + p <- p + ggplot2::geom_point(ggplot2::aes(y = data.link(r2), + col = valid)) + p <- p + ggplot2::scale_y_continuous(sec.axis = + ggplot2::sec_axis(trans = axis_link, name = "R2") + ) + if (any(!slopes$valid)) { + p <- p + ggplot2::scale_colour_manual( + values = c("Red", "Black", "Grey")) + } else { + p <- p + ggplot2::scale_colour_manual(values = c("Black", "Grey")) + } + } + + p <- p + ggplot2::theme_bw() + p <- p + ggplot2::theme(legend.position = "none") + p <- p + ggplot2::labs(y = units::make_unit_label("Slope", slope_unit), + x = "Time") + p <- p + ggplot2::facet_wrap(idprobe ~ ., ncol = 1) + return(p) } #' Plot the metabolic rates #' #' @param input An experiment list with calculated metabolic rates. -#' The output of \code{\link{process_experiment}} or -#' \code{\link{calc_mr}}. +#' The output of \code{\link{process_experiment}} or +#' \code{\link{calc_mr}}. #' @inheritParams plot_meas #' #' @return a ggplot object @@ -388,111 +393,111 @@ plot_slopes <- function(input, cycles, probes, r2 = TRUE, verbose = TRUE) { #' @export #' plot_mr <- function(input, cycles, probes, verbose = TRUE) { - # ggplot variables - date_time <- mr_g <- mr_real <- value <- Method <- NULL - - mr <- input$mr - if ("mr_g" %in% colnames(mr)) { - mr$mr <- mr$mr_g - } else { - mr$mr <- mr$mr_abs - } - - mr$idprobe <- paste0(mr$id, " (", mr$probe, ")") - idprobe_levels <-paste0(input$probe_info$id, - " (", input$probe_info$probe, ")") - mr$idprobe <- factor(mr$idprobe, levels = idprobe_levels) - - if (!missing(cycles)) { - cycles <- check_arg_in_data(cycles, mr$cycle, - "cycles", verbose = verbose) - mr <- mr[mr$cycle %in% cycles, ] - } - - if (!missing(probes)) { - probes <- check_arg_in_data(probes, mr$probe, - "probes", verbose = verbose) - mr <- mr[mr$probe %in% probes, ] - mr$idprobe <- droplevels(mr$idprobe) - } - - p <- ggplot2::ggplot() - - # if (!is.null(input$sim_params)) { - # p <- plot_sim_watermark(p = p, - # x = mr$date_time, - # y = c(mr$mr_real, mr$mr)) - # mr_real_aux <- data.frame(mr_real = input$sim_params$real_mr, - # date_time = input$slopes$date_time) - # p <- p + ggplot2::geom_line(data = mr_real_aux, - # ggplot2::aes(x = date_time, y = mr_real), - # linetype = "dashed") - # p <- p + ggplot2::geom_point(data = mr_real_aux, - # ggplot2::aes(x = date_time, y = mr_real), - # shape = 1) - # } - - p <- p + ggplot2::geom_line(data = mr, - ggplot2::aes(x = date_time, y = mr)) - p <- p + ggplot2::geom_point(data = mr, - ggplot2::aes(x = date_time, y = mr)) - - if (!is.null(input$smr)) { - # if (is.null(input$sim_params)) { - mr_cols <- grepl("_mr", colnames(input$smr)) - smr <- reshape2::melt(input$smr, - id.vars = c("probe", "id"), - measure.vars = colnames(input$smr)[mr_cols]) - smr$Method <- sub("_mr", "", smr$variable) - # } else { - # mr_cols <- colnames(input$smr) %in% c("q0.2_mr", "q0.2_smr_real") - # smr <- reshape2::melt(input$smr, - # id.vars = c("probe", "id"), - # measure.vars = colnames(input$smr)[mr_cols]) - # smr$Method <- sub("_mr", " estimated", smr$variable) - # smr$Method <- sub("_smr_real", " real", smr$Method) - # } - smr$idprobe <- paste0(smr$id, " (", smr$probe, ")") - smr$idprobe <- factor(smr$idprobe, levels = idprobe_levels) - - if (!missing(probes)) { - smr <- smr[smr$probe %in% probes, ] - smr$idprobe <- droplevels(smr$idprobe) - } - - p <- p + ggplot2::geom_hline(data = smr, - ggplot2::aes(yintercept = value, - linetype = Method), - col = "red") - } - - if (!is.null(input$mmr)) { - mmr <- input$mmr - mmr$idprobe <- paste0(mmr$id, " (", mmr$probe, ")") - mmr$idprobe <- factor(mmr$idprobe, levels = idprobe_levels) - - if (!missing(probes)) { - mmr <- mmr[mmr$probe %in% probes, ] - mmr$idprobe <- droplevels(mmr$idprobe) - } - - p <- p + ggplot2::geom_point(data = mmr, - ggplot2::aes(x = date_time, - y = mr), - col = "red", size = 2) - } - - p <- p + ggplot2::theme_bw() - p <- p + ggplot2::facet_wrap(idprobe ~ ., ncol = 1) - p <- p + ggplot2::labs(y = "\u1E40[O[2]]", x = '', linetype = "SMR method") - return(p) + # ggplot variables + date_time <- mr_g <- mr_real <- value <- Method <- NULL + + mr <- input$mr + if ("mr_g" %in% colnames(mr)) { + mr$mr <- mr$mr_g + } else { + mr$mr <- mr$mr_abs + } + + mr$idprobe <- paste0(mr$id, " (", mr$probe, ")") + idprobe_levels <-paste0(input$probe_info$id, + " (", input$probe_info$probe, ")") + mr$idprobe <- factor(mr$idprobe, levels = idprobe_levels) + + if (!missing(cycles)) { + cycles <- check_arg_in_data(cycles, mr$cycle, + "cycles", verbose = verbose) + mr <- mr[mr$cycle %in% cycles, ] + } + + if (!missing(probes)) { + probes <- check_arg_in_data(probes, mr$probe, + "probes", verbose = verbose) + mr <- mr[mr$probe %in% probes, ] + mr$idprobe <- droplevels(mr$idprobe) + } + + p <- ggplot2::ggplot() + + # if (!is.null(input$sim_params)) { + # p <- plot_sim_watermark(p = p, + # x = mr$date_time, + # y = c(mr$mr_real, mr$mr)) + # mr_real_aux <- data.frame(mr_real = input$sim_params$real_mr, + # date_time = input$slopes$date_time) + # p <- p + ggplot2::geom_line(data = mr_real_aux, + # ggplot2::aes(x = date_time, y = mr_real), + # linetype = "dashed") + # p <- p + ggplot2::geom_point(data = mr_real_aux, + # ggplot2::aes(x = date_time, y = mr_real), + # shape = 1) + # } + + p <- p + ggplot2::geom_line(data = mr, + ggplot2::aes(x = date_time, y = mr)) + p <- p + ggplot2::geom_point(data = mr, + ggplot2::aes(x = date_time, y = mr)) + + if (!is.null(input$smr)) { + # if (is.null(input$sim_params)) { + mr_cols <- grepl("_mr", colnames(input$smr)) + smr <- reshape2::melt(input$smr, + id.vars = c("probe", "id"), + measure.vars = colnames(input$smr)[mr_cols]) + smr$Method <- sub("_mr", "", smr$variable) + # } else { + # mr_cols <- colnames(input$smr) %in% c("q0.2_mr", "q0.2_smr_real") + # smr <- reshape2::melt(input$smr, + # id.vars = c("probe", "id"), + # measure.vars = colnames(input$smr)[mr_cols]) + # smr$Method <- sub("_mr", " estimated", smr$variable) + # smr$Method <- sub("_smr_real", " real", smr$Method) + # } + smr$idprobe <- paste0(smr$id, " (", smr$probe, ")") + smr$idprobe <- factor(smr$idprobe, levels = idprobe_levels) + + if (!missing(probes)) { + smr <- smr[smr$probe %in% probes, ] + smr$idprobe <- droplevels(smr$idprobe) + } + + p <- p + ggplot2::geom_hline(data = smr, + ggplot2::aes(yintercept = value, + linetype = Method), + col = "red") + } + + if (!is.null(input$mmr)) { + mmr <- input$mmr + mmr$idprobe <- paste0(mmr$id, " (", mmr$probe, ")") + mmr$idprobe <- factor(mmr$idprobe, levels = idprobe_levels) + + if (!missing(probes)) { + mmr <- mmr[mmr$probe %in% probes, ] + mmr$idprobe <- droplevels(mmr$idprobe) + } + + p <- p + ggplot2::geom_point(data = mmr, + ggplot2::aes(x = date_time, + y = mmr), + col = "red", size = 2) + } + + p <- p + ggplot2::theme_bw() + p <- p + ggplot2::facet_wrap(idprobe ~ ., ncol = 1) + p <- p + ggplot2::labs(y = "\u1E40[O[2]]", x = '', linetype = "SMR method") + return(p) } #' Plot the standard metabolic rates #' #' @param input An experiment list with calculated SMR. -#' The output of \code{\link{process_experiment}} or -#' \code{\link{calc_smr}}. +#' The output of \code{\link{process_experiment}} or +#' \code{\link{calc_smr}}. #' @param probes A string of which probes to plot #' #' @return a ggplot object @@ -500,45 +505,45 @@ plot_mr <- function(input, cycles, probes, verbose = TRUE) { #' @export #' plot_smr <- function(input, probes) { - # ggplot variables - mr_g <- value <- Method <- NULL - - if (!missing(probes)) { - probes <- check_arg_in_data(probes, input$smr$probe, "probes") - input$smr <- input$smr[input$smr$probe %in% probes, ] - input$mr <- input$mr[input$mr$probe %in% probes, ] - } - - mr_cols <- grepl("mr_g", colnames(input$smr)) - aux <- reshape2::melt(input$smr, - id.vars = c("probe", "id"), - measure.vars = colnames(input$smr)[mr_cols]) - aux$Method <- sub("_mr_g", "", aux$variable) - - # change facet labels - input$mr$idprobe <- paste0(input$mr$id, " (", input$mr$probe, ")") - idprobe_levels <- paste0(input$probe_info$id, - " (", input$probe_info$probe, ")") - input$mr$idprobe <- factor(input$mr$idprobe, levels = idprobe_levels) - - aux$idprobe <- paste0(aux$id, " (", aux$probe, ")") - aux$idprobe <- factor(aux$idprobe, levels = idprobe_levels) - - p <- ggplot2::ggplot(data = input$mr) - p <- p + ggplot2::geom_histogram(ggplot2::aes(x = mr_g), - fill = "white", colour = "grey", binwidth = 0.1) - p <- p + ggplot2::geom_vline(data = aux, ggplot2::aes(xintercept = value, - linetype = Method)) - p <- p + ggplot2::labs(x = "\u1E40[O[2]]") - p <- p + ggplot2::facet_grid(idprobe ~ .) - p + # ggplot variables + mr_g <- value <- Method <- NULL + + if (!missing(probes)) { + probes <- check_arg_in_data(probes, input$smr$probe, "probes") + input$smr <- input$smr[input$smr$probe %in% probes, ] + input$mr <- input$mr[input$mr$probe %in% probes, ] + } + + mr_cols <- grepl("mr", colnames(input$smr)) + aux <- reshape2::melt(input$smr, + id.vars = c("probe", "id"), + measure.vars = colnames(input$smr)[mr_cols]) + aux$Method <- sub("_mr", "", aux$variable) + + # change facet labels + input$mr$idprobe <- paste0(input$mr$id, " (", input$mr$probe, ")") + idprobe_levels <- paste0(input$probe_info$id, + " (", input$probe_info$probe, ")") + input$mr$idprobe <- factor(input$mr$idprobe, levels = idprobe_levels) + + aux$idprobe <- paste0(aux$id, " (", aux$probe, ")") + aux$idprobe <- factor(aux$idprobe, levels = idprobe_levels) + + p <- ggplot2::ggplot(data = input$mr) + p <- p + ggplot2::geom_histogram(ggplot2::aes(x = mr_g), + fill = "white", colour = "grey", binwidth = 0.1) + p <- p + ggplot2::geom_vline(data = aux, ggplot2::aes(xintercept = value, + linetype = Method)) + p <- p + ggplot2::labs(x = "\u1E40[O[2]]") + p <- p + ggplot2::facet_grid(idprobe ~ .) + p } #' Plot the standard metabolic rates #' #' @param input An experiment list with calculated rolling MR values. -#' The output of \code{\link{roll_mr}}. +#' The output of \code{\link{roll_mr}}. #' @inheritParams plot_meas #' #' @return a ggplot object @@ -546,88 +551,88 @@ plot_smr <- function(input, probes) { #' @export #' plot_rolling_mr <- function(input, probes, cycles) { - # ggplot variables - phase_time <- date_time <- mr_g <- r2 <- NULL - - mr <- input$rolling_mr$values - max_mr <- input$rolling_mr$max - - if (is.null(mr)) { - stop("Could not find rolling mr values in input") - } - - if (!missing(probes)) { - probes <- check_arg_in_data(probes, mr$probe, - "probes", verbose = TRUE) - mr <- mr[mr$probe %in% probes, ] - max_mr <- max_mr[max_mr$probe %in% probes, ] - } - - if (!missing(cycles)) { - cycles <- check_arg_in_data(cycles, mr$cycle, - "cycles", verbose = TRUE) - mr <- mr[mr$cycle %in% cycles, ] - max_mr <- max_mr[max_mr$cycle %in% cycles, ] - } - - # remove units - o2_unit <- units(mr$mr_g) - mr$mr_g <- as.numeric(mr$mr_g) - max_mr$mr_g <- as.numeric(max_mr$mr_g) - - # make facets - aux <- unique(mr[, c("id", "probe", "cycle")]) - aux <- aux[order(aux$probe), ] - - ipp_levels <- paste0(aux$id, " (", - aux$probe, ") - ", - aux$cycle) - - mr$idprobecycle <- paste0(mr$id, " (", mr$probe, ") - ", mr$cycle) - mr$idprobecycle <- factor(mr$idprobecycle, levels = ipp_levels) - max_mr$idprobecycle <- paste0(max_mr$id, - " (", max_mr$probe, ") - ", max_mr$cycle) - max_mr$idprobecycle <- factor(max_mr$idprobecycle, levels = ipp_levels) - - - p <- ggplot2::ggplot(data = mr, ggplot2::aes(x = phase_time, y = mr_g)) - p <- p + ggplot2::geom_line() - p <- p + ggplot2::geom_point(data = max_mr, col = 'red') - p <- p + ggplot2::theme_bw() - - - rer_mean <- as.numeric(mean(mr$mr_g, na.rm = TRUE)) - r2_mean <- 0.5 - mean_dif <- r2_mean - rer_mean - aux <- as.numeric(range(mr$mr_g, na.rm = TRUE)) - rer_range <- aux[2] - aux[1] - r2_range <- 1 - - if (rer_range == 0 | r2_range == 0) { - compress_factor <- 1 - } else { - compress_factor <- rer_range/r2_range - } - - data_link <- function(x, a = r2_mean, b = rer_mean, - c_factor = compress_factor) { - b + ((x - a) * c_factor) # = y - } - - axis_link <- function(y, a = r2_mean, b = rer_mean, - c_factor = compress_factor) { - (y + a * c_factor - b) / c_factor # = x - } - - p <- p + ggplot2::geom_line(ggplot2::aes(y = data_link(r2)), col = "grey") - p <- p + ggplot2::scale_y_continuous( - sec.axis = ggplot2::sec_axis(trans = axis_link, name = "R2") - ) - - p <- p + ggplot2::labs(y = units::make_unit_label("\u1E40[O[2]]", o2_unit), - x = "Time within phase") - p <- p + ggplot2::facet_wrap(idprobecycle ~ ., ncol = 1) - p + # ggplot variables + phase_time <- date_time <- mr_g <- r2 <- NULL + + mr <- input$rolling_mr$values + max_mr <- input$rolling_mr$max + + if (is.null(mr)) { + stop("Could not find rolling mr values in input") + } + + if (!missing(probes)) { + probes <- check_arg_in_data(probes, mr$probe, + "probes", verbose = TRUE) + mr <- mr[mr$probe %in% probes, ] + max_mr <- max_mr[max_mr$probe %in% probes, ] + } + + if (!missing(cycles)) { + cycles <- check_arg_in_data(cycles, mr$cycle, + "cycles", verbose = TRUE) + mr <- mr[mr$cycle %in% cycles, ] + max_mr <- max_mr[max_mr$cycle %in% cycles, ] + } + + # remove units + o2_unit <- units(mr$mr_g) + mr$mr_g <- as.numeric(mr$mr_g) + max_mr$mr_g <- as.numeric(max_mr$mr_g) + + # make facets + aux <- unique(mr[, c("id", "probe", "cycle")]) + aux <- aux[order(aux$probe), ] + + ipp_levels <- paste0(aux$id, " (", + aux$probe, ") - ", + aux$cycle) + + mr$idprobecycle <- paste0(mr$id, " (", mr$probe, ") - ", mr$cycle) + mr$idprobecycle <- factor(mr$idprobecycle, levels = ipp_levels) + max_mr$idprobecycle <- paste0(max_mr$id, + " (", max_mr$probe, ") - ", max_mr$cycle) + max_mr$idprobecycle <- factor(max_mr$idprobecycle, levels = ipp_levels) + + + p <- ggplot2::ggplot(data = mr, ggplot2::aes(x = phase_time, y = mr_g)) + p <- p + ggplot2::geom_line() + p <- p + ggplot2::geom_point(data = max_mr, col = 'red') + p <- p + ggplot2::theme_bw() + + + rer_mean <- as.numeric(mean(mr$mr_g, na.rm = TRUE)) + r2_mean <- 0.5 + mean_dif <- r2_mean - rer_mean + aux <- as.numeric(range(mr$mr_g, na.rm = TRUE)) + rer_range <- aux[2] - aux[1] + r2_range <- 1 + + if (rer_range == 0 | r2_range == 0) { + compress_factor <- 1 + } else { + compress_factor <- rer_range/r2_range + } + + data_link <- function(x, a = r2_mean, b = rer_mean, + c_factor = compress_factor) { + b + ((x - a) * c_factor) # = y + } + + axis_link <- function(y, a = r2_mean, b = rer_mean, + c_factor = compress_factor) { + (y + a * c_factor - b) / c_factor # = x + } + + p <- p + ggplot2::geom_line(ggplot2::aes(y = data_link(r2)), col = "grey") + p <- p + ggplot2::scale_y_continuous( + sec.axis = ggplot2::sec_axis(trans = axis_link, name = "R2") + ) + + p <- p + ggplot2::labs(y = units::make_unit_label("\u1E40[O[2]]", o2_unit), + x = "Time within phase") + p <- p + ggplot2::facet_wrap(idprobecycle ~ ., ncol = 1) + p } @@ -641,94 +646,94 @@ plot_rolling_mr <- function(input, probes, cycles) { #' @export #' plot_experiment <- function(input, cycles, probe, verbose = FALSE) { - # if (!is.null(input$sim_params) & missing(probe)) { - # probe <- "S1" - # } - - if (length(probe) > 1) { - stop("Please select only one probe at a time.") - } - if (!is.null(input$bg$pre)) { - if (verbose) message('Plotting pre-background') - p_pre <- plot_bg(input$bg$pre, probes = probe) - p_pre <- p_pre + ggplot2::labs(title = 'Pre-background') - } else { - p_pre <- patchwork::wrap_elements( - grid::textGrob('Pre-bg was not included')) - } - - if (!is.null(input$bg$post)) { - if (verbose) message('Plotting post-background') - p_post <- plot_bg(input$bg$post, probes = probe) - p_post <- p_post + ggplot2::labs(title = 'Post-background') - } else { - p_post <- patchwork::wrap_elements( - grid::textGrob('Post-bg was not included')) - } - - if (verbose) message('Plotting measurements') - - # if (is.null(input$sim_params)) { - p_meas <- plot_meas(input, cycles = cycles, probes = probe, - show_temp = TRUE, verbose = verbose) - p_meas <- p_meas + ggplot2::labs(x = '') - x_ruler <- p_meas - # } - - if (verbose) message('Plotting plotting slopes') - - p_slopes <- plot_slopes(input, cycles = cycles, - probes = probe, verbose = FALSE) - p_slopes <- p_slopes + ggplot2::xlab('') - - # if (is.null(input$sim_params)) { - p_slopes <- p_slopes + mimic_x_datetime(x_ruler) - # } else { - # x_ruler <- p_slopes - # } - - if (verbose) message('Plotting metabolic rate') - - probe_check <- any(input$mr$probe == probe) - cycle_check <- (missing(cycles) || any(input$mr$cycle %in% cycles)) - combined_check <- (missing(cycles) || - any(input$mr$probe == probe & input$mr$cycle %in% cycles)) - if (probe_check && cycle_check && combined_check) { - p_mr <- plot_mr(input, cycles = cycles, probes = probe, verbose = FALSE) - p_mr <- p_mr + mimic_x_datetime(x_ruler) - } else { - p_mr <- patchwork::wrap_elements( - grid::textGrob('No valid MO2 values found')) - } - - # if (is.null(input$sim_params)) { - p_final <- p_pre + p_post + p_meas + p_slopes + - p_mr + patchwork::plot_layout(design = 'AB\nCC\nDD\nEE') - # } else { - # p_final <- p_pre + p_post + p_slopes + - # p_mr + patchwork::plot_layout(design = 'AB\nCC\nDD') - # } - return(p_final) + # if (!is.null(input$sim_params) & missing(probe)) { + # probe <- "S1" + # } + + if (length(probe) > 1) { + stop("Please select only one probe at a time.") + } + if (!is.null(input$bg$pre)) { + if (verbose) message('Plotting pre-background') + p_pre <- plot_bg(input$bg$pre, probes = probe) + p_pre <- p_pre + ggplot2::labs(title = 'Pre-background') + } else { + p_pre <- patchwork::wrap_elements( + grid::textGrob('Pre-bg was not included')) + } + + if (!is.null(input$bg$post)) { + if (verbose) message('Plotting post-background') + p_post <- plot_bg(input$bg$post, probes = probe) + p_post <- p_post + ggplot2::labs(title = 'Post-background') + } else { + p_post <- patchwork::wrap_elements( + grid::textGrob('Post-bg was not included')) + } + + if (verbose) message('Plotting measurements') + + # if (is.null(input$sim_params)) { + p_meas <- plot_meas(input, cycles = cycles, probes = probe, + show_temp = TRUE, verbose = verbose) + p_meas <- p_meas + ggplot2::labs(x = '') + x_ruler <- p_meas + # } + + if (verbose) message('Plotting plotting slopes') + + p_slopes <- plot_slopes(input, cycles = cycles, + probes = probe, verbose = FALSE) + p_slopes <- p_slopes + ggplot2::xlab('') + + # if (is.null(input$sim_params)) { + p_slopes <- p_slopes + mimic_x_datetime(x_ruler) + # } else { + # x_ruler <- p_slopes + # } + + if (verbose) message('Plotting metabolic rate') + + probe_check <- any(input$mr$probe == probe) + cycle_check <- (missing(cycles) || any(input$mr$cycle %in% cycles)) + combined_check <- (missing(cycles) || + any(input$mr$probe == probe & input$mr$cycle %in% cycles)) + if (probe_check && cycle_check && combined_check) { + p_mr <- plot_mr(input, cycles = cycles, probes = probe, verbose = FALSE) + p_mr <- p_mr + mimic_x_datetime(x_ruler) + } else { + p_mr <- patchwork::wrap_elements( + grid::textGrob('No valid MO2 values found')) + } + + # if (is.null(input$sim_params)) { + p_final <- p_pre + p_post + p_meas + p_slopes + + p_mr + patchwork::plot_layout(design = 'AB\nCC\nDD\nEE') + # } else { + # p_final <- p_pre + p_post + p_slopes + + # p_mr + patchwork::plot_layout(design = 'AB\nCC\nDD') + # } + return(p_final) } mimic_x_datetime <- function(input) { - ggplot2::xlim(as.POSIXct(ggplot2::layer_scales(input)$x$range$range, - origin = "1970-01-01 00:00.00")) + ggplot2::xlim(as.POSIXct(ggplot2::layer_scales(input)$x$range$range, + origin = "1970-01-01 00:00.00")) } mimic_y_datetime <- function(input) { - ggplot2::ylim(as.POSIXct(ggplot2::layer_scales(input)$y$range$range, - origin = "1970-01-01 00:00.00")) + ggplot2::ylim(as.POSIXct(ggplot2::layer_scales(input)$y$range$range, + origin = "1970-01-01 00:00.00")) } mimic_x <- function(input) { - ggplot2::xlim(ggplot2::layer_scales(input)$x$range$range) + ggplot2::xlim(ggplot2::layer_scales(input)$x$range$range) } mimic_y <- function(input) { - ggplot2::ylim(ggplot2::layer_scales(input)$y$range$range) + ggplot2::ylim(ggplot2::layer_scales(input)$y$range$range) } diff --git a/R/slope_functions.R b/R/slope_functions.R index 2c65d7d..d03934d 100644 --- a/R/slope_functions.R +++ b/R/slope_functions.R @@ -2,18 +2,18 @@ #' #' @param input A dataframe of clean measurements. #' The output of \code{\link{trim_resp}}. -#' @param correct_for_volume Logical: Should the corrected chamber volume be -#' used to remove the volume units from the slope? Defaults to FALSE, which +#' @param correct_for_water_vol Logical: Should water_vol be +#' used to remove the water_vol units from the slope? Defaults to FALSE, which #' is the right choice in most cases. If you are using a chamber to correct -#' the background of other chambers but they don't have the same volume +#' the background of other chambers but they don't have the same water_vol #' (not recommended), then you can set this to TRUE to minimise the impacts -#' of the different volumes on the results. Note that chambers with different -#' formats will have different surface-to-volume relationships, and might +#' of the different water_vols on the results. Note that chambers with different +#' formats will have different surface-to-water_vol relationships, and might #' therefore not be good controls of each other. #' #' @export #' -calc_slopes <- function(input, correct_for_volume = FALSE) { +calc_slopes <- function(input, correct_for_water_vol = FALSE) { # The operation is done by cycle and by probe, # so the dataset is broken twice below by_probe <- split(input$trimmed, input$trimmed$probe) @@ -25,18 +25,18 @@ calc_slopes <- function(input, correct_for_volume = FALSE) { recipient <- lapply(by_cycle, function(the_cycle) { output <- data.frame(probe = the_cycle$probe[1], id = the_cycle$id[1], - mass = NA_real_, - volume = the_cycle$volume[1], + animal_mass = NA_real_, + water_vol = the_cycle$water_vol[1], date_time = the_cycle$date_time[nrow(the_cycle)], phase = the_cycle$phase[1], cycle = the_cycle$cycle[1], temp = mean(the_cycle$temp, na.rm = TRUE)) - if ("mass" %in% colnames(the_cycle)) { - output$mass <- the_cycle$mass[1] + if ("animal_mass" %in% colnames(the_cycle)) { + output$animal_mass <- the_cycle$animal_mass[1] } else { # if mass is not provided, delete column - output$mass <- NULL + output$animal_mass <- NULL } if (all(is.na(the_cycle$o2_delta))) { @@ -67,13 +67,8 @@ calc_slopes <- function(input, correct_for_volume = FALSE) { output <- as.data.frame(data.table::rbindlist(recipient)) - if (correct_for_volume) { - if ("mass" %in% colnames(output)) { - water <- output$volume - conv_w_to_ml(output$mass) - } else { - water <- output$volume - } - output$slope <- output$slope * output$volume + if (correct_for_water_vol) { + output$slope <- output$slope * output$water_vol } if (any(is.na(output$slope))) { @@ -122,7 +117,7 @@ calc_single_slope <- function(input, probe, cycle, the_cycle <- my_probe[my_probe$cycle == cycle, ] - if ((max_duration - skip) < 10) { + if ((max_duration - skip) < 20) { warning("Calculating very short slopes is unadvisable.", immediate. = TRUE, call. = FALSE) } @@ -133,8 +128,8 @@ calc_single_slope <- function(input, probe, cycle, output <- data.frame(probe = the_cycle$probe[1], id = the_cycle$id[1], - mass = the_cycle$mass[1], - volume = the_cycle$volume[1], + animal_mass = the_cycle$animal_mass[1], + water_vol = the_cycle$water_vol[1], date_time = the_cycle$date_time[nrow(the_cycle)], phase = the_cycle$phase[1], cycle = the_cycle$cycle[1], @@ -145,9 +140,8 @@ calc_single_slope <- function(input, probe, cycle, output$se <- NA output$r2 <- NA } else { - # lm and units don't play along well if you intend - # to ask for a summary. Must drop units and reattach - # them later until this is fixed. + # lm and units don't play along well. + # Must drop units and reattach them later until this is fixed. m <- lm(as.numeric(o2_delta) ~ as.numeric(phase_time), data = the_cycle) diff --git a/R/smr_functions.R b/R/smr_functions.R index bfb95d4..bb59a1d 100644 --- a/R/smr_functions.R +++ b/R/smr_functions.R @@ -28,9 +28,9 @@ calc_smr <- function(mr, G = 1:4, q = c(0.2, 0.25), p = 0.1, n = 10){ issue_low_n_warning <- TRUE recipient <- lapply(by_probe, function(the_probe) { - - if ("mass" %in% colnames(the_probe)) { - output <- the_probe[1, c("probe", "id", "mass")] + # cat(the_probe$probe[1], "\n") + if ("animal_mass" %in% colnames(the_probe)) { + output <- the_probe[1, c("probe", "id", "animal_mass")] the_probe$input_mr <- the_probe$mr_g } else { output <- the_probe[1, c("probe", "id")] diff --git a/R/startup.R b/R/startup.R index 2fad152..765d945 100644 --- a/R/startup.R +++ b/R/startup.R @@ -1,7 +1,7 @@ .onAttach <- function(libname, pkgname){ inst.ver <- utils::packageVersion("pyroresp") aux <- as.matrix(data.frame(Package = "pyroresp", LibPath = NA, Version = as.character(utils::packageVersion("pyroresp")), Priority = NA, Built = NA)) - packageStartupMessage("Welcome to pyroresp (", inst.ver, ")!\nRun ?pyroresp for starting tips.") + packageStartupMessage("Welcome to pyroresp (", inst.ver, ")!\nThis package is rapidly evolving.") # new.ver <- tryCatch(old.packages(instPkgs = aux, repos = "https://cloud.r-project.org"), warning = function(w) NULL, error = function(e) NULL) # if (!is.null(new.ver)) { # nocov start # packageStartupMessage(paste0("-------------------------------------------------------------\n!!! A NEW VERSION of pyroresp is available! (v.", inst.ver, " -> v.", new.ver[, "ReposVer"], ")\n!!! You should update pyroresp before continuing.\n!!! You can update your packages by running update.packages()\n-------------------------------------------------------------\n")) diff --git a/R/wrappers.R b/R/wrappers.R index f0d5948..9c65b2f 100644 --- a/R/wrappers.R +++ b/R/wrappers.R @@ -9,20 +9,23 @@ #' @inheritParams load_phases #' @param probe_info a dataframe containing animal information. This #' dataframe must contain the following columns: -#' - id : The ID of the animal -#' - mass : The mass of the animal, in grams -#' - volume : The non-corrected volume of the chamber + tubing +#' - animal_id : The ID of the animal +#' - animal_mass : The mass of the animal, in grams +#' - animal_vol : Optional: The volume of the animal. +#' Assumed to be equal to mass if missing. +#' - chamber_vol : The non-corrected volume of the chamber #' - probe : The device-channel combination for the probe #' - first_cycle : The first cycle of valid data for that animal #' #' @return A list containing a phases dataframe and a pyro list with the #' individual source data frames (in source_data), as well as a single, -#' combined data frame organized by time (in compiled_data). +#' combined data frame organized by time (in compiled_data). #' #' @export #' load_experiment <- function(folder, date_format, tz = Sys.timezone(), - phases, probe_info, encoding = "ISO-8859-1") { + phases, probe_info, encoding = "ISO-8859-1", mass_unit = "g", + vol_unit = "ml") { if (length(folder) == 0 || !dir.exists(folder)) { stop('Could not find target folder') @@ -33,40 +36,25 @@ load_experiment <- function(folder, date_format, tz = Sys.timezone(), } if (!missing(probe_info)) { - if (!"mass" %in% colnames(probe_info)) { - warning("No column 'mass' found in the probe_info.", - " Won't correct chamber volume nor calculate", - " mass-corrected MO2.", - immediate. = TRUE, call. = FALSE) - } - required_cols <- c("id", "volume", "probe") - cols_missing <- !(required_cols %in% colnames(probe_info)) - if (any(cols_missing)) { - stop("The following required columns are missing ", - "from the probe_info input: ", - paste0(required_cols[cols_missing], collapse = ", "), - call. = FALSE) - } + probe_info <- process_probe_info(probe_info, + vol_unit = vol_unit, + mass_unit = mass_unit) + } else { + probe_info <- NULL } if (any(sapply(names(phases), nchar) > 4)) { - warning("Long device names detected in the phases input. Are you sure", - " you appended the device names correctly to the file name?", - " These are the current device names: ", + warning("Long device names detected in the phases input.", + " Confirm these are correct: ", paste(names(phases), collapse = ", "), ".") } - pyro <- load_pyro_data(folder, date_format = date_format, tz = tz, - encoding = encoding) + pyro <- load_pyro_folder(folder, date_format = date_format, + tz = tz, encoding = encoding) - output <- list(phases = phases, pyro = pyro) - - if (!missing(probe_info)) { - output$probe_info <- probe_info - - units(output$probe_info$mass) <- "g" - units(output$probe_info$volume) <- "ml" - } + output <- list(phases = phases, + pyro = pyro, + probe_info = probe_info) return(output) } @@ -81,7 +69,7 @@ load_experiment <- function(folder, date_format, tz = Sys.timezone(), #' #' @export #' -load_pyro_data <- function(folder, date_format, tz, +load_pyro_folder <- function(folder, date_format, tz, type = c("Oxygen", "pH", "Oxygen|pH"), encoding = "ISO-8859-1") { type <- match.arg(type) @@ -115,27 +103,6 @@ load_pyro_data <- function(folder, date_format, tz, #' @export #' compile_sources <- function(input) { - # start_aux <- sapply(input, function(i) { - # as.character(min(i$date_time)) - # }) - # start_aux <- as.POSIXct(start_aux, tz = attributes(input[[1]]$date_time)$tz[1]) - # very_start <- min(start_aux) - - # end_aux <- sapply(input, function(i) { - # as.character(max(i$date_time)) - # }) - # end_aux <- as.POSIXct(end_aux, tz = attributes(input[[1]]$date_time)$tz[1]) - # very_end <- max(end_aux) - - # recipient <- data.frame(date_time = seq(from = very_start, - # to = very_end, by = 1)) - - # for (i in input) { - # new_piece <- i[!duplicated(i$date_time), ] - # recipient <- merge(recipient, new_piece, - # by = 'date_time', all = TRUE) - # } - recipient <- input[[1]] head(recipient) recipient <- recipient[!duplicated(recipient$date_time), ] @@ -161,6 +128,8 @@ compile_sources <- function(input) { #' @param input The output of \code{\link{load_experiment}} #' @param original_o2 The o2 unit the data was captured in. #' @param convert_o2_to The o2 unit desired for the final results. +#' @param conv_with_mean_temp Calculate the mean temperature experienced +#' during a given phase and use that for O2 unit conversion. #' @param min_temp,max_temp #' For temperature ramp experiments. The minimum OR maximum temperatures #' that must be reached before data is considered valid. Discards all phases @@ -172,7 +141,7 @@ compile_sources <- function(input) { #' YYYY-MM-DD HH:MM:SS format. #' @param from_cycle,to_cycle #' Trim the experiment to a specific group of cycles. You may use one or both -#' of these arguments at the same time. Input must be numeric. +#' of these arguments at the same time. Input must be numeric. #' @param verbose Logical. Should steps being taken be detailed with messages. #' Defaults to TRUE. #' @@ -183,6 +152,7 @@ compile_sources <- function(input) { #' process_experiment <- function(input, wait = 0, tail_trim = 0, meas_max = Inf, meas_min = 60, original_o2, convert_o2_to, + conv_with_mean_temp = FALSE, na_action = c("ignore", "linear", "before", "after", "remove"), zero_buffer = 3, first_cycle = 1, min_temp, max_temp, start_time, stop_time, from_cycle, to_cycle, @@ -206,6 +176,30 @@ process_experiment <- function(input, wait = 0, tail_trim = 0, " at a time. See function help for details.") } + if (!missing(start_time)) { + if (verbose) { + message(paste0("M: Discarding readings before ", start_time, ".")) + } + keep <- input$pyro$compiled_data$date_time >= as.POSIXct(start_time) + if (all(!keep)) { + stop ("Data ends before ", start_time, ".") + } else { + input$pyro$compiled_data <- input$pyro$compiled_data[keep, ] + } + } + + if (!missing(stop_time)) { + if (verbose) { + message(paste0("M: Discarding readings after ", stop_time, ".")) + } + keep <- input$pyro$compiled_data$date_time <= as.POSIXct(stop_time) + if (all(!keep)) { + stop ("Data starts after ", stop_time, ".") + } else { + input$pyro$compiled_data <- input$pyro$compiled_data[keep, ] + } + } + if (verbose) { message("M: Merging pyroscience and phases file.") } @@ -237,6 +231,31 @@ process_experiment <- function(input, wait = 0, tail_trim = 0, meas_min = meas_min, first_cycle = first_cycle) + if (conv_with_mean_temp) { + if (verbose) message("M: Calculating mean temp per phase.") + aux <- aggregate(as.numeric(input$trimmed$temp), + list(probe = input$trimmed$probe, + phase = input$trimmed$phase), + mean, na.rm = TRUE) + input_probe_phase <- paste(input$trimmed$probe, input$trimmed$phase) + aux_probe_phase <- paste(aux$probe, aux$phase) + link <- match(input_probe_phase, aux_probe_phase) + input$trimmed$conv_temp <- aux$x[link] + + aux <- aggregate(as.numeric(input$melted$temp), + list(probe = input$melted$probe, + phase = input$melted$phase), + mean, na.rm = TRUE) + input_probe_phase <- paste(input$melted$probe, input$melted$phase) + aux_probe_phase <- paste(aux$probe, aux$phase) + link <- match(input_probe_phase, aux_probe_phase) + input$melted$conv_temp <- aux$x[link] + + } else { + input$trimmed$conv_temp <- as.numeric(input$trimmed$temp) + input$melted$conv_temp <- as.numeric(input$melted$temp) + } + if (verbose) message("M: Calculating air saturation.") o2_conv_cols <- c("o2", "temp", "sal", "pressure") @@ -253,7 +272,7 @@ process_experiment <- function(input, wait = 0, tail_trim = 0, o2 = as.numeric(input$trimmed$o2[not_NA]), from = original_o2, to = "percent_a.s.", - temp = as.numeric(input$trimmed$temp[not_NA]), + temp = as.numeric(input$trimmed$conv_temp[not_NA]), sal = as.numeric(input$trimmed$sal[not_NA]), atm_pres = as.numeric(input$trimmed$pressure[not_NA]) ) @@ -281,7 +300,7 @@ process_experiment <- function(input, wait = 0, tail_trim = 0, o2 = input$trimmed$o2[not_NA], from = original_o2, to = convert_o2_to, - temp = as.numeric(input$trimmed$temp[not_NA]), + temp = as.numeric(input$trimmed$conv_temp[not_NA]), sal = as.numeric(input$trimmed$sal[not_NA]), atm_pres = as.numeric(input$trimmed$pressure[not_NA]) ) @@ -296,13 +315,17 @@ process_experiment <- function(input, wait = 0, tail_trim = 0, o2 = input$melted$o2[not_NA], from = original_o2, to = convert_o2_to, - temp = as.numeric(input$melted$temp[not_NA]), + temp = as.numeric(input$melted$conv_temp[not_NA]), sal = as.numeric(input$melted$sal[not_NA]), atm_pres = as.numeric(input$melted$pressure[not_NA]) ) units(input$melted$o2) <- gsub("_per_", "/", convert_o2_to) } + # remove temporary temperature columns + input$trimmed$conv_temp <- NULL + input$melted$conv_temp <- NULL + if (verbose) message("M: Calculating deltas.") input$trimmed <- calc_delta(input$trimmed, zero_buffer = zero_buffer) @@ -336,40 +359,6 @@ process_experiment <- function(input, wait = 0, tail_trim = 0, input$trimmed <- input$trimmed[input$trimmed$date_time >= time_break, ] } - if (!missing(start_time)) { - if (verbose) { - message(paste0("M: Discarding phases before ", start_time, ".")) - } - the_matches <- which(input$trimmed$date_time >= as.POSIXct(start_time)) - cutoff <- head(the_matches, 1) - if (length(cutoff) == 0) { - stop ("Data ends before ", start_time, ".") - } else { - first_phase <- input$trimmed$phase[cutoff] - the_matches <- which(input$trimmed$phase == first_phase) - first_true <- head(the_matches, 1) - break_ <- input$trimmed$date_time[first_true] - input$trimmed <- input$trimmed[input$trimmed$date_time >= break_, ] - } - } - - if (!missing(stop_time)) { - if (verbose) { - message(paste0("M: Discarding phases after ", stop_time, ".")) - } - the_matches <- which(input$trimmed$date_time <= as.POSIXct(stop_time)) - cutoff <- tail(the_matches, 1) - if (length(cutoff) == 0) { - stop ("Data starts after ", stop_time, ".") - } else { - last_phase <- input$trimmed$phase[cutoff] - the_matches <- which(input$trimmed$phase == last_phase) - last_true <- tail(the_matches, 1) - break_ <- input$trimmed$date_time[last_true] - input$trimmed <- input$trimmed[input$trimmed$date_time <= break_, ] - } - } - if (!missing(from_cycle)) { if (verbose) { message(paste0("M: Discarding cycles prior to cycle ", @@ -401,8 +390,8 @@ process_experiment <- function(input, wait = 0, tail_trim = 0, #' @export #' process_slopes <- function(input, r2 = 0.95, pre, post, method, - correct_for_volume = FALSE) { - input <- calc_slopes(input, correct_for_volume = correct_for_volume) + correct_for_water_vol = FALSE) { + input <- calc_slopes(input, correct_for_water_vol = correct_for_water_vol) input <- subtract_bg(input = input, pre = pre, post = post, method = method) @@ -436,6 +425,10 @@ process_slopes <- function(input, r2 = 0.95, pre, post, method, process_mr <- function(input, G = 1:4, q = c(0.2, 0.25), p = 0.1, n = 10) { + if (is.null(input$good_slopes)) { + stop("input has no good slope data.", + call. = FALSE) + } input$mr <- calc_mr(input$good_slopes) # convert seconds to hours (more common) @@ -454,7 +447,7 @@ process_mr <- function(input, G = 1:4, input$smr <- NULL } else { input$smr <- calc_smr(input$mr, G = G, q = q, p = p, n = n) - keep_these <- !(colnames(input$smr) %in% c("id", "mass", "volume")) + keep_these <- !(colnames(input$smr) %in% c("id", "animal_mass", "water_vol")) smr_aux <- input$smr[, keep_these] input$smr <- merge(input$probe_info, smr_aux, by = "probe", all = TRUE) @@ -464,8 +457,11 @@ process_mr <- function(input, G = 1:4, prefix <- sub("_mr", "", i) if (!is.null(input$bg$pre)) { new_col <- paste0(prefix, "pre_bg_pct") - aux <- input$smr[, i] - aux <- -aux * input$smr$mass / input$smr$volume + aux <- -input$smr[, i] + if ("animal_mass" %in% colnames(input$smr)) { + aux <- aux * input$smr$animal_mass + } + aux <- aux / input$smr$water_vol units(aux) <- units(input$bg$pre$bg$slope) bg_link <- match(input$smr$probe, input$bg$pre$bg$probe) input$smr[, new_col] <- input$bg$pre$bg$slope[bg_link] / aux @@ -474,8 +470,11 @@ process_mr <- function(input, G = 1:4, } if (!is.null(input$bg$post)) { new_col <- paste0(prefix, "post_bg_pct") - aux <- input$smr[, i] - aux <- -aux * input$smr$mass / input$smr$volume + aux <- -input$smr[, i] + if ("animal_mass" %in% colnames(input$smr)) { + aux <- aux * input$smr$animal_mass + } + aux <- aux / input$smr$water_vol units(aux) <- units(input$bg$post$bg$slope) bg_link <- match(input$smr$probe, input$bg$post$bg$probe) input$smr[, new_col] <- input$bg$post$bg$slope[bg_link] / aux @@ -489,26 +488,27 @@ process_mr <- function(input, G = 1:4, input$mmr <- extract_mmr(input$mr) mmr_aux <- input$mmr[, !(colnames(input$mmr) %in% c("id", "mass", "volume"))] input$mmr <- merge(input$probe_info, mmr_aux, - by = "probe", all = TRUE) - - if (!is.null(input$smr) && !is.null(input$bg$pre)) { - aux <- -input$mmr$mr_g * input$mmr$mass / input$smr$volume - units(aux) <- units(input$bg$pre$bg$slope) - bg_link <- match(input$mmr$probe, input$bg$pre$bg$probe) - input$mmr$pre_bg_pct <- input$bg$pre$bg$slope[bg_link] / aux - # units is now "1"; changing to percent automatically multiplies by 100 - units(input$mmr$pre_bg_pct) <- "percent" - } - if (!is.null(input$smr) && !is.null(input$bg$post)) { - aux <- -input$mmr$mr_g * input$mmr$mass / input$smr$volume - units(aux) <- units(input$bg$post$bg$slope) - bg_link <- match(input$mmr$probe, input$bg$post$bg$probe) - input$mmr$post_bg_pct <- input$bg$post$bg$slope[bg_link] / aux - # units is now "1"; changing to percent automatically multiplies by 100 - units(input$mmr$post_bg_pct) <- "percent" - } - # how to make a bg summary when using a reference chamber? + by = "probe", all = TRUE) + + ## calculating bg as a percentage of mmr seems irrelevant. + # if (!is.null(input$smr) && !is.null(input$bg$pre)) { + # aux <- -input$mmr$mr_g * input$mmr$mass / input$smr$volume + # units(aux) <- units(input$bg$pre$bg$slope) + # bg_link <- match(input$mmr$probe, input$bg$pre$bg$probe) + # input$mmr$pre_bg_pct <- input$bg$pre$bg$slope[bg_link] / aux + # # units is now "1"; changing to percent automatically multiplies by 100 + # units(input$mmr$pre_bg_pct) <- "percent" + # } + # if (!is.null(input$smr) && !is.null(input$bg$post)) { + # aux <- -input$mmr$mr_g * input$mmr$mass / input$smr$volume + # units(aux) <- units(input$bg$post$bg$slope) + # bg_link <- match(input$mmr$probe, input$bg$post$bg$probe) + # input$mmr$post_bg_pct <- input$bg$post$bg$slope[bg_link] / aux + # # units is now "1"; changing to percent automatically multiplies by 100 + # units(input$mmr$post_bg_pct) <- "percent" + # } + return(input) } \ No newline at end of file diff --git a/man/calc_mr.Rd b/man/calc_mr.Rd index 6713a26..80afea5 100644 --- a/man/calc_mr.Rd +++ b/man/calc_mr.Rd @@ -4,20 +4,15 @@ \alias{calc_mr} \title{Calculate metabolic rates} \usage{ -calc_mr(slope_data, density = 1) +calc_mr(input) } \arguments{ -\item{slope_data}{a data frame obtained by using either the -function \code{\link{calc_slopes}} or \code{\link{filter_r2}}} - -\item{density}{The density. defaults to 1 (1 g = 1 ml). -Higher densities lead to lower volumes, e.g. if density = 2, -then 1 g = 0.5 ml.} +\item{input}{a data frame with columns slope_cor (mandatory), +water_vol (mandatory), and animal_mass (optional).} } \value{ -A data frame with mr_abs, bg, and mr_cor. +The input with columns mr_abs and mr_g if animal_mass provided. } \description{ -Calculates the absolute and mass-specific metabolic rates, -as well as the relative importance of background (in percentage). +Calculates absolute and mass-specific metabolic rates. } diff --git a/man/calc_slopes.Rd b/man/calc_slopes.Rd index 0db46dc..cd8f0cd 100644 --- a/man/calc_slopes.Rd +++ b/man/calc_slopes.Rd @@ -4,19 +4,19 @@ \alias{calc_slopes} \title{Calculate slopes from clean measurements} \usage{ -calc_slopes(input, correct_for_volume = FALSE) +calc_slopes(input, correct_for_water_vol = FALSE) } \arguments{ \item{input}{A dataframe of clean measurements. The output of \code{\link{trim_resp}}.} -\item{correct_for_volume}{Logical: Should the corrected chamber volume be -used to remove the volume units from the slope? Defaults to FALSE, which +\item{correct_for_water_vol}{Logical: Should water_vol be +used to remove the water_vol units from the slope? Defaults to FALSE, which is the right choice in most cases. If you are using a chamber to correct -the background of other chambers but they don't have the same volume +the background of other chambers but they don't have the same water_vol (not recommended), then you can set this to TRUE to minimise the impacts -of the different volumes on the results. Note that chambers with different -formats will have different surface-to-volume relationships, and might +of the different water_vols on the results. Note that chambers with different +formats will have different surface-to-water_vol relationships, and might therefore not be good controls of each other.} } \description{ diff --git a/man/extract_epoc.Rd b/man/extract_epoc.Rd index 7fa4dd0..5e63538 100644 --- a/man/extract_epoc.Rd +++ b/man/extract_epoc.Rd @@ -4,12 +4,15 @@ \alias{extract_epoc} \title{Extracts the cumulative AUC} \usage{ -extract_epoc(input, max_duration = Inf) +extract_epoc(input, min_duration = 0, max_duration = Inf) } \arguments{ \item{input}{a resp object with an "auc" sub-object. The output of \code{\link{calc_auc}}} +\item{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.} + \item{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.} diff --git a/man/load_experiment.Rd b/man/load_experiment.Rd index 3898fe3..2ff11f2 100644 --- a/man/load_experiment.Rd +++ b/man/load_experiment.Rd @@ -10,7 +10,9 @@ load_experiment( tz = Sys.timezone(), phases, probe_info, - encoding = "ISO-8859-1" + encoding = "ISO-8859-1", + mass_unit = "g", + vol_unit = "ml" ) } \arguments{ @@ -24,9 +26,11 @@ load_experiment( \item{probe_info}{a dataframe containing animal information. This dataframe must contain the following columns: -- id : The ID of the animal -- mass : The mass of the animal, in grams -- volume : The non-corrected volume of the chamber + tubing +- animal_id : The ID of the animal +- animal_mass : The mass of the animal, in grams +- animal_vol : Optional: The volume of the animal. +Assumed to be equal to mass if missing. +- chamber_vol : The non-corrected volume of the chamber - probe : The device-channel combination for the probe - first_cycle : The first cycle of valid data for that animal} diff --git a/man/load_pyro_data.Rd b/man/load_pyro_folder.Rd similarity index 92% rename from man/load_pyro_data.Rd rename to man/load_pyro_folder.Rd index f7eee58..e4f3dcf 100644 --- a/man/load_pyro_data.Rd +++ b/man/load_pyro_folder.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/wrappers.R -\name{load_pyro_data} -\alias{load_pyro_data} +\name{load_pyro_folder} +\alias{load_pyro_folder} \title{Wrapper to scan a pyro folder and load raw data files.} \usage{ -load_pyro_data( +load_pyro_folder( folder, date_format, tz, diff --git a/man/process_probe_info.Rd b/man/process_probe_info.Rd new file mode 100644 index 0000000..aa2fc7e --- /dev/null +++ b/man/process_probe_info.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/aux_functions.R +\name{process_probe_info} +\alias{process_probe_info} +\title{needs a better place in the future. +needs to be exported in the future.} +\usage{ +process_probe_info(input, vol_unit = "ml", mass_unit = "g") +} +\description{ +needs a better place in the future. +needs to be exported in the future. +} +\keyword{internal} diff --git a/man/process_slopes.Rd b/man/process_slopes.Rd index ec366bf..8c47492 100644 --- a/man/process_slopes.Rd +++ b/man/process_slopes.Rd @@ -37,15 +37,6 @@ by cycle. \item "none" - does not perform oxygen consumption subtraction. Not recommended for anything other than checking test data. }} - -\item{correct_for_volume}{Logical: Should the corrected chamber volume be -used to remove the volume units from the slope? Defaults to FALSE, which -is the right choice in most cases. If you are using a chamber to correct -the background of other chambers but they don't have the same volume -(not recommended), then you can set this to TRUE to minimise the impacts -of the different volumes on the results. Note that chambers with different -formats will have different surface-to-volume relationships, and might -therefore not be good controls of each other.} } \value{ An updated experiment list, containing two new objects; one with