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

Grid-search

The best syntax to use function gridcv for LV-based functions (e.g. rr, krr, etc.) is to set parameter lb outside argument pars defining the grid (in general with function mpar). In that case, the computation time is reduced [See the naïve and non-optimized syntax at the end of this script]. This is the same principle when definining the parameter nlv in LV-based functions (eg. plskern, kplsr, lwplsr, etc.).

lb = 10.0.^(-15:.1:3)
model = rr()
rescv = gridcv(model, Xtrain, ytrain; segm, score = rmsep, lb)
@names rescv 
res = rescv.res
res_rep = rescv.res_rep

loglb = round.(log.(10, res.lb), digits = 3)
plotgrid(loglb, res.y1; step = 2, xlabel ="lambda (log scale)", ylabel = "RMSEP-Val").f

f, ax = plotgrid(loglb, res.y1; step = 2, xlabel = "Nb. LVs", ylabel = "RMSEP-CV")
for i = 1:rep, j = 1:K
    zres = res_rep[res_rep.rep .== i .&& res_rep.segm .== j, :]
    lines!(ax, loglb, zres.y1; color = (:grey, .2))
end
lines!(ax, loglb, res.y1; color = :red, linewidth = 1)
f

If other parameters have to be defined in the grid, they have to be set in argument pars built with function mpar, such as in the example below.

pars = mpar(scal = [false; true])
lb = 10.0.^(-15:.1:3)
#lb = logrange(1e-15, 1e3, 50)   # alternative syntax
model = rr()
res = gridcv(model, Xtrain, ytrain; segm, score = rmsep, pars, lb).res

loglb = round.(log.(10, res.lb), digits = 3)
plotgrid(loglb, res.y1, res.scal; step = 2, xlabel ="lambda (log scale)", 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 = rr(nlv = res.lb[u], scal = res.scal[u])
fit!(model, Xtrain, ytrain)
pred = predict(model, Xtest).pred

63×1 Matrix{Float64}:
 27.876444002998458
  3.5904081173640385
  2.32548879045871
  9.169771876981347
 13.228383160517499
 20.88245631069842
 21.875931002536706
 47.68294255512646
  7.591807587883839
  2.0283951227110926
  ⋮
 18.070871456149778
 19.35569672742449
 21.75703948647831
 27.713990081253858
 27.393213664847824
 34.13829277337615
 33.47241171228406
 38.85986950484511
 44.4835756075456

Generalization error

rmsep(pred, ytest)

1×1 Matrix{Float64}:
 2.8631045531296424

Plotting predictions vs. observed data

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

Naïve (not time-efficient) syntax to use gridcv for LV-based functions

Parameter lb can also be set in argument pars (wich is the generic approach to define the grid). This is strictly equivalent (gives the same results) but the computations are slower.

pars = mpar(lb = 10.0.^(-15:.1:3))
model = rr()
gridcv(model, Xtrain, ytrain; segm, score = rmsep, pars).res

pars = mpar(lb = 10.0.^(-15:.1:3), scal = [false; true])
model = rr()
gridcv(model, Xtrain, ytrain; segm, score = rmsep, pars).res

---