Index of functions

summary(object::Cca, X, Y)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.
• Y : The Y-data that was used to fit the model.

summary(object::Ccawold, X, Y)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.
• Y : The Y-data that was used to fit the model.

summary(object::Comdim, Xbl)

Summarize the fitted model.

• object : The fitted model.
• Xbl : The X-data that was used to fit the model.

summary(object::Fda)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.

summary(object::Kpca)

Summarize the fitted model.

• object : The fitted model.

summary(object::Mbpca, Xbl)

Summarize the fitted model.

• object : The fitted model.
• Xbl : The X-data that was used to fit the model.

summary(object::Mbplsr, Xbl)

Summarize the fitted model.

• object : The fitted model.
• Xbl : The X-data that was used to fit the model.

summary(object::Mbplswest, Xbl)

Summarize the fitted model.

• object : The fitted model.
• Xbl : The X-data that was used to fit the model.

summary(object::Pca, X)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.

summary(object::Plscan, X, Y)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.
• Y : The Y-data that was used to fit the model.

summary(object::Plstuck, X, Y)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.
• Y : The Y-data that was used to fit the model.

summary(object::Rasvd, X, Y)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.
• Y : The Y-data that was used to fit the model.

summary(object::Spca, X)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.

summary(object::Union{Pcr, Spcr}, X)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.

summary(object::Union{Plsr, Splsr}, X)

Summarize the fitted model.

• object : The fitted model.
• X : The X-data that was used to fit the model.

aggstat(X, y; algo = mean)
aggstat(X::DataFrame; vary, vargroup, algo = mean)

Compute column-wise statistics by class in a dataset.

• X : Data (n, p).
• y : A categorical variable (n) (class membership).
• algo : Function to compute (default = mean).

Specific for dataframes:

• vary : Vector of the names of the variables to summarize.
• vargroup : Vector of the names of the categorical variables to consider for computations by class.

Variables defined in vary and vargroup must be columns of X.

Return a matrix or, if only argument X::DataFrame is used, a dataframe.

Examples

using Jchemo, DataFrames, Statistics

n, p = 20, 5
X = rand(n, p)
df = DataFrame(X, :auto)
y = rand(1:3, n)
res = aggstat(X, y; algo = sum)
res.X
aggstat(df, y; algo = sum).X

n, p = 20, 5
X = rand(n, p)
df = DataFrame(X, string.("v", 1:p))
df.gr1 = rand(1:2, n)
df.gr2 = rand(["a", "b", "c"], n)
df
aggstat(df; vary = [:v1, :v2], vargroup = [:gr1, :gr2], algo = var)

aggsum(x::Vector, y::Union{Vector, BitVector})

Compute sub-total sums by class of a categorical variable.

• x : A quantitative variable to sum (n)
• y : A categorical variable (n) (class membership).

Return a vector.

Examples

using Jchemo

x = rand(1000)
y = vcat(rand(["a" ; "c"], 900), repeat(["b"], 100))
aggsum(x, y)

aicplsr(X, y; alpha = 2, kwargs...)

Compute Akaike's (AIC) and Mallows's (Cp) criteria for univariate PLSR models.

• X : X-data (n, p).
• y : Univariate Y-data.

Keyword arguments:

• Same arguments as those of function cglsr.
• alpha : Coefficient multiplicating the model complexity (df) to compute AIC.

The function uses function dfplsr_cg.

References

Hansen, P.C., 1998. Rank-Deficient and Discrete Ill-Posed Problems, Mathematical Modeling and Computation. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898719697

Hansen, P.C., 2008. Regularization Tools version 4.0 for Matlab 7.3. Numer Algor 46, 189–194. https://doi.org/10.1007/s11075-007-9136-9

Lesnoff, M., Roger, J.-M., Rutledge, D.N., 2021. Monte Carlo methods for estimating Mallows’s Cp and AIC criteria for PLSR models. Illustration on agronomic spectroscopic NIR data. Journal of Chemometrics n/a, e3369. https://doi.org/10.1002/cem.3369

Examples

using Jchemo, JchemoData, CairoMakie
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlv = 40
res = aicplsr(X, y; nlv) ;
res.crit
res.opt
res.delta

zaic = res.crit.aic
f, ax = plotgrid(0:nlv, zaic; xlabel = "Nb. LVs", ylabel = "AIC")
scatter!(ax, 0:nlv, zaic)
f

aov1(x, Y)
One-factor ANOVA test.

• x : Univariate categorical (factor) data (n).
• Y : Y-data (n, q).

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/iris.jld2") 
@load db dat
@names dat
@head dat.X
x = dat.X[:, 5]
Y = dat.X[:, 1:4]
tab(x) 

res = aov1(x, Y) ;
@names res
res.SSF
res.SSR 
res.F 
res.pval

bias(pred, Y)

Compute the prediction bias, i.e. the opposite of the mean prediction error.

• pred : Predictions.
• Y : Observed data.

Examples

using Jchemo

Xtrain = rand(10, 5) 
Ytrain = rand(10, 2)
ytrain = Ytrain[:, 1]
Xtest = rand(4, 5) 
Ytest = rand(4, 2)
ytest = Ytest[:, 1]

model = plskern(nlv = 2)
fit!(model, Xtrain, Ytrain)
pred = predict(model, Xtest).pred
bias(pred, Ytest)

model = plskern(nlv = 2)
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
bias(pred, ytest)

blockscal(; kwargs...)
blockscal(Xbl; kwargs...)
blockscal(Xbl, weights::Weight; kwargs...)

Scale multiblock X-data.

• Xbl : List of blocks (vector of matrices) of X-data Typically, output of function mblock from (n, p) data.
• weights : Weights (n) of the observations (rows of the blocks). Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• centr : Boolean. If true, each column of blocks in Xbl is centered (before the block scaling).
• scal : Boolean. If true, each column of blocks in Xbl is scaled by its uncorrected standard deviation (before the block scaling).
• bscal : Type of block scaling. Possible values are: :none, :frob, :mfa, :ncol, :sd. See thereafter.

If implemented, the data transformations follow the order: column centering, column scaling and finally block scaling.

Types of block scaling:

• :none : No block scaling.
• :frob : Let D be the diagonal matrix of vector weights.w. Each block X is divided by its Frobenius norm = sqrt(tr(X' * D * X)). After this scaling, tr(X' * D * X) = 1.
• :mfa : Each block X is divided by sv, where sv is the dominant singular value of X (this is the "MFA" approach; "AFM "in French).
• :ncol : Each block X is divided by the nb. of columns of the block.
• :sd : Each block X is divided by sqrt(sum(weighted variances of the block-columns)). After this scaling, sum(weighted variances of the block-columns) = 1.

Examples

using Jchemo
n = 5 ; m = 3 ; p = 10 
X = rand(n, p) 
Xnew = rand(m, p)
listbl = [3:4, 1, [6; 8:10]]
Xbl = mblock(X, listbl) 
Xblnew = mblock(Xnew, listbl) 
@head Xbl[3]

centr = true ; scal = true
bscal = :frob
model = blockscal(; centr, scal, bscal)
fit!(model, Xbl)
## Data transformation
zXbl = transf(model, Xbl) ; 
@head zXbl[3]

zXblnew = transf(model, Xblnew) ; 
zXblnew[3]

calds(; algo = plskern, kwargs...)
calds(X1, X2; algo = plskern, kwargs...)

Direct standardization (DS) for calibration transfer of spectral data.

• X1 : Spectra (n, p) to transfer to the target.
• X2 : Target spectra (n, p).

Keyword arguments:

• algo : Function used as transfer model.
• kwargs : Optional arguments for algo.

X1 and X2 must represent the same n samples ("standards").

The objective is to transform spectra X1 to new spectra as close as possible as the target X2. Method DS fits a model (defined in algo) that predicts X2 from X1.

References

Y. Wang, D. J. Veltkamp, and B. R. Kowalski, “Multivariate Instrument Standardization,” Anal. Chem., vol. 63, no. 23, pp. 2750–2756, 1991, doi: 10.1021/ac00023a016.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/caltransfer.jld2")
@load db dat
@names dat
## Objects X1 and X2 are spectra collected 
## on the same samples. 
## X2 represents the target space. 
## We want to transfer X1 in the same space
## as X2.
## Data to transfer
X1cal = dat.X1cal
X1val = dat.X1val
n = nro(X1cal)
m = nro(X1val)
## Target space
X2cal = dat.X2cal
X2val = dat.X2val

## Fitting the model
fitm = calds(X1cal, X2cal; algo = plskern, nlv = 10) 
#fitm = calds(X1cal, X2cal; algo = mlrpinv)   # less robust 

## Transfer of new spectra X1val 
## expected to be close to X2val
pred = predict(fitm, X1val).pred

i = 1
f = Figure(size = (500, 300))
ax = Axis(f[1, 1])
lines!(X2val[i, :]; label = "x2")
lines!(ax, X1val[i, :]; label = "x1")
lines!(pred[i, :]; linestyle = :dash, label = "x1_corrected")
axislegend(position = :rb, framevisible = false)
f

calpds(; npoint = 5, algo = plskern, kwargs...)
calpds(X1, X2; npoint = 5, algo = plskern, kwargs...)

Piecewise direct standardization (PDS) for calibration transfer of spectral data.

• X1 : Spectra (n, p) to transfer to the target.
• X2 : Target spectra (n, p).

Keyword arguments:

• npoint : Half-window size (nb. points left or right to the given wavelength).
• algo : Function used as transfer model.
• kwargs : Optional arguments for algo.

X1 and X2 must represent the same n standard samples.

The objective is to transform spectra X1 to new spectra as close as possible as the target X2. Method PDS fits models (defined in algo) that predict X2 from X1.

The window used in X1 to predict wavelength "i" in X2 is:

• i - npoint, i - npoint + 1, ..., i, ..., i + npoint - 1, i + npoint

References

Bouveresse, E., Massart, D.L., 1996. Improvement of the piecewise direct targetisation procedure for the transfer of NIR spectra for multivariate calibration. Chemometrics and Intelligent Laboratory Systems 32, 201–213. https://doi.org/10.1016/0169-7439(95)00074-7

Y. Wang, D. J. Veltkamp, and B. R. Kowalski, “Multivariate Instrument Standardization,” Anal. Chem., vol. 63, no. 23, pp. 2750–2756, 1991, doi: 10.1021/ac00023a016.

Wülfert, F., Kok, W.Th., Noord, O.E. de, Smilde, A.K., 2000. Correction of Temperature-Induced Spectral Variation by Continuous Piecewise Direct Standardization. Anal. Chem. 72, 1639–1644. https://doi.org/10.1021/ac9906835

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/caltransfer.jld2")
@load db dat
@names dat
## Objects X1 and X2 are spectra collected 
## on the same samples. 
## X2 represents the target space. 
## We want to transfer X1 in the same space
## as X2.
## Data to transfer
X1cal = dat.X1cal
X1val = dat.X1val
n = nro(X1cal)
m = nro(X1val)
## Target space
X2cal = dat.X2cal
X2val = dat.X2val

## Fitting the model
fitm = calpds(X1cal, X2cal; npoint = 2, algo = plskern, nlv = 2) 

## Transfer of new spectra X1val 
## expected to be close to X2val
pred = predict(fitm, X1val).pred

i = 1
f = Figure(size = (500, 300))
ax = Axis(f[1, 1])
lines!(X2val[i, :]; label = "x2")
lines!(ax, X1val[i, :]; label = "x1")
lines!(pred[i, :]; linestyle = :dash, label = "x1_corrected")
axislegend(position = :rb, framevisible = false)
f

cca(; kwargs...)
cca(X, Y; kwargs...)
cca(X, Y, weights::Weight; kwargs...)
cca!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Canonical correlation Analysis (CCA, RCCA).

• X : First block of data.
• Y : Second block of data.
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. Possible values are: :none, :frob. See functions blockscal.
• tau : Regularization parameter (∊ [0, 1]).
• scal : Boolean. If true, each column of blocks X and Y is scaled by its uncorrected standard deviation (before the block scaling).

This function implements a CCA algorithm using SVD decompositions and presented in Weenink 2003 section 2.

A continuum regularization is available (parameter tau). After block centering and scaling, the function returns block LVs (Tx and Ty) that are proportionnal to the eigenvectors of Projx * Projy and Projy * Projx, respectively, defined as follows:

• Cx = (1 - tau) * X'DX + tau * Ix
• Cy = (1 - tau) * Y'DY + tau * Iy
• Cxy = X'DY
• Projx = sqrt(D) * X * invCx * X' * sqrt(D)
• Projy = sqrt(D) * Y * invCx * Y' * sqrt(D)

where D is the observation (row) metric. Value tau = 0 can generate unstability when inverting the covariance matrices. Often, a better alternative is to use an epsilon value (e.g. tau = 1e-8) to get similar results as with pseudo-inverses.

After normalized (and using uniform weights), the scores returned by the function are expected to be the same as those returned by functions rcc of the R packages CCA (González et al.) and mixOmics (Lê Cao et al.) whith their parameters lambda1 and lambda2 set to:

• lambda1 = lambda2 = tau / (1 - tau) * n / (n - 1)

See function plscan for the details on the summary outputs.

References

González, I., Déjean, S., Martin, P.G.P., Baccini, A., 2008. CCA: An R Package to Extend Canonical Correlation Analysis. Journal of Statistical Software 23, 1-14. https://doi.org/10.18637/jss.v023.i12

Hotelling, H. (1936): “Relations between two sets of variates”, Biometrika 28: pp. 321–377.

Lê Cao, K.-A., Rohart, F., Gonzalez, I., Dejean, S., Abadi, A.J., Gautier, B., Bartolo, F., Monget, P., Coquery, J., Yao, F., Liquet, B., 2022. mixOmics: Omics Data Integration Project. https://doi.org/10.18129/B9.bioc.mixOmics

Weenink, D. 2003. Canonical Correlation Analysis, Institute of Phonetic Sciences, Univ. of Amsterdam, Proceedings 25, 81-99.

Examples

using Jchemo, JchemoData, JLD2
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "linnerud.jld2") 
@load db dat
@names dat
X = dat.X
Y = dat.Y
n, p = size(X)
q = nco(Y)

nlv = 3
bscal = :frob ; tau = 1e-8
model = cca(; nlv, bscal, tau)
fit!(model, X, Y)
@names model
@names model.fitm

@head model.fitm.Tx
@head transfbl(model, X, Y).Tx

@head model.fitm.Ty
@head transfbl(model, X, Y).Ty

res = summary(model, X, Y) ;
@names res
res.cortx2ty
res.rvx2tx
res.rvy2ty
res.rdx2tx
res.rdy2ty
res.corx2tx 
res.cory2ty 

ccawold(; kwargs...)
ccawold(X, Y; kwargs...)
ccawold(X, Y, weights::Weight; kwargs...)
ccawold!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Canonical correlation analysis (CCA, RCCA) - Wold Nipals algorithm.

• X : First block of data.
• Y : Second block of data.
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. Possible values are: :none, :frob. See functions blockscal.
• tau : Regularization parameter (∊ [0, 1]).
• tol : Tolerance value for convergence (Nipals).
• maxit : Maximum number of iterations (Nipals).
• scal : Boolean. If true, each column of blocks X and Y is scaled by its uncorrected standard deviation (before the block scaling).

This function implements the Nipals ccawold algorithm presented by Tenenhaus 1998 p.204 (related to Wold et al. 1984).

In this implementation, after each step of LVs computation, X and Y are deflated relatively to their respective scores (tx and ty).

A continuum regularization is available (parameter tau). After block centering and scaling, the covariances matrices are computed as follows:

• Cx = (1 - tau) * X'DX + tau * Ix
• Cy = (1 - tau) * Y'DY + tau * Iy

where D is the observation (row) metric. Value tau = 0 can generate unstability when inverting the covariance matrices. Often, a better alternative is to use an epsilon value (e.g. tau = 1e-8) to get similar results as with pseudo-inverses.

The normed scores returned by the function are expected (using uniform weights) to be the same as those returned by function rgcca of the R package RGCCA (Tenenhaus & Guillemot 2017, Tenenhaus et al. 2017).

See function plscan for the details on the summary outputs. See function plscan for the details on the summary outputs.

References

Tenenhaus, A., Guillemot, V. 2017. RGCCA: Regularized and Sparse Generalized Canonical Correlation Analysis for Multiblock Data Multiblock data analysis. https://cran.r-project.org/web/packages/RGCCA/index.html

Tenenhaus, M., 1998. La régression PLS: théorie et pratique. Editions Technip, Paris.

Tenenhaus, M., Tenenhaus, A., Groenen, P.J.F., 2017. Regularized Generalized Canonical Correlation Analysis: A Framework for Sequential Multiblock Component Methods. Psychometrika 82, 737–777. https://doi.org/10.1007/s11336-017-9573-x

Wold, S., Ruhe, A., Wold, H., Dunn, III, W.J., 1984. The Collinearity Problem in Linear Regression. The Partial Least Squares (PLS) Approach to Generalized Inverses. SIAM Journal on Scientific and Statistical Computing 5, 735–743. https://doi.org/10.1137/0905052

Examples

using Jchemo, JchemoData, JLD2
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "linnerud.jld2") 
@load db dat
@names dat
X = dat.X
Y = dat.Y
n, p = size(X)
q = nco(Y)

nlv = 2
bscal = :frob ; tau = 1e-4
model = ccawold(; nlv, bscal, tau, tol = 1e-10)
fit!(model, X, Y)
@names model
@names model.fitm

@head model.fitm.Tx
@head transfbl(model, X, Y).Tx

@head model.fitm.Ty
@head transfbl(model, X, Y).Ty

res = summary(model, X, Y) ;
@names res
res.explvarx
res.explvary
res.cortx2ty
res.rvx2tx
res.rvy2ty
res.rdx2tx
res.rdy2ty
res.corx2tx 
res.cory2ty 

center()
center(X)
center(X, weights::Weight)

Column-wise centering of X-data.

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

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

model = center() 
fit!(model, Xtrain)
Xptrain = transf(model, Xtrain)
Xptest = transf(model, Xtest)
colmean(Xptrain)
@head Xptest 
@head Xtest .- colmean(Xtrain)'
plotsp(Xptrain).f
plotsp(Xptest).f

cglsr(; kwargs...)
cglsr(X, y; kwargs...)
cglsr!(X::Matrix, y::Matrix; kwargs...)

Conjugate gradient algorithm for the normal equations (CGLS; Björck 1996).

• X : X-data (n, p).
• y : Univariate Y-data (n).

Keyword arguments:

• nlv : Nb. CG iterations.
• gs : Boolean. If true (default), a Gram-Schmidt orthogonalization of the normal equation residual vectors is done.
• filt : Boolean. If true, CG filter factors are computed (output F). Default = false.
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

CGLS algorithm "7.4.1" Bjorck 1996, p.289. In the present function, the part of the code computing the re-orthogonalization (Hansen 1998) and filter factors (Vogel 1987, Hansen 1998) is a transcription (with few adaptations) of the Matlab function cgls (Saunders et al. https://web.stanford.edu/group/SOL/software/cgls/; Hansen 2008).

References

Björck, A., 1996. Numerical Methods for Least Squares Problems, Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9781611971484

Hansen, P.C., 1998. Rank-Deficient and Discrete Ill-Posed Problems, Mathematical Modeling and Computation. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898719697

Hansen, P.C., 2008. Regularization Tools version 4.0 for Matlab 7.3. Numer Algor 46, 189–194. https://doi.org/10.1007/s11075-007-9136-9

Manne R. Analysis of two partial-least-squares algorithms for multivariate calibration. Chemometrics Intell. Lab. Syst. 1987, 2: 187–197.

Phatak A, De Hoog F. Exploiting the connection between PLS, Lanczos methods and conjugate gradients: alternative proofs of some properties of PLS. J. Chemometrics 2002; 16: 361–367.

Vogel, C. R., "Solving ill-conditioned linear systems using the conjugate gradient method", Report, Dept. of Mathematical Sciences, Montana State University, 1987.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlv = 5 ; scal = true
model = cglsr(; nlv, scal)
fit!(model, Xtrain, ytrain)
@names model.fitm 
@head model.fitm.B
coef(model.fitm).B
coef(model.fitm).int

pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f   

coef(object::Cglsr)

Compute the b-coefficients of a fitted model.

• object : The fitted model.

coef(object::Dkplsr; nlv = nothing)

Compute the b-coefficients of a fitted model.

• object : The fitted model.
• nlv : Nb. LVs, or collection of nb. LVs, to consider.

coef(object::Kplsr; nlv = nothing)

Compute the b-coefficients of a fitted model.

• object : The fitted model.
• nlv : Nb. LVs, or collection of nb. LVs, to consider.

coef(object::Krr; lb = nothing)

Compute the b-coefficients of a fitted model.

• object : The fitted model.
• lb : Ridge regularization parameter "lambda".

coef(object::Pcr; nlv = nothing)

Compute the b-coefficients of a LV model.

• object : The fitted model.
• nlv : Nb. LVs to consider.

For a model fitted from X(n, p) and Y(n, q), the returned object B is a matrix (p, q). If nlv = 0, B is a matrix of zeros. The returned object int is the intercept.

coef(object::Rosaplsr; nlv = nothing)

Compute the X b-coefficients of a model fitted with nlv LVs.

• object : The fitted model.
• nlv : Nb. LVs to consider.

coef(object::Rr; lb = nothing)

Compute the b-coefficients of a fitted model.

• object : The fitted model.
• lb : Ridge regularization parameter "lambda".

coef(object::Mlr)

Compute the coefficients of the fitted model.

• object : The fitted model.

coef(object::Union{Plsr, Pcr, Splsr}; nlv = nothing)

Compute the b-coefficients of a LV model.

• object : The fitted model.
• nlv : Nb. LVs to consider.

For a model fitted from X(n, p) and Y(n, q), the returned object B is a matrix (p, q). If nlv = 0, B is a matrix of zeros. The returned object int is the intercept.

colmad(X)

Compute column-wise median absolute deviations (MAD) of a matrix.

• X : Data (n, p).

Return a vector.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)

colmad(X)

colmean(X)
colmean(X, weights::Weight)

Compute column-wise means of a matrix.

• X : Data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Return a vector.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
w = mweight(rand(n))

colmean(X)
colmean(X, w)

colmed(X)

Compute column-wise medians of a matrix.

• X : Data (n, p).

Return a vector.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)

colmed(X)

colnorm(X)
colnorm(X, weights::Weight)

Compute column-wise norms of a matrix.

• X : Data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

The norm computed for a column x of X is:

• sqrt(x' * x)

The weighted norm is:

• sqrt(x' * D * x), where D is the diagonal matrix of weights.w

Warning: colnorm(X, mweight(ones(n))) = colnorm(X) / sqrt(n).

Return a vector.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
w = mweight(rand(n))

colnorm(X)
colnorm(X, w)

colstd(X)
colstd(X, weights::Weight)

Compute column-wise standard deviations (uncorrected) of a matrix.

• X : Data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Return a vector.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
w = mweight(rand(n))

colstd(X)
colstd(X, w)

colsum(X)
colsum(X, weights::Weight)

Compute column-wise sums of a matrix.

• X : Data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Return a vector.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
w = mweight(rand(n))

colsum(X)
colsum(X, w)

colvar(X)
colvar(X, weights::Weight)

Compute column-wise variances (uncorrected) of a matrix.

• X : Data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Return a vector.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
w = mweight(rand(n))

colvar(X)
colvar(X, w)

comdim(; kwargs...)
comdim(Xbl; kwargs...)
comdim(Xbl, weights::Weight; kwargs...)
comdim!(Xbl::Matrix, weights::Weight; kwargs...)

Common components and specific weights analysis (CCSWA, a.k.a ComDim).

• Xbl : List of blocks (vector of matrices) of X-data. Typically, output of function mblock.
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. global latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. See function blockscal for possible values.
• tol : Tolerance value for convergence (Nipals).
• maxit : Maximum number of iterations (Nipals).
• scal : Boolean. If true, each column of blocks in Xbl is scaled by its uncorrected standard deviation (before the block scaling).

"SVD" algorithm of Hannafi & Qannari 2008 p.84.

The function returns several objects, in particular:

• T : The global LVs (not-normed).
• U : The global LVs (normed).
• W : The block weights (normed).
• Tb : The block LVs (in the metric scale), returned grouped by LV.
• Tbl : The block LVs (in the original scale), returned grouped by block.
• Vbl : The block loadings (normed).
• lb : The block specific weights (saliences) 'lambda'.
• mu : The sum of the block specific weights.

Function summary returns:

• explvarx : Proportion of the total X inertia (squared Frobenious norm) explained by the global LVs.
• explvarxx : Proportion of the total XX' inertia explained by each global score (= indicator "V" in Qannari et al. 2000, Hanafi et al. 2008).
• explxbl : Proportion of the inertia of each block (= Xbl[k]) explained by the global LVs.
• psal2 : Proportion of the squared saliences of each block within each global score.
• contrxbl2t : Contribution of each block to the global LVs.
• rvxbl2t : RV coefficients between each block and the global LVs.
• rdxbl2t : Rd coefficients between each block and the global LVs.
• cortbl2t : Correlations between the block LVs (= Tbl[k]) and the global LVs.
• corx2t : Correlation between the X-variables and the global LVs.

References

Cariou, V., Qannari, E.M., Rutledge, D.N., Vigneau, E., 2018. ComDim: From multiblock data analysis to path modeling. Food Quality and Preference, Sensometrics 2016: Sensometrics-by-the-Sea 67, 27–34. https://doi.org/10.1016/j.foodqual.2017.02.012

Cariou, V., Jouan-Rimbaud Bouveresse, D., Qannari, E.M., Rutledge, D.N., 2019. Chapter 7 - ComDim Methods for the Analysis of Multiblock Data in a Data Fusion Perspective, in: Cocchi, M. (Ed.), Data Handling in Science and Technology, Data Fusion Methodology and Applications. Elsevier, pp. 179–204. https://doi.org/10.1016/B978-0-444-63984-4.00007-7

Ghaziri, A.E., Cariou, V., Rutledge, D.N., Qannari, E.M., 2016. Analysis of multiblock datasets using ComDim: Overview and extension to the analysis of (K + 1) datasets. Journal of Chemometrics 30, 420–429. https://doi.org/10.1002/cem.2810

Hanafi, M., 2008. Nouvelles propriétés de l’analyse en composantes communes et poids spécifiques. Journal de la société française de statistique 149, 75–97.

Qannari, E.M., Wakeling, I., Courcoux, P., MacFie, H.J.H., 2000. Defining the underlying sensory dimensions. Food Quality and Preference 11, 151–154. https://doi.org/10.1016/S0950-3293(99)00069-5

Examples

using Jchemo, JchemoData, JLD2
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "ham.jld2") 
@load db dat
@names dat 
X = dat.X
group = dat.group
listbl = [1:11, 12:19, 20:25]
Xbl = mblock(X[1:6, :], listbl)
Xblnew = mblock(X[7:8, :], listbl)
n = nro(Xbl[1]) 

nlv = 3
bscal = :frob
scal = false
#scal = true
model = comdim(; nlv, bscal, scal)
fit!(model, Xbl)
@names model 
@names model.fitm
## Global scores 
@head model.fitm.T
@head transf(model, Xbl)
transf(model, Xblnew)
## Blocks scores
i = 1
@head model.fitm.Tbl[i]
@head transfbl(model, Xbl)[i]

res = summary(model, Xbl) ;
@names res 
res.explvarx
res.explvarxx
res.psal2 
res.contrxbl2t
res.explxbl   # = model.fitm.lb if bscal = :frob
rowsum(Matrix(res.explxbl))
res.rvxbl2t
res.rdxbl2t
res.cortbl2t
res.corx2t 

conf(pred, y; digits = 1)

Confusion matrix.

• pred : Univariate predictions.
• y : Univariate observed data.

Keyword arguments:

• digits : Nb. digits used to round percentages.

Examples

using Jchemo, CairoMakie

y = ["d"; "c"; "b"; "c"; "a"; "d"; "b"; "d"; 
    "b"; "b"; "a"; "a"; "c"; "d"; "d"]
pred = ["a"; "d"; "b"; "d"; "b"; "d"; "b"; "d"; 
    "b"; "b"; "a"; "a"; "d"; "d"; "d"]
#y = rand(1:10, 200); pred = rand(1:10, 200)

res = conf(pred, y) ;
@names res
res.cnt       # Counts (dataframe built from `A`) 
res.pct       # Row %  (dataframe built from `Apct`))
res.A         
res.Apct
res.diagpct
res.accpct    # Accuracy (% classification successes)
res.lev       # Levels

plotconf(res).f

plotconf(res; cnt = false, ptext = false).f

convertdf(df::DataFrame; miss = nothing, typ)

Convert the columns of a dataframe to given types.

• df : A dataframe.
• miss : The code used in df to identify the data to be declared as missing (of type Missing). See function recod_miss.
• typ : A vector of the targeted types for the columns of the new dataframe.

Examples

using Jchemo, DataFrames

cor2(pred, Y)

Compute the squared linear correlation between data and predictions.

• pred : Predictions.
• Y : Observed data.

Examples

using Jchemo 

Xtrain = rand(10, 5) 
Ytrain = rand(10, 2)
ytrain = Ytrain[:, 1]
Xtest = rand(4, 5) 
Ytest = rand(4, 2)
ytest = Ytest[:, 1]

model = plskern; nlv = 2)
fit!(model, Xtrain, Ytrain)
pred = predict(model, Xtest).pred
cor2(pred, Ytest)

model = plskern; nlv = 2)
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
cor2(pred, ytest)

corm(X) 
corm(X, Y) 
corm(X, weights::Weight)
corm(X, Y, weights::Weight)

Compute a weighted correlation matrix.

• X : Data (n, p).
• Y : Data (n, q).
• weights : Weights (n) of the observations. Object of type Weight (e.g. generated by function mweight).

Uncorrected correlation matrix

• of X-columns : ==> (p, p) matrix
• or between X-columns and Y-columns : ==> (p, q) matrix.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
Y = rand(n, 3)
w = mweight(rand(n))

corm(X, w)
corm(X, Y, w)

corv(x, y)

Compute correlation between two vectors.

• x : vector (n).
• y : vector (n).

References

@Stevengj, https://discourse.julialang.org/t/interesting-post-about-simd-dot-product-and-cosine-similarity/123282.

Examples

using Jchemo

n = 5
x = rand(n)
y = rand(n)

corv(x, y)

cosm(X)
cosm(X, Y)

Compute a cosinus matrix.

• X : Data (n, p).
• Y : Data (n, q).

The function computes the cosinus matrix:

• of the columns of X: ==> (p, p) matrix
• or between columns of X and Y : ==> (p, q) matrix.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
Y = rand(n, 3)

cosm(X)
cosm(X, Y)

cosv(x, y)

Compute cosinus between two vectors.

• x : vector (n).
• y : vector (n).

References

@Stevengj, https://discourse.julialang.org/t/interesting-post-about-simd-dot-product-and-cosine-similarity/123282.

Examples

using Jchemo

n = 5
x = rand(n)
y = rand(n)

cosv(x, y)

covm(X)
covm(X, weights::Weight)
covm(X, Y) 
covm(X, Y, weights::Weight)

Compute a weighted covariance matrix.

• X : Data (n, p).
• Y : Data (n, q).
• weights : Weights (n) of the observations. Object of type Weight (e.g. generated by function mweight).

The function computes the uncorrected covariance matrix:

• of the columns of X: ==> (p, p) matrix
• or between columns of X and Y : ==> (p, q) matrix.

Examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
Y = rand(n, 3)
w = mweight(rand(n))

covm(X, w)
covm(X, Y, w)

cosv(x, y)

Compute uncorrected covariance between two vectors.

• x : vector (n).
• y : vector (n).

References

@Stevengj, https://discourse.julialang.org/t/interesting-post-about-simd-dot-product-and-cosine-similarity/123282.

Examples

using Jchemo

n = 5
x = rand(n)
y = rand(n)

covv(x, y)

cscale()
cscale(X)
cscale(X, weights::Weight)

Column-wise centering and scaling of X-data.

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

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))

db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

model = cscale() 
fit!(model, Xtrain)
Xptrain = transf(model, Xtrain)
Xptest = transf(model, Xtest)
colmean(Xptrain)
colstd(Xptrain)
@head Xptest 
@head (Xtest .- colmean(Xtrain)') ./ colstd(Xtrain)'
plotsp(Xptrain).f
plotsp(Xptest).f

detrend_airpls(; kwargs...)
detrend_airpls(X; kwargs...)

Baseline correction of each row of X-data by adaptive iteratively reweighted penalized least squares algorithm (AIRPLS).

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

Keyword arguments:

• lb : Penalizing (smoothing) parameter "lambda".
• maxit : Maximum number of iterations.
• verbose : If true, nb. iterations are printed.

De-trend transformation: the function fits a baseline by AIRPLS (see Zhang et al. 2010, and Baek et al. 2015 section 2) for each observation and returns the residuals (= signals corrected from the baseline).

References

Baek, S.-J., Park, A., Ahn, Y.-J., Choo, J., 2015. Baseline correction using asymmetrically reweighted penalized least squares smoothing. Analyst 140, 250–257. https://doi.org/10.1039/C4AN01061B

Zhang, Z.-M., Chen, S., Liang, Y.-Z., 2010. Baseline correction using adaptive iteratively reweighted penalized least squares. Analyst 135, 1138–1146. https://doi.org/10.1039/B922045C

https://github.com/zmzhang/airPLS/tree/master

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

## Example on 1 spectrum
i = 2
zX = Matrix(X)[i:i, :]
lb = 1e6
model = detrend_airpls(; lb)
fit!(model, zX)
zXc = transf(model, zX)   # = corrected spectrum 
B = zX - zXc              # = estimated baseline
f, ax = plotsp(zX, wl)
lines!(wl, vec(B); color = :blue)
lines!(wl, vec(zXc); color = :black)
f

detrend_arpls(; kwargs...)
detrend_arpls(X; kwargs...)

Baseline correction of each row of X-data by asymmetrically reweighted penalized least squares smoothing (ARPLS).

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

Keyword arguments:

• lb : Penalizing (smoothness) parameter "lambda".
• tol : Tolerance value for stopping the iterations.
• maxit : Maximum number of iterations.
• verbose : If true, nb. iterations are printed.

De-trend transformation: the function fits a baseline by ARPLS (see Baek et al. 2015 section 3) for each observation and returns the residuals (= signals corrected from the baseline).

References

Baek, S.-J., Park, A., Ahn, Y.-J., Choo, J., 2015. Baseline correction using asymmetrically reweighted penalized least squares smoothing. Analyst 140, 250–257. https://doi.org/10.1039/C4AN01061B

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

## Example on 1 spectrum
i = 2
zX = Matrix(X)[i:i, :]
lb = 1e4
model = detrend_arpls(; lb, p)
fit!(model, zX)
zXc = transf(model, zX)   # = corrected spectrum 
B = zX - zXc              # = estimated baseline
f, ax = plotsp(zX, wl)
lines!(wl, vec(B); color = :blue)
lines!(wl, vec(zXc); color = :black)
f

detrend_asls(; kwargs...)
detrend_asls(X; kwargs...)

Baseline correction of each row of X-data by asymmetric least squares algorithm (ASLS).

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

Keyword arguments:

• lb : Penalizing (smoothness) parameter "lambda".
• p : Asymmetry parameter (0 < p << 1).
• tol : Tolerance value for stopping the iterations.
• maxit : Maximum number of iterations.
• verbose : If true, nb. iterations are printed.

De-trend transformation: the function fits a baseline by ASLS (see Baek et al. 2015 section 2) for each observation and returns the residuals (= signals corrected from the baseline).

Generally 0.001 ≤ p ≤ 0.1 is a good choice (for a signal with positive peaks) and 1e2 ≤ lb ≤ 1e9, but exceptions may occur (Eilers & Boelens 2005).

References

Baek, S.-J., Park, A., Ahn, Y.-J., Choo, J., 2015. Baseline correction using asymmetrically reweighted penalized least squares smoothing. Analyst 140, 250–257. https://doi.org/10.1039/C4AN01061B

Eilers, P. H., & Boelens, H. F. (2005). Baseline correction with asymmetric least squares smoothing. Leiden University Medical Centre Report, 1(1).

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

## Example on 1 spectrum
i = 2
zX = Matrix(X)[i:i, :]
lb = 1e5 ; p = .001
model = detrend_asls(; lb, p)
fit!(model, zX)
zXc = transf(model, zX)   # = corrected spectrum 
B = zX - zXc              # = estimated baseline
f, ax = plotsp(zX, wl)
lines!(wl, vec(B); color = :blue)
lines!(wl, vec(zXc); color = :black)
f

detrend_lo(; kwargs...)
detrend_lo(X; kwargs...)

Baseline correction of each row of X-data by LOESS regression.

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

Keyword arguments:

• span : Window for neighborhood selection (level of smoothing) for the local fitting, typically in 0, 1.
• degree : Polynomial degree for the local fitting.

De-trend transformation: The function fits a baseline by LOESS regression (function loessr) for each observation and returns the residuals (= signals corrected from the baseline).

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

model = detrend_lo(span = .8)
fit!(model, Xtrain)
Xptrain = transf(model, Xtrain)
Xptest = transf(model, Xtest)
plotsp(Xptrain, wl).f
plotsp(Xptest, wl).f

## Example on 1 spectrum
i = 2
zX = Matrix(X)[i:i, :]
model = detrend_lo(span = .75)
fit!(model, zX)
zXc = transf(model, zX)   # = corrected spectrum 
B = zX - zXc            # = estimated baseline
f, ax = plotsp(zX, wl)
lines!(wl, vec(B); color = :blue)
lines!(wl, vec(zXc); color = :black)
f

detrend_pol(; kwargs...)
detrend_pol(X; kwargs...)

Baseline correction of each row of X-data by polynomial linear regression.

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

Keyword arguments:

• degree : Polynom degree.

De-trend transformation: the function fits a baseline by polynomial regression for each observation and returns the residuals (= signals corrected from the baseline).

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

model = detrend_pol(degree = 2)
fit!(model, Xtrain)
Xptrain = transf(model, Xtrain)
Xptest = transf(model, Xtest)
plotsp(Xptrain, wl).f
plotsp(Xptest, wl).f

## Example on 1 spectrum
i = 2
zX = Matrix(X)[i:i, :]
model = detrend_pol(degree = 1)
fit!(model, zX)
zXc = transf(model, zX)   # = corrected spectrum 
B = zX - zXc            # = estimated baseline
f, ax = plotsp(zX, wl)
lines!(wl, vec(B); color = :blue)
lines!(wl, vec(zXc); color = :black)
f

dfplsr_cg(X, y; kwargs...)

Compute the model complexity (df) of PLSR models with the CGLS algorithm.

• X : X-data (n, p).
• y : Univariate Y-data.

Keyword arguments:

• Same as function cglsr.

The number of degrees of freedom (df) of the PLSR model is returned for 0, 1, ..., nlv LVs.

References

Hansen, P.C., 1998. Rank-Deficient and Discrete Ill-Posed Problems, Mathematical Modeling and Computation. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898719697

Hansen, P.C., 2008. Regularization Tools version 4.0 for Matlab 7.3. Numer Algor 46, 189–194. https://doi.org/10.1007/s11075-007-9136-9

Lesnoff, M., Roger, J.-M., Rutledge, D.N., 2021. Monte Carlo methods for estimating Mallows’s Cp and AIC criteria for PLSR models. Illustration on agronomic spectroscopic NIR data. Journal of Chemometrics n/a, e3369. https://doi.org/10.1002/cem.3369

Examples

## The example below reproduces the numerical illustration
## given by Kramer & Sugiyama 2011 on the Ozone data 
## (Fig. 1, center).
## Function "pls.model" used for df calculations
## in the R package "plsdof" v0.2-9 (Kramer & Braun 2019)
## automatically scales the X matrix before PLS.
## The example scales X for consistency with plsdof.

using Jchemo, JchemoData, JLD2, DataFrames, CairoMakie 
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "ozone.jld2") 
@load db dat
@names dat
X = dat.X
dropmissing!(X) 
zX = rmcol(Matrix(X), 4) 
y = X[:, 4] 
## For consistency with plsdof
xscales = colstd(zX)
zXs = fscale(zX, xscales)
## End

nlv = 12 ; gs = true
res = dfplsr_cg(zXs, y; nlv, gs) ;
res.df 
df_kramer = [1.000000, 3.712373, 6.456417, 11.633565, 
    12.156760, 11.715101, 12.349716,
    12.192682, 13.000000, 13.000000, 
    13.000000, 13.000000, 13.000000]
f, ax = plotgrid(0:nlv, df_kramer; step = 2, xlabel = "Nb. LVs", ylabel = "df")
scatter!(ax, 0:nlv, res.df; color = "red")
ablines!(ax, 1, 1; color = :grey, linestyle = :dot)
f

difmean(X1, X2; normx::Bool = false)

Compute a 1-D detrimental matrix by difference of the column-means of two X-datas.

• X1 : Spectra (n1, p).
• X2 : Spectra (n2, p).

Keyword arguments:

• normx : Boolean. If true, the column-means vectors of X1 and X2 are normed before computing their difference.

The function returns a matrix D (1, p) computed by the difference between two mean-spectra, i.e. the column-means of X1 and X2.

D is assumed to contain the detrimental information that can be removed (by orthogonalization) from X1 and X2 for calibration transfer. For instance, D can be used as input of function eposvd.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/caltransfer.jld2")
@load db dat
@names dat
X1cal = dat.X1cal
X1val = dat.X1val
X2cal = dat.X2cal
X2val = dat.X2val

## The objective is to remove a detrimental 
## information (here, D) from spaces X1 and X2
D = difmean(X1cal, X2cal).D
res = eposvd(D; nlv = 1)
## Corrected Val matrices
X1val_c = X1val * res.M
X2val_c = X2val * res.M

i = 1
f = Figure(size = (800, 300))
ax1 = Axis(f[1, 1])
ax2 = Axis(f[1, 2])
lines!(ax1, X1val[i, :]; label = "x1")
lines!(ax1, X2val[i, :]; label = "x2")
axislegend(ax1, position = :cb, framevisible = false)
lines!(ax2, X1val_c[i, :]; label = "x1_correct")
lines!(ax2, X2val_c[i, :]; label = "x2_correct")
axislegend(ax2, position = :cb, framevisible = false)
f

dkplskdeda(; kwargs...)
dkplskdeda(X, y; kwargs...)
dkplskdeda(X, y, weights::Weight; kwargs...)

DKPLS-KDEDA.

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute. Must be >= 1
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• Keyword arguments of function dmkern (bandwidth definition) can also be specified here.
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation in the PLS computation.

Same as function plskdeda (PLS-KDEDA) except that a direct kernel PLSR (function dkplsr), instead of a PLSR (function plskern), is run on the Y-dummy table.

See function dkplslda for examples.

dkplslda(; kwargs...)
dkplslda(X, y; kwargs...)
dkplslda(X, y, weights::Weight; kwargs...)

DKPLS-LDA.

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute. Must be >= 1.
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

Same as function plslda (PLS-LDA) except that a direct kernel PLSR (function dkplsr), instead of a PLSR (function plskern), is run on the Y-dummy table.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlv = 15
gamma = .1
model = dkplslda(; nlv, gamma) 
#model = dkplslda(; nlv, gamma, prior = :prop) 
#model = dkplsqda(; nlv, gamma, alpha = .5) 
#model = dkplskdeda(; nlv, gamma, a = .5) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

embfitm = fitm.fitm.embfitm ;
@head embfitm.T
@head transf(model, Xtrain)
@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

coef(embfitm)

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xtest; nlv = 1:2).pred

dkplsqda(; kwargs...)
dkplsqda(X, y; kwargs...)
dkplsqda(X, y, weights::Weight; kwargs...)

DKPLS-QDA.

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute. Must be >= 1
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• alpha : Scalar (∈ [0, 1]) defining the continuum between QDA (alpha = 0) and LDA (alpha = 1).
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation in the PLS computation.

Same as function plsqda (PLS-QDA) except that a direct kernel PLSR (function dkplsr), instead of a PLSR (function plskern), is run on the Y-dummy table.

See function dkplslda for examples.

dkplsr(; kwargs...)
dkplsr(X, Y; kwargs...)
dkplsr(X, Y, weights::Weight; kwargs...)
dkplsr!(X::Matrix, Y::Union{Matrix, BitMatrix}, weights::Weight; kwargs...)

Direct kernel partial least squares regression (DKPLSR) (Bennett & Embrechts 2003).

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to consider.
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• scal : Boolean. If true, each column of X and Y is scaled by its uncorrected standard deviation.

The method builds kernel Gram matrices and then runs a usual PLSR algorithm on them. This is faster (but not equivalent) to the "true" KPLSR (Nipals) algorithm (function kplsr) described in Rosipal & Trejo (2001).

References

Bennett, K.P., Embrechts, M.J., 2003. An optimization perspective on kernel partial least squares regression, in: Advances in Learning Theory: Methods, Models and Applications, NATO Science Series III: Computer & Systems Sciences. IOS Press Amsterdam, pp. 227-250.

Rosipal, R., Trejo, L.J., 2001. Kernel Partial Least Squares Regression in Reproducing Kernel Hilbert Space. Journal of Machine Learning Research 2, 97-123.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlv = 20
kern = :krbf ; gamma = 1e-1 ; scal = false
#gamma = 1e-4 ; scal = true
model = dkplsr(; nlv, kern, gamma, scal) ;
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
@head model.fitm.T

coef(model)
coef(model; nlv = 3)

@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

res = predict(model, Xtest)
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",  
    ylabel = "Observed").f  

####### Example of fitting the function sinc(x)
####### described in Rosipal & Trejo 2001 p. 105-106 
x = collect(-10:.2:10) 
x[x .== 0] .= 1e-5
n = length(x)
zy = sin.(abs.(x)) ./ abs.(x) 
y = zy + .2 * randn(n) 
nlv = 2
gamma = 1 / 3
model = dkplsr(; nlv, gamma) ;
fit!(model, x, y)
pred = predict(model, x).pred 
f, ax = scatter(x, y) 
lines!(ax, x, zy, label = "True model")
lines!(ax, x, vec(pred), label = "Fitted model")
axislegend("Method")
f

dkplsrda(; kwargs...)
dkplsrda(X, y; kwargs...)
dkplsrda(X, y, weights::Weight; kwargs...)

Discrimination based on direct kernel partial least squares regression (KPLSR-DA).

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute.
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation.

Same as function plsrda (PLSR-DA) except that a direct kernel PLSR (function dkplsr), instead of a PLSR (function plskern), is run on the Y-dummy table.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlv = 15
kern = :krbf ; gamma = .001 
scal = true
model = dkplsrda(; nlv, kern, gamma, scal) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

@head fitm.fitm.T
@head transf(model, Xtrain)
@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

coef(fitm.fitm)

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xtest; nlv = 1:2).pred

dmkern(; kwargs...)
dmkern(X; kwargs...)

Gaussian kernel density estimation (KDE).

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

Keyword arguments:

• h : Define the bandwith, see examples.
• a : Constant for the Scott's rule (default bandwith), see thereafter.

Estimation of the probability density of X (column space) by non parametric Gaussian kernels.

Data X can be univariate (p = 1) or multivariate (p > 1). In the last case, function dmkern computes a multiplicative kernel such as in Scott & Sain 2005 Eq.19, and the internal bandwidth matrix H is diagonal (see the code).

Note: H in the dmkern code is often noted "H^(1/2)" in the litterature (e.g. Wikipedia).

The default bandwith is computed by:

• h = a * n^(-1 / (p + 4)) * colstd(X)

(a = 1 in Scott & Sain 2005).

References

Scott, D.W., Sain, S.R., 2005. 9 - Multidimensional Density Estimation, in: Rao, C.R., Wegman, E.J., Solka, J.L. (Eds.), Handbook of Statistics, Data Mining and Data Visualization. Elsevier, pp. 229–261. https://doi.org/10.1016/S0169-7161(04)24009-3

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "iris.jld2") 
@load db dat
@names dat
X = dat.X[:, 1:4] 
y = dat.X[:, 5]
n = nro(X)
tab(y) 

nlv = 2
model0 = fda(; nlv)
fit!(model0, X, y)
@head T = transf(model0, X)
n, p = size(T)

#### Probability density in the FDA score space (2D)

model = dmkern()
fit!(model, T) 
@names model.fitm
model.fitm.H
u = [1; 4; 150]
predict(model, T[u, :]).pred

h = .3
model = dmkern(; h)
fit!(model, T) 
model.fitm.H
predict(model, T[u, :]).pred

h = [.3; .1]
model = dmkern(; h)
fit!(model, T) 
model.fitm.H
predict(model, T[u, :]).pred

## Bivariate distribution
npoints = 2^7
nlv = 2
lims = [(minimum(T[:, j]), maximum(T[:, j])) for j = 1:nlv]
x1 = LinRange(lims[1][1], lims[1][2], npoints)
x2 = LinRange(lims[2][1], lims[2][2], npoints)
z = mpar(x1 = x1, x2 = x2)
grid = reduce(hcat, z)
m = nro(grid)
model = dmkern() 
#model = dmkern(a = .5) 
#model = dmkern(h = .3) 
fit!(model, T) 

res = predict(model, grid) ;
pred_grid = vec(res.pred)
f = Figure(size = (600, 400))
ax = Axis(f[1, 1];  title = "Density for FDA scores (Iris)", xlabel = "Score 1", 
    ylabel = "Score 2")
co = contour!(ax, grid[:, 1], grid[:, 2], pred_grid; levels = 10, labels = true)
scatter!(ax, T[:, 1], T[:, 2], color = :red, markersize = 5)
#xlims!(ax, -15, 15) ;ylims!(ax, -15, 15)
f

## Univariate distribution
x = T[:, 1]
model = dmkern() 
#model = dmkern(a = .5) 
#model = dmkern(h = .3) 
fit!(model, x) 
pred = predict(model, x).pred 
f = Figure()
ax = Axis(f[1, 1])
hist!(ax, x; bins = 30, normalization = :pdf)  # area = 1
scatter!(ax, x, vec(pred); color = :red)
f

x = T[:, 1]
npoints = 2^8
lims = [minimum(x), maximum(x)]
#delta = 5 ; lims = [minimum(x) - delta, maximum(x) + delta]
grid = LinRange(lims[1], lims[2], npoints)
model = dmkern() 
#model = dmkern(a = .5) 
#model = dmkern(h = .3) 
fit!(model, x) 
pred_grid = predict(model, grid).pred 
f = Figure()
ax = Axis(f[1, 1])
hist!(ax, x; bins = 30, normalization = :pdf)  # area = 1
lines!(ax, grid, vec(pred_grid); color = :red)
f

dmnorm(; kwargs...)
dmnorm(X; kwargs...)
dmnorm!(X::Matrix; kwargs...)
dmnorm(mu, S; kwargs...)
dmnorm!(mu::Vector, S::Matrix; kwargs...)

Normal probability density estimation.

• X : X-data (n, p) used to estimate the mean mu and the covariance matrix S. If X is not given, mu and S must be provided in kwargs.
• mu : Mean vector of the normal distribution.
• S : Covariance matrix of the Normal distribution.

Keyword arguments:

• simpl : Boolean. If true, the constant term and the determinant in the Normal density formula are set to 1.

Data X can be univariate (p = 1) or multivariate (p > 1). See examples.

When simple = true, the determinant of the covariance matrix (object detS) and the constant (2 * pi)^(-p / 2) (object cst) in the density formula are set to 1. The function returns a pseudo density that resumes to exp(-d / 2), where d is the squared Mahalanobis distance to the center mu. This can for instance be useful when the number of columns (p) of X becomes too large, with the possible consequences that:

• detS tends to 0 or, conversely, to infinity;
• cst tends to 0,

which makes impossible to compute the true density.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "iris.jld2") 
@load db dat
@names dat
X = dat.X[:, 1:4] 
y = dat.X[:, 5]
n = nro(X)
tab(y) 

nlv = 2
model0 = fda(; nlv)
fit!(model0, X, y)
@head T = transf(model0, X)
n, p = size(T)

#### Probability density in the FDA score space (2D)
#### Example of class Setosa 
s = y .== "setosa"
zT = T[s, :]
m = nro(zT)

#### Bivariate distribution
model = dmnorm()
fit!(model, zT)
fitm = model.fitm
@names fitm
fitm.Uinv 
fitm.detS
@head pred = predict(model, zT).pred

## Direct syntax
mu = colmean(zT)
S = covm(zT, mweight(ones(m))) * m / (m - 1) # corrected cov. matrix
fitm = dmnorm(mu, S) ; 
@names fitm
fitm.Uinv
fitm.detS

npoints = 2^7
lims = [(minimum(zT[:, j]), maximum(zT[:, j])) for j = 1:nlv]
x1 = LinRange(lims[1][1], lims[1][2], npoints)
x2 = LinRange(lims[2][1], lims[2][2], npoints)
z = mpar(x1 = x1, x2 = x2)
grid = reduce(hcat, z)
model = dmnorm()
fit!(model, zT)
res = predict(model, grid) ;
pred_grid = vec(res.pred)
f = Figure(size = (600, 400))
ax = Axis(f[1, 1];  title = "Density for FDA scores (Iris - Setosa)", 
    xlabel = "Score 1", ylabel = "Score 2")
co = contour!(ax, grid[:, 1], grid[:, 2], pred_grid; levels = 10, labels = true)
scatter!(ax, T[:, 1], T[:, 2], color = :red, markersize = 5)
scatter!(ax, zT[:, 1], zT[:, 2], color = :blue, markersize = 5)
#xlims!(ax, -12, 12) ;ylims!(ax, -12, 12)
f

#### Univariate distribution
j = 1
x = zT[:, j]
model = dmnorm()
fit!(model, x)
pred = predict(model, x).pred 
f = Figure()
ax = Axis(f[1, 1]; xlabel = string("FDA-score ", j))
hist!(ax, x; bins = 30, normalization = :pdf)  # area = 1
scatter!(ax, x, vec(pred); color = :red)
f

x = zT[:, j]
npoints = 2^8
lims = [minimum(x), maximum(x)]
#delta = 5 ; lims = [minimum(x) - delta, maximum(x) + delta]
grid = LinRange(lims[1], lims[2], npoints)
model = dmnorm()
fit!(model, x)
pred_grid = predict(model, grid).pred 
f = Figure()
ax = Axis(f[1, 1]; xlabel = string("FDA-score ", j))
hist!(ax, x; bins = 30, normalization = :pdf)  # area = 1
lines!(ax, grid, vec(pred_grid); color = :red)
f

dmnormlog(; kwargs...)
dmnormlog(X; kwargs...)
dmnormlog!(X::Matrix; kwargs...)
dmnormlog(mu, S; kwargs...)
dmnormlog!(mu::Vector, S::Matrix; kwargs...)

Logarithm of the normal probability density estimation. * X : X-data (n, p) used to estimate the mean mu and the covariance matrix S. If X is not given, mu and S must be provided in kwargs. * mu : Mean vector of the normal distribution. * S : Covariance matrix of the Normal distribution. Keyword arguments: * simpl : Boolean. If true, the constant term and the determinant in the Normal density formula are set to 1.

See the help page of function dmnorm.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "iris.jld2") 
@load db dat
@names dat
X = dat.X[:, 1:4] 
y = dat.X[:, 5]
n = nro(X)
tab(y) 

## Example of class Setosa 
s = y .== "setosa"
zX = X[s, :]

model = dmnormlog()
fit!(model, zX)
fitm = model.fitm
@names fitm
fitm.Uinv 
fitm.logdetS
@head pred = predict(model, zX).pred

## Consistency with dmnorm
model0 = dmnorm()
fit!(model0, zX)
@head pred0 = predict(model0, zX).pred
@head log.(pred0)

dummy(y)

Compute dummy table from a categorical variable.

• y : A categorical variable.

The output Y (dummy table) is a BitMatrix.

Examples

using Jchemo

y = ["d", "a", "b", "c", "b", "c"]
#y =  rand(1:3, 7)
res = dummy(y)
@names res
res.Y

dupl(X; digits = 3)

Find duplicated rows in a dataset.

• X : A dataset.
• digits : Nb. digits used to round X before checking.

Examples

using Jchemo

X = rand(5, 3)
Z = vcat(X, X[1:3, :], X[1:1, :])
dupl(X)
dupl(Z)

M = hcat(X, fill(missing, 5))
Z = vcat(M, M[1:3, :])
dupl(M)
dupl(Z)

ensure_df(X)

Reshape X to a dataframe if necessary.

ensure_mat(X)

Reshape X to a matrix if necessary.

eposvd(D; nlv = 1)

Compute an orthogonalization matrix for calibration transfer of spectral data.

• D : Data (m, p) containing the detrimental information on which spectra (rows of a matrix X) have to be orthogonalized.

Keyword arguments:

• nlv : Nb. of first loadings vectors of D considered for the orthogonalization.

The objective is to remove some detrimental information (e.g. humidity patterns in signals, multiple spectrometers, etc.) from a X-dataset (n, p). The detrimental information is defined by the main row-directions computed from a matrix D (m, p).

Function eposvd returns two objects:

• V (p, nlv) : The matrix of the nlv first loading vectors of the SVD decomposition (non centered PCA) of D.
• M (p, p) : The orthogonalization matrix, used to orthogonolize a given matrix X to directions contained in V.

Any matrix X can then be corrected from D by:

• X_corrected = X * M.

Matrix D can be built from many methods. For instance, two common methods are:

• EPO (Roger et al. 2003, 2018): D is built from a set of differences between spectra collected under different conditions.
• TOP (Andrew & Fearn 2004): Each row of D is the mean spectrum computed for a given spectrometer instrument.

A particular situation is the following. Assume that D is built from some differences between matrices X1 and X2, and that a bilinear model (e.g. PLSR) is fitted on the data {X1corrected, Y} where X1corrected = X1 * M. To predict new data X2new with the fitted model, there is no need to correct X2new.

References

Andrew, A., Fearn, T., 2004. Transfer by orthogonal projection: making near-infrared calibrations robust to between-instrument variation. Chemometrics and Intelligent Laboratory Systems 72, 51–56. https://doi.org/10.1016/j.chemolab.2004.02.004

Roger, J.-M., Chauchard, F., Bellon-Maurel, V., 2003. EPO-PLS external parameter orthogonalisation of PLS application to temperature-independent measurement of sugar content of intact fruits. Chemometrics and Intelligent Laboratory Systems 66, 191-204. https://doi.org/10.1016/S0169-7439(03)00051-0

Roger, J.-M., Boulet, J.-C., 2018. A review of orthogonal projections for calibration. Journal of Chemometrics 32, e3045. https://doi.org/10.1002/cem.3045

Zeaiter, M., Roger, J.M., Bellon-Maurel, V., 2006. Dynamic orthogonal projection. A new method to maintain the on-line robustness of multivariate calibrations. Application to NIR-based monitoring of wine fermentations. Chemometrics and Intelligent Laboratory Systems, 80, 227–235. https://doi.org/10.1016/j.chemolab.2005.06.011

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/caltransfer.jld2")
@load db dat
@names dat
X1cal = dat.X1cal
X1val = dat.X1val
X2cal = dat.X2cal
X2val = dat.X2val

## The objective is to remove a detrimental 
## information (here, D) from spaces X1 and X2
D = X1cal - X2cal
nlv = 2
res = eposvd(D; nlv)
res.M # orthogonalization matrix
res.V # detrimental directions (columns of matrix V = loadings of D)

## Corrected Val matrices
X1val_c = X1val * res.M
X2val_c = X2val * res.M

i = 1
f = Figure(size = (800, 300))
ax1 = Axis(f[1, 1])
ax2 = Axis(f[1, 2])
lines!(ax1, X1val[i, :]; label = "x1")
lines!(ax1, X2val[i, :]; label = "x2")
axislegend(ax1, position = :cb, framevisible = false)
lines!(ax2, X1val_c[i, :]; label = "x1_correct")
lines!(ax2, X2val_c[i, :]; label = "x2_correct")
axislegend(ax2, position = :cb, framevisible = false)
f

errp(pred, y)

Compute the classification error rate (ERRP).

• pred : Predictions.
• y : Observed data (class membership).

Examples

using Jchemo

Xtrain = rand(10, 5) 
ytrain = rand(["a" ; "b"], 10)
Xtest = rand(4, 5) 
ytest = rand(["a" ; "b"], 4)

model = plsrda(; nlv = 2)
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
errp(pred, ytest)

euclsq(X, Y)

Squared Euclidean distances between the rows of X and Y.

• X : Data (n, p).
• Y : Data (m, p).

For X(n, p) and Y (m, p), the function returns an object (n, m) with:

• i, j = distance between row i of X and row j of Y.

Examples

X = rand(5, 3)
Y = rand(2, 3)

euclsq(X, Y)

euclsq(X[1:1, :], Y[1:1, :])

euclsq(X[:, 1], 4)
euclsq(1, 4)

fcenter(X, v)
fcenter!(X::AbstractMatrix, v)

Center each column of a matrix.

• X : Data (n, p).
• v : Centering vector (p).

examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
xmeans = colmean(X)
fcenter(X, xmeans)

fconcat()

Concatenate horizontaly multiblock X-data.

• Xbl : List of blocks (vector of matrices) of X-data Typically, output of function mblock from (n, p) data.

Examples

using Jchemo
n = 5 ; m = 3 ; p = 9 
X = rand(n, p) 
Xnew = rand(m, p)
listbl = [3:4, 1, [6; 8:9]]
Xbl = mblock(X, listbl) 
Xblnew = mblock(Xnew, listbl) 
@head Xbl[3]

fconcat(Xbl)

fcscale(X, u, v)
fcscale!(X, u, v)

Center and fscale each column of a matrix.

• X : Data (n, p).
• u : Centering vector (p).
• v : Scaling vector (p).

examples

using Jchemo

n, p = 5, 6
X = rand(n, p)
xmeans = colmean(X)
xscales = colstd(X)
fcscale(X, xmeans, xscales)

fda(; kwargs...)
fda(X, y; kwargs...)
fda(X, y, weights; kwargs...)
fda!(X::Matrix, y, weights; kwargs...)

Factorial discriminant analysis (FDA).

• X : X-data (n, p).
• y : y-data (n) (class membership).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of discriminant components.
• lb : Ridge regularization parameter "lambda". Can be used when X has collinearities.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

FDA by eigen factorization of Inverse(W) * B, where W is the "Within"-covariance matrix (pooled over the classes), and B the "Between"-covariance matrix.

The function maximizes the consensus:

• p'Bp / p'Wp

i.e. max p'Bp with constraint p'Wp = 1. Vectors p (columns of V) are the linear discrimant coefficients often referred to as "LD".

If X is ill-conditionned, a ridge regularization can be used:

• If lb > 0, W is replaced by W + lb * I, where I is the Idendity matrix.

In the high-level version of the present functions, the observation weights are automatically defined by the given priors (argument prior): the sub-totals by class of the observation weights are set equal to the prior probabilities. The low-level version (argument weights) allows to implement other choices.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie 
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/iris.jld2") 
@load db dat
@names dat
@head dat.X
X = dat.X[:, 1:4]
y = dat.X[:, 5]
n = nro(X)
ntest = 30
s = samprand(n, ntest) 
Xtrain = X[s.train, :]
ytrain = y[s.train]
Xtest = X[s.test, :]
ytest = y[s.test]
tab(ytrain)
tab(ytest)

nlv = 2
model = fda(; nlv)
#model = fdasvd(; nlv)
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
lev = fitm.lev
nlev = length(lev)
aggsum(fitm.weights.w, ytrain)

@head fitm.T 
@head transf(model, Xtrain)
@head transf(model, Xtest)

## X-loadings matrix
## = coefficients of the linear discriminant function
## = "LD" of function lda of the R package MASS
fitm.V
fitm.V' * fitm.V

## Explained variance computed by weighted PCA 
## of the class centers in transformed scale
summary(model).explvarx

## Projections of the class centers 
## to the score space
ct = fitm.Tcenters 
f, ax = plotxy(fitm.T[:, 1], fitm.T[:, 2], ytrain; ellipse = true, title = "FDA",
    xlabel = "Score-1", ylabel = "Score-2")
scatter!(ax, ct[:, 1], ct[:, 2], marker = :star5, markersize = 15, color = :red)  # see available_marker_symbols()
f

fdasvd(; kwargs...)
fdasvd(X, y, weights; kwargs...)
fdasvd!(X::Matrix, y, weights; kwargs...)

Factorial discriminant analysis (FDA).

• X : X-data (n, p).
• y : y-data (n) (class membership).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of discriminant components.
• lb : Ridge regularization parameter "lambda". Can be used when X has collinearities.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

FDA by a weighted SVD factorization of the matrix of the class centers (after spherical transformaton). The function gives the same results as function fda.

See function fda for details and examples.

fdif(; kwargs...)
fdif(X; kwargs...)

Finite differences (discrete derivates) for each row of X-data.

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

Keyword arguments:

• npoint : Nb. points involved in the window for the finite differences. The range of the window (= nb. intervals of two successive colums) is npoint - 1.

The method reduces the column-dimension:

• (n, p) –> (n, p - npoint + 1).

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

model = fdif(npoint = 2) 
fit!(model, Xtrain)
Xptrain = transf(model, Xtrain)
Xptest = transf(model, Xtest)
plotsp(Xptrain).f
plotsp(Xptest).f

findmax_cla(x)
findmax_cla(x, weights::Weight)

Find the most occurent level in x.

• x : A categorical variable.
• weights : Weights (n) of the observations. Object of type Weight (e.g. generated by function mweight).

If ex-aequos, the function returns the first.

Examples

using Jchemo

x = rand(1:3, 10)
tab(x)
findmax_cla(x)

findmiss(X)

Find rows with missing data in a dataset.

• X : A dataset.

For dataframes, see also DataFrames.completecases and DataFrames.dropmissing.

Examples

using Jchemo

X = rand(5, 4)
zX = hcat(rand(2, 3), fill(missing, 2))
Z = vcat(X, zX)
findmiss(X)
findmiss(Z)

frob(X)
frob(X, weights::Weight)
frob2(X)
frob2(X, weights::Weight)

Frobenius norm of a matrix.

• X : A matrix (n, p).
• weights : Weights (n) of the observations. Object of type Weight (e.g. generated by function mweight).

The Frobenius norm of X is:

• sqrt(tr(X' * X)).

The Frobenius weighted norm is:

• sqrt(tr(X' * D * X)), where D is the diagonal matrix of vector w.

Functions frob2 are the squared versions of frob.

References

@Stevengj, https://discourse.julialang.org/t/interesting-post-about-simd-dot-product-and-cosine-similarity/123282.

fscale(X, v)
fscale!(X::AbstractMatrix, v)

Scale each column of a matrix.

• X : Data (n, p).
• v : Scaling vector (p).

Examples

using Jchemo

X = rand(5, 2) 
fscale(X, colstd(X))

fweight(X, v)
fweight!(X::AbstractMatrix, v)

Weight each row of a matrix.

• X : Data (n, p).
• v : A weighting vector (n).

Examples

using Jchemo, LinearAlgebra

X = rand(5, 2) 
w = rand(5) 
fweight(X, w)
diagm(w) * X

fweight!(X, w)
X

getknn(Xtrain, X; metric = :eucl, k = 1)

Return the k nearest neighbors in Xtrain of each row of the query X.

• Xtrain : Training X-data.
• X : Query X-data.

Keyword arguments:

• metric : Type of distance used for the query. Possible values are :eucl (Euclidean), :mah (Mahalanobis), :sam (spectral angular distance), :cor (correlation distance).
• k : Number of neighbors to return.

The distances (not squared) are also returned.

Spectral angular and correlation distances between two vectors x and y:

• Spectral angular distance (x, y) = acos(x'y / normv(x)normv(y)) / pi
• Correlation distance (x, y) = sqrt((1 - cor(x, y)) / 2)

Both distances are bounded within 0 (y = x) and 1 (y = -x).

Examples

using Jchemo
Xtrain = rand(5, 3)
X = rand(2, 3)
x = X[1:1, :]

k = 3
res = getknn(Xtrain, X; k)
res.ind  # indexes
res.d    # distances

res = getknn(Xtrain, x; k)
res.ind

res = getknn(Xtrain, X; metric = :mah, k)
res.ind

gridcv(model, X, Y; segm, score, pars = nothing, nlv = nothing, lb = nothing, 
    verbose = false)

Cross-validation (CV) of a model over a grid of parameters.

• model : Model to evaluate.
• X : Training X-data (n, p).
• Y : Training Y-data (n, q).

Keyword arguments:

• segm : Segments of observations used for the CV (output of functions segmts, segmkf, etc.).
• score : Function computing the prediction score (e.g. rmsep).
• pars : tuple of named vectors of same length defining the parameter combinations (e.g. output of function mpar).
• verbose : If true, predicting information are printed.
• nlv : Value, or vector of values, of the nb. of latent variables (LVs).
• lb : Value, or vector of values, of the ridge regularization parameter "lambda".

The function is used for grid-search: it computed a prediction score (= error rate) for the specified model over the combinations of parameters defined in pars.

For models based on LV or ridge regularization, using arguments nlv and lb allow faster computations than including these parameters in argument `pars. See the examples.

The function returns two outputs:

• res : mean results
• res_p : results per replication.

Examples

######## Regression

using JLD2, CairoMakie, JchemoData
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
model = savgol(npoint = 21, deriv = 2, degree = 2)
fit!(model, X)
Xp = transf(model, X)
s = year .<= 2012
Xtrain = Xp[s, :]
ytrain = y[s]
Xtest = rmrow(Xp, s)
ytest = rmrow(y, s)
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)

## Replicated K-fold CV 
K = 3 ; rep = 10
segm = segmkf(ntrain, K; rep)
## Replicated test-set validation
#m = Int(round(ntrain / 3)) ; rep = 30
#segm = segmts(ntrain, m; rep)

####-- Plsr
model = plskern()
nlv = 0:30
rescv = gridcv(model, Xtrain, ytrain; segm, score = rmsep, nlv) ;
@names rescv
res = rescv.res 
plotgrid(res.nlv, res.y1; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = plskern(; nlv = res.nlv[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    

## Adding pars 
pars = mpar(scal = [false; true])
rescv = gridcv(model, Xtrain, ytrain; segm,  score = rmsep, pars, nlv) ;
res = rescv.res 
typ = res.scal
plotgrid(res.nlv, res.y1, typ; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = plskern(nlv = res.nlv[u], scal = res.scal[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f   

####-- Rr 
lb = (10).^(-8:.1:3)
model = rr() 
rescv = gridcv(model, Xtrain, ytrain; segm, score = rmsep, lb) ;
res = rescv.res 
loglb = log.(10, res.lb)
plotgrid(loglb, res.y1; step = 2, xlabel = "log(lambda)", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = rr(lb = res.lb[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f     
    
## Adding pars 
pars = mpar(scal = [false; true])
rescv = gridcv(model, Xtrain, ytrain; segm, score = rmsep, pars, lb) ;
res = rescv.res 
loglb = log.(10, res.lb)
typ = string.(res.scal)
plotgrid(loglb, res.y1, typ; step = 2, xlabel = "log(lambda)", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = rr(lb = res.lb[u], scal = res.scal[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f   

####-- Kplsr 
model = kplsr()
nlv = 0:30
gamma = (10).^(-5:1.:5)
pars = mpar(gamma = gamma)
rescv = gridcv(model, Xtrain, ytrain; segm,  score = rmsep, pars, nlv) ;
res = rescv.res 
loggamma = round.(log.(10, res.gamma), digits = 1)
plotgrid(res.nlv, res.y1, loggamma; step = 2, xlabel = "Nb. LVs",  ylabel = "RMSEP", 
    leg_title = "Log(gamma)").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = kplsr(nlv = res.nlv[u], gamma = res.gamma[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f   

####-- Knnr 
nlvdis = [15, 25] ; metric = [:mah]
h = [1, 2.5, 5]
k = [1; 5; 10; 20; 50 ; 100] 
pars = mpar(nlvdis = nlvdis, metric = metric, h = h, k = k)
length(pars[1]) 
model = knnr()
rescv = gridcv(model, Xtrain, ytrain; segm, score = rmsep, pars, verbose = true) ;
res = rescv.res 
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = knnr(nlvdis = res.nlvdis[u], metric = res.metric[u], h = res.h[u], 
    k = res.k[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f   

####-- Lwplsr 
nlvdis = 15 ; metric = [:mah]
h = [1, 2, 5] ; k = [200, 350, 500] 
pars = mpar(nlvdis = nlvdis, metric = metric, h = h, k = k)
length(pars[1]) 
nlv = 0:20
model = lwplsr()
rescv = gridcv(model, Xtrain, ytrain; segm, score = rmsep, pars, nlv, verbose = true) ;
res = rescv.res 
group = string.("h=", res.h, " k=", res.k)
plotgrid(res.nlv, res.y1, group; xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
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
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f   

####-- LwplsrAvg 
nlvdis = 15 ; metric = [:mah]
h = [1, 2, 5] ; k = [200, 350, 500] 
nlv = [0:20, 5:20] 
pars = mpar(nlvdis = nlvdis, metric = metric, h = h, k = k, nlv = nlv)
length(pars[1]) 
model = lwplsravg()
rescv = gridcv(model, Xtrain, ytrain; segm, score = rmsep, pars, verbose = true) ;
res = rescv.res 
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = lwplsravg(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
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f   
    
##---- Mbplsr
listbl = [1:525, 526:1050]
Xbltrain = mblock(Xtrain, listbl)
Xbltest = mblock(Xtest, listbl) 
Xblcal = mblock(Xcal, listbl) 
Xblval = mblock(Xval, listbl) 

model = mbplsr()
bscal = [:none, :frob]
pars = mpar(bscal = bscal) 
nlv = 0:30
rescv = gridcv(model, Xbltrain, ytrain; segm,  score = rmsep, pars, nlv) ;
res = rescv.res 
group = res.bscal 
plotgrid(res.nlv, res.y1, group; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = mbplsr(bscal = res.bscal[u], nlv = res.nlv[u])
fit!(model, Xbltrain, ytrain)
pred = predict(model, Xbltest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f   

######## Discrimination
## The principle is the same as for regression

using JLD2, CairoMakie, JchemoData
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
tab(Y.typ)
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)

## Replicated K-fold CV 
K = 3 ; rep = 10
segm = segmkf(ntrain, K; rep)
## Replicated test-set validation
#m = Int(round(ntrain / 3)) ; rep = 30
#segm = segmts(ntrain, m; rep)

####-- Plslda
model = plslda()
nlv = 1:30
prior = [:unif; :prop]
pars = mpar(prior = prior)
rescv = gridcv(model, Xtrain, ytrain; segm, score = errp, pars, nlv)
res = rescv.res
typ = res.prior
plotgrid(res.nlv, res.y1, typ; step = 2, xlabel = "Nb. LVs", ylabel = "ERR").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = plslda(nlv = res.nlv[u], prior = res.prior[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show errp(pred, ytest)
conf(pred, ytest).pct

gridcv_br(X, Y; segm, algo, score, pars, verbose = false)

Working function for gridcv.

See function gridcv for examples.

gridcv_lb(X, Y; segm, algo, score, pars = nothing, lb, verbose = false)

Working function for gridcv.

Specific and faster than gridcv_br for models using ridge regularization (e.g. RR). Argument pars must not contain nlv.

See function gridcv for examples.

gridcv_lv((X, Y; segm, algo, score, pars = nothing, nlv, verbose = false)

Working function for gridcv.

Specific and faster than gridcv_br for models using latent variables (e.g. PLSR). Argument pars must not contain nlv.

See function gridcv for examples.

gridscore(model, Xtrain, Ytrain, X, Y; score, pars = nothing, nlv = nothing, 
    lb = nothing, verbose = false)

Test-set validation of a model over a grid of parameters.

• model : Model to evaluate.
• Xtrain : Training X-data (n, p).
• Ytrain : Training Y-data (n, q).
• X : Validation X-data (m, p).
• Y : Validation Y-data (m, q).

Keyword arguments:

• score : Function computing the prediction score (e.g. rmsep).
• pars : tuple of named vectors of same length defining the parameter combinations (e.g. output of function mpar).
• verbose : If true, predicting information are printed.
• nlv : Value, or vector of values, of the nb. of latent variables (LVs).
• lb : Value, or vector of values, of the ridge regularization parameter "lambda".

The function is used for grid-search: it computed a prediction score (= error rate) for model model over the combinations of parameters defined in pars. The score is computed over sets {X,Y`}.

For models based on LV or ridge regularization, using arguments nlv and lb allow faster computations than including these parameters in argument `pars. See the examples.

Examples

######## Regression 

using JLD2, CairoMakie, JchemoData
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
model = savgol(npoint = 21, deriv = 2, degree = 2)
fit!(model, X)
Xp = transf(model, X)
s = year .<= 2012
Xtrain = Xp[s, :]
ytrain = y[s]
Xtest = rmrow(Xp, s)
ytest = rmrow(y, s)
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)
## Train ==> Cal + Val 
nval = Int(round(.3 * ntrain))
s = samprand(ntrain, nval)
Xcal = Xtrain[s.train, :]
ycal = ytrain[s.train]
Xval = Xtrain[s.test, :]
yval = ytrain[s.test]

##---- Plsr
model = plskern()
nlv = 0:30
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, nlv)
plotgrid(res.nlv, res.y1; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = plskern(nlv = res.nlv[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    

## Adding pars 
pars = mpar(scal = [false; true])
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, pars, nlv)
typ = res.scal
plotgrid(res.nlv, res.y1, typ; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = plskern(nlv = res.nlv[u], scal = res.scal[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    

##---- Rr 
lb = (10).^(-8:.1:3)
model = rr() 
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, lb)
loglb = log.(10, res.lb)
plotgrid(loglb, res.y1; step = 2, xlabel = "log(lambda)", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = rr(lb = res.lb[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    
    
## Adding pars 
pars = mpar(scal = [false; true])
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, pars, lb)
loglb = log.(10, res.lb)
typ = string.(res.scal)
plotgrid(loglb, res.y1, typ; step = 2, xlabel = "log(lambda)", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = rr(lb = res.lb[u], scal = res.scal[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    

##---- Kplsr 
model = kplsr()
nlv = 0:30
gamma = (10).^(-5:1.:5)
pars = mpar(gamma = gamma)
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, pars, nlv)
loggamma = round.(log.(10, res.gamma), digits = 1)
plotgrid(res.nlv, res.y1, loggamma; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP",
    leg_title = "Log(gamma)").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = kplsr(nlv = res.nlv[u], gamma = res.gamma[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    

##---- Knnr 
nlvdis = [15; 25] ; metric = [:mah]
h = [1, 2.5, 5]
k = [1, 5, 10, 20, 50, 100] 
pars = mpar(nlvdis = nlvdis, metric = metric, h = h, k = k)
length(pars[1]) 
model = knnr()
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, pars, verbose = true)
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = knnr(nlvdis = res.nlvdis[u], metric = res.metric[u], h = res.h[u], 
    k = res.k[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    

##---- Lwplsr 
nlvdis = 15 ; metric = [:mah]
h = [1, 2, 5] ; k = [200, 350, 500] 
pars = mpar(nlvdis = nlvdis, metric = metric, h = h, k = k)
length(pars[1]) 
nlv = 0:20
model = lwplsr()
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, pars, nlv, verbose = true)
group = string.("h=", res.h, " k=", res.k)
plotgrid(res.nlv, res.y1, group; xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
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
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    

##---- LwplsrAvg 
nlvdis = 15 ; metric = [:mah]
h = [1, 2, 5] ; k = [200, 350, 500] 
nlv = [0:20, 5:20] 
pars = mpar(nlvdis = nlvdis, metric = metric, h = h, k = k, nlv = nlv)
length(pars[1]) 
model = lwplsravg()
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, pars, verbose = true)
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = lwplsravg(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
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f   

##---- Mbplsr
listbl = [1:525, 526:1050]
Xbltrain = mblock(Xtrain, listbl)
Xbltest = mblock(Xtest, listbl) 
Xblcal = mblock(Xcal, listbl) 
Xblval = mblock(Xval, listbl) 

model = mbplsr()
bscal = [:none, :frob]
pars = mpar(bscal = bscal) 
nlv = 0:30
res = gridscore(model, Xblcal, ycal, Xblval, yval; score = rmsep, pars, nlv)
group = res.bscal 
plotgrid(res.nlv, res.y1, group; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = mbplsr(bscal = res.bscal[u], nlv = res.nlv[u])
fit!(model, Xbltrain, ytrain)
pred = predict(model, Xbltest).pred
@show rmsep(pred, ytest)
plotxy(vec(pred), ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    
    
######## Discrimination
## The principle is the same as for regression

using JLD2, CairoMakie, JchemoData
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
tab(Y.typ)
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)
## Train ==> Cal + Val 
nval = Int(round(.3 * ntrain))
s = samprand(ntrain, nval)
Xcal = Xtrain[s.train, :]
ycal = ytrain[s.train]
Xval = Xtrain[s.test, :]
yval = ytrain[s.test]

##---- Plslda
model = plslda()
nlv = 1:30
prior = [:unif, :prop]
pars = mpar(prior = prior)
res = gridscore(model, Xcal, ycal, Xval, yval; score = errp, pars, nlv)
typ = res.prior
plotgrid(res.nlv, res.y1, typ; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model = plslda(nlv = res.nlv[u], prior = res.prior[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
@show errp(pred, ytest)
conf(pred, ytest).pct

gridscore(model::Pipeline, Xtrain, Ytrain, X, Y; score, pars = nothing, 
    nlv = nothing, lb = nothing, verbose = false)

Test-set validation of a model pipeline over a grid of parameters.

• model : A pipeline of models to evaluate.
• Xtrain : Training X-data (n, p).
• Ytrain : Training Y-data (n, q).
• X : Validation X-data (m, p).
• Y : Validation Y-data (m, q).

Keyword arguments:

• score : Function computing the prediction score (e.g. rmsep).
• pars : tuple of named vectors of same length defining the parameter combinations (e.g. output of function mpar).
• verbose : If true, predicting information are printed.
• nlv : Value, or vector of values, of the nb. of latent variables (LVs).
• lb : Value, or vector of values, of the ridge regularization parameter "lambda".

In the present version of the function, only the last model of the pipeline (= the final predictor) is validated.

For other details, see function gridscore for simple models.

Examples

using JLD2, CairoMakie, JchemoData
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)
## Building Cal and Val 
## within Train
nval = Int(round(.3 * ntrain))
s = samprand(ntrain, nval)
Xcal = Xtrain[s.train, :]
ycal = ytrain[s.train]
Xval = Xtrain[s.test, :]
yval = ytrain[s.test]

####-- Pipeline Snv :> Savgol :> Plsr
## Only the last model is validated
## model1
model1 = snv()
## model2 
npoint = 11 ; deriv = 2 ; degree = 3
model2 = savgol(; npoint, deriv, degree)
## model3
nlv = 0:30
model3 = plskern()
##
model = pip(model1, model2, model3)
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, nlv) ;
plotgrid(res.nlv, res.y1; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP").f
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model3 = plskern(nlv = res.nlv[u])
model = pip(model1, model2, model3)
fit!(model, Xtrain, ytrain)
res = predict(model, Xtest) ; 
@head res.pred 
rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",
      ylabel = "Observed").f

####-- Pipeline Pca :> Svmr
## Only the last model is validated
## model1
nlv = 15 ; scal = true
model1 = pcasvd(; nlv, scal)
## model2
kern = [:krbf]
gamma = (10).^(-5:1.:5)
cost = (10).^(1:3)
epsilon = [.1, .2, .5]
pars = mpar(kern = kern, gamma = gamma, cost = cost, epsilon = epsilon)
model2 = svmr()
##
model = pip(model1, model2)
res = gridscore(model, Xcal, ycal, Xval, yval; score = rmsep, pars, verbose = true)
u = findall(res.y1 .== minimum(res.y1))[1] 
res[u, :]
model2 = svmr(kern = res.kern[u], gamma = res.gamma[u], cost = res.cost[u], epsilon = res.epsilon[u])
model = pip(model1, model2) 
fit!(model, Xtrain, ytrain)
res = predict(model, Xtest) ; 
@head res.pred 
rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",
      ylabel = "Observed").f

gridscore_br(Xtrain, Ytrain, X, Y; algo, score, pars, 
    verbose = false)

Working function for gridscore.

See function gridscore for examples.

gridscore_lb(Xtrain, Ytrain, X, Y; algo, score, pars = nothing, 
    lb, verbose = false)

Working function for gridscore.

Specific and faster than gridscore_br for models using ridge regularization (e.g. RR). Argument pars must not contain lb.

See function gridscore for examples.

gridscore_lv(Xtrain, Ytrain, X, Y; algo, score, pars = nothing, 
    nlv, verbose = false)

Working function for gridscore.

Specific and faster than gridscore_br for models using latent variables (e.g. PLSR). Argument pars must not contain nlv.

See function gridscore for examples.

@head X

Display the first rows of a dataset.

Examples

using Jchemo

X = rand(100, 5)
@head X

interpl(; kwargs...)
interpl(X; kwargs...)

Sampling spectra by interpolation.

• X : Matrix (n, p) of spectra (rows).

Keyword arguments:

• wl : Values representing the column "names" of X. Must be a numeric vector of length p, or an AbstractRange, with growing values.
• wlfin : Final values (within the range of wl) where to interpolate each spectrum. Must be a numeric vector, or an AbstractRange, with growing values.

The function implements a cubic spline interpolation using package DataInterpolations.jl.

References

Package DAtaInterpolations.jl https://github.com/PumasAI/DataInterpolations.jl https://htmlpreview.github.io/?https://github.com/PumasAI/DataInterpolations.jl/blob/v2.0.0/example/DataInterpolations.html

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

wlfin = range(500, 2400, length = 10)
#wlfin = collect(range(500, 2400, length = 10))
model = interpl(; wl, wlfin)
fit!(model, Xtrain)
Xptrain = transf(model, Xtrain)
Xptest = transf(model, Xtest)
plotsp(Xptrain).f
plotsp(Xptest).f

iqrv(x)

Compute the interquartile interval (IQR) of a vector.

• x : A vector (n).

Examples

x = rand(100)
iqrv(x)

isel!(model, X, Y, wl = 1:nco(X); rep = 1, nint = 5, psamp = .3, score = rmsep)

Interval variable selection.

• model : Model to evaluate.
• X : X-data (n, p).
• Y : Y-data (n, q).
• wl : Optional numeric labels (p, 1) of the X-columns.

Keyword arguments:

• rep : Number of replications of the splitting training/test.
• nint : Nb. intervals.
• psamp : Proportion of data used as test set to compute the score.
• score : Function computing the prediction score.

The principle is as follows:

• Data (X, Y) are splitted randomly to a training and a test set.
• Range 1:p in X is segmented to nint intervals, when possible of equal size.
• The model is fitted on the training set and the score (error rate) on the test set, firtsly accounting for all the p variables (reference) and secondly for each of the nint intervals.
• This process is replicated rep times. Average results are provided in the outputs, as well the results per replication.

References

• Nørgaard, L., Saudland, A., Wagner, J., Nielsen, J.V., Munck, L.,

Engelsen, S.B., 2000. Interval Partial Least-Squares Regression (iPLS): A Comparative Chemometric Study with an Example from Near-Infrared Spectroscopy. Appl Spectrosc 54, 413–419. https://doi.org/10.1366/0003702001949500

Examples

using Jchemo, JchemoData, DataFrames, JLD2, CairoMakie
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "tecator.jld2") 
@load db dat
@names dat
X = dat.X
Y = dat.Y 
wl_str = names(X)
wl = parse.(Float64, wl_str) 
ntot, p = size(X)
typ = Y.typ
namy = names(Y)[1:3]
plotsp(X, wl; xlabel = "Wavelength (nm)", ylabel = "Absorbance").f

s = typ .== "train"
Xtrain = X[s, :]
Ytrain = Y[s, namy]
Xtest = rmrow(X, s)
Ytest = rmrow(Y[:, namy], s)
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)

## Work on the j-th y-variable 
j = 2
nam = namy[j]
ytrain = Ytrain[:, nam]
ytest = Ytest[:, nam]

model = plskern(nlv = 5)
nint = 10
res = isel!(model, Xtrain, ytrain, wl; rep = 30, nint) ;
res.res_rep
res.res0_rep
zres = res.res
zres0 = res.res0
f = Figure(size = (650, 300))
ax = Axis(f[1, 1], xlabel = "Wawelength (nm)", ylabel = "RMSEP_Val",
    xticks = zres.lo)
scatter!(ax, zres.mid, zres.y1; color = (:red, .5))
vlines!(ax, zres.lo; color = :grey, linestyle = :dash, linewidth = 1)
hlines!(ax, zres0.y1, linestyle = :dash)
f

kdeda(; kwargs...)
kdeda(X, y; kwargs...)

Discriminant analysis using non-parametric kernel Gaussian density estimation (KDE-DA).

• X : X-data (n, p).
• y : Univariate class membership (n).

Keyword arguments:

• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• Keyword arguments of function dmkern (bandwidth definition) can also be specified here.

The principle is the same as functions qda except that densities by class are estimated from function dmkern instead of function dmnorm.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/iris.jld2")
@load db dat
@names dat
@head dat.X
X = dat.X[:, 1:4]
y = dat.X[:, 5]
n = nro(X)
ntest = 30
s = samprand(n, ntest)
Xtrain = X[s.train, :]
ytrain = y[s.train]
Xtest = X[s.test, :]
ytest = y[s.test]
ntrain = n - ntest
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

prior = :unif
#prior = :prop
model = kdeda(; prior)
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

model = kdeda(; prior, a = .5) 
#model = kdeda(; prior, h = .1) 
fit!(model, Xtrain, ytrain)
model.fitm.fitm[1].H

knnda(; kwargs...)
knnda(X, y; kwargs...)

k-Nearest-Neighbours weighted discrimination (KNN-DA).

• X : X-data (n, p).
• y : Univariate class membership (n).

Keyword arguments:

• metric : Type of dissimilarity used to select the neighbors and to compute the weights. Possible values are: :eucl (Euclidean distance), :mah (Mahalanobis 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.
• scal : Boolean. If true, each column of the global X is scaled by its uncorrected standard deviation before the distance and weight computations.

This function has the same principle as function knnr except that a discrimination replaces the regression. A weighted vote is done over the neighborhood, and the prediction corresponds to the most frequent class.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

metric = :eucl
h = 2 ; k = 10
model = knnda(; metric, h, k) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

res = predict(model, Xtest) ; 
@names res 
res.listnn
res.listd
res.listw
@head res.pred
@show errp(res.pred, ytest)
conf(res.pred, ytest).cnt

## With dimension reduction
model1 = pcasvd(; nlv = 15)
metric = :mah ; h = 1 ; k = 3 
model2 = knnda(; metric, h, k) 
model = pip(model1, model2)
fit!(model, Xtrain, ytrain)
@head pred = predict(model, Xtest).pred 
errp(pred, ytest)

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

k-Nearest-Neighbours weighted regression (KNNR).

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

Keyword arguments:

• metric : Type of dissimilarity used to select the neighbors and to compute the weights. Possible values are: :eucl (Euclidean distance), :mah (Mahalanobis 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.
• scal : Boolean. If true, each column of the global X is scaled by its uncorrected standard deviation before the distance and weight computations.

The general principle of this function is as follows (many other variants of kNNR pipelines can be built): a) For each new observation to predict, the prediction is the weighted mean of y over a selected neighborhood (in X) of size k. b) Within the selected neighborhood, the weights are defined from the dissimilarities between the new observation and the neighborhood, and are computed from function 'winvs'.

In general, for X-data with high dimensions, using the Mahalanobis distance requires a preliminary dimensionality reduction (see examples).

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

h = 1 ; k = 3 
model = knnr(; h, k) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
dump(model.fitm.par)
res = predict(model, Xtest) ; 
@names res 
res.listnn
res.listd
res.listw
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    

## With dimension reduction
model1 = pcasvd(nlv = 15)
metric = :eucl ; h = 1 ; k = 3 
model2 = knnr(; metric, h, k) 
model = pip(model1, model2)
fit!(model, Xtrain, ytrain)
res = predict(model, Xtest) ; 
@head res.pred
@show rmsep(res.pred, ytest)

####### Example of fitting the function sinc(x)
####### described in Rosipal & Trejo 2001 p. 105-106 
x = collect(-10:.2:10) 
x[x .== 0] .= 1e-5
n = length(x)
zy = sin.(abs.(x)) ./ abs.(x) 
y = zy + .2 * randn(n) 
model = knnr(k = 15, h = 5) 
fit!(model, x, y)
pred = predict(model, x).pred 
f, ax = scatter(x, y) 
lines!(ax, x, zy, label = "True model")
lines!(ax, x, vec(pred), label = "Fitted model")
axislegend("Method")
f

kpca(; kwargs...)
kpca(X; kwargs...)
kpca(X, weights::Weight; kwargs...)

Kernel PCA (Scholkopf et al. 1997, Scholkopf & Smola 2002, Tipping 2001).

• X : X-data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. principal components (PCs) to consider.
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

The method is implemented by SVD factorization of the weighted Gram matrix:

• D^(1/2) * Phi(X) * Phi(X)' * D^(1/2)

where X is the cenetred matrix and D is a diagonal matrix of weights (weights.w) of the observations (rows of X).

References

Scholkopf, B., Smola, A., Müller, K.-R., 1997. Kernel principal component analysis, in: Gerstner, W., Germond, A., Hasler, M., Nicoud, J.-D. (Eds.), Artificial Neural Networks, ICANN 97, Lecture Notes in Computer Science. Springer, Berlin, Heidelberg, pp. 583-588. https://doi.org/10.1007/BFb0020217

Scholkopf, B., Smola, A.J., 2002. Learning with kernels: support vector machines, regularization, optimization, and beyond, Adaptive computation and machine learning. MIT Press, Cambridge, Mass.

Tipping, M.E., 2001. Sparse kernel principal component analysis. Advances in neural information processing systems, MIT Press. http://papers.nips.cc/paper/1791-sparse-kernel-principal-component-analysis.pdf

Examples

using Jchemo, JchemoData, JLD2 
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/iris.jld2") 
@load db dat
@names dat
@head dat.X
X = dat.X[:, 1:4]
n = nro(X)
ntest = 30
s = samprand(n, ntest) 
Xtrain = X[s.train, :]
Xtest = X[s.test, :]

nlv = 3
kern = :krbf ; gamma = 1e-4
model = kpca(; nlv, kern, gamma) ;
fit!(model, Xtrain)
@names model.fitm
@head T = model.fitm.T
T' * T
model.fitm.V' * model.fitm.V

@head Ttest = transf(model, Xtest)

res = summary(model) ;
@names res
res.explvarx

kplskdeda(; kwargs...)
kplskdeda(X, y; kwargs...)
kplskdeda(X, y, weights::Weight; kwargs...)

KPLS-KDEDA.

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute. Must be >= 1
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• Keyword arguments of function dmkern (bandwidth definition) can also be specified here.
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation in the PLS computation.

Same as function plskdeda (PLS-KDEDA) except that a kernel PLSR (function kplsr), instead of a PLSR (function plskern), is run on the Y-dummy table.

See function kplslda for examples.

kplslda(; kwargs...)
kplslda(X, y; kwargs...)
kplslda(X, y, weights::Weight; kwargs...)

KPLS-LDA.

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute. Must be >= 1
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation in the PLS computation.

Same as function plslda (PLS-LDA) except that a kernel PLSR (function kplsr), instead of a PLSR (function plskern), is run on the Y-dummy table.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlv = 15
gamma = .1
model = kplslda(; nlv, gamma) 
#model = kplslda(; nlv, gamma, prior = :prop) 
#model = kplsqda(; nlv, gamma, alpha = .5) 
#model = kplskdeda(; nlv, gamma, a = .5) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

embfitm = fitm.fitm.embfitm ;
@head embfitm.T
@head transf(model, Xtrain)
@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

coef(embfitm)

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xtest; nlv = 1:2).pred

kplsqda(; kwargs...)
kplsqda(X, y; kwargs...)
kplsqda(X, y, weights::Weight; kwargs...)

KPLS-QDA.

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute. Must be >= 1
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• alpha : Scalar (∈ [0, 1]) defining the continuum between QDA (alpha = 0) and LDA (alpha = 1).
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation in the PLS computation.

Same as function plsqda (PLS-QDA) except that a kernel PLSR (function kplsr), instead of a PLSR (function plskern), is run on the Y-dummy table.

See function kplslda for examples.

kplsr(; kwargs...)
kplsr(X, Y; kwargs...)
kplsr(X, Y, weights::Weight; kwargs...)
kplsr!(X::Matrix, Y::Union{Matrix, BitMatrix}, weights::Weight; kwargs...)

Kernel partial least squares regression (KPLSR) implemented with a Nipals algorithm (Rosipal & Trejo, 2001).

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to consider.
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• scal : Boolean. If true, each column of X and Y is scaled by its uncorrected standard deviation.

This algorithm becomes slow for n > 1000. Use function dkplsr instead.

References

Rosipal, R., Trejo, L.J., 2001. Kernel Partial Least Squares Regression in Reproducing Kernel Hilbert Space. Journal of Machine Learning Research 2, 97-123.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlv = 20
kern = :krbf ; gamma = 1e-1
model = kplsr(; nlv, kern, gamma) ;
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
@head model.fitm.T

coef(model)
coef(model; nlv = 3)

@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

res = predict(model, Xtest)
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
    ylabel = "Observed").f    

####### Example of fitting the function sinc(x)
####### described in Rosipal & Trejo 2001 p. 105-106 
x = collect(-10:.2:10) 
x[x .== 0] .= 1e-5
n = length(x)
zy = sin.(abs.(x)) ./ abs.(x) 
y = zy + .2 * randn(n) 
nlv = 2
kern = :krbf ; gamma = 1 / 3
model = kplsr(; nlv, kern, gamma) 
fit!(model, x, y)
pred = predict(model, x).pred 
f, ax = scatter(x, y) 
lines!(ax, x, zy, label = "True model")
lines!(ax, x, vec(pred), label = "Fitted model")
axislegend("Method")
f

kplsrda(; kwargs...)
kplsrda(X, y; kwargs...)
kplsrda(X, y, weights::Weight; kwargs...)

Discrimination based on kernel partial least squares regression (KPLSR-DA).

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute.
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation.

Same as function plsrda (PLSR-DA) except that a kernel PLSR (function kplsr), instead of a PLSR (function plskern), is run on the Y-dummy table.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlv = 15
kern = :krbf ; gamma = .001 
scal = true
model = kplsrda(; nlv, kern, gamma, scal) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

@head fitm.fitm.T
@head transf(model, Xtrain)
@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

coef(fitm.fitm)

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xtest; nlv = 1:2).pred

kpol(X, Y; kwargs...)

Compute a polynomial kernel Gram matrix.

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

Keyword arguments:

• gamma : Scale of the polynom.
• coef0 : Offset of the polynom.
• degree : Degree of the polynom.

Given matrices X and Y of sizes (n, p) and (m, p), respectively, the function returns the (n, m) Gram matrix:

• K(X, Y) = Phi(X) * Phi(Y)'.

The polynomial kernel between two vectors x and y is computed by (gamma * (x' * y) + coef0)^degree.

References

Scholkopf, B., Smola, A.J., 2002. Learning with kernels: support vector machines, regularization, optimization, and beyond, Adaptive computation and machine learning. MIT Press, Cambridge, Mass.

Examples

using Jchemo
X = rand(5, 3)
Y = rand(2, 3)
kpol(X, Y; gamma = .1, coef0 = 10, degree = 3)

krbf(X, Y; kwargs...)

Compute a Radial-Basis-Function (RBF) kernel Gram matrix.

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

Keyword arguments:

• gamma : Scale parameter.

Given matrices X and Y of sizes (n, p) and (m, p), respectively, the function returns the (n, m) Gram matrix:

• K(X, Y) = Phi(X) * Phi(Y)'.

The RBF kernel between two vectors x and y is computed by exp(-gamma * ||x - y||^2).

References

Scholkopf, B., Smola, A.J., 2002. Learning with kernels: support vector machines, regularization, optimization, and beyond, Adaptive computation and machine learning. MIT Press, Cambridge, Mass.

Examples

using Jchemo
X = rand(5, 3)
Y = rand(2, 3)
krbf(X, Y; gamma = .1)

krr(; kwargs...)
krr(X, Y; kwargs...)
krr(X, Y, weights::Weight; kwargs...)
krr!(X::Matrix, Y::Union{Matrix, BitMatrix}, weights::Weight; kwargs...)

Kernel ridge regression (KRR) implemented by SVD factorization.

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• lb : Ridge regularization parameter "lambda".
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• scal : Boolean. If true, each column of `X is scaled by its uncorrected standard deviation.

KRR is also referred to as least squared SVM regression (LS-SVMR). The method is close to the particular case of SVM regression where there is no marge excluding the observations (epsilon coefficient set to zero). The difference is that a L2-norm optimization is done, instead of L1 in SVM.

References

Bennett, K.V., Embrechts, M.J., 2003. An optimization perspective on kernel partial least squares regression, in: Advances in Learning Theory: Methods, Models and Applications, NATO Science Series III: Computer & Systems Sciences. IOS Press Amsterdam, pp. 227-250.

Cawley, G.C., Talbot, N.L.C., 2002. Reduced Rank Kernel Ridge Regression. Neural Processing Letters 16, 293-302. https://doi.org/10.1023/A:1021798002258

Krell, M.M., 2018. Generalizing, Decoding, and Optimizing Support Vector Machine Classification. arXiv:1801.04929.

Saunders, C., Gammerman, A., Vovk, V., 1998. Ridge Regression Learning Algorithm in Dual Variables, in: In Proceedings of the 15th International Conference on Machine Learning. Morgan Kaufitmann, pp. 515-521.

Suykens, J.A.K., Lukas, L., Vandewalle, J., 2000. Sparse approximation using least squares support vector machines. 2000 IEEE International Symposium on Circuits and Systems. Emerging Technologies for the 21st Century. Proceedings (IEEE Cat No.00CH36353). https://doi.org/10.1109/ISCAS.2000.856439

Welling, M., n.d. Kernel ridge regression. Department of Computer Science, University of Toronto, Toronto, Canada. https://www.ics.uci.edu/~welling/classnotes/papers_class/Kernel-Ridge.pdf

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

lb = 1e-3
kern = :krbf ; gamma = 1e-1
model = krr(; lb, kern, gamma) ;
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm

coef(model)

res = predict(model, Xtest)
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction", 
   ylabel = "Observed").f    

coef(model; lb = 1e-1)
res = predict(model, Xtest; lb = [.1 ; .01])
@head res.pred[1]
@head res.pred[2]

lb = 1e-3
kern = :kpol ; degree = 1
model = krr(; lb, kern, degree) 
fit!(model, Xtrain, ytrain)
res = predict(model, Xtest)
rmsep(res.pred, ytest)

####### Example of fitting the function sinc(x)
####### described in Rosipal & Trejo 2001 p. 105-106 
x = collect(-10:.2:10) 
x[x .== 0] .= 1e-5
n = length(x)
zy = sin.(abs.(x)) ./ abs.(x) 
y = zy + .2 * randn(n) 
lb = 1e-1
kern = :krbf ; gamma = 1 / 3
model = krr(; lb, kern, gamma) 
fit!(model, x, y)
pred = predict(model, x).pred 
f, ax = scatter(x, y) 
lines!(ax, x, zy, label = "True model")
lines!(ax, x, vec(pred), label = "Fitted model")
axislegend("Method")
f

krrda(; kwargs...)
krrda(X, y; kwargs...)
krrda(X, y, weights::Weight; kwargs...)

Discrimination based on kernel ridge regression (KRR-DA).

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• lb : Ridge regularization parameter "lambda".
• kern : Type of kernel used to compute the Gram matrices. Possible values are: :krbf, :kpol. See respective functions krbf and kpol for their keyword arguments.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

Same as function rrda (RR-DA) except that a kernel RR (function krr), instead of a RR (function rr), is run on the Y-dummy table.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

lb = 1e-5
kern = :krbf ; gamma = .001 
scal = true
model = krrda(; lb, kern, gamma, scal) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

coef(fitm.fitm)

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xtest; lb = [.1, .001]).pred

lda(; kwargs...)
lda(; kwargs...)
lda(X, y; kwargs...)
lda(X, y, weights::Weight; kwargs...)

Linear discriminant analysis (LDA).

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).

In the high-level version of the present functions, the observation weights are automatically defined by the given priors (argument prior): the sub-totals by class of the observation weights are set equal to the prior probabilities. The low-level version (argument weights) allows to implement other choices.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/iris.jld2")
@load db dat
@names dat
@head dat.X
X = dat.X[:, 1:4]
y = dat.X[:, 5]
n = nro(X)
ntest = 30
s = samprand(n, ntest)
Xtrain = X[s.train, :]
ytrain = y[s.train]
Xtest = X[s.test, :]
ytest = y[s.test]
ntrain = n - ntest
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

model = lda()
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni
aggsum(fitm.weights.w, ytrain)

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

list(Q, n::Integer)

Create a Vector {Q}(undef, n).

isassigned(object, i) can be used to check if cell i is empty.

Examples

using Jchemo

list(Float64, 5)
list(Array{Float64}, 5)
list(Matrix{Int}, 5)

list(n::Integer)

Create a Vector{Any}(nothing, n).

isnothing(object, i) can be used to check if cell i is empty.

Examples

using Jchemo

list(5)

locw(Xtrain, Ytrain, X; listnn, listw = nothing, algo, verbose = false, kwargs...)

Compute predictions for a given kNN model.

• Xtrain : Training X-data.
• Ytrain : Training Y-data.
• X : X-data (m observations) to predict.

Keyword arguments:

• listnn : List (vector) of m vectors of indexes.
• listw : List (vector) of m vectors of weights.
• algo : Function computing the model on the m neighborhoods.
• verbose : Boolean. If true, predicting information are printed.
• kwargs : Keywords arguments to pass in function algo. Each argument must have length = 1 (not be a collection).

Each component i of listnn and listw contains the indexes and weights, respectively, of the nearest neighbors of x_i in Xtrain. The sizes of the neighborhood for i = 1,...,m can be different.

locwlv(Xtrain, Ytrain, X; listnn, listw = nothing, algo, nlv, verbose = true, kwargs...)

Compute predictions for a given kNN model.

• Xtrain : Training X-data.
• Ytrain : Training Y-data.
• X : X-data (m observations) to predict.

Keyword arguments:

• listnn : List (vector) of m vectors of indexes.
• listw : List (vector) of m vectors of weights.
• algo : Function computing the model on the m neighborhoods.
• nlv : Nb. or collection of nb. of latent variables (LVs).
• verbose : Boolean. If true, predicting information are printed.
• kwargs : Keywords arguments to pass in function algo. Each argument must have length = 1 (not be a collection).

Same as locw but specific and much faster for LV-based models (e.g. PLSR).

loessr(; kwargs...)
loessr(X, y; kwargs...)

Compute a locally weighted regression model (LOESS).

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

Keyword arguments:

• span : Window for neighborhood selection (level of smoothing) for the local fitting, typically in 0, 1.
• degree : Polynomial degree for the local fitting.
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

The function fits a LOESS model using package `Loess.jl'.

Smaller values of span result in smaller local context in fitting (less smoothing).

References

https://github.com/JuliaStats/Loess.jl

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

Cleveland, W. S., & Grosse, E. (1991). Computational methods for local regression. Statistics and computing, 1(1), 47-62. DOI: 10.1007/BF01890836

Examples

using Jchemo, CairoMakie

####### Example of fitting the function sinc(x)
####### described in Rosipal & Trejo 2001 p. 105-106 
x = collect(-10:.2:10) 
x[x .== 0] .= 1e-5
n = length(x)
zy = sin.(abs.(x)) ./ abs.(x) 
y = zy + .2 * randn(n) 
model = loessr(span = 1 / 3) 
fit!(model, x, y)
pred = predict(model, x).pred 
f = Figure(size = (700, 300))
ax = Axis(f[1, 1], xlabel = "x", ylabel = "y")
scatter!(x, y) 
lines!(ax, x, zy, label = "True model")
lines!(ax, x, vec(pred); label = "Loess")
f[1, 2] = Legend(f, ax, framevisible = false)
f

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

k-Nearest-Neighbours locally weighted multiple linear regression (kNN-LWMLR).

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

Keyword arguments:

• metric : Type of dissimilarity used to select the neighbors and to compute the weights. Possible values are: :eucl (Euclidean distance), :mah (Mahalanobis 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.
• scal : Boolean. If true, each column of the global X is scaled by its uncorrected standard deviation before the distance and weight computations.
• verbose : Boolean. If true, predicting information are printed.

This is the same principle as function lwplsr except that MLR models are fitted on the neighborhoods, instead of PLSR models. The neighborhoods are computed directly on X (there is no preliminary dimension reduction).

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlv = 20
model0 = pcasvd(; nlv) ;
fit!(model0, Xtrain) 
@head Ttrain = model0.fitm.T 
@head Ttest = transf(model0, Xtest)

metric = :eucl 
h = 2 ; k = 100 
model = lwmlr(; metric, h, k) 
fit!(model, Ttrain, ytrain)
@names model
@names model.fitm
dump(model.fitm.par)

res = predict(model, Ttest) ; 
@names res 
res.listnn
res.listd
res.listw
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",  
    ylabel = "Observed").f    

####### Example of fitting the function sinc(x)
####### described in Rosipal & Trejo 2001 p. 105-106 
x = collect(-10:.2:10) 
x[x .== 0] .= 1e-5
n = length(x)
zy = sin.(abs.(x)) ./ abs.(x) 
y = zy + .2 * randn(n) 
model = lwmlr(metric = :eucl, h = 1.5, k = 20) ;
fit!(model, x, y)
pred = predict(model, x).pred 
f, ax = scatter(x, y) 
lines!(ax, x, zy, label = "True model")
lines!(ax, x, vec(pred), label = "Fitted model")
axislegend("Method")
f

lwmlrda(; kwargs...)
lwmlrda(X, y; kwargs...)

k-Nearest-Neighbours locally weighted MLR-based discrimination (kNN-LWMLR-DA).

• X : X-data (n, p).
• y : Univariate class membership (n).

Keyword arguments:

• metric : Type of dissimilarity used to select the neighbors and to compute the weights. Possible values are: :eucl (Euclidean distance), :mah (Mahalanobis 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.
• scal : Boolean. If true, each column of the global X is scaled by its uncorrected standard deviation before the distance and weight computations.
• verbose : Boolean. If true, predicting information are printed.

This is the same principle as function lwmlr except that MLR-DA models, instead of MLR models, are fitted on the neighborhoods.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/iris.jld2")
@load db dat
@names dat
@head dat.X
X = dat.X[:, 1:4]
y = dat.X[:, 5]
n = nro(X)
ntest = 30
s = samprand(n, ntest)
Xtrain = X[s.train, :]
ytrain = y[s.train]
Xtest = X[s.test, :]
ytest = y[s.test]
ntrain = n - ntest
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

metric = :mah
h = 2 ; k = 10
model = lwmlrda(; metric, h, k) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

res = predict(model, Xtest) ; 
@names res 
res.listnn
res.listd
res.listw
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

lwplslda(; kwargs...)
lwplslda(X, y; kwargs...)

kNN-LWPLS-LDA.

• X : X-data (n, p).
• y : Univariate class membership (n).

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. Possible values are: :eucl (Euclidean distance), :mah (Mahalanobis 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.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional).
• scal : Boolean. If true, (a) each column of the global X (and of the global Ydummy 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 Ydummy scaling is also done within each neighborhood (local level) for the weighted PLS.
• verbose : Boolean. If true, predicting information are printed.

This is the same principle as function lwplsr except that a PLS-LDA model, instead of a PLSR model, is fitted on each neighborhoods.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlvdis = 25 ; metric = :mah
h = 2 ; k = 200
nlv = 10
model = lwplslda(; nlvdis, metric, h, k, nlv, prior = :prop) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

res = predict(model, Xtest) ; 
@names res 
res.listnn
res.listd
res.listw
@head res.pred
@show errp(res.pred, ytest)
conf(res.pred, ytest).cnt

lwplsqda(; kwargs...)
lwplsqda(X, y; kwargs...)

kNN-LWPLS-QDA.

• X : X-data (n, p).
• y : Univariate class membership (n).

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. Possible values are: :eucl (Euclidean distance), :mah (Mahalanobis 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.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional).
• alpha : Scalar (∈ [0, 1]) defining the continuum between QDA (alpha = 0) and LDA (alpha = 1).
• scal : Boolean. If true, (a) each column of the global X (and of the global Ydummy 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 Ydummy scaling is also done within each neighborhood (local level) for the weighted PLS.
• verbose : Boolean. If true, predicting information are printed.

This is the same principle as function lwplsr except that a PLS-QDA model, instead of a PLSR model, is fitted on each neighborhoods.

• Warning: The present version of this function can suffer from stops due to non positive definite matrices when doing QDA on neighborhoods. This is due to that some classes within the neighborhood can have very few observations. It is recommended to select a sufficiantly large number of neighbors or/and to use a regularized QDA (alpha > 0).

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlvdis = 25 ; metric = :mah
h = 2 ; k = 200
nlv = 10
model = lwplsqda(; nlvdis, metric, h, k, nlv, prior = :prop, alpha = .5) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

res = predict(model, Xtest) ; 
@names res 
res.listnn
res.listd
res.listw
@head res.pred
@show errp(res.pred, ytest)
conf(res.pred, ytest).cnt

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. Possible values are: :eucl (Euclidean distance), :mah (Mahalanobis 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.

Function lwplsr fits kNN-LWPLSR models such as in Lesnoff et al. 2020. The general principle of the pipeline is as follows (many other variants of pipelines can be built):

LWPLSR is a particular case of weighted PLSR (WPLSR) (e.g. Schaal et al. 2002). In WPLSR, a priori weights, different from the usual 1/n (standard PLSR), are given to the n training observations. These weights are used for calculating (i) the scores and loadings of the WPLS and (ii) the regression model that fits (by weighted least squares) the Y-response(s) to the WPLS scores. The specificity of LWPLSR (compared to WPLSR) is that the weights are computed from dissimilarities (e.g. distances) between the new observation to predict and the training observations ("L" in LWPLSR comes from "localized"). Note that in LWPLSR the weights and therefore the fitted WPLSR model change for each new observation to predict.

In the original LWPLSR, all the n training observations are used for each observation to predict (e.g. Sicard & Sabatier 2006, Kim et al 2011). This can be very time consuming when n is large. A faster (and often more efficient) strategy is to preliminary select, in the training set, a number of k nearest neighbors to the observation to predict (= "weighting 1") and then to apply LWPLSR only to this pre-selected neighborhood (= "weighting 2"). T his strategy corresponds to a kNN-LWPLSR and is the one implemented in function lwplsr.

In lwplsr, the dissimilarities used for weightings 1 and 2 are computed from the raw X-data, or after a dimension reduction, depending on argument nlvdis. In the last case, global PLS2 scores (LVs) are computed from {X, Y} and the dissimilarities are computed over these scores.

In general, for high dimensional X-data, using the Mahalanobis distance requires preliminary dimensionality reduction of the data. In function knnr', the preliminary reduction (argumentnlvdis) is done by PLS on {X,Y`}.

References

Kim, S., Kano, M., Nakagawa, H., Hasebe, S., 2011. Estimation of active pharmaceutical ingredients content using locally weighted partial least squares and statistical wavelength selection. Int. J. Pharm., 421, 269-274.

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, e3209. https://doi.org/10.1002/cem.3209

Schaal, S., Atkeson, C., Vijayamakumar, S. 2002. Scalable techniques from nonparametric statistics for the real time robot learning. Applied Intell., 17, 49-60.

Sicard, E. Sabatier, R., 2006. Theoretical framework for local PLS1 regression and application to a rainfall dataset. Comput. Stat. Data Anal., 51, 1393-1410.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlvdis = 15 ; metric = :mah 
h = 1 ; k = 500 ; nlv = 10
model = lwplsr(; nlvdis, metric, h, k, nlv) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm

res = predict(model, Xtest) ; 
@names res 
res.listnn
res.listd
res.listw
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",  
    ylabel = "Observed").f    

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

Averaging kNN-LWPLSR models with different numbers of latent variables (kNN-LWPLSR-AVG).

• 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. Possible values are: :eucl (Euclidean distance), :mah (Mahalanobis 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 : A range of nb. of latent variables (LVs) to compute 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.

Ensemblist method where the predictions are computed by averaging the predictions of a set of models built with different numbers of LVs, such as in Lesnoff 2023. On each neighborhood, a PLSR-averaging (Lesnoff et al.

1. is done instead of a PLSR.

For instance, if argument nlv is set to nlv = 5:10, the prediction for a new observation is the simple average of the predictions returned by the models with 5 LVs, 6 LVs, ... 10 LVs, respectively.

References

Lesnoff, M., Andueza, D., Barotin, C., Barre, V., Bonnal, L., Fernández Pierna, J.A., Picard, F., Vermeulen, V., Roger, J.-M., 2022. Averaging and Stacking Partial Least Squares Regression Models to Predict the Chemical Compositions and the Nutritive Values of Forages from Spectral Near Infrared Data. Applied Sciences 12, 7850. https://doi.org/10.3390/app12157850

M. Lesnoff, 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 (2023) 105031. https://doi.org/10.1016/j.chemolab.2023.105031.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlvdis = 5 ; metric = :mah 
h = 1 ; k = 200 ; nlv = 4:20
model = lwplsravg(; nlvdis, metric, h, k, nlv) ;
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm

res = predict(model, Xtest) ; 
@names res 
res.listnn
res.listd
res.listw
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",  
    ylabel = "Observed").f  

lwplsrda(; kwargs...)
lwplsrda(X, y; kwargs...)

kNN-LWPLSR-DA.

• X : X-data (n, p).
• y : Univariate class membership (n).

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. Possible values are: :eucl (Euclidean distance), :mah (Mahalanobis 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.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional).
• scal : Boolean. If true, (a) each column of the global X (and of the global Ydummy 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 Ydummy scaling is also done within each neighborhood (local level) for the weighted PLS.
• verbose : Boolean. If true, predicting information are printed.

This is the same principle as function lwplsr except that PLSR-DA models, instead of PLSR models, are fitted on the neighborhoods.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlvdis = 25 ; metric = :mah
h = 2 ; k = 200
nlv = 10
model = lwplsrda(; nlvdis, metric, h, k, nlv) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

res = predict(model, Xtest) ; 
@names res 
res.listnn
res.listd
res.listw
@head res.pred
@show errp(res.pred, ytest)
conf(res.pred, ytest).cnt

madv(x)

Compute the median absolute deviation (MAD) of a vector.

• x : A vector (n).

This is the MAD adjusted by a factor (1.4826) for asymptotically normal consistency.

Examples

using Jchemo

x = rand(100)
madv(x)

mahsq(X, Y)
mahsq(X, Y, Sinv)

Squared Mahalanobis distances between the rows of X and Y.

• X : Data (n, p).
• Y : Data (m, p).
• Sinv : Inverse of a covariance matrix S. If not given, S is computed as the uncorrected covariance matrix of X.

When X and Y are (n, p) and (m, p), repectively, it returns an object (n, m) with:

• i, j = distance between row i of X and row j of Y.

Examples

using StatsBase 

X = rand(5, 3)
Y = rand(2, 3)

mahsq(X, Y)

S = cov(X, corrected = false)
Sinv = inv(S)
mahsq(X, Y, Sinv)
mahsq(X[1:1, :], Y[1:1, :], Sinv)

mahsq(X[:, 1], 4)
mahsq(1, 4, 2.1)

mahsqchol(X, Y)
mahsqchol(X, Y, Uinv)

Compute the squared Mahalanobis distances (with a Cholesky factorization) between the observations (rows) of X and Y.

• X : Data (n, p).
• Y : Data (m, p).
• Uinv : Inverse of the upper matrix of a Cholesky factorization of a covariance matrix S. If not given, the factorization is done on S, the uncorrected covariance matrix of X.

When X and Y are (n, p) and (m, p), repectively, it returns an object (n, m) with:

• i, j = distance between row i of X and row j of Y.

Examples

using LinearAlgebra, StatsBase

X = rand(5, 3)
Y = rand(2, 3)

mahsqchol(X, Y)

S = cov(X, corrected = false)
U = cholesky(Hermitian(S)).U 
Uinv = inv(U)
mahsqchol(X, Y, Uinv)

mahsqchol(X[:, 1], 4)
mahsqchol(1, 4, sqrt(2.1))

matB(X, y, weights::Weight)

Between-class covariance matrix.

• X : X-data (n, p).
• y : A vector (n) defining the class membership.
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Compute the between-class covariance matrix (output B) of X. This is the (non-corrected) covariance matrix of the weighted class centers.

Examples

using Jchemo, StatsBase

n = 20 ; p = 3
X = rand(n, p)
y = rand(1:3, n)
tab(y) 
weights = mweight(ones(n)) 

res = matB(X, y, weights) ;
res.B
res.priors
res.ni
res.lev

res = matW(X, y, weights) ;
res.W
res.Wi

matW(X, y, weights).W + matB(X, y, weights).B
cov(X; corrected = false)

v = mweight(collect(1:n))
matW(X, y, v).priors 
matB(X, y, v).priors 
matW(X, y, v).W + matB(X, y, v).B
covm(X, v)

matW(X, y, weights::Weight)

Within-class covariance matrices.

• X : X-data (n, p).
• y : A vector (n) defing the class membership.
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Compute the (non-corrected) within-class and pooled covariance matrices (outputs Wi and W, respectively) of X.

If class i contains only one observation, Wi is computed by:

• covm(X,weights).

For examples, see function matB.

mavg(; kwargs...)
mavg(X; kwargs...)

Smoothing by moving averages of each row of X-data.

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

Keyword arguments:

• npoint : Nb. points involved in the window.

The function returns a matrix (n, p).

The smoothing is computed by convolution with padding, using function imfilter of package ImageFiltering.jl. The centered kernel is ones(npoint) / npoint. Each returned point is located on the center of the kernel. Assume a signal x of length p (row of X) correponding to a vector wl of p wavelengths (or other indexes).

If npoint = 3, the kernel is kern = [.33, .33, .33], and:

• The output value at index i = 4 is: dot(kern, [x[3], x[4], x[5]]). The output wavelength is: wl[4]
• The output value at index i = 1 is: dot(kern, [x[1], x[1], x[2]]) (padding). The corresponding wavelength is: wl[1].

If npoint = 4, the kernel is kern = [.25, .25, .25, .25], and:

• The output value at index i = 4 is: dot(kern, [x[3], x[4], x[5], x[6]]). The corresponding wavelength is: (wl[4] + wl[5]) / 2.
• The output value at index i = 1 is: dot(kern, x[1], x[1], x[2], x[3]) (padding). The corresponding wavelength is: (wl[1] + wl[2]) / 2.

References

Package ImageFiltering.jl https://github.com/JuliaImages/ImageFiltering.jl

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
year = dat.Y.year
s = year .<= 2012
Xtrain = X[s, :]
Xtest = rmrow(X, s)
wlst = names(dat.X)
wl = parse.(Float64, wlst)
plotsp(X, wl; nsamp = 20).f

model = mavg(npoint = 10) 
fit!(model, Xtrain)
Xptrain = transf(model, Xtrain)
Xptest = transf(model, Xtest)
plotsp(Xptrain).f
plotsp(Xptest).f

mbconcat()
mbconcat(Xbl)

Concatenate horizontaly multiblock X-data.

• Xbl : List of blocks (vector of matrices) of X-data Typically, output of function mblock from (n, p) data.

Examples

using Jchemo
n = 5 ; m = 3 ; p = 9 
X = rand(n, p) 
Xnew = rand(m, p)
listbl = [3:4, 1, [6; 8:9]]
Xbl = mblock(X, listbl) 
Xblnew = mblock(Xnew, listbl) 
@head Xbl[3]

model = mbconcat() 
fit!(model, Xbl)
transf(model, Xbl)
transf(model, Xblnew)

mblock(X, listbl)

Make blocks from a matrix.

• X : X-data (n, p).
• listbl : A vector whose each component defines the colum numbers defining a block in X. The length of listbl is the number of blocks.

The function returns a list (vector) of blocks.

Examples

using Jchemo
n = 5 ; p = 10 
X = rand(n, p) 
listbl = [3:4, 1, [6; 8:10]]

Xbl = mblock(X, listbl)
Xbl[1]
Xbl[2]
Xbl[3]

mbpca(; kwargs...)
mbpca(Xbl; kwargs...)
mbpca(Xbl, weights::Weight; kwargs...)
mbpca!(Xbl::Matrix, weights::Weight; kwargs...)

Consensus principal components analysis (CPCA, a.k.a MBPCA).

• Xbl : List of blocks (vector of matrices) of X-data. Typically, output of function mblock.
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. global latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. See function blockscal for possible values.
• tol : Tolerance value for Nipals convergence.
• maxit : Maximum number of iterations (Nipals).
• scal : Boolean. If true, each column of blocks in Xbl is scaled by its uncorrected standard deviation (before the block scaling).

CPCA algorithm (Westerhuis et a; 1998), a.k.a MBPCA, and reffered to as CPCA-W in Smilde et al. 2003.

Apart eventual block scaling, the MBPCA is equivalent to the PCA of the horizontally concatenated matrix X = X1 X2 ... Xk.

The function returns several objects, in particular:

• T : The global LVs (not-normed).
• U : The global LVs (normed).
• W : The block weights (normed).
• Tb : The block LVs (in the metric scale), returned grouped by LV.
• Tbl : The block LVs (in the original scale), returned grouped by block.
• Vbl : The block loadings (normed).
• lb : The block specific weights ('lambda') for the global LVs.
• mu : The sum of the block specific weights (= eigen values of the global PCA).

Function summary returns:

• explvarx : Proportion of the total X inertia (squared Frobenious norm) explained by the global LVs.
• explxbl : Proportion of the inertia of each block (= Xbl[k]) explained by the global LVs.
• contrxbl2t : Contribution of each block to the global LVs.
• rvxbl2t : RV coefficients between each block and the global LVs.
• rdxbl2t : Rd coefficients between each block and the global LVs.
• cortbl2t : Correlations between the block LVs (= Tbl[k]) and the global LVs.
• corx2t : Correlation between the X-variables and the global LVs.

References

Mangamana, E.T., Cariou, V., Vigneau, E., Glèlè Kakaï, R.L., Qannari, E.M., 2019. Unsupervised multiblock data analysis: A unified approach and extensions. Chemometrics and Intelligent Laboratory Systems 194, 103856. https://doi.org/10.1016/j.chemolab.2019.103856

Smilde, A.K., Westerhuis, J.A., de Jong, S., 2003. A framework for sequential multiblock component methods. Journal of Chemometrics 17, 323–337. https://doi.org/10.1002/cem.811

Westerhuis, J.A., Kourti, T., MacGregor, J.F., 1998. Analysis of multiblock and hierarchical PCA and PLS models. Journal of Chemometrics 12, 301–321. https://doi.org/10.1002/(SICI)1099-128X(199809/10)12:5<301::AID-CEM515>3.0.CO;2-S

Examples

using Jchemo, JchemoData, JLD2
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "ham.jld2") 
@load db dat
@names dat 
X = dat.X
group = dat.group
listbl = [1:11, 12:19, 20:25]
Xbl = mblock(X[1:6, :], listbl)
Xblnew = mblock(X[7:8, :], listbl)
n = nro(Xbl[1]) 

nlv = 3
bscal = :frob
scal = false
#scal = true
model = mbpca(; nlv, bscal, scal)
fit!(model, Xbl)
@names model 
@names model.fitm
## Global scores 
@head model.fitm.T
@head transf(model, Xbl)
transf(model, Xblnew)
## Blocks scores
i = 1
@head model.fitm.Tbl[i]
@head transfbl(model, Xbl)[i]

res = summary(model, Xbl) ;
@names res 
res.explvarx
res.explxbl   # = model.fitm.lb if bscal = :frob
rowsum(Matrix(res.explxbl))
res.contrxbl2t
res.rvxbl2t
res.rdxbl2t
res.cortbl2t
res.corx2t 

mbplskdeda(; kwargs...)
mbplskdeda(Xbl, y; kwargs...)
mbplskdeda(Xbl, y, weights::Weight; kwargs...)

Multiblock PLS-KDEDA.

• Xbl : List of blocks (vector of matrices) of X-data Typically, output of function mblock from (n, p) data.
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. See function blockscal for possible values.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• Keyword arguments of function dmkern (bandwidth definition) can also be specified here.
• scal : Boolean. If true, each column of blocks in Xbl and Ydummy is scaled by its uncorrected standard deviation (before the block scaling) in the MBPLS computation.

The principle is the same as function mbplsqda except that the densities by class are estimated from dmkern instead of dmnorm.

See function mbplslda for examples.

mbplslda(; kwargs...)
mbplslda(Xbl, y; kwargs...)
mbplslda(Xbl, y, weights::Weight; kwargs...)

Multiblock PLS-LDA.

• Xbl : List of blocks (vector of matrices) of X-data Typically, output of function mblock from (n, p) data.
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. See function blockscal for possible values.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of blocks in Xbl and Ydummy is scaled by its uncorrected standard deviation (before the block scaling) in the MBPLS computation.

The method is as follows:

1. The training variable y (univariate class membership) is transformed to a dummy table (Ydummy) containing nlev columns, where nlev is the number of classes present in y. Each column of Ydummy is a dummy (0/1) variable.
2. A multivariate MBPLSR (MBPLSR2) is run on {X, Ydummy}, returning a score matrix T.
3. A LDA is done on {T, y}, returning estimates of posterior probabilities (∊ [0, 1]) of class membership.
4. For a given observation, the final prediction is the class corresponding to the dummy variable for which the probability estimate is the highest.

In the high-level version of the present functions, the observation weights are automatically defined by the given priors (argument prior): the sub-totals by class of the observation weights are set equal to the prior probabilities. The low-level version (argument weights) allows to implement other choices.

Examples

using Jchemo, JLD2, CairoMakie, JchemoData
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
tab(Y.typ)
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)
wlst = names(X)
wl = parse.(Float64, wlst)
#plotsp(X, wl; nsamp = 20).f
##
listbl = [1:350, 351:700]
Xbltrain = mblock(Xtrain, listbl)
Xbltest = mblock(Xtest, listbl) 

nlv = 15
scal = false
#scal = true
bscal = :none
#bscal = :frob
model = mbplslda(; nlv, bscal, scal)
#model = mbplsqda(; nlv, bscal, alpha = .5, scal)
#model = mbplskdeda(; nlv, bscal, scal)
fit!(model, Xbltrain, ytrain) 
@names model 

@head transf(model, Xbltrain)
@head transf(model, Xbltest)

res = predict(model, Xbltest) ; 
@head res.pred 
@show errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xbltest; nlv = 1:2).pred

mbplsqda(; kwargs...)
mbplsqda(Xbl, y; kwargs...)
mbplsqda(Xbl, y, weights::Weight; kwargs...)

Multiblock PLS-QDA.

• Xbl : List of blocks (vector of matrices) of X-data Typically, output of function mblock from (n, p) data.
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. See function blockscal for possible values.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• alpha : Scalar (∈ [0, 1]) defining the continuum between QDA (alpha = 0) and LDA (alpha = 1).
• scal : Boolean. If true, each column of blocks in Xbl and Ydummy is scaled by its uncorrected standard deviation (before the block scaling) in the MBPLS computation.

The method is as follows:

1. The training variable y (univariate class membership) is transformed to a dummy table (Ydummy) containing nlev columns, where nlev is the number of classes present in y. Each column of Ydummy is a dummy (0/1) variable.
2. A multivariate MBPLSR (MBPLSR2) is run on {X, Ydummy}, returning a score matrix T.
3. A QDA (possibly with continuum) is done on {T, y}, returning estimates of posterior probabilities (∊ [0, 1]) of class membership.
4. For a given observation, the final prediction is the class corresponding to the dummy variable for which the probability estimate is the highest.

In the high-level version of the present functions, the observation weights are automatically defined by the given priors (argument prior): the sub-totals by class of the observation weights are set equal to the prior probabilities. The low-level version (argument weights) allows to implement other choices.

See function mbplslda for examples.

mbplsr(; kwargs...)
mbplsr(Xbl, Y; kwargs...)
mbplsr(Xbl, Y, weights::Weight; kwargs...)
mbplsr!(Xbl::Matrix, Y::Union{Matrix, BitMatrix}, weights::Weight; kwargs...)

Multiblock PLSR (MBPLSR).

• Xbl : List of blocks (vector of matrices) of X-data Typically, output of function mblock from (n, p) data.
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. global latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. See function blockscal for possible values.
• scal : Boolean. If true, each column of blocks in Xbl and Y is scaled by its uncorrected standard deviation (before the block scaling).

This function runs a PLSR on {X, Y} where X is the horizontal concatenation of the blocks in Xbl. The function gives the same global LVs and predictions as function mbplswest, but is much faster.

Function summary returns:

• explvarx : Proportion of the total X inertia (squared Frobenious norm) explained by the global LVs.
• rvxbl2t : RV coefficients between each block and the global LVs.
• rdxbl2t : Rd coefficients between each block (= Xbl[k]) and the global LVs.
• corx2t : Correlation between the X-variables and the global LVs.

Examples

using Jchemo, JchemoData, JLD2
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "ham.jld2") 
@load db dat
@names dat 
X = dat.X
Y = dat.Y
y = Y.c1
group = dat.group
listbl = [1:11, 12:19, 20:25]
s = 1:6
Xbltrain = mblock(X[s, :], listbl)
Xbltest = mblock(rmrow(X, s), listbl)
ytrain = y[s]
ytest = rmrow(y, s) 
ntrain = nro(ytrain) 
ntest = nro(ytest) 
ntot = ntrain + ntest 
(ntot = ntot, ntrain , ntest)

nlv = 3
bscal = :frob
scal = false
#scal = true
model = mbplsr(; nlv, bscal, scal)
fit!(model, Xbltrain, ytrain)
@names model 
@names model.fitm
@head model.fitm.T
@head transf(model, Xbltrain)
transf(model, Xbltest)

res = predict(model, Xbltest)
res.pred 
rmsep(res.pred, ytest)

res = summary(model, Xbltrain) ;
@names res 
res.explvarx
res.rvxbl2t
res.rdxbl2t
res.cortbl2t
res.corx2t 

mbplsrda(; kwargs...)
mbplsrda(Xbl, y; kwargs...)
mbplsrda(Xbl, y, weights::Weight; kwargs...)

Discrimination based on multiblock partial least squares regression (MBPLSR-DA).

• Xbl : List of blocks (vector of matrices) of X-data Typically, output of function mblock from (n, p) data.
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. See function blockscal for possible values.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of blocks in Xbl and Ydummy is scaled by its uncorrected standard deviation (before the block scaling).

The method is as follows:

1. The training variable y (univariate class membership) is transformed to a dummy table (Ydummy) containing nlev columns, where nlev is the number of classes present in y. Each column of Ydummy is a dummy (0/1) variable.
2. Then, a multivariate MBPLSR (MBPLSR2) is run on {X, Ydummy}, returning predictions of the dummy variables (= object posterior returned by fuction predict). These predictions can be considered as unbounded estimates (i.e. eventuall outside of [0, 1]) of the class membership probabilities.
3. For a given observation, the final prediction is the class corresponding to the dummy variable for which the probability estimate is the highest.

In the high-level version of the present functions, the observation weights are automatically defined by the given priors (argument prior): the sub-totals by class of the observation weights are set equal to the prior probabilities. The low-level version (argument weights) allows to implement other choices.

Examples

using Jchemo, JLD2, CairoMakie, JchemoData
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
tab(Y.typ)
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)
wlst = names(X)
wl = parse.(Float64, wlst)
#plotsp(X, wl; nsamp = 20).f
##
listbl = [1:350, 351:700]
Xbltrain = mblock(Xtrain, listbl)
Xbltest = mblock(Xtest, listbl) 

nlv = 15
scal = false
#scal = true
bscal = :none
#bscal = :frob
model = mbplsrda(; nlv, bscal, scal)
fit!(model, Xbltrain, ytrain) 
@names model 

@head model.fitm.fitm.T 
@head transf(model, Xbltrain)
@head transf(model, Xbltest)

res = predict(model, Xbltest) ; 
@head res.pred 
@show errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xbltest; nlv = 1:2).pred

mbplswest(; kwargs...)
mbplswest(Xbl, Y; kwargs...)
mbplswest(Xbl, Y, weights::Weight; kwargs...)
mbplswest!(Xbl::Matrix, Y::Matrix, weights::Weight; kwargs...)

Multiblock PLSR (MBPLSR) - Nipals algorithm.

• Xbl : List of blocks (vector of matrices) of X-data Typically, output of function mblock from (n, p) data.
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. global latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. See function blockscal for possible values.
• tol : Tolerance value for convergence (Nipals).
• maxit : Maximum number of iterations (Nipals).
• scal : Boolean. If true, each column of blocks in Xbl and Y is scaled by its uncorrected standard deviation (before the block scaling).

This functions implements the MBPLSR Nipals algorithm such as in Westerhuis et al. 1998. The function gives the same global scores and predictions as function mbplsr.

Function summary returns:

• explvarx : Proportion of the total X inertia (squared Frobenious norm) explained by the global LVs.
• rvxbl2t : RV coefficients between each block and the global LVs.
• rdxbl2t : Rd coefficients between each block (= Xbl[k]) and the global LVs.
• cortbl2t : Correlations between the block LVs (= Tbl[k]) and the global LVs.
• corx2t : Correlation between the X-variables and the global LVs.

References

Westerhuis, J.A., Kourti, T., MacGregor, J.F., 1998. Analysis of multiblock and hierarchical PCA and PLS models. Journal of Chemometrics 12, 301–321. https://doi.org/10.1002/(SICI)1099-128X(199809/10)12:5<301::AID-CEM515>3.0.CO;2-S

Examples

using Jchemo, JchemoData, JLD2
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "ham.jld2") 
@load db dat
@names dat 
X = dat.X
Y = dat.Y
y = Y.c1
group = dat.group
listbl = [1:11, 12:19, 20:25]
s = 1:6
Xbltrain = mblock(X[s, :], listbl)
Xbltest = mblock(rmrow(X, s), listbl)
ytrain = y[s]
ytest = rmrow(y, s) 
ntrain = nro(ytrain) 
ntest = nro(ytest) 
ntot = ntrain + ntest 
(ntot = ntot, ntrain , ntest)

nlv = 3
bscal = :frob
scal = false
#scal = true
model = mbplswest(; nlv, bscal, scal)
fit!(model, Xbltrain, ytrain)
@names model 
@names model.fitm
@head model.fitm.T
@head transf(model, Xbltrain)
transf(model, Xbltest)

res = predict(model, Xbltest)
res.pred 
rmsep(res.pred, ytest)

res = summary(model, Xbltrain) ;
@names res 
res.explvarx
res.rvxbl2t
res.rdxbl2t
res.cortbl2t
res.corx2t 

meanv(x)
meanv(x, weights::Weight)

Compute the mean of a vector.

• x : A vector (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Examples

using Jchemo

n = 100
x = rand(n)
w = mweight(rand(n)) 

meanv(x)
meanv(x, w)

merrp(pred, y)

Compute the mean intra-class classification error rate.

• pred : Predictions.
• y : Observed data (class membership).

ERRP (see function errp) is computed for each class. Function merrp returns the average of these intra-class ERRPs.

Examples

using Jchemo

Xtrain = rand(10, 5) 
ytrain = rand(["a" ; "b"], 10)
Xtest = rand(4, 5) 
ytest = rand(["a" ; "b"], 4)

model = plsrda(; nlv = 2)
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
merrp(pred, ytest)

mlev(x)

Return the sorted levels of a vector or a dataset.

Examples

using Jchemo

x = rand(["a";"b";"c"], 20)
lev = mlev(x)
nlev = length(lev)

X = reshape(x, 5, 4)
mlev(X)

df = DataFrame(g1 = rand(1:2, n), g2 = rand(["a"; "c"], n))
mlev(df)

mlr(; kwargs...)
mlr(X, Y; kwargs...)
mlr(X, Y, weights::Weight; kwargs...)
mlr!(X::Matrix, Y::Union{Matrix, BitMatrix}, weights::Weight; kwargs...)

Compute a mutiple linear regression model (MLR) by using the QR algorithm.

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• noint : Boolean. Define if the model is computed with an intercept or not.

Safe but can be little slower than other methods.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie 
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/iris.jld2") 
@load db dat
@names dat
@head dat.X
X = dat.X[:, 2:4]
y = dat.X[:, 1]
n = nro(X)
ntest = 30
s = samprand(n, ntest) 
Xtrain = X[s.train, :]
ytrain = y[s.train]
Xtest = X[s.test, :]
ytest = y[s.test]

model = mlr()
#model = mlrchol()
#model = mlrpinv()
#model = mlrpinvn() 
fit!(model, Xtrain, ytrain) 
@names model
@names model.fitm
fitm = model.fitm ;
fitm.B
fitm.int 
coef(model) 
res = predict(model, Xtest)
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",  
    ylabel = "Observed").f    

model = mlr(noint = true)
fit!(model, Xtrain, ytrain) 
coef(model) 

model = mlrvec()
fit!(model, Xtrain[:, 1], ytrain) 
coef(model) 

mlrchol()
mlrchol(X, Y)
mlrchol(X, Y, weights::Weight)
mlrchol!mlrchol!(X::Matrix, Y::Matrix, weights::Weight)

Compute a mutiple linear regression model (MLR) using the Normal equations and a Choleski factorization.

• X : X-data, with nb. columns >= 2 (required by function cholesky).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Compute only a model with intercept.

Faster but can be less accurate (based on squared element X'X).

See function mlr for examples.

mlrda(; kwargs...)
mlrda(X, y; kwargs...)
mlrda(X, y, weights::Weight)

Discrimination based on multple linear regression (MLR-DA).

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).

The method is as follows:

1. The training variable y (univariate class membership) is transformed to a dummy table (Ydummy) containing nlev columns, where nlev is the number of classes present in y. Each column of Ydummy is a dummy (0/1) variable.
2. Then, a multiple linear regression (MLR) is run on {X, Ydummy}, returning predictions of the dummy variables (= object posterior returned by fuction predict). These predictions can be considered as unbounded estimates (i.e. eventuall outside of [0, 1]) of the class membership probabilities.
3. For a given observation, the final prediction is the class corresponding to the dummy variable for which the probability estimate is the highest.

In the high-level version of the present functions, the observation weights are automatically defined by the given priors (argument prior): the sub-totals by class of the observation weights are set equal to the prior probabilities. The low-level version (argument weights) allows to implement other choices.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/iris.jld2")
@load db dat
@names dat
@head dat.X
X = dat.X[:, 1:4]
y = dat.X[:, 5]
n = nro(X)
ntest = 30
s = samprand(n, ntest)
Xtrain = X[s.train, :]
ytrain = y[s.train]
Xtest = X[s.test, :]
ytest = y[s.test]
ntrain = n - ntest
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

model = mlrda()
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

mlrpinv()
mlrpinv(X, Y; kwargs...)
mlrpinv(X, Y, weights::Weight; kwargs...)
mlrpinv!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Compute a mutiple linear regression model (MLR) by using a pseudo-inverse.

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• noint : Boolean. Define if the model is computed with an intercept or not.

Safe but can be slower.

See function mlr for examples.

mlrpinvn() 
mlrpinvn(X, Y)
mlrpinvn(X, Y, weights::Weight)
mlrpinvn!mlrchol!(X::Matrix, Y::Matrix, 
    weights::Weight)

Compute a mutiple linear regression model (MLR) by using the Normal equations and a pseudo-inverse.

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Safe and fast for p not too large.

Compute only a model with intercept.

See function mlr for examples.

mlrvec(; kwargs...)
mlrvec(X, Y; kwargs...)
mlrvec(X, Y, weights::Weight; kwargs...)
mlrvec!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Compute a simple (univariate x) linear regression model.

• x : Univariate X-data (n).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• noint : Boolean. Define if the model is computed with an intercept or not.

See function mlr for examples.

mpar(; kwargs...)

Return a tuple with all the combinations of the parameter values defined in kwargs. Keyword arguments:

• kwargs : Vector(s) of the parameter(s) values.

Examples

using Jchemo
nlvdis = 25 ; metric = [:mah] 
h = [1 ; 2 ; Inf] ; k = [500 ; 1000] 
pars = mpar(nlvdis = nlvdis, metric = metric, h = h, k = k) 
length(pars[1])
reduce(hcat, pars)

mse(pred, Y; digits = 3)

Summary of model performance for regression.

• pred : Predictions.
• Y : Observed data.

Examples

using Jchemo

Xtrain = rand(10, 5) 
Ytrain = rand(10, 2)
ytrain = Ytrain[:, 1]
Xtest = rand(4, 5) 
Ytest = rand(4, 2)
ytest = Ytest[:, 1]

model = plskern; nlv = 2)
fit!(model, Xtrain, Ytrain)
pred = predict(model, Xtest).pred
mse(pred, Ytest)

model = plskern; nlv = 2)
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
mse(pred, ytest)

msep(pred, Y)

Compute the mean of the squared prediction errors (MSEP).

• pred : Predictions.
• Y : Observed data.

Examples

using Jchemo 

Xtrain = rand(10, 5) 
Ytrain = rand(10, 2)
ytrain = Ytrain[:, 1]
Xtest = rand(4, 5) 
Ytest = rand(4, 2)
ytest = Ytest[:, 1]

model = plskern; nlv = 2)
fit!(model, Xtrain, Ytrain)
pred = predict(model, Xtest).pred
msep(pred, Ytest)

model = plskern; nlv = 2)
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred
msep(pred, ytest)

mweight(x::Vector)

Return an object of type Weight containing vector w = x / sum(x) (if ad'hoc building, w must sum to 1).

Examples

using Jchemo

x = rand(10)
w = mweight(x)
sum(w.w)

mweightcla(x::AbstractVector; prior::Union{Symbol, Vector} = :unif)
mweightcla(Q::DataType, x::Vector; prior::Union{Symbol, Vector} = :unif)

Compute observation weights for a categorical variable, given specified sub-total weights for the classes.

• x : A categorical variable (n) (class membership).
• Q : A data type (e.g. Float32).

Keyword arguments:

• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (the vector must be sorted in the same order as mlev(x)).

Return an object of type Weight (see function mweight) containing a vector w (n) that sums to 1.

Examples

using Jchemo

x = vcat(rand(["a" ; "c"], 900), repeat(["b"], 100))
tab(x)
weights = mweightcla(x)
#weights = mweightcla(x; prior = :prop)
#weights = mweightcla(x; prior = [.1, .7, .2])
res = aggstat(weights.w, x; algo = sum)
[res.lev res.X]

nco(X)

Return the nb. columns of X.

nipals(X; kwargs...)
nipals(X, UUt, VVt; kwargs...)

Nipals to compute the first score and loading vectors of a matrix.

• X : X-data (n, p).
• UUt : Matrix (n, n) for Gram-Schmidt orthogonalization.
• VVt : Matrix (p, p) for Gram-Schmidt orthogonalization.

Keyword arguments:

• tol : Tolerance value for stopping the iterations.
• maxit : Maximum nb. of iterations.

The function finds:

• {u, v, sv} = argmin(||X - u * sv * v'||)

with the constraints:

• ||u|| = ||v|| = 1

using the alternating least squares algorithm to compute SVD (Gabriel & Zalir 1979).

At the end, X ~ u * sv * v', where:

• u : left singular vector (u * sv = scores)
• v : right singular vector (loadings)
• sv : singular value.

When NIPALS is used on sequentially deflated matrices, vectors u and v can loose orthogonality due to accumulation of rounding errors. Orthogonality can be rebuilt from the Gram-Schmidt method (arguments UUt and VVt).

References

K.R. Gabriel, S. Zamir, Lower rank approximation of matrices by least squares with any choice of weights, Technometrics 21 (1979) 489–498.

Examples

using Jchemo, LinearAlgebra

X = rand(5, 3)

res = nipals(X)
res.niter
res.sv
svd(X).S[1] 
res.v
svd(X).V[:, 1] 
res.u
svd(X).U[:, 1] 

nipalsmiss(X; kwargs...)
nipalsmiss(X, UUt, VVt; kwargs...)

Nipals to compute the first score and loading vectors of a matrix with missing data.

• X : X-data (n, p).
• UUt : Matrix (n, n) for Gram-Schmidt orthogonalization.
• VVt : Matrix (p, p) for Gram-Schmidt orthogonalization.

Keyword arguments:

• tol : Tolerance value for stopping the iterations.
• maxit : Maximum nb. of iterations.

See function nipals.

References

K.R. Gabriel, S. Zamir, Lower rank approximation of matrices by least squares with any choice of weights, Technometrics 21 (1979) 489–498.

Wright, K., 2018. Package nipals: Principal Components Analysis using NIPALS with Gram-Schmidt Orthogonalization. https://cran.r-project.org/

Examples

using Jchemo 

X = [1. 2 missing 4 ; 4 missing 6 7 ; 
    missing 5 6 13 ; missing 18 7 6 ; 
    12 missing 28 7] 

res = nipalsmiss(X)
res.niter
res.sv
res.v
res.u

normv(x)
normv(x, weights::Weight)

Compute the norm of a vector.

• x : A vector (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

The weighted norm of vector x is computed by:

• sqrt(x' * D * x), where D is the diagonal matrix of vector weights.w.

References

@gdkrmr, https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/17

@Stevengj, https://discourse.julialang.org/t/interesting-post-about-simd-dot-product-and-cosine-similarity/123282.

Examples

using Jchemo

n = 1000
x = rand(n)
w = mweight(ones(n))

normv(x)
sqrt(n) * normv(x, w)

nro(X)

Return the nb. rows of X.

occod(; kwargs...)
occod(fitm, X; kwargs...)

One-class classification using PCA/PLS orthognal distance (OD).

• fitm : The preliminary model (e.g. object fitm returned by function pcasvd) that was fitted on the training data assumed to represent the training class.
• X : Training X-data (n, p), on which was fitted the model fitm.

Keyword arguments:

• mcut : Type of cutoff. Possible values are: :mad, :q. See Thereafter.
• cri : When mcut = :mad, a constant. See thereafter.
• risk : When mcut = :q, a risk-I level. See thereafter.

In this method, the outlierness d of an observation is the orthogonal distance (= 'X-residuals') of this observation, ie. the Euclidean distance between the observation and its projection on the score plan defined by the fitted (e.g. PCA) model (e.g. Hubert et al. 2005, Van Branden & Hubert 2005 p. 66, Varmuza & Filzmoser 2009 p. 79).

See function occsd for details on outputs.

References

M. Hubert, V. J. Rousseeuw, K. Vanden Branden (2005). ROBPCA: a new approach to robust principal components analysis. Technometrics, 47, 64-79.

K. Vanden Branden, M. Hubert (2005). Robuts classification in high dimension based on the SIMCA method. Chem. Lab. Int. Syst, 79, 10-21.

K. Varmuza, V. Filzmoser (2009). Introduction to multivariate statistical analysis in chemometrics. CRC Press, Boca Raton.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/challenge2018.jld2") 
@load db dat
@names dat
X = dat.X    
Y = dat.Y
model = savgol(npoint = 21, deriv = 2, degree = 3)
fit!(model, X) 
Xp = transf(model, X) 
s = Bool.(Y.test)
Xtrain = rmrow(Xp, s)
Ytrain = rmrow(Y, s)
Xtest = Xp[s, :]
Ytest = Y[s, :]

## Below, the reference class is "EEH"
cla1 = "EHH" ; cla2 = "PEE" ; cod = "out"   # here cla2 should be detected
#cla1 = "EHH" ; cla2 = "EHH" ; cod = "in"   # here cla2 should not be detected
s1 = Ytrain.typ .== cla1
s2 = Ytest.typ .== cla2
zXtrain = Xtrain[s1, :]    
zXtest = Xtest[s2, :] 
ntrain = nro(zXtrain)
ntest = nro(zXtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)
ytrain = repeat(["in"], ntrain)
ytest = repeat([cod], ntest)

## Group description
model = pcasvd(nlv = 10) 
fit!(model, zXtrain) 
Ttrain = model.fitm.T
Ttest = transf(model, zXtest)
T = vcat(Ttrain, Ttest)
group = vcat(repeat(["1"], ntrain), repeat(["2"], ntest))
i = 1
plotxy(T[:, i], T[:, i + 1], group; leg_title = "Class", xlabel = string("PC", i), 
    ylabel = string("PC", i + 1)).f

#### Occ
## Preliminary PCA fitted model
model0 = pcasvd(nlv = 10) 
fit!(model0, zXtrain)
## Outlierness
model = occod()
#model = occod(mcut = :mad, cri = 4)
#model = occod(mcut = :q, risk = .01)
#model = occsdod()
fit!(model, model0.fitm, zXtrain) 
@names model 
@names model.fitm 
@head d = model.fitm.d
d = d.dstand
f, ax = plotxy(1:length(d), d; size = (500, 300), 
    xlabel = "Obs. index", ylabel = "Standardized distance")
hlines!(ax, 1; linestyle = :dot)
f

res = predict(model, zXtest) ;
@names res
@head res.d
@head res.pred
tab(res.pred)
errp(res.pred, ytest)
conf(res.pred, ytest).cnt
d1 = model.fitm.d.dstand
d2 = res.d.dstand
d = vcat(d1, d2)
f, ax = plotxy(1:length(d), d, group; size = (500, 300), leg_title = "Class", 
    xlabel = "Obs. index", ylabel = "Standardized distance")
hlines!(ax, 1; linestyle = :dot)
f

occsd(; kwargs...)
occsd(fitm; kwargs...)

One-class classification using PCA/PLS score distance (SD).

• fitm : The preliminary model (e.g. object fitm returned by function pcasvd) that was fitted on the training data assumed to represent the training class.

Keyword arguments:

• mcut : Type of cutoff. Possible values are: :mad, :q. See Thereafter.
• cri : When mcut = :mad, a constant. See thereafter.
• risk : When mcut = :q, a risk-I level. See thereafter.

In this method, the outlierness d of an observation is defined by its score distance (SD), ie. the Mahalanobis distance between the projection of the observation on the score plan defined by the fitted (e.g. PCA) model and the center of the score plan.

If a new observation has d higher than a given cutoff, the observation is assumed to not belong to the training (= reference) class. The cutoff is computed with non-parametric heuristics. Noting [d] the vector of outliernesses computed on the training class:

• If mcut = :mad, then cutoff = median([d]) + cri * madv([d]).
• If mcut = :q, then cutoff is estimated from the empirical cumulative density function computed on [d], for a given risk-I (risk).

Alternative approximate cutoffs have been proposed in the literature (e.g.: Nomikos & MacGregor 1995, Hubert et al. 2005, Pomerantsev 2008). Typically, and whatever the approximation method used to compute the cutoff, it is recommended to tune this cutoff depending on the detection objectives.

Outputs

• pval: Estimate of p-value (see functions pval) computed from the training distribution [d].
• dstand: standardized distance defined as d / cutoff. A value dstand > 1 may be considered as extreme compared to the distribution of the training data.
• gh is the Winisi "GH" (usually, GH > 3 is considered as extreme).

Specific for function predict:

• pred: class prediction
  • dstand <= 1 ==> in: the observation is expected to belong to the training class,
  • dstand > 1 ==> out: extreme value, possibly not belonging to the same class as the training.

References

M. Hubert, V. J. Rousseeuw, K. Vanden Branden (2005). ROBPCA: a new approach to robust principal components analysis. Technometrics, 47, 64-79.

Nomikos, V., MacGregor, J.F., 1995. Multivariate SPC Charts for Monitoring Batch Processes. null 37, 41-59. https://doi.org/10.1080/00401706.1995.10485888

Pomerantsev, A.L., 2008. Acceptance areas for multivariate classification derived by projection methods. Journal of Chemometrics 22, 601-609. https://doi.org/10.1002/cem.1147

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/challenge2018.jld2") 
@load db dat
@names dat
X = dat.X    
Y = dat.Y
model = savgol(npoint = 21, deriv = 2, degree = 3)
fit!(model, X) 
Xp = transf(model, X) 
s = Bool.(Y.test)
Xtrain = rmrow(Xp, s)
Ytrain = rmrow(Y, s)
Xtest = Xp[s, :]
Ytest = Y[s, :]

## Below, the reference class is "EEH"
cla1 = "EHH" ; cla2 = "PEE" ; cod = "out"   # here cla2 should be detected
#cla1 = "EHH" ; cla2 = "EHH" ; cod = "in"   # here cla2 should not be detected
s1 = Ytrain.typ .== cla1
s2 = Ytest.typ .== cla2
zXtrain = Xtrain[s1, :]    
zXtest = Xtest[s2, :] 
ntrain = nro(zXtrain)
ntest = nro(zXtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)
ytrain = repeat(["in"], ntrain)
ytest = repeat([cod], ntest)

## Group description
model = pcasvd(nlv = 10) 
fit!(model, zXtrain) 
Ttrain = model.fitm.T
Ttest = transf(model, zXtest)
T = vcat(Ttrain, Ttest)
group = vcat(repeat(["1"], ntrain), repeat(["2"], ntest))
i = 1
plotxy(T[:, i], T[:, i + 1], group; leg_title = "Class", xlabel = string("PC", i),  
    ylabel = string("PC", i + 1)).f

#### Occ
## Preliminary PCA fitted model
model0 = pcasvd(nlv = 30) 
fit!(model0, zXtrain)
## Outlierness
model = occsd()
#model = occsd(mcut = :mad, cri = 4)
#model = occsd(mcut = :q, risk = .01)
fit!(model, model0.fitm) 
@names model 
@names model.fitm 
@head d = model.fitm.d
d = d.dstand
f, ax = plotxy(1:length(d), d; size = (500, 300), xlabel = "Obs. index", 
    ylabel = "Standardized distance")
hlines!(ax, 1; linestyle = :dot)
f

res = predict(model, zXtest) ;
@names res
@head res.d
@head res.pred
tab(res.pred)
errp(res.pred, ytest)
conf(res.pred, ytest).cnt
d1 = model.fitm.d.dstand
d2 = res.d.dstand
d = vcat(d1, d2)
f, ax = plotxy(1:length(d), d, group; size = (500, 300), leg_title = "Class", xlabel = "Obs. index", 
    ylabel = "Standardized distance")
hlines!(ax, 1; linestyle = :dot)
f

occsdod(; kwargs...)
occsdod(object, X; kwargs...)

One-class classification using a consensus between PCA/PLS score and orthogonal (SD and OD) distances.

• fitm : The preliminary model (e.g. object fitm returned by function pcasvd) that was fitted on the training data assumed to represent the training class.
• X : Training X-data (n, p), on which was fitted the model fitm.

Keyword arguments:

• mcut : Type of cutoff. Possible values are: :mad, :q. See Thereafter.
• cri : When mcut = :mad, a constant. See thereafter.
• risk : When mcut = :q, a risk-I level. See thereafter.

In this method, the outlierness d of a given observation is a consensus between the score distance (SD) and the orthogonal distance (OD). The consensus is computed from the standardized distances by:

• dstand = sqrt(dstand_sd * dstand_od).

See functions:

• occsd for details of the outputs,
• and occod for examples.

occstah(; kwargs...)
occstah(X; kwargs...)

One-class classification using the Stahel-Donoho outlierness.

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

Keyword arguments:

• nlv : Nb. dimensions on which X is projected.
• mcut : Type of cutoff. Possible values are: :mad, :q. See Thereafter.
• cri : When mcut = :mad, a constant. See thereafter.
• risk : When mcut = :q, a risk-I level. See thereafter.
• scal : Boolean. If true, each column of X is scaled such as in function outstah.

In this method, the outlierness d of a given observation is the Stahel-Donoho outlierness (see ?outstah).

See function occsd for details on outputs.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/challenge2018.jld2") 
@load db dat
@names dat
X = dat.X    
Y = dat.Y
model = savgol(npoint = 21, deriv = 2, degree = 3)
fit!(model, X) 
Xp = transf(model, X) 
s = Bool.(Y.test)
Xtrain = rmrow(Xp, s)
Ytrain = rmrow(Y, s)
Xtest = Xp[s, :]
Ytest = Y[s, :]

## Below, the reference class is "EEH"
cla1 = "EHH" ; cla2 = "PEE" ; cod = "out"   # here cla2 should be detected
#cla1 = "EHH" ; cla2 = "EHH" ; cod = "in"   # here cla2 should not be detected
s1 = Ytrain.typ .== cla1
s2 = Ytest.typ .== cla2
zXtrain = Xtrain[s1, :]    
zXtest = Xtest[s2, :] 
ntrain = nro(zXtrain)
ntest = nro(zXtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)
ytrain = repeat(["in"], ntrain)
ytest = repeat([cod], ntest)

## Group description
model = pcasvd(nlv = 10) 
fit!(model, zXtrain) 
Ttrain = model.fitm.T
Ttest = transf(model, zXtest)
T = vcat(Ttrain, Ttest)
group = vcat(repeat(["1"], ntrain), repeat(["2"], ntest))
i = 1
plotxy(T[:, i], T[:, i + 1], group; leg_title = "Class", xlabel = string("PC", i),  
    ylabel = string("PC", i + 1)).f

#### Occ
## Preliminary dimension reduction 
## (Not required but often more efficient)
nlv = 50
model0 = pcasvd(; nlv)
fit!(model0, zXtrain)
Ttrain = model0.fitm.T
Ttest = transf(model0, zXtest)
## Outlierness
model = occstah(; nlv, scal = true)
fit!(model, Ttrain) 
@names model 
@names model.fitm 
@head d = model.fitm.d
d = d.dstand
f, ax = plotxy(1:length(d), d; size = (500, 300), xlabel = "Obs. index", 
    ylabel = "Standardized distance")
hlines!(ax, 1; linestyle = :dot)
f

res = predict(model, Ttest) ;
@names res
@head res.d
@head res.pred
tab(res.pred)
errp(res.pred, ytest)
conf(res.pred, ytest).cnt
d1 = model.fitm.d.dstand
d2 = res.d.dstand
d = vcat(d1, d2)
f, ax = plotxy(1:length(d), d, group; size = (500, 300), leg_title = "Class", 
    xlabel = "Obs. index", ylabel = "Standardized distance")
hlines!(ax, 1; linestyle = :dot)
f

out(x)

Return if elements of a vector are strictly outside of a given range.

• x : Univariate data.
• y : Univariate data on which is computed the range (min, max).

Return a BitVector.

Examples

using Jchemo

x = [-200.; -100; -1; 0; 1; 200]
out(x, [-1; .2; 1])
out(x, (-1, 1))

outeucl(X, V; kwargs...)
outeucl!(X::Matrix, V::Matrix; kwargs...)

Compute an outlierness from Euclidean distances to center.

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

Keyword arguments:

• scal : Boolean. If true, each column of X is scaled by its MAD before computing the outlierness.

Outlyingness is calculated by the Euclidean distance between the observation (rows of X) and a robust estimate of the center of the data (in the present function, the spatial median). Such outlyingness was for instance used in the robust PLSR algorithm of Serneels et al. 2005 (PRM).

References

Serneels, S., Croux, C., Filzmoser, V., Van Espen, V.J., 2005. Partial robust M-regression. Chemometrics and Intelligent Laboratory Systems 79, 55-64. https://doi.org/10.1016/j.chemolab.2005.04.007

Examples

using Jchemo, CairoMakie
n = 300 ; p = 700 ; m = 80
ntot = n + m
X1 = randn(n, p)
X2 = randn(m, p) .+ rand(1:3, p)'
X = vcat(X1, X2)

nlv = 10
scal = false
#scal = true
res = outeucl(X; scal) ;
@names res
res.d    # outlierness 
plotxy(1:ntot, res.d).f

outstah(X, V; kwargs...)
outstah!(X::Matrix, V::Matrix; kwargs...)

Compute the Stahel-Donoho outlierness.

• X : X-data (n, p).
• V : A projection matrix (p, nlv) representing the directions of the projection pursuit.

Keyword arguments:

• scal : Boolean. If true, each column of X is scaled by its MAD before computing the outlierness.

See Maronna and Yohai 1995 for details on the outlierness measure.

A projection-pursuit approach is used: given a projection matrix V (p, nlv) (in general built randomly), the observations (rows of X) are projected on the nlv directions and the Stahel-Donoho outlierness is computed for each observation from these projections.

References

Maronna, R.A., Yohai, V.J., 1995. The Behavior of the Stahel-Donoho Robust Multivariate Estimator. Journal of the American Statistical Association 90, 330–341. https://doi.org/10.1080/01621459.1995.10476517

Examples

using Jchemo, CairoMakie

n = 300 ; p = 700 ; m = 80
ntot = n + m
X1 = randn(n, p)
X2 = randn(m, p) .+ rand(1:3, p)'
X = vcat(X1, X2)

nlv = 10
V = rand(0:1, p, nlv)
scal = false
#scal = true
res = outstah(X, V; scal) ;
@names res
res.d    # outlierness 
plotxy(1:ntot, res.d).f

parsemiss(Q, x::Vector{Union{String, Missing}})

Parsing a string vector allowing missing data.

• Q : Type that results from the parsing of type `String'.
• x : A string vector containing missing (of type Missing) observations.

See examples.

Examples

using Jchemo

x = ["1"; "3.2"; missing]
x_p = parsemiss(Float64, x)

pcaeigen(; kwargs...)
pcaeigen(X; kwargs...)
pcaeigen(X, weights::Weight; kwargs...)
pcaeigen!(X::Matrix, weights::Weight; kwargs...)

PCA by Eigen factorization.

• X : X-data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of principal components (PCs).
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

Let us note D the (n, n) diagonal matrix of weights (weights.w) and X the centered matrix in metric D. The function minimizes ||X - T * V'||^2 in metric D, by computing an Eigen factorization of X' * D * X.

See function pcasvd for examples.

pcaeigenk(; kwargs...)
pcaeigenk(X; kwargs...)
pcaeigenk(X, weights::Weight; kwargs...)
pcaeigenk!(X::Matrix, weights::Weight; kwargs...)

PCA by Eigen factorization of the kernel matrix XX'.

• X : X-data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of principal components (PCs).
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

This is the "kernel cross-product" version of the PCA algorithm (e.g. Wu et al. 1997). For wide matrices (n << p, where p is the nb. columns) and n not too large, this algorithm can be much faster than the others.

Let us note D the (n, n) diagonal matrix of weights (weights.w) and X the centered matrix in metric D. The function minimizes ||X - T * V'||^2 in metric D, by computing an Eigen factorization of D^(1/2) * X * X' D^(1/2).

See function pcasvd for examples.

References

Wu, W., Massart, D.L., de Jong, S., 1997. The kernel PCA algorithms for wide data. Part I: Theory and algorithms. Chemometrics and Intelligent Laboratory Systems 36, 165-172. https://doi.org/10.1016/S0169-7439(97)00010-5

pcanipals(; kwargs...)
pcanipals(X; kwargs...)
pcanipals(X, weights::Weight; kwargs...)
pcanipals!(X::Matrix, weights::Weight; kwargs...)

PCA by NIPALS algorithm.

• X : X-data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of principal components (PCs).
• gs : Boolean. If true (default), a Gram-Schmidt orthogonalization of the scores and loadings is done before each X-deflation.
• tol : Tolerance value for stopping the iterations.
• maxit : Maximum nb. of iterations.
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

Let us note D the (n, n) diagonal matrix of weights (weights.w) and X the centered matrix in metric D. The function minimizes ||X - T * V'||^2 in metric D by NIPALS.

See function pcasvd for examples.

References

Andrecut, M., 2009. Parallel GPU Implementation of Iterative PCA Algorithms. Journal of Computational Biology 16, 1593-1599. https://doi.org/10.1089/cmb.2008.0221

K.R. Gabriel, S. Zamir, Lower rank approximation of matrices by least squares with any choice of weights, Technometrics 21 (1979) 489–498.

Gabriel, R. K., 2002. Le biplot - Outil d'exploration de données multidimensionnelles. Journal de la Société Française de la Statistique, 143, 5-55.

Lingen, F.J., 2000. Efficient Gram-Schmidt orthonormalisation on parallel computers. Communications in Numerical Methods in Engineering 16, 57-66. https://doi.org/10.1002/(SICI)1099-0887(200001)16:1<57::AID-CNM320>3.0.CO;2-I

Tenenhaus, M., 1998. La régression PLS: théorie et pratique. Editions Technip, Paris, France.

Wright, K., 2018. Package nipals: Principal Components Analysis using NIPALS with Gram-Schmidt Orthogonalization. https://cran.r-project.org/

pcanipalsmiss(; kwargs...)
pcanipals(X; kwargs...)
pcanipals(X, weights::Weight; kwargs...)
pcanipals!(X::Matrix, weights::Weight; kwargs...)

PCA by NIPALS algorithm allowing missing data.

• X : X-data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of principal components (PCs).
• gs : Boolean. If true (default), a Gram-Schmidt orthogonalization of the scores and loadings is done before each X-deflation.
• tol : Tolerance value for stopping the iterations.
• maxit : Maximum nb. of iterations.
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

References

Wright, K., 2018. Package nipals: Principal Components Analysis using NIPALS with Gram-Schmidt Orthogonalization. https://cran.r-project.org/

Examples

X = [1 2. missing 4 ; 4 missing 6 7 ; 
    missing 5 6 13 ; missing 18 7 6 ; 
    12 missing 28 7] 

nlv = 3 
tol = 1e-15
scal = false
#scal = true
gs = false
#gs = true
model = pcanipalsmiss(; nlv, tol, gs, maxit = 500, scal)
fit!(model, X)
@names model 
@names model.fitm
fitm = model.fitm ;
fitm.niter
fitm.sv
fitm.V
fitm.T
## Orthogonality 
## only if gs = true
fitm.T' * fitm.T
fitm.V' * fitm.V

## Impute missing data in X
model = pcanipalsmiss(; nlv = 2, gs = true) ;
fit!(model, X)
Xfit = xfit(model.fitm)
s = ismissing.(X)
X_imp = copy(X)
X_imp[s] .= Xfit[s]
X_imp

pcaout(; kwargs...)
pcaout(X; kwargs...)
pcaout(X, weights::Weight; kwargs...)
pcaout!(X::Matrix, weights::Weight; kwargs...)

Robust PCA using outlierness.

• X : X-data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of principal components (PCs).
• prm : Proportion of the data removed (hard rejection of outliers) for each outlierness measure.
• scal : Boolean. If true, each column of X is scaled by its MAD when computing the outlierness and by its uncorrected standard deviation when computing weighted PCA.

Robust PCA combining outlyingness measures and weighted PCA (WPCA).

The objective is to remove the effect of multivariate X-outliers that have potentially bad leverages. Observations (X-rows) receive weights depending on two outlyingness indicators:

1. The Stahel-Donoho outlyingness (Maronna and Yohai, 1995) is computed (function outstah) on X. The proportion prm of the observations with the highest outlyingness values receive a weight w1 = 0 (the other receive a weight w1 = 1).
2. An outlyingness based on the Euclidean distance to center (function outstah) is computed. The proportion prm of the observations with the highest outlyingness values receive a weight w2 = 0 (the other receive a weight w2 = 1).

The final weights of the observations are computed by weights.w * w1 * w2 that is used in a weighted PCA.

By default, the function uses prm = .3 (such as in the ROBPCA algorithm of Hubert et al. 2005, 2009).

References

Hubert, M., Rousseeuw, V.J., Vanden Branden, K., 2005. ROBPCA: A New Approach to Robust Principal Component Analysis. Technometrics 47, 64-79. https://doi.org/10.1198/004017004000000563

Hubert, M., Rousseeuw, V., Verdonck, T., 2009. Robust PCA for skewed data and its outlier map. Computational Statistics & Data Analysis 53, 2264-2274. https://doi.org/10.1016/j.csda.2008.05.027

Maronna, R.A., Yohai, V.J., 1995. The Behavior of the Stahel-Donoho Robust Multivariate Estimator. Journal of the American Statistical Association 90, 330–341. https://doi.org/10.1080/01621459.1995.10476517

Examples

using Jchemo, JchemoData, JLD2, CairoMakie 
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "octane.jld2") 
@load db dat
@names dat
X = dat.X 
wlst = names(X)
wl = parse.(Float64, wlst)
n = nro(X)

nlv = 3
model = pcaout(; nlv)  
#model = pcasvd(; nlv) 
fit!(model, X)
@names model
@names model.fitm
@head T = model.fitm.T
## Same as:
transf(model, X)

i = 1
plotxy(T[:, i], T[:, i + 1]; zeros = true, xlabel = string("PC", i), 
    ylabel = string("PC", i + 1)).f

pcapp(; kwargs...)
pcapp(X; kwargs...)
pcapp!(X::Matrix; kwargs...)

Robust PCA by projection pursuit.

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

Keyword arguments:

• nlv : Nb. of principal components (PCs).
• nsim : Nb. of additional (to X-rows) simulated directions for the projection pursuit.
• scal : Boolean. If true, each column of X is scaled by its MAD.

For nsim = 0, this is the Croux & Ruiz-Gazen (C-R, 2005) PCA algorithm that uses a projection pursuit (PP) method. Data X are robustly centered by the spatial median, and the observations are projected to the "PP" directions defined by the observations (rows of X) after they are normed. The first PCA loading vector is the direction (within the PP directions) that maximizes a given 'projection index', here the median absolute deviation (MAD). Then, X is deflated to this loading vector, and the process is re-run to define the next loading vector. And so on.

A possible extension of this algorithm is to randomly simulate additionnal candidate PP directions to the n row observations. If nsim > 0, the function simulates nsim additional PP directions to the n initial ones, as proposed in Hubert et al. (2005): random couples of observations are sampled in X and, for each couple, the direction passes through the two observations of the couple (see function simpphub).

References

Croux, C., Ruiz-Gazen, A., 2005. High breakdown estimators for principal components: the projection-pursuit approach revisited. Journal of Multivariate Analysis 95, 206–226. https://doi.org/10.1016/j.jmva.2004.08.002

Hubert, M., Rousseeuw, V.J., Vanden Branden, K., 2005. ROBPCA: A New Approach to Robust Principal Component Analysis. Technometrics 47, 64-79. https://doi.org/10.1198/004017004000000563

Examples

using Jchemo, JchemoData, JLD2, CairoMakie 
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "octane.jld2") 
@load db dat
@names dat
X = dat.X 
wlst = names(X)
wl = parse.(Float64, wlst)
n = nro(X)

nlv = 3
model = pcapp(; nlv, nsim = 2000)  
#model = pcasvd(; nlv) 
fit!(model, X)
@names model
@names model.fitm
@head T = model.fitm.T
## Same as:
@head transf(model, X)

i = 1
plotxy(T[:, i], T[:, i + 1]; zeros = true, xlabel = string("PC", i), 
    ylabel = string("PC", i + 1)).f

pcasph(; kwargs...)
pcasph(X; kwargs...)
pcasph(X, weights::Weight; kwargs...)
pcasph!(X::Matrix, weights::Weight; kwargs...)

Spherical PCA.

• X : X-data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of principal components (PCs).
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

Spherical PCA (Locantore et al. 1990, Maronna 2005, Daszykowski et al. 2007). Matrix X is centered by the spatial median computed by function Jchemo.colmedspa.

References

Daszykowski, M., Kaczmarek, K., Vander Heyden, Y., Walczak, B., 2007. Robust statistics in data analysis - A review. Chemometrics and Intelligent Laboratory Systems 85, 203-219. https://doi.org/10.1016/j.chemolab.2006.06.016

Locantore N., Marron J.S., Simpson D.G., Tripoli N., Zhang J.T., Cohen K.L. Robust principal component analysis for functional data, Test 8 (1999) 1–7

Maronna, R., 2005. Principal components and orthogonal regression based on robust scales, Technometrics, 47:3, 264-273, DOI: 10.1198/004017005000000166

Examples

using Jchemo, JchemoData, JLD2, CairoMakie 
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "octane.jld2") 
@load db dat
@names dat
X = dat.X 
wlst = names(X)
wl = parse.(Float64, wlst)
n = nro(X)

nlv = 3
model = pcasph(; nlv)  
#model = pcasvd(; nlv) 
fit!(model, X)
@names model
@names model.fitm
@head T = model.fitm.T
## Same as:
transf(model, X)

i = 1
plotxy(T[:, i], T[:, i + 1]; zeros = true, xlabel = string("PC", i), 
    ylabel = string("PC", i + 1)).f

pcasvd(; kwargs...)
pcasvd(X; kwargs...)
pcasvd(X, weights::Weight; kwargs...)
pcasvd!(X::Matrix, weights::Weight; kwargs...)

PCA by SVD factorization.

• X : X-data (n, p).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of principal components (PCs).
• scal : Boolean. If true, each column of X is scaled by its uncorrected standard deviation.

Let us note D the (n, n) diagonal matrix of weights (weights.w) and X the centered matrix in metric D. The function minimizes ||X - T * V'||^2 in metric D, by computing a SVD factorization of sqrt(D) * X:

• sqrt(D) * X ~ U * S * V'

Outputs are:

• T = D^(-1/2) * U * S
• V = V
• The diagonal of S

Examples

using Jchemo, JchemoData, JLD2, CairoMakie 
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/iris.jld2") 
@load db dat
@names dat
@head dat.X
X = dat.X[:, 1:4]
n = nro(X)
ntest = 30
s = samprand(n, ntest) 
@head Xtrain = X[s.train, :]
@head Xtest = X[s.test, :]

nlv = 3
model = pcasvd(; nlv)
#model = pcaeigen(; nlv)
#model = pcaeigenk(; nlv)
#model = pcanipals(; nlv)
fit!(model, Xtrain)
@names model
@names model.fitm
@head T = model.fitm.T
## Same as:
@head transf(model, X)
T' * T
@head V = model.fitm.V
V' * V

@head Ttest = transf(model, Xtest)

res = summary(model, Xtrain) ;
@names res
res.explvarx
res.contr_var
res.coord_var
res.cor_circle

pcr(; kwargs...)
pcr(X, Y; kwargs...)
pcr(X, Y, weights::Weight; kwargs...)
pcr!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Principal component regression (PCR) with a SVD factorization.

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• Same as function pcasvd

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlv = 15
model = pcr(; nlv) ;
fit!(model, Xtrain, ytrain)
@names model
fitm = model.fitm ;
@names fitm
@names fitm.fitm

@head fitm.fitm.T
@head transf(model, X)

coef(model)
coef(model; nlv = 3)

@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

res = predict(model, Xtest)
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",  
    ylabel = "Observed").f    

res = predict(model, Xtest; nlv = 1:2)
@head res.pred[1]
@head res.pred[2]

res = summary(model, Xtrain) ;
@names res
z = res.explvarx
plotgrid(z.nlv, z.cumpvar; step = 2, xlabel = "Nb. LVs", ylabel = "Prop. Explained X-Variance").f

pip(args...)

Build a pipeline of models.

• args... : Succesive models, see examples.

Examples

using JLD2, CairoMakie, JchemoData
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)

## Pipeline Snv :> Savgol :> Pls :> Svmr

model1 = snv()
model2 = savgol(npoint = 11, deriv = 2, degree = 3)
model3 = plskern(nlv = 15)
model4 = svmr(gamma = 1e3, cost = 1000, epsilon = .1)
model = pip(model1, model2, model3, model4)
fit!(model, Xtrain, ytrain)
res = predict(model, Xtest) ; 
@head res.pred 
rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",
      ylabel = "Observed").f

plotconf(object; size = (500, 400), cnt = true, ptext = true, 
    fontsize = 15, coldiag = :red, )

Plot a confusion matrix.

• object : Output of function conf.

Keyword arguments:

• size : Size (horizontal, vertical) of the figure.
• cnt : Boolean. If true, plot the occurrences, else plot the row %s.
• ptext : Boolean. If true, display the value in each cell.
• fontsize : Font size when ptext = true.
• coldiag : Font color when ptext = true.

See examples in help page of function conf. ```

plotgrid(indx::AbstractVector, r; size = (500, 300), step = 5, 
    color = nothing, kwargs...)
plotgrid(indx::AbstractVector, r, group; size = (700, 350), 
    step = 5, color = nothing, leg = true, leg_title = "Group", kwargs...)

Plot error/performance rates of a model.

• indx : A numeric variable representing the grid of model parameters, e.g. the nb. LVs if PLSR models.
• r : The error/performance rate.

Keyword arguments:

• group : Categorical variable defining groups. A separate line is plotted for each level of group.
• size : Size (horizontal, vertical) of the figure.
• step : Step used for defining the xticks.
• color : Set color. If group if used, must be a vector of same length as the number of levels in group.
• leg : Boolean. If group is used, display a legend or not.
• leg_title : Title of the legend.
• kwargs : Optional arguments to pass in Axis of CairoMakie.

To use plotgrid, a backend (e.g. CairoMakie) has to be specified.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

model = plskern() 
nlv = 0:20
res = gridscore(model, Xtrain, ytrain, 
    Xtest, ytest; score = rmsep, nlv)
plotgrid(res.nlv, res.y1; xlabel = "Nb. LVs", ylabel = "RMSEP").f

model = lwplsr() 
nlvdis = 15 ; metric = [:mah]
h = [1 ; 2.5 ; 5] ; k = [50 ; 100] 
pars = mpar(nlvdis = nlvdis, metric = metric, 
    h = h, k = k)
nlv = 0:20
res = gridscore(model, Xtrain, ytrain, Xtest, ytest; score = rmsep, 
    pars, nlv)
group = string.("h=", res.h, " k=", res.k)
plotgrid(res.nlv, res.y1, group; xlabel = "Nb. LVs", ylabel = "RMSECV").f

plotsp(X, wl = 1:nco(X); size = (500, 300), color = nothing, nsamp = nothing, 
    kwargs...)

Plotting spectra.

• X : X-data (n, p).
• wl : Column names of X. Must be numeric.

Keyword arguments:

• size : Size (horizontal, vertical) of the figure.
• color : Set a unique color (and eventually transparency) to the spectra.
• nsamp : Nb. spectra (X-rows) to plot. If nothing, all spectra are plotted.
• kwargs : Optional arguments to pass in Axis of CairoMakie.

The function plots the rows of X.

To use plotxy, a backend (e.g. CairoMakie) has to be specified.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X
wlst = names(X)
wl = parse.(Float64, wlst) 

plotsp(X).f
plotsp(X; color = (:red, .2)).f
plotsp(X, wl; xlabel = "Wavelength (nm)", ylabel = "Absorbance").f

tck = collect(wl[1]:200:wl[end]) ;
plotsp(X, wl; xlabel = "Wavelength (nm)", ylabel = "Absorbance", xticks = tck).f

f, ax = plotsp(X, wl; color = (:red, .2))
xmeans = colmean(X)
lines!(ax, wl, xmeans; color = :black, linewidth = 2)
vlines!(ax, 1200)
f

plotxy(x, y; size = (500, 300), color = nothing, ellipse::Bool = false, 
    prob = .95, circle::Bool = false, bisect::Bool = false, zeros::Bool = false,
    xlabel = "", ylabel = "", title = "", kwargs...)
plotxy(x, y, group; size = (600, 350), color = nothing, ellipse::Bool = false, 
    prob = .95, circle::Bool = false, bisect::Bool = false, zeros::Bool = false,
    xlabel = "", ylabel = "", title = "", leg::Bool = true, leg_title = "Group", 
    kwargs...)

Scatter plot of (x, y) data

• x : A x-vector (n).
• y : A y-vector (n).
• group : Categorical variable defining groups (n).

Keyword arguments:

• size : Size (horizontal, vertical) of the figure.
• color : Set color(s). If group if used, color must be a vector of same length as the number of levels in group.
• ellipse : Boolean. Draw an ellipse of confidence, assuming a Ch-square distribution with df = 2. If group is used, one ellipse is drawn per group.
• prob : Probability for the ellipse of confidence.
• bisect : Boolean. Draw a bisector.
• zeros : Boolean. Draw horizontal and vertical axes passing through origin (0, 0).
• xlabel : Label for the x-axis.
• ylabel : Label for the y-axis.
• title : Title of the graphic.
• leg : Boolean. If group is used, display a legend or not.
• leg_title : Title of the legend.
• kwargs : Optional arguments to pass in function scatter of Makie.

To use plotxy, a backend (e.g. CairoMakie) has to be specified.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
lev = mlev(year)
nlev = length(lev)

model = pcasvd(nlv = 5)  
fit!(model, X) 
@head T = model.fitm.T

plotxy(T[:, 1], T[:, 2]; color = (:red, .5)).f

plotxy(T[:, 1], T[:, 2], year; ellipse = true, xlabel = "PC1", ylabel = "PC2").f

i = 2
colm = cgrad(:Dark2_5, nlev; categorical = true)
plotxy(T[:, i], T[:, i + 1], year; color = colm, xlabel = string("PC", i), 
    ylabel = string("PC", i + 1), zeros = true, ellipse = true).f

plotxy(T[:, 1], T[:, 2], year).lev

plotxy(1:5, 1:5).f

y = reshape(rand(5), 5, 1)
plotxy(1:5, y).f

## Several layers can be added
## (same syntax as in Makie)
A = rand(50, 2)
f, ax = plotxy(A[:, 1], A[:, 2]; xlabel = "x1", ylabel = "x2")
ylims!(ax, -1, 2)
hlines!(ax, 0.5; color = :red, linestyle = :dot)
f

plscan(; kwargs...)
plscan(X, Y; kwargs...)
plscan(X, Y, weights::Weight; kwargs...)
plscan!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Canonical partial least squares regression (Canonical PLS).

• X : First block of data.
• Y : Second block of data.
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. Possible values are: :none, :frob. See functions blockscal.
• scal : Boolean. If true, each column of blocks X and Y is scaled by its uncorrected standard deviation (before the block scaling).

Canonical PLS with the Nipals algorithm (Wold 1984, Tenenhaus 1998 chap.11), referred to as PLS-W2A (i.e. Wold PLS mode A) in Wegelin 2000. The two blocks X and Y play a symmetric role. After each step of scores computation, X and Y are deflated by the x- and y-scores, respectively.

Function summary returns:

• cortx2ty: Correlations between the X- and Y-LVs.

and for block X:

• explvarx : Proportion of the block inertia (squared Frobenious norm) explained by the block LVs (Tx).
• rvx2tx : RV coefficients between the block and the block LVs.
• rdx2tx : Rd coefficients between the block and the block LVs.
• corx2tx : Correlation between the block variables and the block LVs.

The same is returned for block Y.

References

Tenenhaus, M., 1998. La régression PLS: théorie et pratique. Editions Technip, Paris.

Wegelin, J.A., 2000. A Survey of Partial Least Squares (PLS) Methods, with Emphasis on the Two-Block Case (No. 371). University of Washington, Seattle, Washington, USA.

Wold, S., Ruhe, A., Wold, H., Dunn, III, W.J., 1984. The Collinearity Problem in Linear Regression. The Partial Least Squares (PLS) Approach to Generalized Inverses. SIAM Journal on Scientific and Statistical Computing 5, 735–743. https://doi.org/10.1137/0905052

Examples

using Jchemo, JchemoData, JLD2
mypath = dirname(dirname(pathof(JchemoData)))
db = joinpath(mypath, "data", "linnerud.jld2") 
@load db dat
@names dat
X = dat.X
Y = dat.Y
n, p = size(X)
q = nco(Y)

nlv = 2
bscal = :frob
model = plscan(; nlv, bscal)
fit!(model, X, Y)
@names model
@names model.fitm

fitm = model.fitm
@head fitm.Tx
@head transfbl(model, X, Y).Tx

@head fitm.Ty
@head transfbl(model, X, Y).Ty

res = summary(model, X, Y) ;
@names res
res.explvarx
res.explvary
res.cortx2ty
res.rvx2tx
res.rvy2ty
res.rdx2tx
res.rdy2ty
res.corx2tx 
res.cory2ty 

plskdeda(; kwargs...)
plskdeda(X, y; kwargs...)
plskdeda(X, y, weights::Weight; kwargs...)

KDE-DA on PLS latent variables (PLS-KDEDA).

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute. Must be >= 1.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation in the PLS computation.
• Keyword arguments of function dmkern (bandwidth definition) can also be specified here.

The principle is the same as function plsqda except that the densities by class are estimated from dmkern instead of dmnorm.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlv = 15
model = plskdeda(; nlv) 
#model = plskdeda(; nlv, a = .5)
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

embfitm = fitm.fitm.embfitm ;
@head embfitm.T
@head transf(model, Xtrain)
@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

coef(embfitm)

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xtest; nlv = 1:2).pred
summary(embfitm, Xtrain)

plskern(; kwargs...)
plskern(X, Y; kwargs...)
plskern(X, Y, weights::Weight; kwargs...)
plskern!(X::Matrix, Y::Union{Matrix, BitMatrix}, weights::Weight; kwargs...)

Partial least squares regression (PLSR) with the "improved kernel algorithm #1" (Dayal & McGegor, 1997).

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute.
• scal : Boolean. If true, each column of X and Y is scaled by its uncorrected standard deviation.

About the row-weighting in PLS algorithms (weights): See in particular Schaal et al. 2002, Siccard & Sabatier 2006, Kim et al. 2011, and Lesnoff et al. 2020.

References

Dayal, B.S., MacGregor, J.F., 1997. Improved PLS algorithms. Journal of Chemometrics 11, 73-85.

Kim, S., Kano, M., Nakagawa, H., Hasebe, S., 2011. Estimation of active pharmaceutical ingredients content using locally weighted partial least squares and statistical wavelength selection. Int. J. Pharm., 421, 269-274.

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. e3209. https://onlinelibrary.wiley.com/doi/abs/10.1002/cem.3209

Schaal, S., Atkeson, C., Vijayamakumar, S. 2002. Scalable techniques from nonparametric statistics for the real time robot learning. Applied Intell., 17, 49-60.

Sicard, E. Sabatier, R., 2006. Theoretical framework for local PLS1 regression and application to a rainfall dataset. Comput. Stat. Data Anal., 51, 1393-1410.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlv = 15
model = plskern(; nlv) ;
#model = plsnipals(; nlv) ;
#model = plswold(; nlv) ;
#model = plsrosa(; nlv) ;
#model = plssimp(; nlv) ;
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
@head model.fitm.T

coef(model)
coef(model; nlv = 3)

@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

res = predict(model, Xtest)
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",  
    ylabel = "Observed").f    

res = predict(model, Xtest; nlv = 1:2)
@head res.pred[1]
@head res.pred[2]

res = summary(model, Xtrain) ;
@names res
z = res.explvarx
plotgrid(z.nlv, z.cumpvar; step = 2, xlabel = "Nb. LVs", 
    ylabel = "Prop. Explained X-Variance").f

plslda(; kwargs...)
plslda(X, y; kwargs...)
plslda(X, y, weights::Weight; kwargs...)

LDA on PLS latent variables (PLS-LDA).

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute. Must be >= 1.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation in the PLS computation.

LDA on PLS latent variables. The method is as follows:

1. The training variable y (univariate class membership) is transformed to a dummy table (Ydummy) containing nlev columns, where nlev is the number of classes present in y. Each column of Ydummy is a dummy (0/1) variable.
2. A multivariate PLSR (PLSR2) is run on {X, Ydummy}, returning a score matrix T.
3. A LDA is done on {T, y}, returning estimates of posterior probabilities (∊ [0, 1]) of class membership.
4. For a given observation, the final prediction is the class corresponding to the dummy variable for which the probability estimate is the highest.

In the high-level version of the present functions, the observation weights are automatically defined by the given priors (argument prior): the sub-totals by class of the observation weights are set equal to the prior probabilities. The low-level version (argument weights) allows to implement other choices.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlv = 15
model = plslda(; nlv) 
#model = plslda(; nlv, prior = :prop) 
#model = plsqda(; nlv, alpha = .1) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni

embfitm = fitm.fitm.embfitm ;
@head embfitm.T
@head transf(model, Xtrain)
@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

coef(embfitm)

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xtest; nlv = 1:2).pred
summary(embfitm, Xtrain)

plsnipals(; kwargs...)
plsnipals(X, Y; kwargs...)
plsnipals(X, Y, weights::Weight; kwargs...)
plsnipals!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Partial Least Squares Regression (PLSR) with the Nipals algorithm.

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute.
• scal : Boolean. If true, each column of X and Y is scaled by its uncorrected standard deviation.

In this function, for PLS2 (multivariate Y), the Nipals iterations are replaced by a direct computation of the PLS weights (w) by SVD decomposition of matrix X'Y (Hoskuldsson 1988 p.213).

See function plskern for examples.

References

Hoskuldsson, A., 1988. PLS regression methods. Journal of Chemometrics 2, 211-228.https://doi.org/10.1002/cem.1180020306

Tenenhaus, M., 1998. La régression PLS: théorie et pratique. Editions Technip, Paris, France.

Wold, S., Sjostrom, M., Eriksson, l., 2001. PLS-regression: a basic tool for chemometrics. Chem. Int. Lab. Syst., 58, 109-130.

plsqda(; kwargs...)
plsqda(X, y; kwargs...)
plsqda(X, y, weights::Weight; kwargs...)

QDA on PLS latent variables (PLS-QDA) with continuum.

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute. Must be >= 1.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• alpha : Scalar (∈ [0, 1]) defining the continuum between QDA (alpha = 0) and LDA (alpha = 1).
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation in the PLS computation.

QDA on PLS latent variables. The method is as follows:

1. The training variable y (univariate class membership) is transformed to a dummy table (Ydummy) containing nlev columns, where nlev is the number of classes present in y. Each column of Ydummy is a dummy (0/1) variable.
2. A multivariate PLSR (PLSR2) is run on {X, Ydummy}, returning a score matrix T.
3. A QDA (possibly with continuum) is done on {T, y}, returning estimates of posterior probabilities (∊ [0, 1]) of class membership.
4. For a given observation, the final prediction is the class corresponding to the dummy variable for which the probability estimate is the highest.

In the high-level version of the present functions, the observation weights are automatically defined by the given priors (argument prior): the sub-totals by class of the observation weights are set equal to the prior probabilities. The low-level version (argument weights) allows to implement other choices.

See function plslda for examples.

plsravg(; kwargs...)
plsravg(X, Y; kwargs...)
plsravg(X, Y, weights::Weight; kwargs...)
plsravg!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Averaging PLSR models with different numbers of latent variables (PLSR-AVG).

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : A range of nb. of latent variables (LVs) to compute.
• scal : Boolean. If true, each column of X and Y is scaled by its uncorrected standard deviation.

Ensemblist method where the predictions are computed by averaging the predictions of a set of models built with different numbers of LVs.

For instance, if argument nlv is set to nlv = 5:10, the prediction for a new observation is the simple average of the predictions returned by the models with 5 LVs, 6 LVs, ... 10 LVs, respectively.

References

Lesnoff, M., Andueza, D., Barotin, C., Barre, V., Bonnal, L., Fernández Pierna, J.A., Picard, F., Vermeulen, V., Roger, J.-M., 2022. Averaging and Stacking Partial Least Squares Regression Models to Predict the Chemical Compositions and the Nutritive Values of Forages from Spectral Near Infrared Data. Applied Sciences 12, 7850. https://doi.org/10.3390/app12157850

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2") 
@load db dat
@names dat
X = dat.X 
Y = dat.Y
@head Y
y = Y.ndf
#y = Y.dm
n = nro(X)
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(y, s)
Xtest = X[s, :]
ytest = y[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)

nlv = 0:30
#nlv = 5:20
#nlv = 25
model = plsravg(; nlv) ;
fit!(model, Xtrain, ytrain)

res = predict(model, Xtest)
@head res.pred
res.predlv   # predictions for each nb. of LVs 
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",  
    ylabel = "Observed").f    

plsrda(; kwargs...)
plsrda(X, y; kwargs...)
plsrda(X, y, weights::Weight; kwargs...)

Discrimination based on partial least squares regression (PLSR-DA).

• X : X-data (n, p).
• y : Univariate class membership (n).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute.
• prior : Type of prior probabilities for class membership. Possible values are: :unif (uniform), :prop (proportional), or a vector (of length equal to the number of classes) giving the prior weight for each class (in case of vector, it must be sorted in the same order as mlev(y)).
• scal : Boolean. If true, each column of X and Ydummy is scaled by its uncorrected standard deviation.

This is the usual "PLSDA". The method is as follows:

1. The training variable y (univariate class membership) is transformed to a dummy table (Ydummy) containing nlev columns, where nlev is the number of classes present in y. Each column of Ydummy is a dummy (0/1) variable.
2. Then, a multivariate PLSR (PLSR2) is run on {X, Ydummy}, returning predictions of the dummy variables (= object posterior returned by fuction predict). These predictions can be considered as unbounded estimates (i.e. eventuall outside of [0, 1]) of the class membership probabilities.
3. For a given observation, the final prediction is the class corresponding to the dummy variable for which the probability estimate is the highest.

In the high-level version of the present functions, the observation weights are automatically defined by the given priors (argument prior):

• the sub-totals by class of the observation weights are set equal to the prior probabilities.

The low-level version (argument weights) allows to implement other choices.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/forages2.jld2")
@load db dat
@names dat
X = dat.X
Y = dat.Y
n = nro(X) 
s = Bool.(Y.test)
Xtrain = rmrow(X, s)
ytrain = rmrow(Y.typ, s)
Xtest = X[s, :]
ytest = Y.typ[s]
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = n, ntrain, ntest)
tab(ytrain)
tab(ytest)

nlv = 15
model = plsrda(; nlv) 
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
fitm = model.fitm ;
fitm.lev
fitm.ni
@names fitm.fitm
aggsum(fitm.fitm.weights.w, ytrain)

@head fitm.fitm.T
@head transf(model, Xtrain)
@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

coef(fitm.fitm)

res = predict(model, Xtest) ;
@names res
@head res.posterior
@head res.pred
errp(res.pred, ytest)
conf(res.pred, ytest).cnt

predict(model, Xtest; nlv = 1:2).pred
summary(fitm.fitm, Xtrain)

plsrosa(; kwargs...)
plsrosa(X, Y; kwargs...)
plsrosa(X, Y, weights::Weight; kwargs...)
plsrosa!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Partial Least Squares Regression (PLSR) with the ROSA algorithm (Liland et al. 2016).

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute.
• scal : Boolean. If true, each column of X and Y is scaled by its uncorrected standard deviation.

Note: The function has the following differences with the original algorithm of Liland et al. (2016):

• Scores T (LVs) are not normed.
• Multivariate Y is allowed.

See function plskern for examples.

References

Liland, K.H., Næs, T., Indahl, U.G., 2016. ROSA—a fast extension of partial least squares regression for multiblock data analysis. Journal of Chemometrics 30, 651–662. https://doi.org/10.1002/cem.2824

plsrout(; kwargs...)
plsrout(X, Y; kwargs...)
plsrout(X, Y, weights::Weight; kwargs...)
pcaout!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Robust PLSR using outlierness.

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. of latent variables (LVs).
• prm : Proportion of the data removed (hard rejection of outliers) for each outlierness measure.
• scal : Boolean. If true, each column of X is scaled by its MAD when computing the outlierness and by its uncorrected standard deviation when computing weighted PCA.

Robust PLSR combining outlyingness measures and weighted PLSR (WPLSR). This is the same principle as function pcaout (see the help page) but the final step is a weighted PLSR instead of a weighted PCA.

Examples

using Jchemo, JchemoData, JLD2, CairoMakie
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/cassav.jld2") 
@load db dat
@names dat
X = dat.X 
y = dat.Y.tbc
year = dat.Y.year
tab(year)
s = year .<= 2012
Xtrain = X[s, :]
ytrain = y[s]
Xtest = rmrow(X, s)
ytest = rmrow(y, s)

nlv = 15
model = plsrout(; nlv)
fit!(model, Xtrain, ytrain)
@names model
@names model.fitm
@head model.fitm.T

coef(model)
coef(model; nlv = 3)

@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)

res = predict(model, Xtest)
@head res.pred
@show rmsep(res.pred, ytest)
plotxy(res.pred, ytest; color = (:red, .5), bisect = true, xlabel = "Prediction",  
    ylabel = "Observed").f    

res = predict(model, Xtest; nlv = 1:2)
@head res.pred[1]
@head res.pred[2]

plssimp(; kwargs...)
plssimp(X, Y; kwargs...)
plssimp(X, Y, weights::Weight; kwargs...)
plssimp!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Partial Least Squares Regression (PLSR) with the SIMPLS algorithm (de Jong 1993).

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute.
• scal : Boolean. If true, each column of X and Y is scaled by its uncorrected standard deviation.

Note: In this function, scores T (LVs) are not normed, conversely to the original algorithm of de Jong (2013).

See function plskern for examples.

References

de Jong, S., 1993. SIMPLS: An alternative approach to partial least squares regression. Chemometrics and Intelligent Laboratory Systems 18, 251–263. https://doi.org/10.1016/0169-7439(93)85002-X

plstuck(; kwargs...)
plstuck(X, Y; kwargs...)
plstuck(X, Y, weights::Weight; kwargs...)
plstuck!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Tucker's inter-battery method of factor analysis

• X : First block of data.
• Y : Second block of data.
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs = scores) to compute.
• bscal : Type of block scaling. Possible values are: :none, :frob. See functions blockscal.
• scal : Boolean. If true, each column of blocks X and Y is scaled by its uncorrected standard deviation (before the block scaling).

Inter-battery method of factor analysis (Tucker 1958, Tenenhaus 1998 chap.3). The two blocks X and X play a symmetric role. This method is referred to as PLS-SVD in Wegelin 2000. The method factorizes the covariance matrix X'Y by SVD.

See function plscan for the details on the summary outputs.

References

Tenenhaus, M., 1998. La régression PLS: théorie et pratique. Editions Technip, Paris.

Tishler, A., Lipovetsky, S., 2000. Modelling and forecasting with robust canonical analysis: method and application. Computers & Operations Research 27, 217–232. https://doi.org/10.1016/S0305-0548(99)00014-3

Tucker, L.R., 1958. An inter-battery method of factor analysis. Psychometrika 23, 111–136. https://doi.org/10.1007/BF02289009

Wegelin, J.A., 2000. A Survey of Partial Least Squares (PLS) Methods, with Emphasis on the Two-Block Case (No. 371). University of Washington, Seattle, Washington, USA.

Examples

using Jchemo, JchemoData, JLD2
path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/linnerud.jld2") 
@load db dat
@names dat
X = dat.X 
Y = dat.Y

model = plstuck(nlv = 3)
fit!(model, X, Y) 
@names model
@names model.fitm

fitm = model.fitm
@head fitm.Tx
@head transfbl(model, X, Y).Tx

@head fitm.Ty
@head transfbl(model, X, Y).Ty

res = summary(model, X, Y) ;
@names res
res.explvarx
res.explvary
res.cortx2ty
res.rvx2tx
res.rvy2ty
res.rdx2tx
res.rdy2ty
res.corx2tx 
res.cory2ty 

plswold(; kwargs...)
plswold(X, Y; kwargs...)
plswold(X, Y, weights::Weight; kwargs...)
plswold!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)

Partial Least Squares Regression (PLSR) with the Wold algorithm

• X : X-data (n, p).
• Y : Y-data (n, q).
• weights : Weights (n) of the observations. Must be of type Weight (see e.g. function mweight).

Keyword arguments:

• nlv : Nb. latent variables (LVs) to compute.
• tol : Tolerance for the Nipals algorithm.
• maxit : Maximum number of iterations for the Nipals algorithm.
• scal : Boolean. If true, each column of X and Y is scaled by its uncorrected standard deviation.

Wold Nipals PLSR algorithm: Tenenhaus 1998 p.204.

See function plskern for examples.

References

Tenenhaus, M., 1998. La régression PLS: théorie et pratique. Editions Technip, Paris, France.

Wold, S., Ruhe, A., Wold, H., Dunn, III, W.J., 1984. The Collinearity Problem in Linear Regression. The Partial Least Squares (PLS). Approach to Generalized Inverses. SIAM Journal on Scientific and Statistical Computing 5, 735–743. https://doi.org/10.1137/0905052

predict(object::Calds, X; kwargs...)

Compute predictions from a fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Calpds, X; kwargs...)

Compute predictions from a fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Cglsr, X; nlv = nothing)

Compute Y-predictions from a fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.
• nlv : Nb. iterations, or collection of nb. iterations, to consider.

predict(object::Dkplsr, X; nlv = nothing)

Compute Y-predictions from a fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.
• nlv : Nb. LVs, or collection of nb. LVs, to consider.

predict(object::Dmkern, x)

Compute predictions from a fitted model.

• object : The fitted model.
• x : Data (vector) for which predictions are computed.

predict(object::Dmnorm, X)

Compute predictions from a fitted model.

• object : The fitted model.
• X : Data (vector) for which predictions are computed.

predict(object::Knnda1, X)

Compute the y-predictions from the fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Knnr, X)

Compute the Y-predictions from the fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Kplsr, X; nlv = nothing)

Compute Y-predictions from a fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.
• nlv : Nb. LVs, or collection of nb. LVs, to consider.

If nothing, it is the maximum nb. LVs.

predict(object::Krr, X; lb = nothing)

Compute Y-predictions from a fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.
• lb : Regularization parameter, or collection of regularization parameters, "lambda" to consider.

predict(object::Loessr, X)

Compute y-predictions from a fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Lwmlr, X)

Compute the Y-predictions from the fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Lwmlrda, X)

Compute y-predictions from the fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Lwplslda, X; nlv = nothing)

Compute the y-predictions from the fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Lwplsqda, X; nlv = nothing)

Compute the y-predictions from the fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Lwplsr, X; nlv = nothing)

Compute the Y-predictions from the fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Lwplsravg, X)

Compute the Y-predictions from the fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Lwplsrda, X; nlv = nothing)

Compute the y-predictions from the fitted model.

• object : The fitted model.
• X : X-data for which predictions are computed.

predict(object::Mbplsprobda, Xbl; nlv = nothing)

Compute Y-predictions from a fitted model.

• object : The fitted model.
• Xbl : A list of blocks (vector of matrices) of X-data for which predictions are computed.
• nlv : Nb. LVs, or collection of nb. LVs, to consider.

predict(object::Mbplsrda, Xbl; nlv = nothing)

Compute Y-predictions from a fitted model.

• object : The fitted model.
• Xbl : A list of blocks (vector of matrices) of X-data for which predictions are computed.
• nlv : Nb. LVs, or collection of nb. LVs, to consider.

predict(object::Mlrda, X)