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)  # the classes 
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)

Since X-columns are highly colinear, the Fda requires a regularization. It can be reached for instance 1) by preliminary dimension reduction, or 2) by using a ridge regularization.

1) Fda on preliminary computed Pca scores

model0 = pcasvd(nlv = 10)
fit!(model0, Xtrain)
@names model0
@names model0.fitm

(:T, :V, :sv, :xmeans, :xscales, :weights, :niter, :par)

@head Ttrain_pca = model0.fitm.T

3×10 Matrix{Float64}:
 -0.00669306  -0.00650046  -0.00570717  …  -0.000492699   0.000134327
 -0.00387969  -0.0068879   -0.00728328      0.0003361    -0.000581481
  0.0022825    0.00305238   0.00109557     -0.00164726   -0.0014889
... (323, 10)

model = fda(nlv = 2)
#model = fdasvd(nlv = 2)     # alternative algorithm (same result)
fit!(model, Ttrain_pca, ytrain) 
fitm = model.fitm
@names fitm

(:T, :V, :Tcenters, :eig, :sstot, :W, :xmeans, :xscales, :weights, :lev, :ni, :par)

@head Ttrain = model.fitm.T

3×2 Matrix{Float64}:
  2.41503  -1.46642
  1.41666  -1.4342
 -1.54614  -0.344225
... (323, 2)

lev = fitm.lev
nlev = length(lev)

3

ct = fitm.Tcenters

3×2 Matrix{Float64}:
 -2.1618    -0.223723
  1.49851   -1.67024
  0.791993   0.694047

f, ax = plotxy(Ttrain[:, 1], Ttrain[:, 2], ytrain; title = "Fda", ellipse = true)
scatter!(ax, ct[:, 1], ct[:, 2], markersize = 10, color = :red)
f

color = cgrad(:lightrainbow, nlev; categorical = true, alpha = .5)
f, ax = plotxy(Ttrain[:, 1], Ttrain[:, 2], ytrain; color, title = "Fda", ellipse = true)
scatter!(ax, ct[:, 1], ct[:, 2]; markersize = 10, color = :red)
f

Projection of Xtest

@head Ttest_pca = transf(model0, Xtest)

3×10 Matrix{Float64}:
 -0.00881891  -0.00355513    0.00143197   …   0.000140379  -0.000265216
  0.00647425   0.00296032   -0.000761146     -0.00113056   -0.00160408
  0.00351213   0.000177464  -0.000182268     -0.00140921   -0.00232458
... (162, 10)

@head Ttest = transf(model, Ttest_pca)

3×2 Matrix{Float64}:
  1.76619  -0.105807
 -2.41847  -0.524465
 -2.57702  -0.907449
... (162, 2)

i = 1  # class "i" in test
color = cgrad(:lightrainbow, nlev; categorical = true, alpha = .5)
f, ax = plotxy(Ttrain[:, 1], Ttrain[:, 2], ytrain; color, title = "Fda", ellipse = true)
scatter!(ax, ct[:, 1], ct[:, 2]; markersize = 10, color = :red)
s = ytest .== lev[i]
zT = Ttest[s, :]
scatter!(ax, zT[:, 1], zT[:, 2]; markersize = 10, color = :blue)
f

2) Fda directly on X but using a ridge regularization

Note: If the ridge parameter lb is too small (e.g. here from around 1e-10), the model can be overfitted and predictions of new observations are misleading.

lb = 1e-5
model = fda(; nlv = 2, lb)
#model = fdasvd(; nlv = 2, lb)
fit!(model, Xtrain, ytrain)
fitm = model.fitm

Jchemo.Fda([1.7829806394884713 -0.8923344270339282; 1.3203919349887605 -1.087112071286717; … ; -0.7042782276724117 0.523450171819936; -0.4631097344316628 1.0465334079722217], [5.171712666455356 -9.508299170939283; 5.153707863984817 -9.06078711679522; … ; 3.032245770823203 21.347000304208176; -2.8066707021640216 21.53541187066701], [-1.3213215617847378 -0.25815326934600685; 1.2872849511152809 -0.8352132611860519; 0.3595461012935211 0.43465430874861954], [0.8946608889469787, 0.2392545735461798, 8.113420527620786e-17, 7.522381965283366e-17, 1.2357532069459468e-17, 1.2357532069459468e-17, -5.640384069100635e-17, -5.640384069100635e-17, 3.953552944990857e-17, 3.953552944990857e-17  …  1.7581207434812442e-20, 2.101291977425629e-20, -1.8276419310707133e-20, 1.1596325005033801e-20, -3.75017290833722e-21, -3.75017290833722e-21, -9.096083385580036e-21, 8.960398722442791e-21, 5.4950099649394124e-21, 5.4950099649394124e-21], 1.1339154624931584, [1.0069729150340446e-5 6.666992825610627e-8 … 2.0918098869706948e-8 2.1162650478714127e-8; 6.666992825610627e-8 1.0063873539470998e-5 … 2.008643294430349e-8 2.0247774844963594e-8; … ; 2.0918098869706948e-8 2.008643294430349e-8 … 1.009178260196924e-5 9.092935251310007e-8; 2.1162650478714127e-8 2.0247774844963594e-8 … 9.092935251310007e-8 1.0095765380348832e-5], [-0.00019330623941241694, -0.0001415769474021643, -5.770661245116858e-5, 4.723558902506759e-5, 0.00014758462542250662, 0.00021809097323352518, 0.0002341317111996633, 0.000252799982295447, 0.0002796860235471441, 0.0003040280867910269  …  -0.0018442000247695933, -0.0019172866156863932, -0.0019496441042941074, -0.0019473165883640694, -0.0018947013042759696, -0.0016716820169876765, -0.0013441905794144805, -0.0009918653793508437, -0.0006005900823885629, -0.00025528401464108526], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], Jchemo.Weight{Float64}([0.0030959752321981413, 0.0030959752321981417, 0.0030959752321981417, 0.0030959752321981413, 0.0030959752321981413, 0.0030959752321981417, 0.0030959752321981417, 0.0030959752321981417, 0.0030959752321981417, 0.0030959752321981417  …  0.0030959752321981417, 0.0030959752321981417, 0.0030959752321981417, 0.0030959752321981417, 0.0030959752321981417, 0.0030959752321981413, 0.0030959752321981413, 0.0030959752321981413, 0.0030959752321981413, 0.0030959752321981413]), ["Cereal and grass forages", "Forage trees", "Legume forages"], [100, 56, 167], Jchemo.ParFda(2, 1.0e-5, :prop, false))

lev = fitm.lev
nlev = length(lev)

3

ct = fitm.Tcenters

3×2 Matrix{Float64}:
 -1.32132   -0.258153
  1.28728   -0.835213
  0.359546   0.434654

@head Ttrain = fitm.T

3×2 Matrix{Float64}:
  1.78298   -0.892334
  1.32039   -1.08711
 -0.881565  -0.00540248
... (323, 2)

f, ax = plotxy(Ttrain[:, 1], Ttrain[:, 2], ytrain; title = "Fda", ellipse = true)
scatter!(ax, ct[:, 1], ct[:, 2]; markersize = 10, color = :red)
f

color = cgrad(:lightrainbow, nlev; categorical = true, alpha = .7)
f, ax = plotxy(Ttrain[:, 1], Ttrain[:, 2], ytrain; color, title = "Fda", ellipse = true)
scatter!(ax, ct[:, 1], ct[:, 2]; markersize = 10, color = :red)
f

Projection of Xtest

@head Ttest = transf(model, Xtest)

3×2 Matrix{Float64}:
  1.47243    0.224802
 -1.39627   -0.2853
 -0.995825  -0.334509
... (162, 2)

i = 1  # class "i" in test
color = cgrad(:lightrainbow, nlev; categorical = true, alpha = .7)
f, ax = plotxy(Ttrain[:, 1], Ttrain[:, 2], ytrain; color, title = "Projection test-class $(lev[i]) (blue points)", 
    ellipse = true)
scatter!(ax, ct[:, 1], ct[:, 2]; markersize = 10, color = :red)
s = ytest .== lev[i]
zT = Ttest[s, :]
scatter!(ax, zT[:, 1], zT[:, 2]; markersize = 10, color = :blue)
f

---