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   # response variable (class membership)
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

Model fitting

'PLSDA' is a family of methods, for instance PLS-MLR-DA (or PLSR-DA), PLS-LDA, PLS-QDA or PLS-KDE-DA.

Note:

• In the functions below, the default value for argument prior is set to :prop. If the classes to discriminate are highly unbalanced, it is recommended to set prior = :unif to avoid bias in the predictions (especially when using function plsrda that is very sensitive to unbalanced classes). See the help of the respective functions for more details.

• For unbalanced classes, it is also recommended to use merrp instead of errp (see the respectve help pages) to compute the overall prediction-error rates in CV processes and/or on test sets.

nlv = 15
model = plsrda(; nlv) 
#model = plslda(; nlv) 
#model = plsqda(; nlv) 
#model = plsqda(; nlv, alpha = 0.5)   # 'alpha' = continuum parameter
#model = plskdeda(; nlv) 
#model = plskdeda(; nlv, a = .5)    # 'a' = bandwidth parameter (see also parameter 'h')

Jchemo.JchemoModel{typeof(plsrda), Base.Pairs{Symbol, Int64, Tuple{Symbol}, @NamedTuple{nlv::Int64}}}(Jchemo.plsrda, nothing, Base.Pairs(:nlv => 15))

fit!(model, Xtrain, ytrain)
@names model

(:algo, :fitm, :kwargs)

fitm = model.fitm
@names fitm

(:fitm, :lev, :ni, :par)

typeof(fitm.fitm)

Jchemo.Plsr

@names fitm.fitm

(:T, :V, :R, :W, :C, :TT, :xmeans, :xscales, :ymeans, :yscales, :weights, :niter, :par)

res = predict(model, Xtest)
@names res

(:pred, :posterior)

@head pred = res.pred

3×1 Matrix{String}:
 "Legume forages"
 "Cereal and grass forages"
 "Cereal and grass forages"
... (162, 1)

@head res.posterior   # predicted posterior probabilities

3×3 Matrix{Float64}:
 -0.00450279   0.0717177    0.932785
  0.912764    -0.0344078    0.121644
  1.07725     -0.00150365  -0.075744
... (162, 3)

@head predict(model, Xtest; nlv = 2).pred

3×1 Matrix{String}:
 "Legume forages"
 "Cereal and grass forages"
 "Cereal and grass forages"
... (162, 1)

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

3-element Vector{Matrix{String}}:
 ["Legume forages"; "Legume forages"; … ; "Legume forages"; "Legume forages";;]
 ["Legume forages"; "Cereal and grass forages"; … ; "Cereal and grass forages"; "Cereal and grass forages";;]
 ["Legume forages"; "Cereal and grass forages"; … ; "Cereal and grass forages"; "Legume forages";;]

errp(pred, ytest)

1×1 Matrix{Float64}:
 0.1111111111111111

merrp(pred, ytest)

1×1 Matrix{Float64}:
 0.11910331384015593

Confusion matrix

cf = conf(pred, ytest)
@names cf

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

cf.cnt

cf.pct

cf.diagpct

cf.accpct

plotconf(cf).f

plotconf(cf; cnt = false).f

plotconf(cf; ptext = false).f

---