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.

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))
#> 15473 injections, 6 cylinders, 5 peaks;  P_ref = 1005.4 hPa

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.)

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 = -6.902e-04  (per hPa)

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.

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.

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):

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))
#>            peak RSD_raw_pct RSD_corrected_pct
#> CH4         CH4       0.653             0.142
#> C2H6       C2H6       0.772             0.145
#> C3H8       C3H8       0.791             0.144
#> i.C4H10 i.C4H10       0.817             0.162
#> n.C4H10 n.C4H10       0.810             0.158

Reduced chi-square ≈ 1 means the corrected areas are down to the short-term measurement-noise floor:

goodness_of_fit(fit)[[1]]
#>         reduced.chisq    chisq  dof      p.value
#> CH4          1.113141 6100.011 5480 5.358257e-09
#> C2H6         1.196429 6556.433 5480 1.778801e-22
#> C3H8         1.128207 6182.574 5480 5.585007e-11
#> i.C4H10      1.286267 7048.742 5480 2.197287e-43
#> n.C4H10      1.251244 6856.819 5480 1.674147e-34

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.

plot_residuals(fit, sample = 1, peak = "n.C4H10")

A compact summary of the fitted object (print() adds the per-sample goodness of fit):

fit
#> <Appac> atmospheric-pressure peak-area correction
#>   samples: 6 (CTL.1, CTL.2, CTL.3, CTL.4, CTL.5, CTL.6)
#>   peaks  : 5 (CH4, C2H6, C3H8, i.C4H10, n.C4H10)
#>   data   : 15473 injections, 2018-12-20 to 2024-02-16
#>   kappa[pressure]: -0.0006902  (P_ref = 1005.44)

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.

get_changepoints(fit@samples)            # episode level breakpoints
#>  [1] "2019-05-11" "2019-09-27" "2020-01-06" "2020-01-24" "2020-02-01"
#>  [6] "2020-05-10" "2020-06-19" "2021-04-09" "2021-04-13" "2021-11-27"
#> [11] "2022-07-26" "2022-09-03" "2023-05-01" "2023-05-28" "2023-06-25"
#> [16] "2023-07-21" "2023-09-26" "2023-11-21"

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.