using Jchemo, JchemoData
using JLD2, CairoMakie

Data importation

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

(:X, :Y)

X = dat.X
@head X

... (178, 100)

Y = dat.Y 
@head Y

... (178, 4)

typ = Y.typ
tab(typ)

OrderedCollections.OrderedDict{String, Int64} with 3 entries:
  "test"  => 31
  "train" => 115
  "val"   => 32

wlst = names(X)
wl = parse.(Float64, wlst) 
#plotsp(X, wl; xlabel = "Wavelength (nm)", ylabel = "Absorbance").f

100-element Vector{Float64}:
  850.0
  852.0
  854.0
  856.0
  858.0
  860.0
  862.0
  864.0
  866.0
  868.0
    ⋮
 1032.0
 1034.0
 1036.0
 1038.0
 1040.0
 1042.0
 1044.0
 1046.0
 1048.0

Preprocessing

model1 = snv()
model2 = savgol(npoint = 15, deriv = 2, degree = 3)
model = pip(model1, model2)
fit!(model, X)
@head Xp = transf(model, X) 
#plotsp(Xp, wl; xlabel = "Wavelength (nm)", ylabel = "Absorbance").f

3×100 Matrix{Float64}:
 0.000397076  0.000495203  0.000598623  …  0.00827138  0.00917311  0.00946072
 0.00242055   0.00244366   0.00234233      0.00580631  0.00689249  0.00749316
 0.0011927    0.00122721   0.00120098      0.0101019   0.0108142   0.0108444
... (178, 100)

Split Tot to Train/Test

s = typ .== "train"
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)

(ntot = 178, ntrain = 115, ntest = 63)

Working response y

namy = names(Y)[1:3]
j = 2  
nam = namy[j]    # work on the j-th y-variable
ytrain = Ytrain[:, nam]
ytest = Ytest[:, nam]

63-element Vector{Float64}:
 29.8
  1.4
  4.6
 11.0
 17.0
 22.4
 27.9
 46.5
  6.1
  2.0
  ⋮
 18.1
 19.4
 24.8
 27.2
 28.4
 31.3
 33.8
 35.5
 42.5

CV-Segments for model tuning

Replicated K-fold CV

K = 3     # nb. folds (segments)
rep = 10   # nb. replications
segm = segmkf(ntrain, K; rep = rep)

10-element Vector{Vector{Vector{Int64}}}:
 [[2, 6, 7, 9, 16, 18, 19, 21, 30, 32  …  94, 96, 98, 99, 100, 107, 110, 111, 113, 114], [1, 4, 11, 17, 26, 28, 29, 37, 40, 41  …  93, 95, 101, 102, 104, 105, 106, 109, 112, 115], [3, 5, 8, 10, 12, 13, 14, 15, 20, 22  …  61, 71, 74, 76, 80, 81, 92, 97, 103, 108]]
 [[2, 3, 9, 10, 15, 18, 21, 23, 24, 25  …  88, 93, 94, 96, 98, 102, 104, 112, 114, 115], [7, 8, 12, 13, 14, 16, 26, 27, 29, 30  …  99, 100, 101, 103, 105, 106, 107, 109, 110, 113], [1, 4, 5, 6, 11, 17, 19, 20, 22, 32  …  79, 81, 83, 84, 85, 90, 91, 95, 108, 111]]
 [[2, 6, 11, 16, 18, 20, 22, 24, 25, 29  …  90, 97, 101, 103, 106, 108, 109, 110, 111, 113], [5, 8, 12, 13, 15, 21, 23, 26, 27, 35  …  82, 91, 94, 95, 98, 99, 102, 104, 114, 115], [1, 3, 4, 7, 9, 10, 14, 17, 19, 28  …  87, 88, 89, 92, 93, 96, 100, 105, 107, 112]]
 [[1, 5, 7, 10, 18, 21, 22, 23, 26, 29  …  86, 87, 94, 95, 97, 103, 106, 109, 111, 112], [2, 4, 8, 11, 12, 13, 15, 17, 19, 24  …  82, 89, 91, 92, 98, 99, 100, 107, 113, 114], [3, 6, 9, 14, 16, 20, 25, 28, 33, 34  …  90, 93, 96, 101, 102, 104, 105, 108, 110, 115]]
 [[3, 7, 13, 16, 21, 23, 26, 28, 31, 34  …  85, 86, 93, 98, 101, 104, 105, 106, 107, 114], [1, 2, 4, 5, 6, 8, 9, 11, 12, 14  …  83, 88, 91, 92, 96, 108, 109, 110, 111, 113], [10, 15, 17, 19, 20, 24, 29, 30, 33, 36  …  90, 94, 95, 97, 99, 100, 102, 103, 112, 115]]
 [[4, 5, 7, 14, 15, 16, 18, 19, 22, 23  …  83, 84, 89, 99, 100, 104, 105, 106, 107, 109], [1, 6, 9, 10, 11, 17, 28, 31, 33, 35  …  88, 90, 93, 94, 96, 97, 101, 102, 110, 112], [2, 3, 8, 12, 13, 20, 21, 24, 25, 26  …  91, 92, 95, 98, 103, 108, 111, 113, 114, 115]]
 [[2, 7, 15, 16, 20, 22, 23, 27, 29, 30  …  91, 92, 94, 95, 98, 101, 102, 106, 108, 115], [4, 11, 13, 14, 18, 19, 21, 24, 32, 33  …  84, 86, 96, 97, 103, 107, 110, 111, 112, 114], [1, 3, 5, 6, 8, 9, 10, 12, 17, 25  …  87, 89, 90, 93, 99, 100, 104, 105, 109, 113]]
 [[2, 7, 10, 12, 14, 16, 22, 24, 29, 30  …  96, 99, 100, 101, 104, 105, 106, 111, 112, 115], [1, 3, 4, 5, 11, 13, 15, 18, 21, 26  …  91, 92, 93, 94, 95, 97, 102, 103, 110, 113], [6, 8, 9, 17, 19, 20, 23, 25, 27, 28  …  77, 78, 84, 88, 89, 98, 107, 108, 109, 114]]
 [[2, 3, 6, 8, 13, 15, 16, 20, 24, 25  …  81, 90, 91, 92, 94, 104, 106, 109, 110, 115], [1, 4, 9, 10, 11, 12, 14, 17, 28, 30  …  95, 98, 99, 100, 102, 105, 111, 112, 113, 114], [5, 7, 18, 19, 21, 22, 23, 26, 31, 32  …  80, 85, 86, 87, 96, 97, 101, 103, 107, 108]]
 [[4, 8, 10, 14, 22, 24, 30, 33, 36, 37  …  95, 96, 98, 100, 104, 107, 108, 110, 112, 113], [3, 5, 6, 9, 12, 13, 17, 18, 19, 20  …  73, 78, 85, 88, 92, 93, 97, 103, 114, 115], [1, 2, 7, 11, 15, 16, 21, 23, 28, 29  …  86, 87, 91, 99, 101, 102, 105, 106, 109, 111]]

Grid-search

nlvdis = [5; 10; 15]; metric = [:mah]   # here only Mahalanobis distance is considered
h = [1; 2; 4; 6; Inf]; k = [30; 50; 100]  
nlv = 0:15
pars = mpar(nlvdis = nlvdis, metric = metric, h = h, k = k) 
length(pars[1])

45

model = lwplsr()
rescv = gridcv(model, Xtrain, ytrain; segm, score = rmsep, pars, nlv)
@names rescv 
res = rescv.res
res_rep = rescv.res_rep

group = string.("nvldis=", res.nlvdis, " h=", res.h, " k=", res.k)
plotgrid(res.nlv, res.y1, group; step = 2, xlabel ="Nb. LVs", ylabel = "RMSEP-Val").f

Selection of the best parameter combination

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

Final prediction (Test) using the optimal model

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

63×1 Matrix{Float64}:
 29.233919470654392
  1.3973931780971185
  4.794708197493841
 11.214014955076081
 15.40594497215125
 20.431938204221815
 27.011452750213174
 47.98443636636166
  5.256365240224742
  1.869392428285745
  ⋮
 17.70770463351476
 19.115522708954746
 21.597608245182037
 27.614157967846275
 27.667973750722332
 29.66398089648304
 33.922975891981565
 36.54196651755939
 41.810372978453216

Generalization error

rmsep(pred, ytest)

1×1 Matrix{Float64}:
 1.0303218875423956

Plotting predictions vs. observed data

plotxy(pred, ytest; size = (500, 400), color = (:red, .5), bisect = true, title = string("Test set - variable ", nam), 
    xlabel = "Prediction", ylabel = "Observed").f

---