Package 'specProc'

Description

This function generates the adjusted boxplot, which is a robust graphical method for visualizing skewed data distributions. It provides a more accurate representation of the data's spread and skewness compared to standard boxplot, especially in the presence of outliers.

Usage

adjusted_boxplot(
  x,
  plot = TRUE,
  id = NULL,
  group = NULL,
  scales = c("free_y", "fixed"),
  points = c("outliers", "all", "none"),
  label_outliers = FALSE,
  show_n = TRUE,
  show_mean = FALSE,
  annotate = NULL,
  horizontal = FALSE,
  log = FALSE,
  fill = "grey85",
  xlab = NULL,
  ylab = NULL,
  title = NULL,
  caption = TRUE,
  base_size = 11,
  x_labels_angle = 0,
  box_width = 0.5,
  notch = FALSE,
  notch_width = 0.5,
  staple_width = 0.5,
  xlabels.angle = deprecated(),
  xlabels.vjust = deprecated(),
  xlabels.hjust = deprecated(),
  box.width = deprecated(),
  notchwidth = deprecated(),
  staplewidth = deprecated()
)

Arguments

Details

The function is based on the medcouple (MC) measure computed on the data and which robustly measures skewness. This measure is bounded between -1 and 1. The medcouple is equal to zero when the observed distribution is symmetric, whereas a positive (resp. negative) value of MC corresponds to a right (resp. left) tailed distribution. The fences are, for $\text{MC} \geq 0$,

$[Q_1 - 1.5 e^{-4\,\text{MC}}\,\text{IQR};\; Q_3 + 1.5 e^{3\,\text{MC}}\,\text{IQR}]$

and, for $\text{MC} < 0$,

$[Q_1 - 1.5 e^{-3\,\text{MC}}\,\text{IQR};\; Q_3 + 1.5 e^{4\,\text{MC}}\,\text{IQR}],$

calibrated by Hubert and Vandervieren (2008) so that about 0.7% of the observations of a clean, moderately skewed distribution are flagged, as with Tukey's boxplot for normal data. The whiskers end at the most extreme observations within the fences, and the observations outside are drawn as points. The statistics are computed by robustbase::adjboxStats(), whose quartiles are Tukey's hinges, as in boxplot(). The method suits distributions that are not excessively skewed, i.e. with $|\text{MC}| \leq 0.6$.

Layout. Each variable is drawn in its own panel, with its own value axis (scales = "free_y"), so that variables of different scales or units stay readable; scales = "fixed" draws them in one panel, on a common axis. With group, the boxes of the groups are side by side in the panel of each variable. Below each box, n is the number of non-missing values.

Points. By default, the outlying observations are drawn, slightly spread sideways so that equal values do not hide each other. points = "all" draws every observation, the outlying ones in black, as some journals require; label_outliers = TRUE names the outlying observations by id (else by row number), and the outliers table of plot = FALSE gives them all.

Annotations. annotate adds, for each box:

• "shape": the medcouple (MC, marked ⁠(> 0.6)⁠ beyond the range of the method), or the g and h of the generalized boxplot, which explain the asymmetry of the whiskers;

• "outliers": the number of flagged values, with the number expected in clean data (about 0.7% for the adjusted boxplot, which Hubert and Vandervieren calibrated as Tukey's boxplot for normal data, and alpha for the generalized one): many more flagged values than expected point to real outliers;

• "fences": the fences beyond which values are flagged, as dashed marks (on the sides with outlying values), the whiskers ending at the last values inside;

• "median": the median with its distribution-free 95% confidence interval, from the order statistics (Md); the notches are made to compare two medians, and are not a confidence interval of one;

• "spread": the interquartile range and the robust coefficient of variation $\text{rCV} = \text{IQR} / (1.349\,\text{median})$, in percent, read as a relative standard deviation;

• "location": the biweight location (biweight_location()), a robust mean, marked by a cross (the mean, with show_mean, is pulled by the outlying values);

• "test": with group, the p-value of the Kruskal-Wallis test between the groups, below the variable name. With several variables, adjust the p-values for multiple testing (see p.adjust() and the tests table);

• "missing": the number of missing values, next to n.

annotate = "all" adds them all (the test only with group). The caption gives their key. With several boxes in a panel, each statistic takes two short lines: for many groups, make the figure wider or choose fewer annotations. All the values are in the stats and tests tables of plot = FALSE.

Publication figures. Give the units of the variables in their names (such as "Ca (%)", with check.names = FALSE in data.frame()) or in ylab. Make the plot at its printed size, with base_size set to the text size of the journal (often 7 to 9 points), and save it with ggplot2::ggsave(), for example ggsave("figure.pdf", p, width = 85, height = 70, units = "mm") for a single column. The caption says which boxplot is drawn, because readers assume Tukey's whiskers (1.5 IQR) otherwise; use caption = FALSE and state it in the figure legend instead, for example: "Adjusted boxplots (Hubert and Vandervieren, 2008): the boxes show the quartiles and the median, the whiskers extend to the most extreme values within the medcouple-adjusted fences, and the points are the values outside."

Value

• If plot = TRUE, a ggplot2 object.

• If plot = FALSE, a list of tibbles: stats, with one row per variable (and group): the numbers of values n and of missing values n_missing, the whisker ends lower and upper (the most extreme values within the fences), the quartiles q1 and q3, the median, the fences, the notch limits, the medcouple, the 95% confidence interval of the median (median_lower, median_upper), the iqr, the robust coefficient of variation rcv (%), the biweight location, the mean, and the numbers of outlying values n_outliers and of expected ones in clean data expected_outliers; outliers, with the outlying values, their row in x, their id and their tail out ("lower" or "upper"); and, with group, tests, the Kruskal-Wallis test between the groups of each variable (statistic, df, p_value).

Author(s)

Christian L. Goueguel

References

The adjusted boxplot is based on the methodology described in:

• Brys, G., Hubert, M., Struyf, A., (2004). A Robust Measure of Skewness. Journal of Computational and Graphical Statistics, 13(4):996-1017

• Hubert, M., Vandervieren, E., (2008). An adjusted boxplot for skewed distributions. Computational Statistics and Data Analysis, 52(12):5186-5201

• McGill, R., Tukey, J. W., Larsen, W. A. (1978). Variations of box plots. The American Statistician, 32(1):12-16 (notches)

See Also

generalized_boxplot() for skewed and heavy-tailed distributions.

Examples

data(forageLIBS)
minerals <- forageLIBS[c("Measurement", "Ca", "Mg", "P", "K", "S")]

# mineral contents, each on its own axis, the outlying samples named: a
# figure of a journal page width
p <- adjusted_boxplot(minerals, id = Measurement, label_outliers = TRUE,
                      ylab = "Content (%)", base_size = 8)
p
# ggplot2::ggsave("minerals.pdf", p, width = 175, height = 70, units = "mm")

# the statistics and the outlying samples
res <- adjusted_boxplot(minerals, id = Measurement, plot = FALSE)
res$stats
res$outliers

# potassium and phosphorus by calcium level, with notches and means: a
# figure of a journal column
minerals$Ca_level <- cut(minerals$Ca, quantile(minerals$Ca, 0:3 / 3),
                         labels = c("Low Ca", "Mid Ca", "High Ca"),
                         include.lowest = TRUE)
p <- adjusted_boxplot(minerals[c("K", "P", "Ca_level")], group = Ca_level,
                      notch = TRUE, show_mean = TRUE, ylab = "Content (%)",
                      title = "Potassium and phosphorus by calcium level",
                      base_size = 8)
p
# ggplot2::ggsave("by_calcium.pdf", p, width = 85, height = 75, units = "mm")

# trace elements on a logarithmic axis, horizontal, with all the samples
traces <- forageLIBS[c("Fe", "Mn", "Zn")]
adjusted_boxplot(traces, log = TRUE, horizontal = TRUE, points = "all",
                 ylab = "Content (mg/kg)")

Description

Relabels the wavelengths of spectra with the correction fitted by wavelength_calibration(). The intensities are not changed.

Usage

apply_calibration(spectra, calibration)

Arguments

Details

The spectra must have the wavelength axis of the calibration (the same channels, in the same order), so that each channel gets the correction of its detector segment. The corrected wavelengths are rounded to 4 decimals (0.1 pm) for the column names.

Value

spectra, with the corrected wavelengths as names.

See Also

wavelength_calibration()

Examples

data(forageLIBS)
spectra <- forageLIBS[-(1:14)]
wl <- as.numeric(names(spectra))
first_detector <- spectra[seq_len(which(diff(wl) < 0)[1])]
reference <- c(`Mg I` = 285.213, `Ca II` = 317.933, `Ca II` = 393.366, `K I` = 404.414,
               `Ca I` = 422.673, `Mg I` = 518.360, `Na I` = 588.995, `Ca I` = 643.907)
cal <- wavelength_calibration(first_detector, reference)
corrected <- apply_calibration(first_detector, cal)
head(names(corrected))

Description

This function efficiently averages the samples spectra in a large dataset. It provides flexibility by computing either the overall mean across all spectra or group-wise means based on the values of a specified grouping column.

Usage

average(x, .group_by = NULL)

Arguments

Details

The function leverages the power of Rcpp to perform the mean calculations in C++. The underlying C++ implementation has a time complexity of O(n x m), where n is the number of rows and m is the number of columns in the data. Missing values are ignored in the computation of each mean.

Value

• If .group_by = NULL, a one-row tibble containing the mean of each column.

• If .group_by is provided, a tibble with one row per group: the first column holds the group labels and the remaining columns hold the group means.

Author(s)

Christian L. Goueguel

Examples

data(forageLIBS)
# one mean spectrum per sample: three samples were measured twice
means <- average(forageLIBS[-c(1, 3:14)], Sample)
dim(means)

Description

Baseline correction based on asymmetric least squares (ALS) algorithm as proposed by Eilers et al. (2005).

Usage

baseline_als(x, lambda = 1000, p = 0.001, max.iter = 10)

Arguments

Details

The function applies Eilers' method based on a Whittaker filter. The algorithm estimates a baseline curve by minimizing the asymmetric least squares criterion, which allows for different weights for positive and negative residuals. The resulting baseline curve is subtracted from the input data, providing a baseline-corrected version.

The penalized system is pentadiagonal and is solved in C++ with a banded Cholesky decomposition, so the cost grows linearly with the number of spectral channels. Negative values in the corrected spectra are kept as they are (they typically reflect noise around the baseline).

Value

A list containing two tibbles:

• correction: The baseline-corrected spectral matrix.

• background: The fitted background emission.

Author(s)

Christian L. Goueguel

References

• Eilers, P.H.C., Boelens, H.F.M., (2005). Baseline correction with asymmetric least squares smoothing. Leiden University Medical Centre report.

Examples

data(forageLIBS)
spectrum <- forageLIBS[1, -(1:14)]
wl <- as.numeric(names(spectrum))
in_region <- wl > 240 & wl < 300
region <- spectrum[in_region]
res <- baseline_als(region, lambda = 1e6)
oldpar <- par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))
# the spectrum and its fitted baseline
plot(wl[in_region], unlist(region), type = "l", col = "grey40", ylim = c(800, 3000),
     xlab = "Wavelength (nm)", ylab = "Counts", main = "Spectrum and baseline")
lines(wl[in_region], unlist(res$background), col = "red")
# the corrected spectrum: the background is now around zero
plot(wl[in_region], unlist(res$correction), type = "l", col = "grey40",
     ylim = c(-200, 2000), xlab = "Wavelength (nm)", ylab = "Counts",
     main = "Corrected spectrum")
abline(h = 0, col = "red", lty = 2)
par(oldpar)

Description

Baseline correction based on asymmetrically reweighted penalized least squares smoothing algorithm as proposed by Baek et al. (2015).

Usage

baseline_arpls(x, lambda = 1000, ratio = 0.05, max.iter = 10)

Arguments

Details

The algorithm estimates a baseline curve by iteratively updating weights based on the residuals and minimizing a penalized least squares criterion. The resulting baseline curve is subtracted from the input data, providing a baseline-corrected version.

The penalized system is pentadiagonal and is solved in C++ with a banded Cholesky decomposition, so the cost grows linearly with the number of spectral channels.

Value

A list containing two tibbles:

• correction: The baseline-corrected spectral matrix.

• background: The fitted background emission.

Author(s)

Christian L. Goueguel

References

• Baek, S.-J., Park, A., Ahn, Y.-J., Choo, J., (2015). Baseline correction using asymmetrically reweighted penalized least squares smoothing. Analyst, 140(1):250–257.

Examples

data(forageLIBS)
spectrum <- forageLIBS[1, -(1:14)]
wl <- as.numeric(names(spectrum))
in_region <- wl > 240 & wl < 300
region <- spectrum[in_region]
res <- baseline_arpls(region, lambda = 1e5)
oldpar <- par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))
# the spectrum and its fitted baseline
plot(wl[in_region], unlist(region), type = "l", col = "grey40", ylim = c(800, 3000),
     xlab = "Wavelength (nm)", ylab = "Counts", main = "Spectrum and baseline")
lines(wl[in_region], unlist(res$background), col = "red")
# the corrected spectrum: the background is now around zero
plot(wl[in_region], unlist(res$correction), type = "l", col = "grey40",
     ylim = c(-200, 2000), xlab = "Wavelength (nm)", ylab = "Counts",
     main = "Corrected spectrum")
abline(h = 0, col = "red", lty = 2)
par(oldpar)

Description

A dials parameter for the lambda argument of step_baseline(), on a log10 scale.

Usage

baseline_lambda(range = c(2, 8), trans = scales::transform_log10())

Arguments

Value

A quant_param object.

Examples

baseline_lambda()

Description

This function performs baseline correction on the spectral matrix by estimating and removing the continuous background emission using least-squares polynomial curve fitting approach.

Usage

baseline_lsp(
  x,
  degree = 4,
  tol = 0.001,
  max.iter = 100,
  method = c("imodpoly", "modpoly")
)

Arguments

Details

Two iterative polynomial fits are available:

• "modpoly": the modified polynomial fit of Lieber and Mahadevan-Jansen (2003). A polynomial is fitted to the spectrum, the points of the spectrum above the fit are replaced by the fit, and the fit is repeated, so that the peaks are progressively removed. The first fits are pulled up by strong peaks, and the ends of the spectrum, where nothing pulls the polynomial back, can then fall far below the background.

• "imodpoly" (default): the improved modified polynomial fit of Zhao et al. (2007). The points more than one standard deviation of the residuals above the first fit (the peaks) are left out of the following fits, and the remaining points are clipped at the fit plus the standard deviation of the residuals, rather than at the fit. The iterations stop when this standard deviation changes by less than tol (relative). The baseline follows the background up to the ends of the spectrum and passes through its noise.

Value

A list with two elements:

• correction: The baseline-corrected spectral matrix.

• background: The fitted background emission.

Author(s)

Christian L. Goueguel

References

• Lieber, C.A., Mahadevan-Jansen, A., (2003). Automated method for subtraction of fluorescence from biological Raman spectra. Applied Spectroscopy, 57(11):1363-1367

• Zhao, J., Lui, H., McLean, D.I., Zeng, H., (2007). Automated autofluorescence background subtraction algorithm for biomedical Raman spectroscopy. Applied Spectroscopy, 61(11):1225-1232.

Examples

data(forageLIBS)
spectrum <- forageLIBS[1, -(1:14)]
wl <- as.numeric(names(spectrum))
in_region <- wl > 240 & wl < 300
region <- spectrum[in_region]
res <- baseline_lsp(region, degree = 4)
oldpar <- par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))
# the spectrum and its fitted baseline
plot(wl[in_region], unlist(region), type = "l", col = "grey40", ylim = c(800, 3000),
     xlab = "Wavelength (nm)", ylab = "Counts", main = "Spectrum and baseline")
lines(wl[in_region], unlist(res$background), col = "red")
# the corrected spectrum: the background is now around zero
plot(wl[in_region], unlist(res$correction), type = "l", col = "grey40",
     ylim = c(-200, 2000), xlab = "Wavelength (nm)", ylab = "Counts",
     main = "Corrected spectrum")
abline(h = 0, col = "red", lty = 2)
par(oldpar)

Description

Reduces the number of channels of spectra by replacing each group of width adjacent channels by their mean.

Usage

bin_spectra(x, width = 3, segments = TRUE)

Arguments

Details

Averaging keeps the intensities on their scale and divides the noise of independent channels by about sqrt(width), at the cost of spectral resolution: keep width smaller than the width of the emission lines (a few channels in LIBS spectra) to keep them resolved. Fewer channels also make the models faster to fit.

The channels are grouped from the first one; when their number is not a multiple of width, the last group is smaller. When segments = TRUE, the spectrum is split between the detectors of a multi-spectrometer system (where the wavelengths step back, or jump by more than 5 times the median spacing), and each detector is binned on its own, so that no group straddles two detectors. The binned channels are named by their mean wavelength (to 10 significant digits), or, when the names are not wavelengths, by the name of their first channel; a channel left alone in its group keeps its name.

Value

The binned spectra, of the same type as x (a tibble for a data frame), with ceiling(n / width) channels for each detector segment of n channels.

Author(s)

Christian L. Goueguel

References

• Vrábel, J., Képeš, E., Duponchel, L., et al. (2020). Classification of challenging laser-induced breakdown spectroscopy soil sample data - EMSLIBS contest. Spectrochimica Acta Part B, 169:105872.

See Also

step_bin_spectra(), gap_derivative(), average()

Examples

data(forageLIBS)
spectra <- forageLIBS[1:3, -(1:14)]
binned <- bin_spectra(spectra, width = 3)
ncol(spectra)
ncol(binned)
names(binned)[1:3]

Description

This function computes the biweight location, a robust measure of central tendency for a numeric vector. The biweight location is less sensitive to outliers than the sample mean.

Usage

biweight_location(
  x,
  loc = NULL,
  c = 6,
  tol = 1e-06,
  max_iter = 50,
  drop.na = FALSE
)

Arguments

Details

Starting from the initial location $M$ (the median by default), the biweight location is

$\zeta = M + \frac{\sum_{|u_i|<1} (x_i - M)(1 - u_i^2)^2}{\sum_{|u_i|<1} (1 - u_i^2)^2}, \quad u_i = \frac{x_i - M}{c \cdot \mathrm{MAD}}$

where MAD is the (unscaled) median absolute deviation. The estimate is iterated, replacing $M$ by $\zeta$, until the change is smaller than tol or max_iter iterations are reached.

Value

The biweight location of x.

Author(s)

Christian L. Goueguel

References

• Mosteller, F., and Tukey, J. W. (1977). Data Analysis and Regression: A Second Course in Statistics. Addison-Wesley, pp. 203-209.

• Beers, T.C., Flynn, K., Gebhardt, K., (1990). Measures of location and scale for velocities in clusters of galaxies - A robust approach. The Astronomical Journal, 100:32-46.

Examples

data(forageLIBS)
# iron contents (mg/kg): a few samples are far above the others
fe <- forageLIBS$Fe[!is.na(forageLIBS$Fe)]
c(mean = mean(fe), median = stats::median(fe), biweight = biweight_location(fe))

Description

This function computes the biweight midcorrelation between two numeric vectors. The biweight midcorrelation is a robust measure of correlation that is less sensitive to outliers than the traditional Pearson's correlation coefficient.

Usage

biweight_midcorrelation(x, y)

Arguments

Details

The biweight midcorrelation is calculated using the biweight midvariances and biweight midcovariance, as described by Wilcox (1994). It is bounded between -1 and 1.

Value

The biweight midcorrelation between x and y.

Author(s)

Christian L. Goueguel

References

• Wilcox, R. R. (1994). The Biweight Midcorrelation: A Robust Correlation Technique for Two Samples. Journal of Statistical Computation and Simulation, 48(2):103-110.

Examples

data(forageLIBS)
ok <- !is.na(forageLIBS$Fe)
c(pearson = stats::cor(forageLIBS$Fe[ok], forageLIBS$Mn[ok]),
  biweight = biweight_midcorrelation(forageLIBS$Fe[ok], forageLIBS$Mn[ok]))

Description

This function computes the biweight midcovariance, a robust measure of covariance between two numerical vectors. The biweight midcovariance is less sensitive to outliers than the traditional covariance.

Usage

biweight_midcovariance(x, y)

Arguments

Value

The biweight midcovariance between x and y.

Author(s)

Christian L. Goueguel

References

• Wilcox, R., (1997). Introduction to Robust Estimation and Hypothesis Testing. Academic Press

Examples

data(forageLIBS)
# iron and manganese contents (mg/kg), with a few iron-rich samples
ok <- !is.na(forageLIBS$Fe)
c(covariance = stats::cov(forageLIBS$Fe[ok], forageLIBS$Mn[ok]),
  biweight = biweight_midcovariance(forageLIBS$Fe[ok], forageLIBS$Mn[ok]))

Description

This function calculates the biweight midvariance of a numeric vector, which is a robust measure of scale that can be used to estimate the variability of the data while being resistant to the influence of outliers.

Usage

biweight_midvariance(x, drop.na = FALSE)

Arguments

Details

For scale estimators, the standard deviation (or variance) is the optimal estimator for Gaussian data. However, it is not resistant and it does not have robustness of efficiency. In robust statistics, the median absolute deviation (MAD) is a resistant estimate, but it has only modest robustness of efficiency, while the biweight midvariance estimator is both resistant and robust of efficiency.

Value

The biweight midvariance of the input vector.

Author(s)

Christian L. Goueguel

References

• Wilcox, R., (1997). Introduction to Robust Estimation and Hypothesis Testing. Academic Press

Examples

data(forageLIBS)
# iron contents (mg/kg): a few samples are far above the others
fe <- forageLIBS$Fe[!is.na(forageLIBS$Fe)]
c(variance = stats::var(fe), biweight = biweight_midvariance(fe))

Description

This function computes the biweight scale, a robust measure of scale or dispersion for a numeric vector. The biweight scale is less sensitive to outliers than the sample standard deviation.

Usage

biweight_scale(x, loc = NULL, c = 9, reduced = FALSE, drop.na = FALSE)

Arguments

Details

The biweight scale is the square root of the biweight midvariance:

$\zeta = \frac{\sqrt{n \sum_{|u_i|<1} (x_i - M)^2 (1 - u_i^2)^4}} {\left| \sum_{|u_i|<1} (1 - u_i^2)(1 - 5u_i^2) \right|}, \quad u_i = \frac{x_i - M}{c \cdot \mathrm{MAD}}$

where $M$ is the location (the median by default) and MAD is the (unscaled) median absolute deviation. For normally distributed data the biweight scale is a consistent estimator of the standard deviation.

Value

The biweight scale of x.

Author(s)

Christian L. Goueguel

References

• Mosteller, F., and Tukey, J. W. (1977). Data Analysis and Regression: A Second Course in Statistics. Addison-Wesley, pp. 203-209.

• Beers, T.C., Flynn, K., Gebhardt, K., (1990). Measures of location and scale for velocities in clusters of galaxies - A robust approach. The Astronomical Journal, 100:32-46.

Examples

data(forageLIBS)
# iron contents (mg/kg): a few samples are far above the others
fe <- forageLIBS$Fe[!is.na(forageLIBS$Fe)]
c(sd = stats::sd(fe), mad = stats::mad(fe), biweight = biweight_scale(fe))

Description

Estimates the excitation temperature of a plasma in local thermodynamic equilibrium (LTE) from the intensities of several lines of the same species, by a Boltzmann plot.

Usage

boltzmann(lines, units = "energy")

Arguments

Details

In LTE, the integrated intensity of a line emitted from an upper level of energy $E_k$ and degeneracy $g_k$, with transition probability $A_{ki}$ and wavelength $\lambda$, satisfies

$\ln\frac{I\lambda}{g_k A_{ki}} = -\frac{E_k}{k_B T} + C$

for intensities in energy units, or $\ln(I / g_k A_{ki})$ for photon counts (units = "photons"). A least-squares line through the points gives the temperature from its slope, $T = -1/(k_B\,\mathrm{slope})$, and its standard error from that of the slope.

The intensities must be corrected for the spectral response of the instrument, and the lines must be optically thin (see self_absorption()). A wide range of upper-level energies gives a more precise temperature. Atomic data ($A_{ki}$, $g_k$, $E_k$) can be taken from the NIST Atomic Spectra Database (https://physics.nist.gov/PhysRefData/ASD/lines_form.html).

Value

An object of class specproc_boltzmann, a list with

• temperature and temperature_se: the temperature and its standard error, in K;

• r_squared: the coefficient of determination of the fit;

• points: a tibble with the abscissa x (eV) and ordinate y of each line;

• fit: the lm fit of y on x.

Use plot_boltzmann() to draw the plot.

Author(s)

Christian L. Goueguel

References

• Cristoforetti, G., De Giacomo, A., Dell'Aglio, M., Legnaioli, S., Tognoni, E., Palleschi, V., Omenetto, N. (2010). Local thermodynamic equilibrium in laser-induced breakdown spectroscopy: beyond the McWhirter criterion. Spectrochimica Acta Part B, 65(1):86-95.

See Also

saha_boltzmann(), plot_boltzmann(), mcwhirter_criterion()

Examples

data(forageLIBS)
mean_spectrum <- colMeans(forageLIBS[-(1:14)])
# calcium lines, with their atomic data from the NIST database
atomic <- data.frame(
  species = c(rep("Ca I", 6), rep("Ca II", 4)),
  stage = c(rep(1, 6), rep(2, 4)),
  wavelength = c(428.301, 430.253, 431.865, 445.478, 612.222, 616.217,
                 315.887, 317.933, 370.603, 373.690),
  Aki = c(4.34e7, 1.36e8, 7.40e7, 8.70e7, 2.87e7, 4.77e7, 3.10e8, 3.60e8, 8.80e7, 1.70e8),
  gk = c(5, 5, 3, 7, 3, 3, 4, 6, 2, 2),
  Ek = c(4.780, 4.780, 4.769, 4.681, 3.910, 3.910, 7.047, 7.050, 6.468, 6.468)
)
lines <- line_intensities(mean_spectrum, atomic, baseline = TRUE)
fit <- boltzmann(lines[lines$stage == 1, ])
fit
plot_boltzmann(fit)

Description

Fits a univariate calibration curve (signal against concentration), with its figures of merit, limits of detection and quantification, and tests of linearity. The curve predicts the concentration of new samples from their signal, with confidence intervals.

Usage

calibration_curve(
  data,
  signal,
  concentration,
  model = "linear",
  weights = NULL,
  blank = NULL,
  lod_method = NULL,
  level = 0.95
)

Arguments

Details

The signal, such as a line intensity from line_intensities() or a line ratio, is modeled as a straight line or a second-degree polynomial of the concentration, fitted by (weighted) least squares. Weights "1/x" or "1/x2" suit signals whose variance grows with concentration, which is common in LIBS.

Figures of merit. The sensitivity is the slope of the curve at zero concentration. The limits of detection and quantification are $LOD = 3.3\,\sigma / S$ and $LOQ = 10\,\sigma / S$ (ICH, 2005), where $S$ is the sensitivity and $\sigma$ the standard deviation of the signal of a blank, estimated from

• "residual": the residual standard deviation of the fit (default without blank);

• "intercept": the standard error of the intercept;

• "blank": the standard deviation of replicate blank signals (default when given).

With weights, the residual standard deviation depends on the scale of the weights, so prefer "blank" or "intercept".

Linearity. For a straight line, Mandel's test compares its residual variance with that of a quadratic fit; a small p-value means that the curvature is significant (as with self-absorption or detector saturation). When concentrations are replicated, the lack-of-fit test compares the residuals with the pure error of the replicates.

Intervals. The confidence band of the curve shows where the mean signal lies at each concentration; the prediction band, wider, where a single new measurement is expected to fall, since it adds the noise of the measurement. plot_calibration() draws either or both, and predict() computes them for given concentrations (type = "signal").

Inverse prediction. The concentration of a sample is the solution of the calibration equation for its signal (for a quadratic curve, the root within or nearest to the calibration range). Its standard error, by the delta method, combines the uncertainty of the curve with that of the new signal, the mean of replicates measurements.

Value

An object of class specproc_calibration, a list with

• coefficients: a tibble of the estimates, standard errors and confidence intervals (lower, upper) at level;

• figures_of_merit: a tibble with the number of standards n, the sensitivity, r_squared, the residual standard deviation sigma, the blank standard deviation used (sigma_blank), lod and loq (in concentration units);

• linearity: a tibble with the statistic, degrees of freedom and p-value of Mandel's test and of the lack-of-fit test (when they apply);

• fit: the lm fit, and data, the calibration data.

Use predict() for new samples and plot_calibration() to draw the curve.

Author(s)

Christian L. Goueguel

References

• ICH (2005). Validation of analytical procedures: text and methodology Q2(R1). International Conference on Harmonisation.

• Mandel, J. (1964). The Statistical Analysis of Experimental Data. Interscience, New York.

• Miller, J.N., Miller, J.C. (2018). Statistics and Chemometrics for Analytical Chemistry, 7th ed. Pearson, Harlow.

See Also

predict.specproc_calibration(), plot_calibration(), line_intensities(), nas()

Examples

data(forageLIBS)
# the K I 769.90 nm line, normalized to the C I 247.86 nm line of the matrix
lines <- line_intensities(forageLIBS[-(1:14)], c(C = 247.856, K = 769.896), baseline = TRUE)
standards <- data.frame(
  K = forageLIBS$K,
  signal = lines$intensity[lines$line == "K"] / lines$intensity[lines$line == "C"]
)
cal <- calibration_curve(standards[1:300, ], signal, K)
cal
# the concentrations of new samples, from their signals
head(predict(cal, standards$signal[301:368]))

Description

Robust PCA that handles outlying observations (casewise outliers), outlying cells (cellwise outliers) and missing values by minimizing a single objective function: the cellPCA method of Centofanti, Hubert and Rousseeuw. Each cell and each observation gets a weight between 0 and 1 that reflects its outlyingness, and regular cells and observations are not downweighted, which makes cellPCA more efficient than macropca(). The iterations are computed in C++. The results have the same form as those of robpca(), so that predict(), plot_outlier_map() and plot_cell_map() work the same way.

Usage

cellpca(
  x,
  k = NULL,
  alpha = 0.5,
  kmax = 10,
  var_explained = 0.8,
  scale = FALSE,
  ndir = 250,
  maxiter = 1000,
  tol = 1e-06,
  max_col_frac = 0.5,
  od_cutoff = c("simulated", "chisq")
)

Arguments

Details

cellPCA approximates the data by a fit $\hat{X} = 1_n \mu^T + U V^T$ of rank k that minimizes

$\hat\sigma_2^2 \frac{1}{n} \sum_{i=1}^n m_i \rho_2\left(\frac{t_i}{\hat\sigma_2}\right), \quad t_i = \sqrt{\frac{1}{m_i} \sum_{j=1}^p m_{ij} \hat\sigma_{1,j}^2 \rho_1\left(\frac{x_{ij} - \hat{x}_{ij}}{\hat\sigma_{1,j}}\right)}$

where $m_{ij}$ is 0 for a missing cell and 1 otherwise and $m_i$ is the number of observed cells of observation $i$. The bounded function $\rho_1$ limits the effect of the outlying cells (the residuals of variable $j$ divided by their scale $\hat\sigma_{1,j}$), and $\rho_2$ that of the outlying observations (the casewise total deviations $t_i$ divided by their scale $\hat\sigma_2$). Both are hyperbolic tangent functions (Hampel et al., 1981): $\rho_1$ with $b = 1.5$ and $c = 4$, and $\rho_2$ with $b$ and $c$ the 0.70 and 0.99 quantiles of the standardized total deviations of simulated Gaussian residuals. Their weights $\psi(z)/z$ are 1 in the central region and 0 beyond $c$.

The algorithm follows the reference code of the authors:

1. Initial fit. A macropca() fit (with alpha, kmax, var_explained, scale and ndir), which also chooses k when it is NULL, with its scores computed from the data whose flagged and missing cells are imputed by the fit.

2. Scales. $\hat\sigma_{1,j}$ is the M-scale of the residuals of variable $j$, and $\hat\sigma_2$ that of the total deviations (with $\rho_{1.5,4}$: 50% breakdown, consistent at the normal distribution). They are kept fixed.

3. Iteratively reweighted least squares. The loadings (one variable at a time, with the casewise and cellwise weights) and the scores (one observation at a time, with the cellwise weights) are updated by weighted least squares, the loadings are orthonormalized, the center is updated, and so are the weights, until the fit $U V^T$ changes by less than tol (relative), at most maxiter times. Each iteration decreases the objective. If more than max_col_frac of the cells of a variable get a zero weight, the previous iteration is kept.

4. Principal directions. The center and the eigenvectors within the subspace are estimated by the deterministic MCD of the scores of the observations with a non-zero casewise weight (the exact MCD when k = 1), and the sign of each loading vector is set so that its largest element is positive.

The residuals of each variable are then standardized by their M-scale, and the cells beyond $\sqrt{\chi^2_{1, 0.99}}$ are flagged. As in the enhanced outlier map of the paper, od is the norm of the standardized residuals of each observation, and sd the score distance of its projection on the subspace (of its robust scores when it has missing cells). The cut-off of sd is $\sqrt{\chi^2_{k, 0.99}}$. That of od is, by default (od_cutoff = "simulated"), the 0.99 quantile of the od of a cellPCA fit to clean data simulated from the fit, as in the reference code (this needs a second fit); "chisq" uses $\sqrt{\chi^2_{p, 0.99}}$ instead.

New observations, with missing or outlying cells, are projected by the robust regression of their observed cells on the loadings, also in C++.

Value

An object of class specproc_cellpca (inheriting from specproc_robpca), with the components described in robpca() and macropca(), and:

• cell_weights: the cellwise weights (0 for missing cells).

• case_weights: the casewise weights.

• deviation: the standardized casewise total deviations.

• fitted: the fitted values $\hat{X}$.

• imputed: the data in which the outlying cells are moved toward the fit in proportion to their weights, and the missing cells are replaced by the fit, so that the projection of each observation on the subspace is its fitted value.

• sigma1, sigma2: the scales of the cellwise residuals and of the casewise total deviations; resid_scale: the scales of the final residuals.

• objective: the objective at each iteration; iterations and converged.

Author(s)

Christian L. Goueguel

References

• Centofanti, F., Hubert, M., Rousseeuw, P.J. (2026). Robust principal components by casewise and cellwise weighting. Technometrics. doi:10.1080/00401706.2026.2643216

• Hampel, F.R., Rousseeuw, P.J., Ronchetti, E. (1981). The change-of- variance curve and optimal redescending M-estimators. Journal of the American Statistical Association, 76(375):643-648.

• Hubert, M., Rousseeuw, P.J., Van den Bossche, W. (2019). MacroPCA: an all-in-one PCA method allowing for missing values as well as cellwise and rowwise outliers. Technometrics, 61(4):459-473.

See Also

macropca(), robpca(), step_cellpca(), plot_outlier_map(), plot_cell_map()

Examples

data(forageLIBS)
# the K I resonance lines
wl <- suppressWarnings(as.numeric(names(forageLIBS)))
spectra <- forageLIBS[which(wl > 760 & wl < 780)]
set.seed(1)
fit <- cellpca(spectra, k = 2, od_cutoff = "chisq")
fit
plot_outlier_map(fit, shade = TRUE, relative = TRUE, log = TRUE, colour_by = "distance")
plot_cell_map(fit, order = "od")

Description

Function to perform mean-centering or median-centering on a numeric matrix or data frame.

Usage

center(x, method = "mean", drop.na = FALSE)

Arguments

Details

Mean-centering calculates the mean of each column and subtracts this from the column. Median-centering is very similar to mean-centering except that the reference point is the median of each column rather than the mean.

Value

A numeric matrix with the same dimensions as the input, with columns centered according to the specified method. The column centers are stored in the "center" attribute.

Author(s)

Christian L. Goueguel

Examples

data(forageLIBS)
spectra <- forageLIBS[-(1:14)]  # the spectral channels
center(spectra)[1:3, 1:4]
center(spectra, method = "median")[1:3, 1:4]

Description

Estimates the elemental composition of a sample from the intensities of its emission lines, without calibration standards, by the calibration-free LIBS (CF-LIBS) method of Ciucci et al. (1999).

Usage

cf_libs(
  lines,
  electron_density = NULL,
  temperature = NULL,
  method = NULL,
  partition = NULL,
  ionization_energy = NULL,
  reference = NULL,
  units = "energy",
  timeout = 120
)

Arguments

Details

In a plasma in local thermodynamic equilibrium (LTE), optically thin and of the same composition as the sample, the intensity of a line of species $s$ satisfies

$\ln\frac{I\lambda}{g_k A_{ki}} = -\frac{E_k}{k_B T} + \ln\frac{F C_s}{U_s(T)}$

where $C_s$ is the concentration of the species, $U_s(T)$ its partition function and $F$ an experimental factor common to all lines. The Boltzmann plots of all species are therefore parallel lines:

1. The temperature is estimated from their common slope, by a least-squares fit with one intercept per plot (unless temperature is given).

2. The relative number density of each species follows from the intercept $q_s$ of its plot, $F C_s = U_s(T) e^{q_s}$.

3. The density of an ionization stage that is not observed is computed with the Saha equation, from the electron density (see electron_density()) and the ionization energy.

4. The factor $F$ is eliminated by closure (the concentrations of the elements sum to 1), or, with reference, by the known mass fraction of one element (internal reference), which is useful when some elements of the sample are not measured.

Two methods place the lines on the plots:

• "boltzmann" (Ciucci et al., 1999): one Boltzmann plot per species, with its own intercept. The temperature rests on the range of upper-level energies within each species, often narrow (1 to 2 eV), so it can be imprecise.

• "saha-boltzmann" (the default when electron_density is given): one Saha-Boltzmann plot per element, on which the lines of the ion are placed with the Saha equation, as in saha_boltzmann(). The energy range is extended by the ionization energy, which gives a much more precise temperature, and the ionization balance of each element follows from the Saha equation. It needs the ionization energies of all elements.

At least one plot needs two lines or more of different upper-level energies.

When the partition functions are computed from energy levels, the levels of the neutral atoms above their ionization energy, lowered by the Debye-Hückel correction when electron_density is given, are left out (see partition_function()). This needs the ionization energies of all elements observed as neutral atoms.

The mass fractions are computed from the atomic fractions with the standard atomic weights. The result is only as good as its assumptions: all major elements must be measured (with closure), the lines must be optically thin (see self_absorption()), the intensities corrected for the spectral response of the instrument, and the atomic data accurate (prefer lines of accuracy B or better in nist_lines()).

Value

An object of class specproc_cflibs, a list with

• composition: a tibble with, for each element, the atomic_fraction and mass_fraction, and the stages observed (a stage computed by the Saha equation is noted "Saha");

• species: a tibble with, for each species, the number of lines, the intercept of its own Boltzmann plot, its partition function, its relative number density and whether it was observed;

• temperature, temperature_se (K) and electron_density;

• points: the points of the plots (for the Saha-Boltzmann method, the abscissa of ions includes the ionization energy and the ordinate the Saha correction), intercept, the intercept of each plot, and fit, the lm fit;

• method and the number of iterations of the Saha-Boltzmann fit.

Use plot_boltzmann() to draw the Boltzmann plots.

Author(s)

Christian L. Goueguel

References

• Ciucci, A., Corsi, M., Palleschi, V., Rastelli, S., Salvetti, A., Tognoni, E. (1999). New procedure for quantitative elemental analysis by laser-induced plasma spectroscopy. Applied Spectroscopy, 53(8):960-964.

• Tognoni, E., Cristoforetti, G., Legnaioli, S., Palleschi, V. (2010). Calibration-free laser-induced breakdown spectroscopy: state of the art. Spectrochimica Acta Part B, 65(1):1-14.

See Also

saha_boltzmann(), electron_density(), nist_lines(), nist_levels(), plot_boltzmann()

Examples

data(forageLIBS)
mean_spectrum <- colMeans(forageLIBS[-(1:14)])
# calcium and magnesium lines, with their atomic data from the NIST database
atomic <- data.frame(
  species = c("Ca I", "Ca I", "Ca I", "Ca II", "Ca II", "Mg I", "Mg I", "Mg I", "Mg II", "Mg II"),
  wavelength = c(428.301, 430.253, 445.478, 315.887, 317.933,
                 516.732, 517.268, 518.360, 279.078, 279.800),
  Aki = c(4.34e7, 1.36e8, 8.70e7, 3.10e8, 3.60e8, 1.13e7, 3.37e7, 5.61e7, 4.01e8, 4.79e8),
  gk = c(5, 5, 7, 4, 6, 3, 3, 3, 4, 6),
  Ek = c(4.780, 4.780, 4.681, 7.047, 7.050, 5.108, 5.108, 5.108, 8.864, 8.864)
)
# partition functions, interpolated from a table (see partition_function())
partition <- function(species, temperature) {
  grid <- c(6000, 8000, 10000, 12000)
  table <- rbind(`Ca I` = c(1.407, 2.401, 4.499, 8.356), `Ca II` = c(2.389, 2.917, 3.560, 4.262),
                 `Mg I` = c(1.049, 1.207, 1.597, 2.463), `Mg II` = c(2.001, 2.010, 2.036, 2.087))
  vapply(species, function(s) exp(stats::approx(grid, log(table[s, ]), xout = temperature)$y),
         numeric(1))
}
fit <- cf_libs(line_intensities(mean_spectrum, atomic, baseline = TRUE),
               electron_density = 1.9e17, partition = partition,
               ionization_energy = c(Ca = 6.113, Mg = 7.646))
fit
plot_boltzmann(fit)

Description

Computes how much each variable (wavelength) contributes to the Q residual or to Hotelling's $T^2$ of each sample, keeping the sign of the deviation, to find the variables that make a sample an outlier. Relative contributions, computed against reference samples, show what differs between a sample and the reference ones.

Usage

contributions(
  model,
  k = NULL,
  statistic = c("q", "t2"),
  data = NULL,
  samples = NULL,
  reference = NULL
)

Arguments

Details

For a sample $x$ (centered and scaled as in the model), with scores $t = xP$ on the first $k$ components (loadings $P$, variances of the scores $\lambda$):

• the Q contributions are the residuals $e = x - tP^T$, whose sum of squares is Q (the squared orthogonal distance of robust fits);

• the $T^2$ contributions are $t \Lambda^{-1/2} P^T$, the scaled scores projected back onto the variables, whose sum of squares is $T^2$ (the squared score distance of robust fits). For rospca() fits, whose loadings are not exactly orthogonal, the sum of squares is close to, but not exactly, $T^2$.

Q contributions show which variables the model does not describe for a sample: large contributions grouped on a few emission lines point to a systematic deviation (a contamination, a matrix effect, a saturated or shifted line), small ones spread over all channels to random noise. $T^2$ contributions show which variables place a sample far from the center in the score space.

Relative contributions. Normally, contributions are relative to the model (its center, for $T^2$). With reference, the mean contribution of the reference samples is subtracted from the contribution of each sample, which shows what differs between them: for example, whether two samples have large Q residuals for the same reason, or what moves a sample from the regular samples to where it lies in the score space. reference = "regular" uses all the regular samples: of the robust fit, or at the highest confidence level of q_residuals() for a stats::prcomp() fit.

The contributions are in the units of the preprocessed data of the model: centered (and scaled, if the model scaled the variables).

Value

A tibble of class specproc_contributions with one row per sample: sample (the row name, or number) and one column per variable. The attribute total holds the statistic of each sample (Q or $T^2$, before any reference is subtracted). Draw it with plot_contributions().

Author(s)

Christian L. Goueguel

References

• Westerhuis, J.A., Gurden, S.P., Smilde, A.K. (2000). Generalized contribution plots in multivariate statistical process monitoring. Chemometrics and Intelligent Laboratory Systems, 51(1):95-114.

See Also

plot_contributions(), q_residuals(), plot_outlier_map()

Examples

spectra_id <- names(forageLIBS)[1:2]
minerals <- names(forageLIBS)[3:14]
spectra <- forageLIBS |> dplyr::select(-dplyr::all_of(c(spectra_id, minerals)))
set.seed(1)
fit <- robpca(spectra, k = 3)

# Q contributions of two outlying samples, relative to the regular ones
q <- contributions(fit, data = spectra, samples = c(49, 127), reference = "regular")
q[, 1:5]

Description

Corrects the intensities of emission lines for self-absorption by the internal reference method of Sun and Yu (2009): for each species, the self-absorption coefficient of every line is measured against a line of the same species that is little self-absorbed.

Usage

correct_self_absorption(
  lines,
  temperature = NULL,
  electron_density = NULL,
  ionization_energy = NULL,
  reference = NULL,
  units = "energy",
  tol = 1e-08,
  timeout = 120
)

Arguments

Details

In an optically thin plasma in LTE, the intensity ratio of two lines of the same species depends only on their atomic data and on the temperature. With a reference line $r$ assumed optically thin, the self-absorption coefficient of a line $j$ is

$SA_j = \frac{I_j}{I_r} \frac{g_r A_r \lambda_j}{g_j A_j \lambda_r} \exp\left(\frac{E_j - E_r}{k_B T}\right)$

(for intensities in energy units; without the wavelengths for photon counts), where $g$, $A$ and $E$ are the statistical weight, transition probability and energy of the upper level. The corrected intensity is $I_j / SA_j$.

Self-absorption grows with the optical depth of the line, roughly proportional to $g_k A_{ki} \lambda^4 \exp(-E_i / k_B T)$, where $E_i$ is the energy of the lower level: resonance lines ($E_i = 0$) and strong lines are the most absorbed. By default, the reference of each species is the line of smallest optical depth, which needs the lower-level energies (column Ei, as returned by nist_lines()); otherwise, give the reference lines with reference. The reference line must be well measured: its noise propagates to every coefficient.

The corrected lines of a species lie, by construction, on the Boltzmann plot of its reference line at the temperature used, so they cannot give the temperature by themselves. Give temperature, measured on optically thin lines (for example with saha_boltzmann() or cf_libs()), or give electron_density: the temperature is then the one at which the corrected lines of the neutral atom and of the ion of each element lie on a common Saha-Boltzmann plot (as in cf_libs()), which needs lines of both stages for at least one element. This is the iterative scheme of Sun and Yu (2009), solved directly for the temperature.

Coefficients above 1 mean that the line is stronger than expected from the reference, because of noise, errors of the atomic data, or a self-absorbed reference; they are kept, with a warning.

Value

lines, as a tibble, with the corrected intensity, the measured_intensity, the self-absorption coefficient SA and a logical reference, ready for boltzmann() or cf_libs(). The temperature used is in the attribute temperature.

Author(s)

Christian L. Goueguel

References

• Sun, L., Yu, H. (2009). Correction of self-absorption effect in calibration-free laser-induced breakdown spectroscopy by an internal reference method. Talanta, 79(2):388-395.

See Also

self_absorption() for the coefficient from line widths, cf_libs(), line_intensities()

Examples

data(forageLIBS)
mean_spectrum <- colMeans(forageLIBS[-(1:14)])
# potassium lines: the resonance doublet at 766.49 and 769.90 nm ends on
# the ground state (Ei = 0) and is strongly reabsorbed
k_lines <- data.frame(
  species = "K I", wavelength = c(404.414, 404.721, 691.108, 693.877, 766.490, 769.896),
  Aki = c(1.150e6, 1.070e6, 2.500e6, 4.956e6, 3.779e7, 3.734e7),
  gk = c(4, 2, 2, 2, 4, 2),
  Ek = c(3.065, 3.063, 3.403, 3.403, 1.617, 1.610),
  Ei = c(0, 0, 1.610, 1.617, 0, 0)
)
corrected <- correct_self_absorption(line_intensities(mean_spectrum, k_lines), temperature = 8000)
corrected[c("wavelength", "measured_intensity", "SA", "reference")]

Description

Computes various correlation coefficients between one or several response variables and each of the remaining variables in a given data frame or tibble. The available correlation methods are Pearson's product-moment correlation (parametric), Spearman's rank correlation, Kendall's tau correlation (non-parametric), Chatterjee's new correlation coefficient, and the biweight midcorrelation (a robust correlation measure).

Usage

correlation(
  x,
  var,
  method = "pearson",
  plot = FALSE,
  color = c("#1f4e79", "#c0392b"),
  interactive = FALSE,
  top = NULL,
  cluster = FALSE
)

Arguments

Details

The Pearson correlation coefficient measures the linear relationship between two continuous variables and is suitable when the data follows a bivariate normal distribution. The Spearman and Kendall correlations are non-parametric measures of monotonic association, making them suitable for non-linear relationships and when the data deviates from normality. The Chatterjee correlation coefficient $\xi_n(X, Y)$ measures how much the response var is a (possibly non-monotonic) function of each other variable; it lies between 0 and 1 (asymptotically) and is not symmetric. The biweight midcorrelation is a robust correlation measure that downweights the influence of outliers and is recommended when the data contains extreme values or deviates significantly from normality.

Missing values are handled pairwise: each coefficient uses the observations where both the response and the other variable are available. With several responses, each one uses its own observations, so a response with many missing values does not reduce the data of the others.

Value

• If plot = FALSE, a tibble with columns variable, .correlation and method, sorted by decreasing correlation. With several responses, it starts with a column outcome, and is sorted within each response.

• If plot = TRUE, a list containing the tibble (correlation) and a ggplot2 object (plot; a patchwork object with cluster = TRUE), and with cluster = TRUE the clustering (clustering, an stats::hclust() object).

• If plot = TRUE and interactive = TRUE, a plotly object.

Plots

When the variables are named by wavelengths (the channels of spectra), the plot is a correlation spectrum: the correlation at each wavelength, with the thresholds of significance at the 5% level (for the Pearson, Spearman and biweight coefficients, from the t distribution; not corrected for multiple testing). Otherwise, it is a chart of the correlation of each variable, sorted, colored by sign and labeled with its value. The interactive versions show the variable and its correlation on hover; the correlation spectrum uses WebGL, so it stays fast with thousands of channels.

With several responses, the plot is a heatmap: one row per response, in the order given, and one column per wavelength (each tile as wide as the spacing of the channels, with gaps between detectors) or per variable. The color scale is fixed, from -1 to 1 (0 to 1 for Chatterjee's coefficient), so that a color means the same correlation in every row and in every plot.

With cluster = TRUE, the responses are clustered on their correlation profiles, the correlations with all the variables (as a function of wavelength, for spectra): the distance between two responses is $1 - r$, where $r$ is the Pearson correlation between their profiles, and the clusters are merged by average linkage. Responses whose correlations rise and fall at the same wavelengths, such as elements with lines in the same regions or that vary together in the samples, are then adjacent, and the dendrogram on the left shows how similar they are. The static plot then needs the patchwork package; the interactive one is only reordered.

Author(s)

Christian L. Goueguel

References

• Chatterjee, S. (2021). A new coefficient of correlation. Journal of the American Statistical Association, 116(536):2009-2022.

• Wilcox, R. (2012). Introduction to robust estimation and hypothesis testing (3rd ed.). Academic Press. (ISBN 978-0123869838).

Examples

# LIBS spectra of forage samples and their mineral contents
data(forageLIBS)
spectra_id <- names(forageLIBS)[1:2]
minerals <- names(forageLIBS)[3:14]

# potassium against the other mineral contents
correlation(forageLIBS[minerals], K)
correlation(forageLIBS[minerals], "K", method = "bicor")
correlation(forageLIBS[minerals], K, plot = TRUE)$plot

spectra <- forageLIBS[setdiff(names(forageLIBS), spectra_id)]

# a correlation spectrum: potassium against every channel
k <- correlation(spectra[setdiff(names(spectra), setdiff(minerals, "K"))], K, plot = TRUE)
k$plot


# a heatmap: every mineral against every channel
all_minerals <- correlation(spectra, dplyr::all_of(minerals), plot = TRUE)
all_minerals$plot

# the same, with the minerals clustered by their correlation profiles
if (requireNamespace("patchwork", quietly = TRUE)) {
  correlation(spectra, dplyr::all_of(minerals), plot = TRUE, cluster = TRUE)$plot
}

Description

This function implements the Direct Orthogonalization (DO) algorithm, as proposed by Andersson (1999), to filter out variation from predictor variables, $\textbf{X}$, that is orthogonal to the response variable(s), $\textbf{Y}$.

Usage

direct_orthogonal(x, y, ncomp = 2, center = TRUE, scale = FALSE)

Arguments

Details

Contrary to the Orthogonal Signal Correction (OSC) algorithm, Wold et al. (1998), which uses inverse Partial Least Squares (PLS) regression to filter out the orthogonal signal, DO filters out the orthogonal signal directly by orthogonalization of the $\textbf{X}$ matrix:

$\textbf{Z} = \textbf{X} - \textbf{Y}(\textbf{Y}^T\textbf{Y})^{-1}\textbf{Y}^T\textbf{X}$

Principal Components Analysis (PCA) is performed on $\textbf{Z}$ to obtain the loadings $\textbf{P}$, and the corrected matrix is $\textbf{X}_{DO} = \textbf{X} - \textbf{XPP}^T$. Direct Orthogonalization is typically simpler and faster than OSC.

To correct new data, preprocess it with the returned center and scale vectors and apply $\textbf{X}_{new} - \textbf{X}_{new}\textbf{PP}^T$.

Value

An object of class specproc_direct_orthogonal (a list), which predict() applies to new spectra, with the following components:

• correction: The corrected matrix.

• loading: The loadings matrix $\textbf{P}$.

• score: The scores matrix $\textbf{XP}$.

• center, scale: The column centers and scales applied to x.

Author(s)

Christian L. Goueguel

References

• Andersson, C.A., (1999). Direct orthogonalization. Chemometrics Intell. Lab. Syst., 47(1):51-63

• Pierna, J.A.F., Massart, D.L., de Noord, O.E., Ricoux, P., (2001). Direct orthogonalization: some case studies. Chemometrics Intell. Lab. Syst., 55(1-2):101-108

• Wold, S., Antti, H., Lindgren, F., Ohman, J. (1998). Orthogonal signal correction of near-infrared spectra. Chemometrics Intell. Lab. Syst., 44(1):175-185.

See Also

predict.specproc_filter() to correct new spectra with the fitted filter, step_direct_orthogonal() to use it in a tidymodels recipe. Before specProc 0.4.0, nas() returned the same correction; it now computes the net analyte signal and figures of merit.

Examples

data(forageLIBS)
wl <- suppressWarnings(as.numeric(names(forageLIBS)))
spectra <- forageLIBS[which(wl > 760 & wl < 780)]  # the K I resonance lines
res <- direct_orthogonal(spectra, forageLIBS$K, ncomp = 2)
dim(res$correction)

Description

This function implements the Direct Orthogonal Signal Correction (DOSC) algorithm, as proposed by Westerhuis et al. (2001), to remove systematic variation from predictor variables, $\textbf{X}$, that is orthogonal to the response variable(s), $\textbf{Y}$.

Usage

direct_osc(x, y, ncomp = 10, center = TRUE, scale = FALSE, tol = 0.001)

Arguments

Details

Different from the Orthogonal Signal Correction (OSC) algorithm, Wold et al. (1998), the DOSC algorithm is non-iterative:

1. $\textbf{Y}$ is projected onto the column space of $\textbf{X}$: $\hat{\textbf{Y}} = \textbf{XX}^{+}\textbf{Y}$.

2. $\textbf{X}$ is orthogonalized with respect to $\hat{\textbf{Y}}$: $\textbf{Z} = \textbf{X} - \hat{\textbf{Y}}\hat{\textbf{Y}}^{+}\textbf{X}$.

3. PCA of $\textbf{Z}$ gives the orthogonal scores $\textbf{T}$.

4. The weights $\textbf{W} = \textbf{X}^{+}\textbf{T}$ express the scores as a linear combination of $\textbf{X}$, the loadings are $\textbf{P} = \textbf{X}^T\textbf{T}(\textbf{T}^T\textbf{T})^{-1}$, and the corrected matrix is $\textbf{X} - \textbf{TP}^T$.

New data are corrected with $\textbf{X}_{new} - \textbf{X}_{new}\textbf{WP}^T$ after applying the returned center and scale.

Value

An object of class specproc_direct_osc (a list), which predict() applies to new spectra, with the following components:

• correction: The corrected matrix.

• loading: The loadings matrix $\textbf{P}$.

• score: The scores matrix $\textbf{T}$.

• weight: The weights matrix $\textbf{W}$.

• center, scale: The column centers and scales applied to x.

Author(s)

Christian L. Goueguel

References

• Westerhuis, J.A., Jong, S.D., Smilde, A.K., (2001). Direct orthogonal signal correction. Chemometrics Intell. Lab. Syst., 56(1):13-25

• Wold, S., Antti, H., Lindgren, F., Ohman, J. (1998). Orthogonal signal correction of near-infrared spectra. Chemometrics Intell. Lab. Syst., 44(1):175-185.

See Also

predict.specproc_filter() to correct new spectra with the fitted filter, and step_direct_osc() to use it in a tidymodels recipe.

Examples

data(forageLIBS)
wl <- suppressWarnings(as.numeric(names(forageLIBS)))
spectra <- forageLIBS[which(wl > 760 & wl < 780)]  # the K I resonance lines
res <- direct_osc(spectra, forageLIBS$K, ncomp = 2)
dim(res$correction)

Description

This function computes the directional outlyingness of a numeric vector, as proposed by Rousseeuw et al. (2018), which is a measure of outlyingness of data point that takes the skewness of the underlying distribution into account.

Usage

directional_outlyingness(
  x,
  cutoff.quantile = 0.995,
  rmZeroes = FALSE,
  maxRatio = NULL,
  precScale = 1e-10
)

Arguments

Details

Directional outlyingness takes the potential skewness of the underlying distribution into account, by splitting the univariate dataset in two half samples around the median. And then apply one-step M-estimator with Huber $\rho$-function for scaling each part.

Value

A tibble with columns:

• data: The original numeric values.

• score: The calculated outlyingness score.

• flag: A logical vector indicating whether each value is a potential outlier or not.

Author(s)

Christian L. Goueguel

References

• Rousseeuw, P.J., Raymaekers, J., Hubert, M., (2018). A Measure of Directional Outlyingness With Applications to Image Data and Video. Journal of Computational and Graphical Statistics, 27(2):345–359.

Examples

data(forageLIBS)
# iron contents (mg/kg): a few samples are far above the others
fe <- forageLIBS$Fe[!is.na(forageLIBS$Fe)]
head(directional_outlyingness(fe))

Description

Computes the distance of each sample to a principal component analysis (PCA) model in the space of the variables (DModX): the residual standard deviation of the sample, normalized by that of the calibration samples, with its limits at one or more confidence levels, and Hotelling's $T^2$ for the influence plot.

Usage

dmodx(
  model,
  k,
  newdata = NULL,
  conf_level = 0.975,
  normalized = TRUE,
  df = "effective",
  t2_method = "f",
  center = TRUE,
  scale = FALSE
)

Arguments

Details

With the residuals $e_{ij}$ of the $K$ variables after k components, for $N$ calibration samples (and $A_0 = 1$ if the data are centered),

$s_i = \sqrt{\frac{\sum_j e_{ij}^2}{K - k}}\,\sqrt{\frac{N}{N - k - A_0}}, \qquad s_0 = \sqrt{\frac{\sum_i \sum_j e_{ij}^2}{(N - k - A_0)(K - k)}}$

(the correction factor $\sqrt{N / (N - k - A_0)}$ applies to the calibration samples only). DModX is $s_i / s_0$ (normalized, default) or $s_i$. Its squared ratio follows an F distribution with $\nu$ and $(N - k - A_0)\nu$ degrees of freedom, so the limit of the normalized DModX is $\sqrt{F_{1-\alpha}}$ (times $s_0$ for the absolute DModX).

The nominal degrees of freedom are $\nu = K - k$ (df = "nominal"), as if the residuals of the variables were independent. For spectra, with far more (correlated) channels than samples, this gives a limit close to 1 that flags a large share of ordinary samples. With df = "effective" (default), $\nu$ is the effective number of residual dimensions, $\theta_1^2 / \theta_2$ from the eigenvalues $\lambda_a$ of the components left out ($\theta_j = \sum_{a > k} \lambda_a^j$, as in Box's approximation of Q), which is $K - k$ when the residuals are independent with equal variances. DModX is then consistent with the Q residuals of q_residuals(), of which it is a scaled square root.

Value

A tibble of class specproc_influence, with one row per sample: sample, t2 and its limits (as for q_residuals()), dmodx, its limits at each confidence level (dmodx_limit_97.5, ...), and outlier, the type of the sample at the highest confidence level. Draw it with plot_influence().

Author(s)

Christian L. Goueguel

References

• Eriksson, L., Johansson, E., Kettaneh-Wold, N., Trygg, J., Wikström, C., Wold, S. (2006). Multi- and Megavariate Data Analysis, Part I, 2nd ed. Umetrics Academy, Umeå.

• Box, G.E.P. (1954). Some theorems on quadratic forms applied in the study of analysis of variance problems, I. Annals of Mathematical Statistics, 25(2):290-302.

See Also

q_residuals(), plot_influence(), hotelling_t2()

Examples

if (rlang::is_installed("HotellingEllipse", version = "1.3.0")) {
  data(forageLIBS)
  # the 380-430 nm window (Ca II H and K lines)
  wl <- suppressWarnings(as.numeric(names(forageLIBS)))
  spectra <- forageLIBS[which(wl > 380 & wl < 430)]
  pca <- stats::prcomp(spectra)
  d <- dmodx(pca, k = 3)
  d[d$outlier != "regular", ]
  plot_influence(d, label = forageLIBS$Measurement)
}

Description

Estimates the electron density of a plasma from the Stark width of an emission line: from the Stark broadening parameters of the line (STARK-B, or a reference width), or from the width of the hydrogen H$\alpha$ line at 656.28 nm.

Usage

electron_density(
  width,
  method = "stark",
  stark = NULL,
  wavelength = NULL,
  temperature = NULL,
  reference_width = NULL,
  reference_density = 1e+17,
  tolerance = 0.5
)

Arguments

Details

width is the Stark full width at half maximum, in nm, with the instrumental and Doppler broadening removed. The line must be optically thin: self-absorption broadens lines and overestimates the density (see self_absorption()).

• For isolated lines (method = "stark"), the Stark profile is Lorentzian: use the Lorentzian width wL of a Voigt fit with peak_fit() or multipeak_fit(), whose Gaussian part absorbs the instrumental and Doppler broadening.

• The Stark profile of H$\alpha$ (method = "halpha") is not Lorentzian: use the FWHM of the whole line, for example voigt_fwhm() of the fitted widths, corrected for the instrumental width when it is not negligible (for Gaussian broadening, approximately $\sqrt{w^2 - w_{inst}^2}$).

• method = "stark": in the impact approximation the Stark width is proportional to the electron density,

  $N_e = N_{ref} \frac{w}{w_{ref}(T)}$

  where $w_{ref}(T)$ is the electron-impact width at density $N_{ref}$. It is taken from stark (a table returned by starkb_lines(), interpolated at temperature with stark_width()), or given directly as reference_width (nm) at reference_density. Tables from other sources can be built with stark_table() or read from saved STARK-B files with read_starkb(). Ion broadening is neglected, which is usually justified for lines of ions and at LIBS densities.

• method = "halpha": the computer-simulation fit of Gigosos, González and Cardeñoso (2003) for the H$\alpha$ line,

  $w = 1.098\,(N_e / 10^{17})^{0.67823}\ \mathrm{nm},$

  which includes the ion dynamics and depends only weakly on temperature.

Value

A numeric vector of electron densities, in cm$^{-3}$.

References

• Gigosos, M.A., González, M.Á., Cardeñoso, V. (2003). Computer simulated Balmer-alpha, -beta and -gamma Stark line profiles for non-equilibrium plasmas diagnostics. Spectrochimica Acta Part B, 58(8):1489-1504.

• Konjević, N. (1999). Plasma broadening and shifting of non-hydrogenic spectral lines: present status and applications. Physics Reports, 316(6):339-401.

See Also

starkb_lines(), stark_width(), peak_fit(), voigt_fwhm()

Examples

# From the H-alpha line: a Stark FWHM of 0.5 nm
electron_density(0.5, method = "halpha")

# From a line with known Stark width 0.012 nm at 1e17 cm-3
electron_density(c(0.006, 0.012, 0.024), reference_width = 0.012)

Description

This function implements extended multiplicative signal correction (EMSC), as proposed by Martens and Stark (1991). EMSC extends msc() by modeling, in addition to the multiplicative and additive effects, a smooth polynomial baseline and, optionally, known interferent spectra.

Usage

emsc(
  x,
  xref = NULL,
  degree = 2,
  interferents = NULL,
  wavelength = NULL,
  robust = TRUE
)

Arguments

Details

Each spectrum $\textbf{x}_i$ is modeled as a linear combination of a reference spectrum $\textbf{m}$, polynomials of the wavelength $\lambda$ and interferent spectra $\textbf{k}_j$:

$\textbf{x}_i = b_i\textbf{m} + a_i + d_{i1}\lambda + d_{i2}\lambda^2 + \dots + \sum_j h_{ij}\textbf{k}_j + \textbf{e}_i$

The coefficients are estimated by least squares, and the corrected spectrum is

$\textbf{x}_i^{corr} = (\textbf{x}_i - a_i - d_{i1}\lambda - \dots - \sum_j h_{ij}\textbf{k}_j) / b_i$

With degree = 0 and no interferents, EMSC is identical to msc() with drop.offset = TRUE.

The wavelengths are rescaled to $[-1, 1]$ before the polynomials are formed, for numerical stability. Interferent spectra describe variation that should be removed, for example the spectrum of a contaminant or the dominant directions of repeated measurements of the same samples (see epo()).

The reference spectrum is estimated from x unless xref is given, so new spectra should be corrected with predict(), which uses the reference and design of the calibration data.

Value

An object of class specproc_emsc (a list), which predict() applies to new spectra, with the following components:

• correction: The corrected spectra.

• coefficients: A tibble with the estimated coefficients of each spectrum: slope ($b_i$), offset ($a_i$), poly1, ... ($d_{ik}$) and interferent1, ... ($h_{ij}$).

• reference: The reference spectrum.

• design: The design matrix (one column per model term).

Author(s)

Christian L. Goueguel

References

• Martens, H., Stark, E. (1991). Extended multiplicative signal correction and spectral interference subtraction: new preprocessing methods for near infrared spectroscopy. Journal of Pharmaceutical and Biomedical Analysis, 9(8):625-635.

• Afseth, N.K., Kohler, A. (2012). Extended multiplicative signal correction in vibrational spectroscopy, a tutorial. Chemometrics and Intelligent Laboratory Systems, 117:92-99.

See Also

msc(), step_emsc() to use EMSC in a tidymodels recipe.

Examples

data(forageLIBS)
spectra <- forageLIBS[-(1:14)]  # the spectral channels
fit <- emsc(spectra[1:300, ], degree = 2)
head(fit$coefficients)
# new spectra are corrected with the calibration reference
corrected <- predict(fit, spectra[301:368, ])
dim(corrected)

Description

This function implements the External Parameter Orthogonalization (EPO) algorithm as proposed by Roger et al. (2003). The EPO algorithm aims to remove interferences or clutter present in the spectral matrix, effectively separating the signal of interest from unwanted perturbations.

Usage

epo(x, ncomp = 2, clutter = NULL)

Arguments

Details

The EPO algorithm works by detecting the dominant directions (right singular vectors) of a clutter matrix $\textbf{D}$ that describe the external variations, e.g. the differences between spectra of the same samples measured under different temperatures or moisture levels. It then projects the spectral matrix onto the subspace orthogonal to these variations, effectively removing the unwanted perturbations and extracting the signal of interest.

Let $\textbf{X}$ be the spectral matrix. The EPO algorithm aims to split $\textbf{X}$ into:

$\textbf{X} = \textbf{XP} + \textbf{XQ} + \textbf{R}$

where $\textbf{P}$ and $\textbf{Q}$ are, respectively, the projection matrices of $\textbf{X}$ onto the useful and perturbation (clutter) subspaces. $\textbf{R}$ is the residual matrix. With $\textbf{V}$ the first ncomp right singular vectors of $\textbf{D}$, $\textbf{Q} = \textbf{VV}^T$ and the corrected matrix is $\textbf{X} - \textbf{XVV}^T$.

The singular value decomposition is computed in C++ (Eigen), and the $p \times p$ projection matrix is never formed explicitly.

Value

An object of class specproc_epo (a list), which predict() applies to new spectra, with four components:

• correction: The orthogonalized matrix, representing the signal of interest.

• clutter: The clutter part of x, $\textbf{XVV}^T$.

• loadings: The singular vectors $\textbf{V}$ used for the orthogonalization.

• singular_values: The corresponding singular values of the clutter matrix.

Author(s)

Christian L. Goueguel

References

• Roger, J.-M., Chauchard, F., Bellon-Maurel, V. (2003). EPO-PLS external parameter orthogonalization of PLS application to temperature-independent measurement of sugar content of intact fruits. Chemometrics and Intelligent Laboratory Systems, 66(2):191-204.

See Also

predict.specproc_filter() to correct new spectra with the fitted filter, and step_epo() to use it in a tidymodels recipe.

Examples

data(forageLIBS)
spectra <- forageLIBS[-(1:14)]  # the spectral channels
# the three samples measured twice
twice <- forageLIBS$Sample[duplicated(forageLIBS$Sample)]
first <- match(twice, forageLIBS$Sample)
second <- vapply(twice, function(s) max(which(forageLIBS$Sample == s)), integer(1))
# their differences describe the variation between repeated measurements
clutter <- as.matrix(spectra[second, ]) - as.matrix(spectra[first, ])
res <- epo(spectra, ncomp = 1, clutter = clutter)
dim(res$correction)