using Jchemo, JchemoData
using JLD2, CairoMakie
using FreqTables

RR-DA is a alternative to PLSR-DA, where ridge regression is used instead of PLSR to predict the Y-dummy table.

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

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

model = rrda(; lb = 1e-5)
fit!(model, Xtrain, ytrain)
@names model

(:algo, :fitm, :kwargs)

fitm = model.fitm
@names fitm

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

typeof(fitm.fitm)

Jchemo.Rr

@names fitm.fitm

(:V, :TtDY, :sv, :xmeans, :xscales, :ymeans, :weights, :par)

res = predict(model, Xtest)
@names res
@head pred = res.pred
@head res.posterior   # predicted posterior probabilities

3×1 Matrix{String}:
 "Legume forages"
 "Cereal and grass forages"
 "Cereal and grass forages"
... (162, 1)
 
3×3 Matrix{Float64}:
 0.277612  -0.212171   0.934559
 0.647358   0.072594   0.280048
 0.611353  -0.0854871  0.474134
... (162, 3)

@head predict(model, Xtest).pred

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

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

2-element Vector{Matrix{String}}:
 ["Legume forages"; "Legume forages"; … ; "Legume forages"; "Legume forages";;]
 ["Legume forages"; "Legume forages"; … ; "Legume forages"; "Legume forages";;]

errp(pred, ytest)

1×1 Matrix{Float64}:
 0.08641975308641975

merrp(pred, ytest)

1×1 Matrix{Float64}:
 0.096588693957115

Confusion matrix

cf = conf(pred, ytest)
@names cf

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

cf.cnt

cf.pct

plotconf(cf).f

---