The present note introduces how to implement the efficient and versatile kNN-LWPLSR algorithm with the Julia package Jchemo.The use of the algorithm is illustrated on the dataset challenge2008 about near-infrared (NIR) spectroscopy . Note: almost the entire content of this note can be directly transposed to the kNN-LWPLSR algorithm provided in the R package rchemo.

kNN-LWPLSR combines nearest neighborhood selection (kNN) and partial least squares regression. kNN-LWPLSR is well suited when the predictive variables (X) are multicollinear and when heterogeneity in the data generates non-linear relations between X and the response variables to predict (Y).

More generally, kNN-LWPLS is family of algorithms that can also be used for discrimination problems (i.e. when the response is categorical).

1. Preliminaries

1.1 NIRS regression

NIR spectrometry is a fast and nondestructive analytical method used in many contexts, for instance in agronomy to evaluate the nutritive quality of forages. Basically, spectral data X (matrix of n observations and p columns representing wavelengths) are collected on samples of the material to study (e.g. forages) using a spectrometer. In parallel, a set of variables of interest (the 'response' variables) Y = { y1, …, yq } (q vectors of n observations) (e.g. chemical compositions) are measured precisely in laboratory. Regression models (or discrimination models when Y represent class memberships) are then fitted on the data {X, Y} and used to predict the response variables from new spectral observations.

Spectral data have the characteristic to have highly collinear columns and, in general, matrix X is ill-conditioned. To be solved, the regression problem requires regularization methods. A very popular approach used for NIR data is the partial least squares regression (PLSR). The first step of PLSR (regularization step) is to reduce the dimension (nb. columns) of X to a limited number a << p of orthogonal vectors n x 1, maximizing the squared covariance with Y. These vectors are referred to as the PLS scores, or latent variables (LVs), often noted t and gathered in matrix T (n x a). The second step (regression) is to fit a multiple linear regression (MLR) that predicts Y from T.

PLSR is in general very efficient when the relationship between X and Y is linear. The method is fast (in particular with the Dayal & McGregor kernel #1 algorithm), even for large data. The parameter to tune is the number of scores (nb. LVs) considered in the regression model.

However, and especially in agronomy, new developments in data acquisition have resulted in increasingly large and complex datasets, posing new challenges for PLSR. Many modern datasets contain clusters due to variations in the collected products (categories, period, and areas of collection, etc.), often leading to non-linear relationships between X and Y. Since PLS relies on the assumption of linearity, its performance tends to decrease when applied to this type of data.

To address this challenge, the family of Local PLSR methods extends the standard PLS by adapting to local patterns, making it more effective in handling complex datasets. The general principle is, for each new observation to predict (xnew), to do a pre-selection of the k nearest neighbors of xnew (kNN selection step) and then to fit a PLSR model to the neighborhood. The PLSR model is used to predict the response Ynew. Two illustrations of neighborhood selection are presented below

The local PLS family includes many variants, depending essentially on

• (a) how are selected the neighborhood, and

• (b) the type of PLSR model implemented on the neighborhoods.

One of the variants is the kNN-LWPLSR approach described in section 1.2. This approach is implemented in the Jchemo function lwplsr, which is used in the present note. kNN-LWPLSR has the advantages to be fast, easy to tune and efficient for many types of data.

1.2 Theory (summary)

kNN-LWPLSR runs locally weighted PLSR (LWPLSR) on each neighborhood, instead of standard PLSR as in usual simpler local PLSR methods. LWPLSR is a particular case of weighted PLSR (WPLSR):

• In WPLSR, a n x 1 vector of weights w = ( w[1], w[2], … w[n] ) is embedded into the PLSR equations. The PLS scores are computed by maximizing w-weighted squared covariances between the scores and the response variables Y. The MLR prediction equation is computed by regressing Y on the scores using w-weighted least-squares. The w-weighting is also embedded in the centering and the eventual scaling of the data. Note: in standard PLSR, a uniform weight, 1 / n, is given to all the training observations and therefore w can be removed from the equations (incidentally, this is a particular case of WPLSR).

• Specifically in LWPLSR, the weight vector w is computed from a decreasing function, say f, of the dissimilarities (e.g., distances) between the n training observations and xnew, the observation to predict. Closer is xi to xnew, higher is the weight w[i] in the PLSR equations and therefore its importance to the prediction. This is the same distance-based principle as in the well-known locally weighted regression algorithm LOESS.

Compared to LWPLSR, kNN-LWPLSR simply adds a preliminary step: a neighborhood is selected in the training around xnew and then, for prediction, LWPLSR is applied to this neighborhood.

Concerning the prediction results, kNN-LWPLSR is equivalent to run LWPLSR on the overall training data but with a double weighting: a first binary weighting (0: xi is not a neighbor, 1: xi is a neighbor) and a second weighting, intra-neighborhood, defined by function f. In practice and for large datasets, however, kNN-LWPLSR is much faster than to compute LWPLSR on the all training set.

2. Function lwplsr

This section details the kNN-LWPLSR pipeline as defined in function lwplsr of package Jchemo.

2.1 Keyword parameters

using Jchemo   # if not loaded before

Function lwplsr has several keyword parameters that can be specified and presented in the help-page of the function

?lwplsr

  lwplsr(; kwargs...)
  lwplsr(X, Y; kwargs...)

  k-Nearest-Neighbours locally weighted partial least squares regression (kNN-LWPLSR).

    •  X : X-data (n, p).

    •  Y : Y-data (n, q).

  Keyword arguments:

    •  nlvdis : Number of latent variables (LVs) to consider in the global PLS used for the dimension reduction before computing the dissimilarities. If nlvdis = 0, there is no dimension reduction.

    •  metric : Type of dissimilarity used to select the neighbors and to compute the weights (see function getknn). Possible values are: :eucl (Euclidean), :mah (Mahalanobis), :sam (spectral angular
       distance), :cor (correlation distance).

    •  h : A scalar defining the shape of the weight function computed by function winvs. Lower is h, sharper is the function. See function winvs for details (keyword arguments criw and squared of winvs
       can also be specified here).

    •  k : The number of nearest neighbors to select for each observation to predict.

    •  tolw : For stabilization when very close neighbors.

    •  nlv : Nb. latent variables (LVs) for the local (i.e. inside each neighborhood) models.

    •  scal : Boolean. If true, (a) each column of the global X (and of the global Y if there is a preliminary PLS reduction dimension) is scaled by its uncorrected standard deviation before to compute
       the distances and the weights, and (b) the X and Y scaling is also done within each neighborhood (local level) for the weighted PLSR.

    •  verbose : Boolean. If true, predicting information are printed.

  [...]

The default values of the parameters can be displayed by

@pars lwplsr

Jchemo.ParLwplsr
  nlvdis: Int64 0
  metric: Symbol eucl
  h: Float64 Inf
  k: Int64 1
  criw: Float64 4.0
  squared: Bool false
  tolw: Float64 0.0001
  nlv: Int64 1
  scal: Bool false
  verbose: Bool false

The five main parameters to consider are: nlvdis, metric, h, k and nlv. They are are detailed below in the below sections 2.2 and 2.3.

2.2 Neighborhood computation

A first step is to choose if the dissimilarities between observations are computed after a dimension reduction of X or not. This is managed by parameter nlvdis

• If nlvdis = 0, there is not dimension reduction.

• If nlvdis > 0, a preliminary global PLS with nlvdis LVs is done on the entire dataset {X, Y} and the dissimilarities are computed on the resulting score matrix T (n x nlvdis). Note: This option only concerns the dissimilarities computation. The regression steps (LWPLSR) are always computed on the original X data.

Then, the type of dissimilarities has to be chosen, with parameter metric. The available metrics are those proposed in function getknn:

•  metric : Type of distance used for the query. Possible values are :eucl (Euclidean), :mah (Mahalanobis), 
    :sam (spectral angular distance), :cos (cosine distance), :cor (correlation distance).

To compute Mahalanobis distances on on 15-25 global nlvdis scores is often a good choice. But the best choice of metric is often dataset-dependent and no general rule can be recommanded.

Note: If X has collinear columns (which is the case for NIRS data), the use of Mahalanobis distance requires a preliminary dimension reduction since the inverse of the covariance matrix cov(X) can not be computed with stability.

2.2 Weighting function

The next paramater to set is parameter h, which determines the shape of the weight function f for LWPLSR models. Function f has a negative exponential shape whose h defines its sharpness: lower is h, sharper is function f and therefore more the closest neighbors of xnew have importance in the LWPLSR fit. The case h = Inf is the unweighted situation (all the components of w are equal), which corresponds to a more usual kNN-PLSR.

In function f, weights are computed by exp(-d / (h * MAD(d))) and are set to 0 for extreme (potentially outlier) distances such as d > Median(d) + 4 * MAD(d). Finally, weights are standardized to their maximal value. An illustration of the effect of h is given below

Many alternatives could have been considered to define the weight functions f (e.g. bicube or tricube functions). The actual negative exponential f chosen in Jchemo is versatile and easily tunable (many shapes can be represented from the variation of a single parameter only).

2.3 Dimensions of the neighborhood and the local models

The two final parameters to set are

• k: the number of observations defining the neighborhood for each observation xnew to predict.

• nlv: the number of LVs considered in the LWPLSR model fitted on the neighborhood.

In this version of the algorithm, k and nlv are the same for all the observations to predict. Note that if k is larger than the number of training observations, kNN-LWPLSR reduces to LWLSR.

3. Case study

Dataset challenge2008 was built for the prediction-challenge organized in 2018 at congress Chemometrics2018 in Paris. It consists of 4,075 NIR spectra collected on various materials typically analysed in agronomic laboratories. These materials are animal feed, rapeseed, corn gluten, grass silage, wheat, full-fat soya, milk powder and whey, maize, sunflower seed (grounded), and soya meal. The univariate response y to predict was the protein concentration.

3.1 Data importation

The dataset contains

• Object X (4075 x 680): The spectra, with wavelengths of 1120-2478 nm and a 2-nm step.

• Object Y (4075 x 4): Variable conc (protein concentration) and other meta-data.

## Preliminary loading of packages
using Jchemo       # if not loaded before
using JchemoData   # a library of various benchmark datasets
using JLD2         # package for loading/saving JLD2 data  
using CairoMakie   # plotting backend 
using FreqTables   # utilities for frequency tables

path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/challenge2018.jld2") 
@load db dat
@names dat

(:X, :Y)

X = dat.X
@head X

... (4075, 680)

Y = dat.Y
@head Y

... (4075, 4)

y = Y.conc        # variable to predict (protein concentration)

4075-element Vector{Float64}:
 12.7399998
 35.721199
 12.0
 13.8449764
 19.2999992
  7.6191959
 24.332531
 36.5717583
 55.9027786
 11.2599993
  ⋮
 65.1012344
 56.1860695
 50.8610802
  8.1399994
 13.3100004
 11.6800003
 18.2399998
 67.6700745
 25.0300007

wlst = names(X)     # wavelengths
wl = parse.(Float64, wlst)

680-element Vector{Float64}:
 1120.0
 1122.0
 1124.0
 1126.0
 1128.0
 1130.0
 1132.0
 1134.0
 1136.0
 1138.0
    ⋮
 2462.0
 2464.0
 2466.0
 2468.0
 2470.0
 2472.0
 2474.0
 2476.0
 2478.0

ntot, p = size(X)

(4075, 680)

freqtable(string.(Y.typ, " - ", Y.label))

10-element Named Vector{Int64}
Dim1                      │ 
──────────────────────────┼────
ANF - animal feed         │ 391
CLZ - rapeseed(ung)       │ 420
CNG - corn gluten         │ 395
EHH - grass silage        │ 422
FFS - full fat soya       │ 432
FRG - wheat (ung)         │ 411
MPW - milk powder & whey  │ 410
PEE - maize wp            │ 407
SFG - sun flower seed(gr) │ 281
TTS - soya meal           │ 506

The spectra can be plotted by

## For illustration purpose, only 30 spectra (randomly chosen) are plotted 
plotsp(X, wl; size = (500, 300), nsamp = 30, xlabel = "Wavelength (nm)", ylabel = "Reflectance").f

3.2 Data preprocessing

Two preprocessing steps are implemented to remove eventual non-informative physical effects in the spectra: a standard normal variation transformation (SNV), followed by a 2nd-order Savitsky-Golay derivation.

model1 = snv()
model2 = savgol(npoint = 21, deriv = 2, degree = 3)
model = pip(model1, model2)
fit!(model, X)
@head Xp = transf(model, X)

3×680 Matrix{Float64}:
 -0.00393533  -0.00441755  -0.00477681  …  0.000514478  0.000481081
 -0.00121436  -0.0013095   -0.00135921     0.000996321  0.000930884
 -0.00355712  -0.00399785  -0.00434838     0.00044084   0.000421751
... (4075, 680)

The observation of the resulting spectra clearly indicates the presence of high heterogeneity in the data

plotsp(Xp, wl; size = (500, 300), nsamp = 30, xlabel = "Wavelength (nm)").f

which is confirmed by the highly clustered pattern observed in the PCA score space (related to the various types of materials analyzed)

3.3 Split training vs. test sets

In this illustration, the split of the total data is the one given at the Paris prediction-challenge (variable Y.test in the dataset), but any other splitting could be chosen by ad'hoc sampling.

freqtable(Y.test)

2-element Named Vector{Int64}
Dim1  │ 
──────┼─────
0     │ 3701
1     │  374

The final data are given by

s = Bool.(Y.test)  # same as: s = Y.test .== 1
Xtrain = rmrow(Xp, s)
Ytrain = rmrow(Y, s)
ytrain = rmrow(y, s)
typtrain = rmrow(Y.typ, s)
Xtest = Xp[s, :]
Ytest = Y[s, :]
ytest = y[s]
typtest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = ntot, ntrain, ntest)

(ntot = 4075, ntrain = 3701, ntest = 374)

It is convenient to check that the test set is well represented or not by the training set. A first look can for instance be made by projecting the test spectra in a PCA score space built from the training spectra.

But such a 2-D (or 3-D) representation can be misleading since differences can exist in higher dimensions. A better approach is to build the plot of the score (SD) vs. orthogonal (OD) distances. It has the advantage to jointly account for all the dimensions of the PCA model.

It is also usefull to check the representativity of the y-variable.

summ(y, Y.test)

Class: 0
1×7 DataFrame
 Row │ variable  mean     std      min      max      n      nmissing
     │ Symbol    Float64  Float64  Float64  Float64  Int64  Int64
─────┼───────────────────────────────────────────────────────────────
   1 │ x1         31.894   20.297    3.061   76.604   3701         0


Class: 1
1×7 DataFrame
 Row │ variable  mean     std      min      max      n      nmissing
     │ Symbol    Float64  Float64  Float64  Float64  Int64  Int64
─────┼───────────────────────────────────────────────────────────────
   1 │ x1         32.288   20.874    2.766  75.8559    374         0

3.4 Model tuning

A usual approach of model tuning in machine learning is to split the training dataset to calibration sets vs. validation sets, and then do a grid search:

• the parameter space is reduced to a discrete grid of parameter combinations on which is explored the model performance on the validation sets. The combination showing the better performances is can generally be considered as the best model based on the available data.

The popular K-fold cross-validation (CV) is such a strategy. Nevertheless, K-fold CV requires to predict all the observations of the training set (in this illustration, ntrain = 3,701 obs.), this for each parameter combination of the grid. This can be too time consuming for local PLSR if the datasets and/or the grid size are large. A slighter approach is to do single split 'calibration/validation' (CAL and VAL sets, respectively) from the training data and to run the grid search on this single split. This strategy is used in the present note.

Below, the VAL set is selected by systematic sampling along the data but other sampling designs (e.g. random) could be chosen.

nval = 300
s = sampsys(1:ntrain, nval)
Xcal = Xtrain[s.train, :]
ycal = ytrain[s.train]
Xval = Xtrain[s.test, :]
yval = ytrain[s.test]
ncal = ntrain - nval
(ntot = ntot, ntrain, ntest, ncal, nval)

(ntot = 4075, ntrain = 3701, ntest = 374, ncal = 3401, nval = 300)

Then the grid of parameters is built by

## For this illustration, the grid has been built with a relatively low number of parameter combinations.
## More extended combinations could be considered. 
nlvdis = [15]; metric = [:mah] 
h = [1; 2; 4; 6; Inf]
k = [200; 350; 500; 1000]  
nlv = 0:15 
pars = mpar(nlvdis = nlvdis, metric = metric, h = h, k = k)  # the grid
length(pars[1])  # nb. parameter combinations considered

20

The grid search can easily be run with function gridscore, which is a generic function working that can be used by all the predictive models of package Jchemo. Note: The equivalent function for K-fold CV is gridcv.

model = lwplsr()
res = gridscore(model, Xcal, ycal, Xval, yval; 
    score = rmsep,  # performance criterion computed on VAL 
    pars,           # defined grid of parameters 
    nlv             # parameter 'nlv' has been set out of 'pars' to decrease the computation time (see gridscore help-page)  
    )
@head res   # first rows of the result table

... (320, 6)

which gives graphically

group = string.("nlvdis=", res.nlvdis, ", h=", res.h, ", k=", res.k)
plotgrid(res.nlv, res.y1, group; step = 1, xlabel ="Nb. LVs", ylabel = "RMSEP (Validation set)").f

3.5 Final predictions

The final model can be defined by selecting the best parameter combination.

u = findall(res.y1 .== minimum(res.y1))[1]
res[u, :]

and the the final prediction of the test set is given by

model = lwplsr(nlvdis = res.nlvdis[u], metric = res.metric[u], h = res.h[u], k = res.k[u], nlv = res.nlv[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@head pred

3×1 Matrix{Float64}:
 19.76328396932776
  6.200351184655297
 11.381457704579901
... (374, 1)

Summary of the predictions

mse(pred, ytest)
rmsep(pred, ytest)    # estimate of generalization error

1×1 Matrix{Float64}:
 0.7157501840557478

plotxy(pred, ytest; color = (:red, .5), bisect = true, xlabel = "Predictions (Test)", 
    ylabel = "Observed data (Test)", title = "Protein concentration (%)").f

4. Related references (not exhaustive)

• Andersson M. A comparison of nine PLS1 algorithms. J Chemom. 2009;23(10):518-529. doi:10.1002/cem.1248.

• Cleveland, W. S. (1979). Robust locally weighted regression and smoothing scatterplots. Journal of the American statistical association, 74(368), 829-836. DOI: 10.1080/01621459.1979.10481038

• Cleveland, W. S., & Devlin, S. J. (1988). Locally weighted regression: an approach to regression analysis by local fitting. Journal of the American statistical association, 83(403), 596-610. DOI:10.1080/01621459.1988.10478639

• Davrieux F, Dufour D, Dardenne P, et al. LOCAL regression algorithm improves near infrared spectroscopy predictions when the target constituent evolves in breeding populations. J Infrared Spectrosc. 2016;24(2):109. doi:10.1255/jnirs.1213

• Davrieux F, Dufour D, Dardenne P, et al. LOCAL regression algorithm improves near infrared spectroscopy predictions when the target constituent evolves in breeding populations. J Infrared Spectrosc. 2016;24(2):109. doi:10.1255/jnirs.1213

• Dayal BS, MacGregor JF. Improved PLS algorithms. J Chemom. 1997;11(1):73-85. doi:10.1002/(SICI)1099-128X(199701)11:1<73::AID-CEM435>3.0.CO;2-#.

• Höskuldsson A. PLS regression methods. J Chemom. 1988;2(3):211-228. doi:10.1002/cem.1180020306.

• Kim S, Kano M, Nakagawa H, Hasebe S. Estimation of active pharmaceutical ingredients content using locally weighted partial least squares and statistical wavelength selection. Int J Pharm. 2011;421(2):269-274. doi:10.1016/j.ijpharm.2011.10.007

• Lesnoff, M., Metz, M., Roger, J.-M., 2020. Comparison of locally weighted PLS strategies for regression and discrimination on agronomic NIR data. Journal of Chemometrics n/a, e3209. https://doi.org/10.1002/cem.3209

• Lesnoff, M., 2024. Averaging a local PLSR pipeline to predict chemical compositions and nutritive values of forages and feed from spectral near infrared data. Chemometrics and Intelligent Laboratory Systems 244, 105031. https://doi.org/10.1016/j.chemolab.2023.105031

• Manne R. Analysis of two partial-least-squares algorithms for multivariate calibration. Chemom Intell Lab Syst. 1987;2(1-3):187-197. doi:10.1016/0169-7439(87)80096-5.

• Schaal, S., Atkeson, C.G., Vijayakumar, S., 2002. Scalable Techniques from Nonparametric Statistics for Real Time Robot Learning. Applied Intelligence 17, 49–60. https://doi.org/10.1023/A:1015727715131

• Shenk J, Westerhaus M, Berzaghi P. Investigation of a LOCAL calibration procedure for near infrared instruments. J Infrared Spectrosc. 1997;5(1):223. doi:10.1255/jnirs.115

• Sicard E, Sabatier R. Theoretical framework for local PLS1 regression, and application to a rainfall data set. Comput Stat Data Anal. 2006;51(2):1393-1410. doi:10.1016/j.csda.2006.05.002.

• Wold H. Nonlinear iterative partial least squares (NIPALS) modeling: some current developments. In: Multivariate Analysis II. Wright State University, Dayton, Ohio, USA. June 19–24, 1972. New York : Academic Press: Krishnaiah , P. R.; 1973:383 – 407.

• Wold S, Sjöström M, Eriksson L. PLS-regression: a basic tool of chemometrics. Chemom Intell Lab Syst. 2001;58(2):109-130. doi:10.1016/S0169-7439(01)00155-1.

• Yoshizaki R, Kano M, Tanabe S, Miyano T. Process Parameter Optimization based on LW-PLS in Pharmaceutical Granulation Process - This work was partially supported by Japan Society for the Promotion of Science (JSPS), Grant-in-Aid for Scientific Research (C) 24560940. IFAC-Pap. 2015;48(8):303-308. doi:10.1016/j.ifacol.2015.08.198.

• Zhang X, Kano M, Li Y. Locally weighted kernel partial least squares regression based on sparse nonlinear features for virtual sensing of nonlinear time-varying processes. Comput Chem Eng. 2017;104:164-171. doi:10.1016/j.compchemeng.2017.04.014.

---