using Jchemo, JchemoData
using JLD2, CairoMakie
using FreqTables

Data importation

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

(:X, :Y)

X = dat.X 
@head X

... (485, 700)

Y = dat.Y
@head Y

... (485, 4)

y = Y.typ
test = Y.test
tab(y)

OrderedCollections.OrderedDict{String, Int64} with 3 entries:
  "Cereal and grass forages" => 160
  "Forage trees"             => 101
  "Legume forages"           => 224

freqtable(y, test)

3×2 Named Matrix{Int64}
             Dim1 ╲ Dim2 │   0    1
─────────────────────────┼─────────
Cereal and grass forages │ 100   60
Forage trees             │  56   45
Legume forages           │ 167   57

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

700-element Vector{Int64}:
 1100
 1102
 1104
 1106
 1108
 1110
 1112
 1114
 1116
 1118
    ⋮
 2482
 2484
 2486
 2488
 2490
 2492
 2494
 2496
 2498

Note:: X-data are already preprocessed (SNV + Savitsky-Golay 2nd deriv).

Split Tot to Train/Test

The model is fitted on Train, and the generalization error is estimated on Test. In this example, Train is already defined in variable typ of the dataset, and Test is defined by the remaining samples. But Tot could also be split a posteriori, for instance by sampling (random, systematic or any other designs). See for instance functions samprand, sampsys, etc.

s = Bool.(test)
Xtrain = rmrow(X, s)
ytrain = rmrow(y, s)
Xtest = X[s, :]
ytest = y[s]
ntot = nro(X)
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = ntot, ntrain, ntest)

(ntot = 485, ntrain = 323, ntest = 162)

tab(ytrain)

OrderedCollections.OrderedDict{String, Int64} with 3 entries:
  "Cereal and grass forages" => 100
  "Forage trees"             => 56
  "Legume forages"           => 167

tab(ytest)

OrderedCollections.OrderedDict{String, Int64} with 3 entries:
  "Cereal and grass forages" => 60
  "Forage trees"             => 45
  "Legume forages"           => 57

K-fold CV

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

1-element Vector{Vector{Vector{Int64}}}:
 [[3, 5, 7, 9, 12, 14, 15, 16, 17, 21  …  286, 290, 291, 292, 294, 296, 299, 300, 317, 322], [1, 8, 10, 13, 20, 22, 23, 26, 27, 29  …  308, 310, 311, 312, 314, 316, 318, 319, 320, 321], [2, 4, 6, 11, 18, 19, 25, 31, 33, 36  …  289, 293, 298, 303, 304, 307, 309, 313, 315, 323]]

nlvdis = [15; 25]; metric = [:mah]
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])

30

model = lwplsrda()
res = gridcv(model, Xtrain, ytrain; segm, score = merrp, pars, nlv).res

group = string.("nvldis=", res.nlvdis, " h=", res.h, " k=", res.k)
plotgrid(res.nlv, res.y1, group; step = 2, xlabel = "Nb. LVs", ylabel = "ERRP-CV", leg_title = "Continuum").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 = lwplsrda(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

162×1 Matrix{String}:
 "Legume forages"
 "Cereal and grass forages"
 "Cereal and grass forages"
 "Legume forages"
 "Cereal and grass forages"
 "Cereal and grass forages"
 "Legume forages"
 "Forage trees"
 "Forage trees"
 "Forage trees"
 ⋮
 "Cereal and grass forages"
 "Cereal and grass forages"
 "Forage trees"
 "Forage trees"
 "Cereal and grass forages"
 "Legume forages"
 "Legume forages"
 "Legume forages"
 "Legume forages"

Generalization error

errp(pred, ytest)

1×1 Matrix{Float64}:
 0.08641975308641975

merrp(pred, ytest)

1×1 Matrix{Float64}:
 0.08849902534113059

cf = conf(pred, ytest)
@names cf

(:cnt, :pct, :A, :Apct, :diagpct, :accpct, :lev)

cf.cnt

cf.pct

---