--- title: "APPAC: correcting GC peak areas for atmospheric pressure" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{APPAC: correcting GC peak areas for atmospheric pressure} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} %\VignetteDepends{ggplot2, patchwork} --- ```{r 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) ``` **APPAC** (*Atmospheric Pressure Peak Area Correction*) corrects gas-chromatography peak areas for the influence of ambient **air pressure** on standard detectors open to the atmosphere — chiefly the flame-ionization detector (FID). This vignette first shows the everyday **usage** on the bundled `PLOT_FID` data, then explains the **method**: how a single common pressure sensitivity, per-cylinder drift and measurement noise are separated by a principal-component decomposition. # Usage APPAC ships the `PLOT_FID` example dataset — six control cylinders, five hydrocarbon peaks and roughly five years of real FID injections with expert-annotated peak integration. The plots use the suggested `ggplot2` / `patchwork`. ```{r packages} library(appac) ``` ## Load the data `PLOT_FID` is long-format (one row per injection × peak). `check_cols()` maps the raw column names to the canonical ones and cleans the sample/peak names. The reference pressure `P_ref` is the mid-range of the observed air pressure. ```{r 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 the correction Two passes: `appac()` estimates the common pressure sensitivity **kappa** and the drift model; `debias_ct()` then finds the centre that minimises the chi-square of the residuals (the bias–variance trade-off), and a second `appac()` runs with that de-biased centre. *(`npt = 7` keeps the vignette quick; production uses `npt = 21`.)* ```{r 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))) ``` ## Inspect the fit `plot_area_pressure_fit()` shows exactly what enters the kappa estimate: the reference-scaled correlated area for every (cylinder, peak), binned by pressure deviation. All series share the single common slope kappa. ```{r kappa-fit, eval = have_plots} plot_area_pressure_fit(fit) ``` Raw areas scatter with the day-to-day pressure swings; the corrected areas collapse onto the reference. Detected episode breakpoints are marked. ```{r area-date, eval = have_plots} plot_area_date(fit, sample = 1, peak = "n.C4H10", show_changepoints = TRUE) ``` ## Run-to-run scatter and diagnostics The headline metric — relative standard deviation of the areas before and after correction (cylinder 1): ```{r 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)) ``` Reduced chi-square ≈ 1 means the corrected areas are down to the short-term measurement-noise floor: ```{r gof} goodness_of_fit(fit)[[1]] ``` `plot_residuals()` panels the residual diagnostics — a (leptokurtic, heavy-tailed) histogram with a fitted normal, the normal Q–Q, and residual vs. date / vs. pressure. Both scatter panels are flat, confirming the structured pressure and drift artifacts are gone and only the measurement noise remains. ```{r residuals, fig.width = 10, fig.height = 8, eval = have_plots} plot_residuals(fit, sample = 1, peak = "n.C4H10") ``` A compact summary of the fitted object (`print()` adds the per-sample goodness of fit): ```{r show} fit ``` # Method: the decomposition ## The physical effect and the forward model Boček, Novák & Janák (1969) showed experimentally that the FID's ionization efficiency is strongly pressure-dependent: normal atmospheric swings already change the sensitivity by up to ±5 % relative to the value at a fixed pressure. In routine data this pressure signal is **masked by slow instrument drift**, so it must be *separated* from the drift before it can be corrected. APPAC models the raw area of peak \(j\) in cylinder \(i\) at injection \(t\) as a multiplicative combination of a reference area, a pressure factor, a slow drift / daily factor and noise: \[ A_{ijt} \;=\; \mu_{ij}\,\bigl(1 + \kappa\,(P_t - P_\mathrm{ref})\bigr)\, (1 + \delta_{it})\,(1 + \varepsilon_{ijt}), \] where \(\mu_{ij}\) is the per-cylinder, per-peak reference (centre), \(\kappa\) the **common** pressure sensitivity shared by all peaks and cylinders, \(P_t\) the ambient pressure, \(\delta_{it}\) a slow drift plus a daily factor, and \(\varepsilon_{ijt}\) heavy-tailed measurement noise. ## Separating the components by PCA Per cylinder, each peak is standardised and a one-component principal-component analysis is run on the (dates × peaks) matrix. The components carry distinct physics: * **PC1 — correlated component:** the dominant common mode across peaks — the pressure response together with the smooth drift. * **PC2 — uncorrelated component:** the second mode, carrying the per-cylinder bias and *abrupt* level shifts (the change-point signal). * **PC3 and higher — noise:** the per-peak, injection-to-injection measurement noise. Because the noise is multiplicative — the per-peak standard deviation grows with the level, \(\mathrm{sd} \approx a\,\mu + b\,\mu^2\) — a large peak would otherwise dominate the PCA. Each peak is therefore referenced to its centre and **dispersion-normalised**: divided by its modelled dispersion so that every peak enters the PCA with comparable variance. This is the variance-stabilising (diagonal-scaling) half of a *whitening* transform; the PCA then supplies the decorrelation, so standardise-then-PCA together sphere the per-peak signals. The whitening is *across peaks* (amplitude), not *in time*: the injection-to-injection noise itself stays reddened (brown, AR(1)-like), as the `Synth_data` fixture deliberately illustrates. ## Estimating the common kappa Kappa is estimated on a **drift-reduced** signal: a robust loess low-pass removes the slow drift on a *copy* of the areas, so \(\kappa\) is not confounded with the drift. The correlated component is then regressed on the pressure deviation with a heavy-tail-robust generalized linear model (`robustbase::glmrob`), which down-weights the tails instead of being dragged by them; the single common slope yields \(\kappa = 1/\sum_k \beta_k\). Robust regression is the appropriate tool here because the corrected-area distribution is markedly leptokurtic. ## Drift and the daily factor (NA-tolerant imputation) The uncorrelated (PC2) daily factors are combined **across cylinders** to estimate the instrument-wide drift and bias. Because the cylinders are not all measured on the same days, the cross-cylinder daily-factor matrix has gaps. APPAC fills them by reconstruction rather than dropping incomplete days: a base-R NA-tolerant PCA (unit-variance scaling, gaps seeded at the mean, then an iterated rank-\(k\) SVD reconstruction — the svdImpute / EM scheme of probabilistic PCA; Tipping & Bishop 1999; Roweis 1998; Dempster, Laird & Rubin 1977). The first two components give the common trend and per-cylinder bias; a parallel decomposition of the correlated part gives the common daily factor. The assembled multiplier (kappa pressure factor + drift + daily factor) is then divided out of the raw areas. ## Refining the centre: chi-square minimisation The correction references each peak to a centre \(\mu_{ij}\); the whole-series mean is a convenient first estimate but a biased one — it is itself pulled by the residual structure the correction has not yet removed. `debias_ct()` refines it by **minimising the residual chi-square**. For a small grid of scale factors \(c\) around 1, it re-runs `appac()` with the centre scaled to \(c\,\mu_{ij}\) and records, per peak, the chi-square of the corrected-area residuals about that centre. That statistic is very nearly quadratic in \(c\), so its minimum is the closed-form parabola vertex \[ \chi^2(c) \;\approx\; b_0 + b_1 c + b_2 c^2 \ (b_2 > 0), \qquad c^{*} \;=\; -\,\frac{b_1}{2\,b_2}, \] with \(b_0, b_1, b_2\) fitted by ordinary least squares over the sweep. The de-biased centre \(c^{*}\mu_{ij}\) trades a little variance for reduced bias and is fed back into a second `appac()` pass. (The closed form replaces a root-find on the derivative, which failed whenever the optimum fell outside the sweep — precisely the large-bias case the refinement exists for.) ## Change points: episode and precision breaks The PC2 drift signal carries **real structural breaks** — cylinder changes, detector maintenance, column trims. `get_changepoints()` detects these episode boundaries with a deterministic structural-break model (an OLS-MOSUM fluctuation test gating BIC-optimal breakpoint dating from **strucchange**; Zeileis et al. 2002). Its second-moment counterpart, `get_variance_changepoints()`, detects changes in measurement **precision** (variance) on the noise-energy signal — a distinct instrument event from a level shift. ```{r changepoints} get_changepoints(fit@samples) # episode level breakpoints ``` Because these breaks are what separates one calibration regime from the next, they are central to the correction rather than an afterthought — and, unlike the single-calibration-period pressure compensation of Ayers & Clardy (1985), let APPAC correct a multi-year record without per-period recalibration. # References * Boček P., Novák J., Janák J. (1969). Effect of pressure on the performance of the flame ionization detector. *J. Chromatogr.* **43**, 431–436. * Ayers B. O., Clardy D. E. (1985). *Barometric pressure correction for a chromatograph.* US Patent 4,512,181. * Tipping M. E., Bishop C. M. (1999). Probabilistic principal component analysis. *J. R. Stat. Soc. B* **61**(3), 611–622. * Roweis S. (1998). EM algorithms for PCA and SPCA. *Adv. Neural Inf. Process. Syst.* **10**. * Dempster A. P., Laird N. M., Rubin D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. *J. R. Stat. Soc. B* **39**(1), 1–38. * Zeileis A., Leisch F., Hornik K., Kleiber C. (2002). strucchange: An R package for testing for structural change in linear regression models. *J. Stat. Softw.* **7**(2).