diff --git a/.gitignore b/.gitignore
index 09eeec8..93d31df 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,6 @@
 # OS files
 .DS_Store
 inst/doc
+
+# qaqc outputs
+*.html
diff --git a/NAMESPACE b/NAMESPACE
index b789f43..87a6f79 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/models.R b/R/models.R
index e8ec099..d7fa9ef 100644
--- a/R/models.R
+++ b/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))
diff --git a/R/qaqc.R b/R/qaqc.R
new file mode 100644
index 0000000..3395c7c
--- /dev/null
+++ b/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)
+}
diff --git a/R/zzz.R b/R/zzz.R
index 3fc5fb1..ffbda40 100644
--- a/R/zzz.R
+++ b/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(...)
diff --git a/inst/qaqc.Rmd b/inst/qaqc.Rmd
new file mode 100644
index 0000000..ac98ff2
--- /dev/null
+++ b/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)
+```
diff --git a/man/wtf_isquiet.Rd b/man/wtf_isquiet.Rd
new file mode 100644
index 0000000..ea7ed3b
--- /dev/null
+++ b/man/wtf_isquiet.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/zzz.R
+\name{wtf_isquiet}
+\alias{wtf_isquiet}
+\title{Check if whattheflux.quiet option is TRUE}
+\usage{
+wtf_isquiet()
+}
+\value{
+TRUE or FALSE
+}
+\description{
+Check if whattheflux.quiet option is TRUE
+}
+\keyword{internal}
diff --git a/man/wtf_normalize_time.Rd b/man/wtf_normalize_time.Rd
index 95fe28d..dc3bb26 100644
--- a/man/wtf_normalize_time.Rd
+++ b/man/wtf_normalize_time.Rd
@@ -18,7 +18,4 @@ TRUE) or the original vector if not.
 \description{
 Normalize a vector of times
 }
-\examples{
-wtf_normalize_time(2:4) # returns 0:2
-wtf_normalize_time(2:4, FALSE) # returns 2:4
-}
+\keyword{internal}
diff --git a/man/wtf_qaqc.Rd b/man/wtf_qaqc.Rd
new file mode 100644
index 0000000..30eaf83
--- /dev/null
+++ b/man/wtf_qaqc.Rd
@@ -0,0 +1,40 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/qaqc.R
+\name{wtf_qaqc}
+\alias{wtf_qaqc}
+\title{Generate a QA/QC document}
+\usage{
+wtf_qaqc(
+  flux_data,
+  group_column,
+  output_file = "qaqc.html",
+  output_dir = getwd(),
+  open_output = TRUE
+)
+}
+\arguments{
+\item{flux_data}{A data frame from \code{\link{wtf_compute_fluxes}} or similar}
+
+\item{group_column}{Name of the grouping label column in \code{flux_data},
+character; pass NULL to run with no grouping}
+
+\item{output_file}{Name of the output file}
+
+\item{output_dir}{Name of the output directory; default is current working directory}
+
+\item{open_output}{Automatically open the output HTML file?}
+}
+\value{
+The path of the output file
+}
+\description{
+Generate a QA/QC document
+}
+\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
+}