Initial pass at whattheflux QA/QC (#35)

* Initial structure for a QAQC Rmarkdown doc

* Open output HTML file

* Rearrange things

* Pass data to Rmd; first, initial, QA/QC plots

* Update qaqc.Rmd

.gitignore

@@ -6,3 +6,6 @@
 # OS files
 .DS_Store
 inst/doc
+
+# qaqc outputs
+*.html

NAMESPACE

@@ -3,8 +3,8 @@
 export(wtf_compute_fluxes)
 export(wtf_fit_models)
 export(wtf_metadata_match)
-export(wtf_normalize_time)
 export(wtf_ppm_to_umol)
+export(wtf_qaqc)
 export(wtf_read_LGR915)
 export(wtf_read_LI7810)
 export(wtf_read_LI7820)
@@ -21,6 +21,7 @@ importFrom(lubridate,ymd_hms)
 importFrom(stats,coefficients)
 importFrom(stats,lm)
 importFrom(stats,na.omit)
+importFrom(utils,browseURL)
 importFrom(utils,head)
 importFrom(utils,read.csv)
 importFrom(utils,read.table)

R/models.R

@@ -77,10 +77,7 @@ wtf_fit_models <- function(time, conc, area, volume) {
 #' @param normalize Normalize the values so that first is zero? Logical
 #' @return A numeric vector of normalized values (if \code{normalize_time} is
 #' TRUE) or the original vector if not.
-#' @export
-#' @examples
-#' wtf_normalize_time(2:4) # returns 0:2
-#' wtf_normalize_time(2:4, FALSE) # returns 2:4
+#' @keywords internal
 wtf_normalize_time <- function(time, normalize = TRUE) {
   if(normalize) {
     as.numeric(time) - as.numeric(min(time, na.rm = TRUE))

R/qaqc.R

@@ -0,0 +1,48 @@
+# qaqc.R
+
+
+#' Generate a QA/QC document
+#'
+#' @param flux_data A data frame from \code{\link{wtf_compute_fluxes}} or similar
+#' @param group_column Name of the grouping label column in \code{flux_data},
+#' character; pass NULL to run with no grouping
+#' @param output_file Name of the output file
+#' @param output_dir Name of the output directory; default is current working directory
+#' @param open_output Automatically open the output HTML file?
+#' @importFrom utils browseURL
+#' @return The path of the output file
+#' @export
+#'
+#' @examples
+#' # Toy data
+#' cars$Plot <- c("A", "B")
+#' fd <- wtf_compute_fluxes(cars, "Plot", "speed", "dist")
+#' x <- wtf_qaqc(fd, group_column = "Plot")
+#' file.remove(x) # clean up
+#' # See the introductory vignette for a fully-worked example with real data
+wtf_qaqc <- function(flux_data,
+                     group_column,
+                     output_file = "qaqc.html",
+                     output_dir = getwd(),
+                     open_output = TRUE) {
+  if(!requireNamespace("rmarkdown", quietly = TRUE)) {
+    stop("To run this function, please install the rmarkdown package")
+  }
+
+  # Save the flux data into a temporary file so as to pass the
+  # fully-qualified filename as a parameter to our Rmarkdown file
+  f <- system.file("qaqc.Rmd", package = "whattheflux")
+  td <- tempdir()
+  tf_flux_data <- file.path(td, "flux_data")
+  saveRDS(flux_data, tf_flux_data)
+
+  # Render
+  fout <- rmarkdown::render(f,
+                            output_file = output_file,
+                            output_dir = output_dir,
+                            quiet = wtf_isquiet(),
+                            params = list(flux_data = tf_flux_data,
+                                          group_column = group_column))
+  if(open_output) browseURL(paste0('file://', fout))
+  invisible(fout)
+}

R/zzz.R

@@ -1,7 +1,14 @@
 # zzz.R - miscellany
 
+#' Check if whattheflux.quiet option is TRUE
+#' @return TRUE or FALSE
+#' @keywords internal
+wtf_isquiet <- function() {
+  getOption("whattheflux.quiet", default = FALSE)
+}
+
 wtf_message <- function(...) {
-  if (getOption("whattheflux.quiet", default = FALSE)) {
+  if (wtf_isquiet()) {
     return()
   }
   message(...)

inst/qaqc.Rmd

@@ -0,0 +1,77 @@
+---
+title: "whattheflux QA/QC"
+date: "`r Sys.Date()`"
+output:
+  html_document:
+    code_folding: hide
+params:
+  flux_data: ""
+  group_column: NULL
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+
+library(ggplot2)
+theme_set(theme_bw())
+
+# Utility function: add labels to the lowest X%
+add_labels <- function(x, col1, col2, level, gc = params$group_column) {
+  reldiff <- abs((x[[col1]] - x[[col2]]) / x[[col1]])
+  outthere <- reldiff > quantile(reldiff, probs = level)
+  x$label <- NA
+  if(is.null(gc)) {
+    x$label[outthere] <- seq_len(nrow(x))[outthere]
+  } else {
+    x$label[outthere] <- x[[gc]][outthere]
+  }
+  return(x)
+}
+
+# Read the flux data
+fd <- readRDS(params$flux_data)
+```
+
+Rows of flux data: `r nrow(fd)`
+
+## Flux distribution
+
+```{r, flux-distribution}
+p <- ggplot(fd, aes(slope_estimate)) +
+  geom_histogram(bins = 30)
+print(p)
+```
+
+## Linear model versus robust linear model
+
+Fluxes that depart from the 1:1 may have influential outliers in the underlying data.
+
+```{r, linear-versus-robust}
+fd <- add_labels(fd, "slope_estimate", "slope_estimate_robust", level = 0.75)
+p <- ggplot(fd, aes(slope_estimate, slope_estimate_robust, label = label)) +
+  geom_point() +
+  geom_text(size = 5, nudge_y = 0.1, na.rm = TRUE) +
+  geom_abline() + 
+  theme(legend.position = "none")
+
+#if(!is.null(group_column)) p <- p + aes(color = .data[[group_column]])
+
+print(p)
+```
+
+## Linear model versus polynomial model
+
+Fluxes that depart from the 1:1 may have nonlinearity issues in the underlying data.
+
+```{r, linear-versus-polynomial}
+fd <- add_labels(fd, "adj.r.squared", "r.squared_poly", level = 0.75)
+p <- ggplot(fd, aes(adj.r.squared, r.squared_poly, label = label)) +
+  geom_point() +
+  geom_text(size = 5, nudge_y = 0.01, na.rm = TRUE) +
+  geom_abline() + 
+  theme(legend.position = "none")
+
+#if(!is.null(group_column)) p <- p + aes(color = .data[[group_column]])
+
+print(p)
+```