using Jchemo, JchemoData
using JLD2, DataFrames, 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)

X data

wlst = names(X) 
wl = parse.(Int, wlst)

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

The X-data are already preprocessed by SNV and second derivation.

plotsp(X, wl; xlabel = "Wavelength (nm)", ylabel = "Absorbance").f

Y data

typ = Y.typ
tab(typ)

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

test = Y.test  # training/test (0/1) observations
tab(test)
freqtable(Y.typ, test)

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

namy = names(Y)[1:2]
summ(Y[:, namy]).res

summ(Y[:, namy], test)

Class: 0
2×7 DataFrame
 Row │ variable  mean     std      min      max      n      nmissing
     │ Symbol    Float64  Float64  Float64  Float64  Int64  Int64
─────┼───────────────────────────────────────────────────────────────
   1 │ dm         92.022    1.992   84.25   96.27      323         0
   2 │ ndf        50.988   13.036   18.425  83.3012    323         0


Class: 1
2×7 DataFrame
 Row │ variable  mean     std      min      max      n      nmissing
     │ Symbol    Float64  Float64  Float64  Float64  Int64  Int64
─────┼───────────────────────────────────────────────────────────────
   1 │ dm         92.257    1.976   80.903  96.18      162         0
   2 │ ndf        50.233   13.56    22.959  85.7235    162         0

A given variable

j = 2
nam = namy[2]
y = Float64.(Y[:, nam])
s = test .== 0
ytrain = y[s] 
ytest = rmrow(y, s)

162-element Vector{Float64}:
 37.579963135639154
 56.99851095511593
 47.21391972672929
 25.309323906319047
 54.28931875525652
 80.29120067524794
 57.55875502858353
 32.66967330326697
 61.936710276193665
 46.46783498421683
  ⋮
 48.90913691567336
 69.05194981637266
 45.154438712841326
 41.50911821540868
 52.42695116229647
 22.959435724677174
 35.71112405578161
 55.53583902046241
 34.64568924391239

f = Figure(size = (400, 300))
ax = Axis(f[1, 1]; xlabel = uppercase(nam), ylabel = "Nb. observations")
hist!(ax, ytrain; bins = 50, label = "Train")
hist!(ax, ytest; bins = 50, label = "Test")
axislegend(position = :rt)
f

f = Figure(size = (500, 400))
offs = [20; 0]
ax = Axis(f[1, 1]; xlabel = uppercase(nam),  ylabel = "Nb. observations", yticks = (offs, ["Train" ; "Test"]))
hist!(ax, ytrain; offset = offs[1], bins = 50)
hist!(ax, ytest; offset = offs[2], bins = 50)
f

f = Figure(size = (400, 300))
ax = Axis(f[1, 1]; xlabel = uppercase(nam), ylabel = "Density")
bdw = 1
density!(ax, ytrain; bandwidth = bdw, color = :blue, label = "Train")
density!(ax, ytest; bandwidth = bdw, color = (:red, .5), label = "Test")
axislegend(position = :rt)  
f

f = Figure(size = (500, 400))
offs = [.1; 0]
ax = Axis(f[1, 1]; xlabel = uppercase(nam), ylabel = "Density", yticks = (offs, ["Train" ; "Test"]))
bdw = 1
density!(ax, ytrain; offset = offs[1], color = (:slategray, 0.5), bandwidth = bdw)
density!(ax, ytest; offset = offs[2], color = (:slategray, 0.5), bandwidth = bdw)
f

f = Figure(size = (400, 300))
ax = Axis(f[1, 1]; xticks = (0:1, ["Train", "Test"]), xlabel = "Group", ylabel = uppercase(nam))
boxplot!(ax, test, y; width = .3, show_notch = true)
f

lev = mlev(test)
nlev = length(lev)
tsp = .5
#colm = [(:blue, tsp), (:orange, tsp)]
colm = cgrad(:Dark2_5; categorical = true, alpha = .8)[1:nlev]
cols = colm[indexin(test, unique(test))]
#cols = (:red, .5)
f = Figure(size = (600, 300))
ax = Axis(f[1, 1]; xticks = (0:1, ["Train", "Test"]), ylabel = uppercase(nam))
rainclouds!(ax, test, y; clouds = hist, markersize = 5, color = cols, gap = .3)
f

---