---
title: "3. Calibration curves and figures of merit"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{3. Calibration curves and figures of merit}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4,
  out.width = "100%",
  dpi = 96
)
```

A calibration relates a signal to the concentration of an element, from
samples of known composition, and predicts the concentration of new
samples. This vignette builds univariate calibration curves from line
intensities, with their limits of detection and quantification, and
compares them with a multivariate model of the whole spectrum, whose
figures of merit come from the net analyte signal.

```{r setup, message = FALSE}
library(specProc)
library(dplyr)
library(tidyr)
library(ggplot2)

data("forageLIBS")
spectra_id <- forageLIBS |> select(1:2) |> names()
minerals <- forageLIBS |> select(3:14) |> names()
spectra <- forageLIBS |> select(-all_of(c(spectra_id, minerals)))
```

## Line intensities

We remove the baseline and measure, in every spectrum, the intensities of
a potassium line, an iron line and the carbon line C I 247.86 nm. Carbon
is the main element of the forage matrix, and dividing by its line
compensates for the changes of ablated mass and plasma from spectrum to
spectrum (the internal standard of `vignette("preprocessing")`):

```{r intensities}
corrected <- baseline_arpls(spectra, lambda = 1e5, max.iter = 20)$correction

lines <- c(C = 247.856, K = 769.896, Fe = 259.940)
intensities <- line_intensities(corrected, lines) |>
  select(spectrum, line, intensity) |>
  pivot_wider(names_from = line, values_from = intensity)

data <- forageLIBS |>
  select(all_of(spectra_id), K, Fe) |>
  mutate(K_signal = intensities$K / intensities$C,
         Fe_signal = intensities$Fe / intensities$C)
```

We set aside a quarter of the samples to test the calibrations:

```{r split}
set.seed(1)
test_samples <- sample(unique(data$Sample), round(0.25 * n_distinct(data$Sample)))
calibration <- filter(data, !Sample %in% test_samples)
test <- filter(data, Sample %in% test_samples)
c(calibration = nrow(calibration), test = nrow(test))
```

## A calibration curve for potassium

`calibration_curve()` fits the signal as a function of the concentration,
and returns the coefficients, the figures of merit and tests of linearity:

```{r k-curve}
k_curve <- calibration_curve(calibration, K_signal, K)
k_curve
```

```{r k-numbers, include = FALSE}
k_fom <- k_curve$figures_of_merit
mandel <- k_curve$linearity |> filter(test == "Mandel")
```

-   **The sensitivity** is the slope: the change of the signal for 1% of K.
-   **The limits of detection and quantification**, \(3.3\,\sigma/S\) and
    \(10\,\sigma/S\), use here the residual standard deviation \(\sigma\) of
    the fit. It includes every source of scatter around the curve, not only
    the noise of a blank, so the limits, `r signif(k_fom$lod, 2)`% and
    `r signif(k_fom$loq, 2)`% K, describe this calibration rather than the
    instrument.
-   **Mandel's test** compares the straight line with a quadratic curve.
    Its p-value, `r signif(mandel$p_value, 2)`, points to a slight
    curvature, as expected for a resonance line, which the plasma partly
    reabsorbs at high concentration.

`plot_calibration()` draws the curve with its confidence band (where the
mean signal lies) and its prediction band (where a single new measurement
is expected):

```{r k-plot, fig.height = 4.5}
plot_calibration(k_curve)
```

`predict()` inverts the curve: it gives the concentration of new samples
from their signals, with a prediction interval, and flags the values below
the limits of detection and quantification:

```{r k-predict}
k_predicted <- predict(k_curve, test$K_signal)
head(k_predicted)
sqrt(mean((k_predicted$concentration - test$K)^2))   # RMSEP, % K
```

The limit of quantification of this curve, `r signif(k_fom$loq, 2)`% K,
is above the potassium content of almost every sample, so
`r round(100 * mean(k_predicted$below_loq))`% of the predictions are
flagged `below_loq`: a single line ranks the samples, but does not
quantify them individually with the usual \(10\,\sigma\) margin.

## Trace elements: weights and limits of detection

Iron ranges from `r min(data$Fe, na.rm = TRUE)` to
`r max(data$Fe, na.rm = TRUE)` ppm, and its signal is noisier at high
concentration. An unweighted fit is then dominated by the few rich samples,
and so is its residual standard deviation, which inflates the limit of
detection. Weights `"1/x"` give each sample a weight inversely proportional
to its concentration. With weights, the residual standard deviation depends
on their scale, so the blank standard deviation is taken from the standard
error of the intercept (`lod_method = "intercept"`):

```{r fe-curves}
fe_calibration <- filter(calibration, !is.na(Fe))
fe_unweighted <- calibration_curve(fe_calibration, Fe_signal, Fe)
fe_weighted <- calibration_curve(fe_calibration, Fe_signal, Fe, weights = "1/x",
                                 lod_method = "intercept")

bind_rows(unweighted = fe_unweighted$figures_of_merit,
          weighted = fe_weighted$figures_of_merit, .id = "fit") |>
  select(fit, n, r_squared, lod, loq)
```

```{r fe-plot, fig.height = 4.5}
plot_calibration(fe_weighted, interval = "prediction")
```

The weighted fit gives a limit of detection of
`r round(fe_weighted$figures_of_merit$lod)` ppm, instead of
`r round(fe_unweighted$figures_of_merit$lod)` ppm. Its R² is lower,
because it no longer favors the rich samples, which dominate R².

```{r fe-predict}
fe_test <- filter(test, !is.na(Fe))
fe_predicted <- predict(fe_weighted, fe_test$Fe_signal)
count(fe_predicted, below_loq)
median(abs(fe_predicted$concentration - fe_test$Fe) / fe_test$Fe)   # relative error
```

## A multivariate model and its figures of merit

A single line uses a small part of the spectrum. Partial least squares
(PLS) regression uses all of it: we divide the baseline-corrected spectra by
the carbon line, and choose the number of components by cross-validation
on the calibration samples, with the pls package:

```{r pls}
x <- as.matrix(corrected) / intensities$C
is_calibration <- !data$Sample %in% test_samples
x_calibration <- x[is_calibration, ]

cv <- pls::plsr(K ~ x, data = data.frame(K = data$K[is_calibration], x = I(x_calibration)),
                ncomp = 15, validation = "CV", segments = 10)
ncomp <- pls::selectNcomp(cv, method = "onesigma")
ncomp
```

`nas()` fits the same model and computes its net analyte signal (NAS): the
part of each spectrum that is unique to potassium, orthogonal to the
contributions of the other constituents. The figures of merit follow from
it:

-   **Sensitivity**: the NAS produced by 1% of K.
-   **Selectivity**: the fraction of each spectrum used for prediction.
-   **Limits of detection and quantification**, \(3.3\,\sigma_x/\mathrm{SEN}\)
    and \(10\,\sigma_x/\mathrm{SEN}\), where \(\sigma_x\) is the noise of the
    spectra (`noise`).

We estimate the noise from the three samples measured twice, as the
standard deviation of their difference spectra divided by \(\sqrt{2}\):

```{r nas}
repeated <- data$Sample[duplicated(data$Sample)]
differences <- sapply(repeated, \(s) {
  rows <- which(data$Sample == s)
  x[rows[1], ] - x[rows[2], ]
})
noise <- sd(differences) / sqrt(2)

k_nas <- nas(x_calibration, data$K[is_calibration], ncomp = ncomp, noise = noise)
k_nas$figures_of_merit
```

```{r nas-predict}
nas_predicted <- predict(k_nas, x[!is_calibration, ])
sqrt(mean((nas_predicted$predicted - test$K)^2))   # RMSEP, % K
summary(nas_predicted$selectivity)
```

```{r nas-numbers, include = FALSE}
fom <- k_nas$figures_of_merit
rmsep_pls <- sqrt(mean((nas_predicted$predicted - test$K)^2))
rmsep_line <- sqrt(mean((k_predicted$concentration - test$K)^2))
```

```{r comparison-plot, fig.height = 4, fig.width = 7}
bind_rows(
  `K I 769.90 nm curve` = tibble(laboratory = test$K, predicted = k_predicted$concentration),
  `PLS (NAS)` = tibble(laboratory = test$K, predicted = nas_predicted$predicted),
  .id = "model"
) |>
  ggplot(aes(laboratory, predicted)) +
  geom_abline(colour = "grey50", linetype = "dashed") +
  geom_point(alpha = 0.6) +
  facet_wrap(~ model) +
  coord_equal() +
  labs(x = "Laboratory K (%)", y = "Predicted K (%)") +
  theme_bw()
```

The PLS model predicts the test samples with an error of
`r sprintf("%.2f", rmsep_pls)`% K, against `r sprintf("%.2f", rmsep_line)`%
for the calibration curve. Its limit of detection, `r signif(fom[["lod"]], 2)`% K,
accounts for the noise of the spectra only, so it is a lower bound, far
below the limit of the calibration curve, which includes all the scatter of
the samples. Its mean selectivity, `r round(fom[["selectivity"]], 2)`,
means that only about a tenth of each spectrum is specific to potassium:
the rest is shared with the other constituents and the plasma.

## Summary

-   **Normalize line intensities** to an internal standard, such as the
    carbon line of an organic matrix.
-   **Check the linearity** with Mandel's test and the lack-of-fit test of
    `calibration_curve()`; resonance lines often curve at high
    concentration.
-   **Weight the fit** when the signal is noisier at high concentration, and
    take the limit of detection from the intercept or from blanks.
-   **Report how the limits were estimated**: from the residuals of a curve,
    they include the scatter of the samples; from the noise of the spectra
    (`nas()`), they are lower bounds.
-   **Validate on samples left out** of the calibration: the prediction
    error on test samples is the most honest figure of merit.

The last vignette, `vignette("plasma-diagnostics")`, measures the plasma
itself: electron density, temperature and self-absorption.