rsimpls(): robust PLS regression by the RSIMPLS algorithm of Hubert and
Vanden Branden (2003), the SIMPLS algorithm computed from a ROBPCA
covariance matrix of the spectra and the responses, followed by a robust
regression on the scores. Outlying spectra and wrong reference values have
little influence on the model. One or several responses; the response is
block-scaled before ROBPCA so that wrong reference values are detected
whatever its units. The fit gives the score, residual and orthogonal
distances of each observation, and the robust R² of the model.
predict() gives the predictions or scores of new spectra, with their
score and orthogonal distances, and plot_outlier_map() draws the
regression outlier map (vertical outliers, good and bad leverage points)
or, with map = "score", the score outlier map of the spectra.
step_rsimpls() is its recipe step, the robust counterpart of
recipes::step_pls(): it replaces the spectra by their robust PLS scores,
and optionally their score and orthogonal distances."rsimpls" engine of parsnip::pls() fits rsimpls() models in
tidymodels, with one or several outcomes; the other arguments of
rsimpls() are engine arguments (set_engine("rsimpls", kmax = , alpha = )). multi_predict() gives the predictions of the models with fewer
components from the same fit, so tune::tune_grid() fits one model per
resample for a grid of num_comp. predict() of an rsimpls() model
gains ncomp, for the models with 1 to kmax components, which the fit
now keeps (models).rpcr(): robust principal component regression (RPCR) of Hubert and
Verboven (2003), the robust PCA of the spectra (robpca()) followed by
the LTS regression of the response on the scores, or the MCD regression of
several responses. It has the outputs, predict() method and outlier maps
of rsimpls().robust_rmsecv(): the robust RMSECV of rsimpls() or rpcr() models
with 1 to kmax components, which leaves the outliers out of the
cross-validated error, and the robust component selection (RCS) criterion
of Engelen and Hubert (2005), to choose the number of components. The
folds are fitted in parallel when a future::plan() is set."lts" engine of parsnip::linear_reg() fits the reweighted LTS
regression (lts_fit(), from robustbase::ltsReg()); after
step_robpca() in a workflow, it gives robust PCR in tidymodels.rsimca(): robust SIMCA classification of Vanden Branden and Hubert
(2005), a robust PCA model of each class and the assignment to the closest
class from the scaled score and orthogonal distances, for whole spectra.
predict() gives the classes, or the distances to the classes and
whether a spectrum is outlying for all of them, and plot_coomans() draws
the distances to two classes against each other (Coomans plot), with the
class and classification boundaries. simca() is a new
parsnip model of SIMCA classification, with rsimca() as its "rsimca"
engine, to fit and tune it in tidymodels: num_comp, and the rule of
classification, gamma and squared, with the new dials parameters
simca_gamma() and simca_squared().robust_da(): robust linear and quadratic discriminant analysis of Hubert
and Van Driessen (2004), from the MCD estimates of the classes, for a few
variables such as line intensities or robust principal component scores.
predict() gives the classes or posterior probabilities. It is the
"mcd" engine of parsnip::discrim_linear() and
parsnip::discrim_quad().select_wavelengths(): selection of the informative wavelengths for PLS
regression, by the variable importance in projection (VIP), the
selectivity ratio (SR), backward elimination on either, or forward
interval PLS (iPLS). With robust = TRUE, the selection is made from an
rsimpls() model, without the regression outliers. The candidate
intervals of iPLS are evaluated in parallel when a future::plan() is
set. plot_wavelength_selection() shows the selected regions on the mean
spectrum with the importance of each variable.plot_shots(): views of the laser shots and of their rejection by
reject_shots(), to judge the rejection before averaging the shots: a map
of the samples by the shot number, colored by the robust z-scores in
classes up to the cutoff, the rejected shots outlined ("heatmap", the
default; acquisition orders the samples as measured, arrange = "extreme" puts the most extreme first, and smooth gives running medians
along the acquisition, which show the drift of the shots during a
session); the robust z-scores of the shots with the cutoff
("criteria"); the shots of a sample with the rejected ones in color
("spectra"); the intensity, shape and rejected share along the shot
number, which shows cleaning shots and the drift of the ablation
("order"); and the rejected shots per sample with the RSD of the
samples before and after the rejection ("samples"). plot = FALSE
returns the table of each view.forageShots: the single laser shots (8 per measurement, in firing order)
of 20 of the measurements of forageLIBS, in the windows 380-430 nm and
760-780 nm, to study the shot-to-shot variability.step_select_wavelengths(): the selection as a recipe step, repeated on
every resample; num_terms, num_intervals (new dials parameter
num_intervals()) and num_comp can be tuned.gap_derivative(): the gap-segment derivative of Norris and Williams
(1984), the first or second derivative computed from the means of
segments of segment channels separated by gap channels, so that the
smoothing and the distance of the difference are set separately (the
window of savitzky_golay() sets both). As in savitzky_golay(), every
channel is kept and the detectors of a multi-spectrometer system are
filtered separately. step_gap_derivative() is its recipe step;
derivative, gap and segment can be tuned.bin_spectra(): the means of groups of width adjacent channels, within
each detector, which reduce the number of channels and the noise at the
cost of resolution. step_bin_spectra() is its recipe step; width can
be tuned. With step_spectral_norm(method = "area") and
step_gap_derivative(), it gives the preprocessing of the winning team of
the EMSLIBS 2019 classification contest (Vrábel et al., 2020).The labels of the plots no longer overlap: they are placed by ggrepel
(now imported), away from each other and from the points, with a segment
to their point when moved far. This applies to the outlying samples of
plot_outlier_map(), plot_influence(), plot_coomans(),
plot_outliers() and plot_embedding() (and its biplot loadings), the
peaks of plot_loadings(), plot_contributions(), plot_cell_map() and
plot_fit(), the outliers of adjusted_boxplot() and
generalized_boxplot() (where many crowd together, some are left
unlabeled), the rejected shots of plot_shots(), the spectra of
plot_spectra(label_spectra = TRUE), and the lines of
plot_wavelength_calibration(), which are now all labeled.
robpca() and rospca() are faster on wide data (up to about 4 times on
full LIBS spectra): the reduction to the subspace spanned by the
observations is now computed in C++ from the smaller cross-product matrix
instead of an SVD. The results are unchanged up to numerical precision.
plot_spectra() has new arguments: panel, a grouping variable that
splits the spectra into panels, arranged one above the other or side by
side (layout = "vertical" or "horizontal"), with the offset restarting
in each panel; grid, to draw grid lines; and color_as = "fill", to show
the colors of colvar (or id) as the fill of the area under each
spectrum instead of the line color. Spectra stacked with a vertical offset
are then filled with opaque colors, the front ones hiding the back ones, as
in a waterfall plot.
plot_spectra() makes figures ready for publication:
legend shows the legend of colvar (a color bar) or id, on any side
or inside the plot, titled by legend_title (such as "K (%)");xlab, ylab and title set the axis titles and the plot title (bold,
split into a title and a subtitle when long), and the intensities are
written in plain notation ("100,000" rather than "1e+05");lines, a line list from libs_lines(), marks the emission lines with
dashed lines, labeled with their species and wavelength above the plot
(close lines share a label);palette sets the colors: a viridis palette ("viridis", the default,
"magma", "cividis", ...) or any colors;base_size and linewidth size the text and the lines for the printed
figure (such as 8 points and 0.3 for a journal column);xlim limits the wavelength range, and scales = "free_y" fits the
intensity axis to each panel;panel_tags tags the panels (a), (b), (c), ...;label_spectra labels each spectrum at its right end, with its id or
value of colvar;summary = "mean" or "median" draws the mean spectrum of each group
(panel and id) with a band of one standard deviation, or the median
with the quartiles;rasterize draws the spectra as an image within an otherwise vector
plot, for small PDF files of many long spectra (needs the ragg package).The help page gives the sizes for journal figures. The colors of colvar
now go from dark purple to yellow-green (viridis), instead of from blue to
red, readable by colorblind readers and in grayscale; use
palette = c("blue", "red") for the former colors. The axis lines are
slightly thinner.
adjusted_boxplot() and generalized_boxplot() make figures ready for
publication:
scales = "free_y", the default), so that variables of different
scales stay readable; scales = "fixed" keeps a common axis;group draws the boxes of groups (sites, treatments, ...) side by side,
and id names the observations: the outliers table gives the row
and id of every outlying value, and label_outliers labels them;show_n), and
show_mean marks the means;points = "all" draws every observation ("outliers" by default,
spread sideways so that equal values do not hide each other, or
"none");caption); the help pages give a figure legend;xlab, ylab, title (bold), base_size, horizontal, log (a
logarithmic axis) and fill (one color, or one per group or variable);stats table gives n, the fences, the notch limits, the mean
and n_outliers of each variable (and group);annotate writes statistics with each box: the shape measured by the
method (medcouple, or g and h), the number of flagged values with the
number expected in clean data, the median with its distribution-free
95% confidence interval, the IQR and the robust coefficient of variation;
it draws the fences and the biweight location, gives the p-value of a
Kruskal-Wallis test between the groups, and the number of missing
values. The stats table gains these values (n_missing,
median_lower, median_upper, iqr, rcv, biweight,
expected_outliers), and a tests table gives the tests.plot_outliers() makes figures ready for publication: readable variable
names, a "Robust z-score" axis, outliers in red among grey samples
("Regular" and "Outlier"), dashed univariate limits (z = -2.5 and 2.5),
and a caption saying how the outliers are flagged and how to read the
plot. type = "distance" draws the robust distance of each sample with
the cutoffs. id names the samples, in the table and with
label_outliers; color_by ("outlier", "distance" or "both"),
plot = FALSE (the table: row, id, distance, cutoff, outlier, weight and
robust z-scores), title, ylab and base_size. Rows with missing values
are left out, with a message, instead of stopping.
plot_fit() gives, above each panel, the center, FWHM and area of each
peak with their standard errors (rounded to two significant digits of
them), and R^2 and RMSE. The FWHM of a Voigt profile combines its widths
(Olivero and Longbothum, 1977), with a standard error from the covariance
of the estimates. show_fwhm draws the FWHM at half maximum; the baseline
is dotted, the peaks of a multi-peak fit are numbered, the residuals have
dashed limits at plus or minus twice the RMSE, and a caption names the
profiles. The data are small open circles, with a legend; xlab, ylab,
base_size, ncol (several spectra: a pair of panels each) and caption
are added.
plot_embedding() makes figures ready for publication: plain axis numbers
("-100,000" rather than "-1e+05"); base_size, which scales the text,
points, lines and labels; a caption saying what the ellipses and limits
are; legend_title; a point shape per group (shapes), for grayscale
printing; palette; aspect_ratio = "equal"; legend (position) and
panel = "white". The labels of the flagged samples alternate above and
below the points, and the short loadings of a biplot are drawn without
labels, which overlapped near the origin.
hotelling_t2() and plot_embedding() give the limit of a new sample,
method = "new" (the F limit times (n + 1) / n). The docs no longer call
the F limit the limit of a new sample.
The axis titles of plot_spectra() and plot_fit() give the units in
parentheses, as in the other plots: "Wavelength (nm)" and
"Intensity (arb. units)".
specProc now requires ggplot2 3.5.0 or later.
generalized_boxplot() flags fewer observations: its default detection
rate alpha is now about 0.7% (2 * pnorm(-4 * qnorm(0.75))), that of
Tukey's boxplot for normal data, as in the authors' implementation (the
Stata command robbox), instead of 5%. With 5%, about 18 of the 368
forageLIBS samples were flagged per mineral even in clean data; use
alpha = 0.05 for the former fences (see also the bug fixes).adjusted_boxplot() and generalized_boxplot() draw one panel per
variable, boxes in light grey (the colors of each variable added
nothing, and depended on whether ggsci was installed) and horizontal
labels. scales = "fixed" and fill give the former layout and colors.xlabels.angle, box.width, notchwidth and
staplewidth of adjusted_boxplot() and generalized_boxplot() are
renamed x_labels_angle, box_width, notch_width and staple_width,
and xlabels.hjust and xlabels.vjust are deprecated (the justification
follows the angle). The former names still work, with a warning.outliers table of adjusted_boxplot() has the columns row and
out (the tail), as that of generalized_boxplot().reject_shots() (and step_reject_shots()) rejects far fewer regular
shots. With the few shots of a sample, the MAD of the sample is unstable:
shots that differed by little from very similar shots got large
z-scores, and the correlations near 1, bounded and skewed, made small
differences look extreme. On the 2944 shots of the forage data set, 235
shots were rejected, many with a correlation above 0.99 with the median
spectrum of their sample; in simulations, 2.5% to 5% of regular shots
were. The z-scores of the correlation are now those of Fisher's z, the
intensity and the distance are on a log scale, and the scale of a sample
is its MAD but not less than the pooled MAD of all the samples (scale = "floor", the default; "sample" and "pooled" are the others): 10
shots of the forage data set are rejected (weak or strong plasmas of
other shapes), and 0.1% to 0.5% of regular shots in simulations, while
the shots of a plasma weakened by a quarter or more are still detected.
The result also gives the shot number .shot (new argument shot), the
raw criteria .intensity (relative to the median shot of the sample),
.correlation and .distance, and its settings.plot_fit() are renamed
point_size, point_color, line_width and fit_color (from pt.size,
pt.colour, line.size and line.colour); pt.shape, pt.fill,
linetype and the resid.* arguments are deprecated and ignored.plot_outliers() flags the samples beyond the adaptive cutoff of
Filzmoser, Garrett and Reimann (2005), as aq.plot() of the mvoutlier
package, instead of the fixed 97.5% chi-squared quantile, which flags
about 2.5% of clean data (the adaptive part had no effect: the cutoff was
the smaller of the two). In forageLIBS, 29 of the 368 samples are
flagged on Ca, Mg, P and K instead of 38. cutoff = "quantile" gives the
former flags. show.outlier and show.mahal are deprecated in favor of
color_by and plot = FALSE, and the table gains the columns row,
cutoff and weight.robpca() now have a fixed sign, their largest element
positive, as those of macropca() and cellpca(). The signs of loadings
and scores may differ from those of earlier versions.generalized_boxplot() now follows Bruffaerts et al. (2014) and the
authors' implementation: the transformed data are standardized by their
interquartile range divided by diff(qnorm(c(0.25, 0.75))) = 1.349 (was
1.3426); a negative tail parameter h (tails lighter than normal) is
kept instead of being set to 0, which widened the fences of such
variables (in forageLIBS, at alpha = 0.7%, the upper fence of Ca was
1.25 instead of 1.19); and the fences are never inside the box.notch = TRUE failed in adjusted_boxplot() and generalized_boxplot()
("object 'notchlower' not found"); the notches span the median
plus or minus 1.58 IQR / sqrt(n) (McGill et al., 1978).plot_spectra() drew straight lines across the gaps between detectors
(in forageLIBS, 781.5 to 789.2 nm and 800.7 to 813.9 nm), and joined the
channels of overlapping detectors in a zigzag (near 766 nm, on the K I
766.49 nm line). Each detector is now drawn as its own line.robpca() and rospca() failed, or flagged most observations as
orthogonal outliers, when k equaled the rank of the data (for example
k = 5 with 5 variables, k = n - 1 with more variables than
observations, or a k chosen by var_explained for data with few
variables): the orthogonal distances to a subspace that is the whole space
are rounding errors, and their cut-off was zero. As in Hubert, Rousseeuw
and Vanden Branden (2005), ROBPCA then reduces to the reweighted MCD of the
data: the reweighting on the orthogonal distances is skipped, and H1 is the
subset H0 of the least outlying observations (as in the rospca package).
The orthogonal distances and their cut-off are now zero, with no
orthogonal outliers, in robpca(), rospca() and macropca(), and for
new observations of predict() that lie in the subspace.robpca(), rospca() and macropca() no longer depend on the units of
the data: their numerical tolerances (singular subsets of FAST-MCD, zero
scales of the outlyingness, the cut-off of the orthogonal distances, the
rank and the C-steps of macropca(), and the constant variables of the
detection of deviating cells) are relative to the scale of the data
instead of absolute. Robust PCA failed on data with a small scale ("No
non-singular subset found" for values around 1e-8, or a zero cut-off of
the orthogonal distances when they were all below 1.5e-8), and
macropca() gave different results. With the same seed,
robpca(x * c) now gives the same subsets and outlier types as
robpca(x) for any c > 0. The final step of macropca() now always
uses the deterministic MCD when it succeeds: cellWise compares its
log-determinant with the determinant of the C-steps, a comparison that
depends on the units, and keeps it in the usual case.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.