## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>",
  fig.width = 8, fig.height = 5, dpi = 96,
  message = FALSE, warning = FALSE
)
## the plotting functions need the suggested ggplot2 / patchwork; skip those
## chunks gracefully if they are not installed.
have_plots <- requireNamespace("ggplot2", quietly = TRUE) &&
              requireNamespace("patchwork", quietly = TRUE)

## ----packages-----------------------------------------------------------------
library(appac)

## ----load---------------------------------------------------------------------
acn <- list(sample_col = "sample.name", peak_col = "peak.name",
            date_col = "injection.date", pressure_col = "air.pressure",
            area_col = "raw.area")
data  <- check_cols(PLOT_FID, acn)
ap    <- as.numeric(data[, "Air_Pressure"])
P_ref <- mean(range(ap, na.rm = TRUE))
cat(sprintf("%d injections, %d cylinders, %d peaks;  P_ref = %.1f hPa\n",
            nrow(data) / length(unique(data$Peak_Name)),
            length(unique(data$Sample_Name)),
            length(unique(data$Peak_Name)), P_ref))

## ----fit----------------------------------------------------------------------
fit1 <- appac(data = data, P_ref = P_ref)
ct   <- debias_ct(fit1, data = data, P_ref = P_ref, npt = 7, quiet = TRUE)
fit  <- appac(data = data, ct = ct, P_ref = P_ref)
cat(sprintf("kappa = %.3e  (per hPa)\n", unlist(fit@correction@coefficients)))

## ----kappa-fit, eval = have_plots---------------------------------------------
plot_area_pressure_fit(fit)

## ----area-date, eval = have_plots---------------------------------------------
plot_area_date(fit, sample = 1, peak = "n.C4H10", show_changepoints = TRUE)

## ----rsd----------------------------------------------------------------------
rsd <- function(x) sd(x) / mean(x) * 100
raw <- fit@samples[[1]]$raw.area
cor <- fit@samples[[1]]$corrected.area
data.frame(peak              = colnames(raw),
           RSD_raw_pct       = round(apply(raw, 2, rsd), 3),
           RSD_corrected_pct = round(apply(cor, 2, rsd), 3))

## ----gof----------------------------------------------------------------------
goodness_of_fit(fit)[[1]]

## ----residuals, fig.width = 10, fig.height = 8, eval = have_plots-------------
plot_residuals(fit, sample = 1, peak = "n.C4H10")

## ----show---------------------------------------------------------------------
fit

## ----changepoints-------------------------------------------------------------
get_changepoints(fit@samples)            # episode level breakpoints