cellpca(): robust PCA by casewise and cellwise weighting (cellPCA,
Centofanti, Hubert and Rousseeuw). It handles outlying observations,
outlying cells and missing values by minimizing a single objective, with
weights between 0 and 1 for every cell and observation. The iteratively
reweighted least squares and the predictions of new data are computed in
C++; the results agree with the reference code of the authors to
numerical precision. predict(), plot_outlier_map()
(enhanced outlier map), plot_cell_map() and flagged_regions() accept
its fits, and step_cellpca() is its recipe step.macropca() and
cellpca() now follows the DDC algorithm of the cellWise package: 1-step
M estimators for the standardization and the residual scales, the robust
Gnanadesikan-Kettenring correlations for the neighbors (wrapped
correlations, found exactly, for more than 750 variables), robust slopes
refined by least squares, each cell included in its own prediction, and
the standalone variables. It flags the same cells as cellWise, and more
of the outlying cells with far fewer false flags than before.macropca() now follows the MacroPCA algorithm of Hubert, Rousseeuw and
Van den Bossche (2019), built on ROBPCA, as implemented in cellWise, whose
results it reproduces (with scale = FALSE): projection pursuit on the
data where only the observations with the fewest flagged cells are
cell-imputed, iterative subspace estimation on a fixed set of
observations, one reweighting step on the orthogonal distances,
concentration steps and the deterministic MCD within the subspace, and
distances from the data with only the missing cells imputed. New data are
analyzed as by MacroPCApredict (with the DDC model of the fit). Its
results change. tol is now the largest angle between successive
subspaces as a fraction of a right angle (default 0.005, as in the
paper), the cut-offs of the score and orthogonal distances are at the 99%
level as in the paper, the residual scales are 1-step M scales, and the
fit returns the observations set aside by DDC (flagged_rows).baseline_lsp() uses the improved modified polynomial fit (IModPoly) of
Zhao et al. (2007) by default (method = "imodpoly"). The modified
polynomial fit of Lieber and Mahadevan-Jansen (2003), still available with
method = "modpoly", is pulled up by strong lines, and its baseline could
then fall far below the background at the ends of the spectrum (for
example above 290 nm in the 240-300 nm window of the forageLIBS
spectra). The default max.iter is now 100.ggsci, prospectr and XICOR move from Imports to Suggests.
generalized_boxplot() uses the default ggplot2 palette when ggsci is
not installed, and correlation(method = "chatterjee") asks for XICOR.
y_gradient_glsw() computes its Savitzky-Golay derivatives with the
package's own filter (same results).mt is no longer suggested.step_*() function checks that recipes is installed.step_*() functions,
predict() methods, reject_shots(), pds(), robust_bcyj(), umad(),
the wavelength calibration helpers and the tuning parameters).rospca() example runs on a spectral window, and the plot_loadings()
example that queries the NIST database is run with \donttest and try().q_residuals(), and the
Q limits of plot_influence()) was wrong when its exponent h0 is close to
zero or negative, which happens for spectra (a few large eigenvalues of
the components left out, then a long tail of small ones): it fell below
the mean of Q, so that nearly every sample was flagged. For example, with
3 components on the forageLIBS spectra, the 99% limit was 1.8e8 instead
of 2.0e9, and 355 of the 368 spectra were flagged. The sign of h0 is now
kept, as its derivation requires, and the limit as h0 tends to zero is
used near zero.opls() has the arguments of the other orthogonalization methods:
opls(x, y, ncomp = NULL, center = TRUE, scale = FALSE, crossval = 7, permutation = 0). ncomp (the number of orthogonal components, NULL
for the automatic choice) replaces ncomp.ortho, and the logical center
and scale replace the character scale; for Pareto scaling, apply
pareto_scale() to x first (or step_pareto_scale() before
step_opls()), which gives the same filter. The former arguments still
work, with a deprecation warning. The permutation test is now skipped by
default (permutation = 0), as it refits the model as many times.projected_osc()'s ncomp is now the number of orthogonal components
removed, as in opls(), osc() and direct_osc(), instead of the number
of PLS components (one more): replace ncomp = k + 1 with ncomp = k.
The default is 4 (formerly 5, the same filter). step_projected_osc() is
unchanged.macropca() and robust_bcyj() (and step_macropca(),
step_robust_bcyj()) are implemented natively, after the published
algorithms, and no longer depend on cellWise, which moves to Suggests.
macropca() detects the deviating cells with DDC (Rousseeuw and Van
den Bossche, 2018), whose neighbor search and cell predictions are
computed in C++ by blocks of variables: on the 368 x 7152 forageLIBS
spectra, DDC takes about 5 s instead of 13 s, and macropca() about
12 s instead of 16 s. The results are close to, but not the same as,
those of cellWise::MacroPCA() (orthogonal distances correlated at
0.99 on forageLIBS, nearly the same subspace); the score and
orthogonal distance cut-offs are those of robpca(). In simulations
with known rowwise and cellwise outliers and missing values, it
estimates the PCA subspace more accurately than cellWise::MacroPCA(),
imputes the missing values as accurately, and flags far fewer regular
observations as outliers (cellWise compares the distances computed with
the outlying cells to a cut-off estimated with them imputed). Its arguments
scale, ndir, maxiter and tol replace the ... passed to
cellWise::MacroPCA(), the result has no fit element any more,
and imputed holds the data with the missing values imputed.robust_bcyj() fits the transformations of Raymaekers and Rousseeuw
(2021), with nearly the same parameters as cellWise::transfo(), and
leaves the variables that cannot be transformed unchanged (method
"none").boltzmann_plot() and saha_boltzmann_plot() are renamed boltzmann()
and saha_boltzmann(), so that they are not confused with
plot_boltzmann(), which draws their result. The old names still work,
with a deprecation warning.dmodx(), the nominal degrees of freedom (K - k) are now
selected with df = "nominal".soilLIBS data set is removed, to keep the package data under the
5 MB CRAN guideline: data(soilLIBS) no longer works. The examples use
forageLIBS instead; the example of reject_shots(), which relied on
the soil replicates, is removed. The raw data and the script that built it remain in
data-raw/ of the source repository.som() fits a self-organizing map of spectra (batch SOM, in C++), which
maps similar spectra onto the same or neighboring units of a hexagonal or
rectangular grid and follows nonlinear relations (matrix effects,
self-absorption, plasma changes) that PCA can miss. It starts from the
plane of the first two principal components (deterministic training),
chooses the grid from the number of spectra, and reports the quantization
and topographic errors. robust = TRUE down-weights outlying spectra by
Huber weights of their quantization errors. predict() places new
spectra on their best-matching unit and flags as novel those beyond a
robust cut-off of the training quantization errors. plot_som() draws
the counts, the U-matrix, the quantization errors, the component planes
(one map per wavelength or line), the spectra on the map and the
prototype spectra, and som_stability() measures how consistently the
spectra keep their neighbors over resampled fits. On the 368 x 7152
forageLIBS spectra, a map takes about 2 s.plot_spectra() has an offset argument, to shift successive spectra
vertically and horizontally (offset = c(horizontal, vertical), or a
single number for a vertical offset), for example to stack them.plot_embedding() recognizes UMAP maps (axes named UMAP1, UMAP2, ...,
as by embed::step_umap()): it leaves out the lines through the origin,
whose position means nothing on such a map, and ignores hotelling
("all" or "group"), with a warning, since its T-squared limits
assume linear scores. On other maps, the lines through the origin are light grey
instead of black.correlation() takes several responses (var = c(K, Ca), or
dplyr::all_of(minerals)), each with its own observations, so that a
response with many missing values does not reduce the data of the others.
The result gains an outcome column, and plot = TRUE draws a heatmap of
the responses against the wavelengths (or the other variables), on a
fixed scale from -1 to 1. Pearson and Spearman correlations of spectra are
computed in one step: the 12 minerals of forageLIBS against its 7152
channels take under a second.
With cluster = TRUE, the rows are ordered by a hierarchical clustering
of the responses on their correlation profiles (1 - r, average linkage),
so that responses whose correlations rise and fall at the same
wavelengths are adjacent, with the dendrogram on the left.plot_calibration()'s "Linear calibration" above
its R², LOD and LOQ. This applies to every plot with a title, including
the interactive (plotly) ones and plot_cell_map(). plot_fit() now
shows its title as a title rather than a subtitle.forageLIBS spectra and mineral contents instead
of simulated data, wherever the function applies to them (some taken
from the vignettes). Examples of pure helper functions (line profiles,
plasma criteria, tuning parameters) keep simple inputs, and the examples
of pds() and step_reject_shots(), for which forageLIBS has no
suitable data, are removed. plot_influence() gains an example.dplyr::, so that it runs
without dplyr attached.forageLIBS spectra: fitting emission lines
(vignette("peak-fitting"), which replaces "line-fitting"),
preprocessing with recipe steps (vignette("preprocessing")),
calibration curves and figures of merit (vignette("calibration")) and
plasma diagnostics (vignette("plasma-diagnostics")). The
orthogonalization vignette is folded into the preprocessing one.reject_shots(), msc() and
emsc() with robust = TRUE, center(method = "median"), robpca()):
reject_shots() is about 20 times faster on 400 spectra of 7152 channels.q_residuals() and dmodx() compute the residuals of the calibration
samples from the scores of the components left out, without
reconstructing the data (about 100 times faster), and
line_intensities() measures all the spectra of a line at once (about 3
times faster). The results are unchanged.opls() is implemented natively and no longer needs the Bioconductor
package ropls. It reproduces ropls::opls() (scores, loadings, R2X, R2Y,
Q2, RMSEE, VIP and permutation p-values), but returns an object of class
specproc_opls instead of a list holding the ropls model: the model
element is removed, and the model has new elements correction (the
OPLS-filtered data), fitted, coefficients, vip, ortho_vip,
components (R2X, R2Y and Q2 of each component) and permutation.
x and y can be matrices or vectors as well as data frames. When
ncomp.ortho = NA and only the predictive component is significant, the
model has no orthogonal component (with a warning) instead of failing.
predict() filters new data, or predicts the response or the scores
(type), and crossval = 0 skips the cross-validation.specLIBS is now soilLIBS, and fourrage
is now forageLIBS. Replace data(specLIBS) with data(soilLIBS) and
data(fourrage) with data(forageLIBS).macropca() (and step_macropca()) now only centers the variables by
default, instead of also scaling them as cellWise::MacroPCA() does:
scaling spectra gives noise and continuum channels as much weight as
emission lines. Use macropca(x, scale = TRUE) for the previous
behavior.plot_cell_map() is redrawn with ggplot2 and patchwork instead of
cellWise::cellMap(), for spectra: the map fills the plot, with a
wavelength axis, cells combined into at most resolution blocks (within
detector segments), flagged cells colored by the sign and size of their
residual, a strip of the outlier type of each observation, rows
optionally sorted by orthogonal distance (order = "od"), and a top
panel of the share of flagged cells at each wavelength (profile). The
arguments nrowsinblock, ncolumnsinblock and ... are replaced by
resolution. The profile shows the mean spectrum in grey (of the
imputed data, or of spectra) and labels the most flagged regions, with
the emission lines they match (lines), and order = "cluster" groups
the observations with similar flagged cells.flagged_regions() lists the wavelength regions where many observations
have cells flagged by macropca() (at least threshold of them), with
their extent, share, direction and matching emission lines.plot_influence() is drawn like plot_outlier_map(): points outlined in
black and filled by type (in the colors of the outlier map), dashed limits
at the highest confidence level (dotted, ... at the others, with a legend
only for several levels), a compact legend at the bottom, and labels on
the labels (default 3) most outlying samples only, instead of all the
flagged ones. It gains the options of plot_outlier_map(): relative
(distances divided by their limits), shade (the outlying regions),
log, colour_by = "distance" and point styles in ..., including a
size per sample. The arguments after label changed order: name log
and title. In both plots, labels near the right edge are placed on the
left of their point, so that they are not cut off.
step_opls() removes the orthogonal components of an opls() model in a
recipe, with num_comp tunable. It gives the same filtered data as
step_projected_osc(), and also offers Pareto scaling.
step_o2pls() removes the outcome-orthogonal components of an o2pls()
model in a recipe, estimated against one or several outcomes (recipe(K + Ca ~ ., ...)), with num_comp and joint_comp tunable for a single
outcome. Only the predictors are filtered. Its help page shows a workflow
predicting several outcomes and, since tune_grid() does not support
several outcomes, how to tune one workflow per outcome with a workflow set
(workflow_map("tune_grid", ...)), each filter still estimated against
all the outcomes.
plot_embedding() is drawn with a grey panel outside the outermost
ellipse (of T-squared, or of each group) and white
inside, with no grid, a fixed aspect_ratio (0.7) and lines through the
origin. The axis titles give the explained variance of the components of
PCA fits, and the T-squared ellipses of the groups (hotelling = "group") are filled with their color. size can be a variable (such as
a concentration) that sets the size of the points, and biplot = TRUE
draws the loadings of a PCA over the scores, with the biplot_top most
important variables (emission lines, for spectra) as labeled arrows.
With hotelling = "all", the T-squared of a PCA model is that of the
model (its eigenvalues, centered at 0): for a robust fit, the ellipses are
no longer inflated or tilted by the outlying samples, and the flagged
samples are the leverage points of plot_outlier_map(). Without groups,
ellipse = TRUE with hotelling = "all" draws only the T-squared
ellipse, with a warning, instead of two ellipses of all the samples. The
samples beyond the T-squared limit are counted in the subtitle, circled
in red only with flag = TRUE, and labeled only with label (TRUE, or
a column or vector of labels), independently of each other.
plot_outlier_map() takes a size per sample (size), and
colour_by = "distance" colors its points by their reduced distance from
the origin (the larger of SD and OD over their cut-offs), on a rainbow
from dark red to blue; with relative = TRUE, where both cut-offs are at
1, yellow marks the cut-offs, between the regular observations (dark red
to orange) and the outlying ones (green to blue).
The legends of plot_embedding() and plot_outlier_map() are more
compact: smaller text and keys, three sizes, the T-squared levels written
on their ellipses instead of in a legend, and no key for the kind of
samples without new samples.
contributions() computes the contributions of each variable to the Q
residual or to Hotelling's T-squared of samples (whose squares add up to
Q and T-squared), for a prcomp()
fit or a robust fit (where they decompose the orthogonal and score
distances), and relative contributions against reference samples
(reference, or "regular" for all regular samples).
plot_contributions() draws them against wavelength with the largest
labeled and matched to emission lines.
plot_loadings() plots the loadings of a PCA (prcomp(), robpca(),
rospca() or macropca()) against wavelength, one panel per component,
and labels the wavelengths that contribute most to each component, with
the emission lines they match in a line list from libs_lines(). The
mean spectrum can be drawn behind the loadings, type = "contribution"
plots the variance of each wavelength explained by the components, and
interactive = TRUE makes a plotly figure. loading_peaks() returns the
labeled peaks as a table, with their candidate lines and a flag for
derivative shapes (line shifts).
plot_outlier_map() gains relative, to plot the reduced distances
(divided by their cut-offs) so that maps of different models share one
scale, shade, to shade the good leverage, orthogonal outlier and bad
leverage regions in greys of increasing darkness and name them, and
log, for logarithmic axes. The points are outlined in black, and
... passes styling such as alpha or size to geom_point().
dmodx() computes the distance to the model in the space of the
variables (DModX) of a PCA model, normalized or absolute,
for the calibration samples or new ones. Its limits use the effective
number of residual dimensions by default (df = "effective"), because
the nominal degrees of freedom (df = "nominal") flag a large share of
ordinary samples when spectra have far more channels than samples.
plot_influence() draws DModX or Q against T-squared.
hotelling_t2(), q_residuals(), dmodx(), plot_influence() and the
T-squared ellipses of plot_embedding() take any confidence level(s)
(conf_level, 0.975 by default; in plot_embedding(), one conf_level
sets both its confidence and T-squared ellipses),
and the Beta distribution of T-squared for the samples of the model
(method = "beta"). hotelling_t2() and the T-squared ellipses of
plot_embedding() now use the HotellingEllipse package (1.3.0 or later).
wavelength_calibration() fits a correction of the wavelength axis from
reference lines (for example, from the NIST database), located to a
fraction of a channel by Gaussian interpolation, with one correction per
detector segment (constant, linear or quadratic) and automatic rejection
of edge, weak, saturated and outlying lines. apply_calibration()
relabels the wavelengths of spectra without changing their intensities,
predict() corrects any wavelength, and plot_wavelength_calibration()
draws the offsets and the fit.
savitzky_golay() also splits the axis where overlapping detectors make
the wavelengths step back.
The periodic table of line_finder() can be hidden (button Hide
table, or show_table = FALSE at start) to give the spectrum the whole
height of the window; the selected elements stay listed in the panel
header.
savitzky_golay() and step_savgol(): Savitzky-Golay smoothing and first
or second derivatives, keeping every channel (polynomial fits at the
edges) and filtering the segments between detector gaps separately. The
window, polynomial order and derivative are tunable (new dials parameter
savgol_derivative()).
normalize() gains the L1 ("l1"), L2 or vector ("l2") and maximum
("max") norms, and step_spectral_norm() normalizes each spectrum in a
recipe by its L1 norm, total area, L2 norm or maximum, with the method
tunable (new dials parameter spectral_norm_method()).
wavelet_features() and step_wavelet(): discrete wavelet transform of
spectra (Haar, Daubechies d4/d6/d8 and least asymmetric la8 wavelets) as
compressed features: the approximation at a level or all coefficients,
optionally the coefficients of largest variance in the training data. The
level and number of coefficients are tunable (new dials parameter
wavelet_level()).
q_residuals() computes Hotelling's T-squared and the Q residual (SPE) of
each sample for a PCA model, for the calibration samples or new ones,
with their 95% and 99% limits (Jackson-Mudholkar or Box for Q), and
classifies the samples as regular, extreme, residual or both.
plot_influence() draws Q against T-squared with the limits.
calibration_curve() gains prediction intervals: plot_calibration()
draws the confidence band, the prediction band or both (interval), and
shows new samples (newdata) at their predicted concentrations with
their intervals; predict(type = "signal") gives the expected signal at
given concentrations with a confidence or prediction interval; the
coefficients have confidence intervals.
The plots of correlation() are redesigned. For spectra (variables named
by wavelength), the plot is a correlation spectrum with the 5%
significance thresholds; otherwise, a sorted chart colored by sign and
labeled with the values, with top to show only the strongest
correlations. The interactive versions are built directly with plotly
(WebGL for spectra) instead of converted from ggplot2. color now takes
one color or two (positive and negative correlations).
plot_embedding() draws two-dimensional embeddings (UMAP from
embed::step_umap(), principal components from recipes::step_pca(),
prcomp() or robpca()), colored by a variable, with optional
confidence ellipses per group (classical or robust, normal or Hotelling
quantiles) from the ConfidenceEllipse package, and Hotelling's T-squared
95% and 99% ellipses with the outlying samples labeled (hotelling, for
all samples or within groups, on k components). The preprocessing
vignette uses it to compare PCA and UMAP maps of the samples.
hotelling_t2() computes Hotelling's T-squared of each sample on k
components of an embedding, with its 95% and 99% limits, for all samples
or within groups.
macropca() failed when k was not given: cellWise::MacroPCA() with
k = 0 only reports the explained variance. k is now chosen from the
new arguments var_explained (0.8) and kmax (10), as in robpca(), and
MacroPCA's message and scree plot are no longer shown.cf_libs(): calibration-free LIBS quantification (Ciucci et al., 1999).
The temperature is estimated from the common slope of parallel Boltzmann
plots (one per species) or, when the electron density is known,
Saha-Boltzmann plots (one per element, much more precise). The densities
follow from the intercepts and the partition functions, and the
ionization stages not observed from the Saha equation. The composition
(atomic and mass fractions) is normalized by closure or by the known
concentration of one element. plot_boltzmann() draws the plots of the
fit. The plasma diagnostics vignette applies it to the forage samples and
compares it with their laboratory values.nist_levels() retrieves energy levels from the NIST Atomic Spectra
Database, and partition_function() computes partition functions from
them, optionally truncated at a maximum energy (such as the lowered
ionization energy).reject_shots() flags the outlying laser shots of each sample (total
intensity, correlation with or distance to the median spectrum, by robust
z-scores), before they are averaged with average().line_intensities() measures the area (or height, or Voigt-fitted area)
of emission lines in spectra, searching each peak near its tabulated
wavelength, with a signal-to-noise ratio and saturation check. It keeps
the columns of a table of lines, so its result feeds boltzmann_plot(),
cf_libs() and calibration_curve(). step_line_intensities() does the
same in a recipe.correct_self_absorption() corrects line intensities for
self-absorption by the internal reference method of Sun and Yu (2009),
with the temperature given or estimated from a Saha-Boltzmann plot.calibration_curve() fits univariate (linear or quadratic, optionally
weighted) calibration curves, with the sensitivity, LOD and LOQ, Mandel's
and lack-of-fit tests of linearity, inverse prediction of concentrations
with confidence intervals (predict()) and plot_calibration().step_reject_shots() removes the outlying shots from
the training data (skipped on new data), and step_line_ratio()
normalizes each spectrum to the area or height of a reference line
(internal standard).line_finder(): a Shiny app to identify the emission lines of LIBS
spectra. Elements are selected on a periodic table, and the strongest lines
of their ionization stages, from the NIST Atomic Spectra Database, are
overlaid on a spectrum in an interactive plotly graph (one color per
stage), with controls for the temperature, wavelength range, number of
lines, wavelength shift and spectrum (mean, group mean or single spectrum).
Fetched lines can be saved and reloaded for offline use.libs_lines() lists the lines of selected species expected to be the
strongest in a plasma at a given temperature (relative intensities in
LTE), and plot_lines() overlays them on a spectrum (plotly or ggplot2).This release adds plasma diagnostics for LIBS: electron density from Stark broadening, excitation temperature from Boltzmann and Saha-Boltzmann plots, and checks of LTE, self-absorption and detector saturation, with atomic data from the NIST Atomic Spectra Database and Stark parameters from STARK-B.
starkb_lines() retrieves Stark widths and shifts of the lines of an
atom or ion from the STARK-B database (Sahal-Bréchot, Dimitrijević and
Moreau) through its VAMDC service, on demand; read_starkb() reads
STARK-B data saved as XSAMS files, for offline and reproducible work;
stark_table() builds the same table from user-supplied widths or from
fitted temperature laws. stark_width() interpolates the width at a
given temperature and electron density, and scales multiplet data to a
line of the multiplet (lambda-squared rule).electron_density() estimates the electron density from the Stark
(Lorentzian) width of a line, with STARK-B data or a reference width, or
from the H-alpha line (Gigosos et al., 2003).boltzmann_plot() and saha_boltzmann_plot() estimate the excitation
temperature from line intensities and atomic data, and
plot_boltzmann() draws the plots.mcwhirter_criterion() checks the McWhirter criterion for LTE, and
self_absorption() computes self-absorption coefficients from line
widths (El Sherbini et al., 2005).saturation_summary() finds channels at the saturation limit of the
detector.nist_lines() and nist_ionization_energy() retrieve transition
probabilities, level energies, statistical weights and ionization
energies from the NIST Atomic Spectra Database, on demand.voigt_fwhm() computes the full width at half maximum of a Voigt profile
(Olivero and Longbothum, 1977).specLIBS and fourrage spectra.This release adds robust PCA, a proper net analyte signal with figures of merit, and rewrites the documentation around tidymodels.
nas() now computes the net analyte signal (Lorber, Faber and Kowalski,
1997; Faber, 1998) of an inverse PLS or PCR calibration model, with its
figures of merit: sensitivity, selectivity and, given the noise level,
analytical sensitivity, limits of detection and quantification, and
signal-to-noise ratios. A predict() method gives the NAS, selectivity and
predicted concentration of new samples. Previously, nas() returned
spectra filtered by direct orthogonalization, identical to
direct_orthogonal(), which remains available for that purpose.robpca() (ROBPCA, Hubert, Rousseeuw and Vanden Branden,
2005) and rospca() (robust sparse PCA, Hubert, Reynkens, Schmitt and
Verdonck, 2016), implemented from the published algorithms with C++
kernels for the Stahel-Donoho outlyingness, FAST-MCD and the grid-search
sparse PCA; and macropca(), a wrapper around cellWise::MacroPCA() for
cellwise outliers and missing values. All three return score and
orthogonal distances with cut-offs, an outlier classification, and a
predict() method for new observations.plot_outlier_map() draws the outlier map (score distance against
orthogonal distance) of any of the three, optionally with new
observations; plot_cell_map() draws the cell map of a MacroPCA fit.step_robpca(), step_rospca() and step_macropca(),
robust counterparts of recipes::step_pca() that can also add the score
and orthogonal distances as columns, and step_robust_bcyj(), a robust
counterpart of recipes::step_BoxCox() and recipes::step_YeoJohnson().rsample::nested_cv() and tune::tune_grid(), and
the orthogonalization vignette tunes every filter together with the PLS
model. Vignettes that need suggested packages stop with a message when
those are not installed.This release connects specProc to the tidymodels framework. Preprocessing and
orthogonalization can now be estimated on calibration data and applied to new
spectra, either directly with predict() or as recipe steps that are
re-estimated on every resample and tuned together with the model.
epo(), osc(), direct_orthogonal(), direct_osc(), projected_osc()
and o2pls() return classed objects with a predict() method, which
corrects new data with the centers, scales and components estimated on the
calibration data. Existing fields are unchanged.step_epo(), step_glsw(), step_osc(),
step_direct_orthogonal(), step_direct_osc(), step_projected_osc() and
step_y_gradient_glsw(), with tidy() and tunable() methods, and the
dials parameter glsw_alpha(). In a workflow, the filters are re-estimated
on every resample and can be tuned with the model. The GLSW steps store the
filter in factored form instead of a p x p matrix. recipes and dials are
suggested, not required.emsc(): extended multiplicative signal correction (Martens and Stark,
1991), with a polynomial baseline and optional interferent spectra, and a
predict() method for new spectra. With degree = 0 it equals msc().step_baseline() (arPLS, ALS
or LSP), step_snv(), step_msc(), step_emsc(), step_pareto_scale()
and step_poisson_scale(). MSC and EMSC estimate the reference spectrum,
and the scaling steps the column scales, on the training data only. The
baseline lambda or degree and the EMSC degree are tunable; the new
dials parameter baseline_lambda() covers lambda.fourrage: LIBS spectra of 365 forage samples with reference
contents of 12 elements, and a vignette based on it, "Removing unwanted
variation: a comparison of orthogonalization methods" (what EPO, GLSW, the
OSC family, DO/NAS, DOSC, POSC/OPLS and y-gradient GLSW remove, whether it
improves potassium predictions, and which methods are equivalent).This release fixes a large number of bugs so that every exported function now runs and returns numerically correct results. It adds a test suite (400+ expectations), and several results are checked against independent reference implementations.
Functions were renamed to a consistent snake_case scheme. The old names still
work but give a deprecation warning (see ?"specProc-deprecated"):
whittaker() -> baseline_als(), lorentzian() -> lorentzian_profile(),
pseudo_voigt() -> pseudo_voigt_profile(), peakfit() -> peak_fit(),
multipeakfit() -> multipeak_fit(), plotfit() -> plot_fit(),
plotSpec() -> plot_spectra(), outlierplot() -> plot_outliers(),
directOutlyingness() -> directional_outlyingness(),
iqrMethod() -> iqr_outliers(), robustBCYJ() -> robust_bcyj(),
rousseeuwCroux() -> rousseeuw_croux(), summaryStats() -> summary_stats(),
tukeyGH() -> tukey_gh(), yGradientglsw() -> y_gradient_glsw(),
pareto() -> pareto_scale().
gaussian() is renamed to gaussian_profile() without an alias, because
exporting gaussian() masked stats::gaussian(), which breaks
glm(family = gaussian) when specProc is loaded.
baseline_arpls(), baseline_lsp() and whittaker() return tibbles (they
previously returned plain lists of columns). They no longer replace negative
corrected values with 1.
average() with .group_by now returns the group labels as the first column.
biweight_location() and biweight_scale() use the unscaled MAD, following
the standard definition (Beers et al., 1990), so that biweight_scale()^2
equals biweight_midvariance(). biweight_scale() loses its unused tol
and max_iter arguments. loc = NULL now means "median of x".
direct_orthogonal() loses its unused max_iter and tol arguments, and
direct_osc() loses max_iter. tol now sets the pseudo-inverse tolerance.
epo() gains a clutter argument (the matrix describing the external
variation). Its default ncomp is now 2, as documented. It now removes the
dominant singular directions; before, it removed the smallest ones.
o2pls() has a new interface: o2pls(x, y, ncomp, nx, ny, center, scale).
pds() returns transfer_matrix (was transfert_matrix), alpha now blends
the local and global models, and progress is no longer printed.
projected_osc() always returns the corrected training data. New data are
preprocessed with the training centers and scales.
quantile_weight() uses whole-distribution quantiles, as in Brys et al.
(2006), and defaults to p = 0.125, q = 0.875.
summaryStats(robust = TRUE) drops the rsd column. It applied the 1.4826
consistency factor twice, and mad already holds the correct value.
opls() now requires the Bioconductor package ropls to be installed
(it moved to Suggests). It no longer prints output or draws plots.
R >= 4.1.0 is required.
baseline_arpls(), whittaker(), baseline_lsp() (data
frames), correlation(method = "bicor"), directOutlyingness(maxRatio = ),
direct_osc(), msc(), multipeakfit(), nas(), o2pls(), osc()
("wold", "fearn"), pds(), peakfit(), projected_osc() and
yGradientglsw().glsw() and yGradientglsw() squared the eigenvalues of the covariance
matrix when building the filter.msc() regressed each wavelength instead of each spectrum on the reference.generalized_boxplot() transformed the ranks instead of the data. This made
the g-and-h estimates meaningless. It also failed when the two tails had
different numbers of outliers.pseudo_voigt() counted the area and the offset twice.biweight_scale() used a wrong exponent in the denominator. It also ignored
reduced = TRUE.biweight_location() did not iterate.average() did not ignore missing values. It also read out of bounds for
factors with unused levels or missing group values.normalize(method = "background") rejected data frames for bkg.plotSpec() connected all spectra into a single line when id was not given.outlierplot() no longer changes the user's random number generator state.iqrMethod() no longer changes global options. Missing values no longer
cause errors.Three vignettes based on specLIBS: "Preprocessing LIBS spectra"
(baseline, normalization judged by replicate RSD and intraclass
correlation, robust screening of shots), "Fitting emission lines"
(profile choice, sources of uncertainty, identifiability) and
"Predicting soil clay content from LIBS spectra" (a compositional
log-ratio PLS model with nested, repeated cross-validation by sample, and
the optimism of common shortcuts).
plot_fit() draws the fitted profiles on a fine wavelength grid.
Functions taking a response now accept 1-d arrays, such as the output of
tapply().
New voigt_profile(): the exact Voigt profile, computed in C++ from the
Faddeeva function with Weideman's (1994) rational approximation (relative
error below 1e-10). peak_fit() and multipeak_fit() now use it for
profile = "voigt"; the pseudo-Voigt approximation is available as
profile = "pseudo_voigt".
New data set specLIBS: 400 LIBS spectra (50 soil samples x 8 locations,
7152 channels, 199-822 nm) with the clay, sand and silt content of each
sample.
generalized_boxplot() whiskers now end at the most extreme observations
within the fences, as in standard boxplots. The fences are returned as
lower_fence and upper_fence.
The README has been rewritten around a complete worked example on
specLIBS.
peakfit() and multipeakfit() estimate starting values from the data when
wL, wG or A are not given. multipeakfit() fits all lines
simultaneously and returns each line's contribution. A spectrum that fails
to fit produces a warning instead of stopping the whole batch.NEWS.md file to track changes to the package.