using Jchemo, JchemoData
using JLD2, CairoMakie

Data importation

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

(:X, :Y)

X = dat.X
@head X

... (178, 100)

Y = dat.Y 
@head Y

... (178, 4)

typ = Y.typ
tab(typ)

OrderedCollections.OrderedDict{String, Int64} with 3 entries:
  "test"  => 31
  "train" => 115
  "val"   => 32

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

100-element Vector{Float64}:
  850.0
  852.0
  854.0
  856.0
  858.0
  860.0
  862.0
  864.0
  866.0
  868.0
    ⋮
 1032.0
 1034.0
 1036.0
 1038.0
 1040.0
 1042.0
 1044.0
 1046.0
 1048.0

Preprocessing

model1 = snv()
model2 = savgol(npoint = 15, deriv = 2, degree = 3)
model = pip(model1, model2)
fit!(model, X)
@head Xp = transf(model, X) 
#plotsp(Xp, wl; xlabel = "Wavelength (nm)", ylabel = "Absorbance").f

3×100 Matrix{Float64}:
 0.000397076  0.000495203  0.000598623  …  0.00827138  0.00917311  0.00946072
 0.00242055   0.00244366   0.00234233      0.00580631  0.00689249  0.00749316
 0.0011927    0.00122721   0.00120098      0.0101019   0.0108142   0.0108444
... (178, 100)

Split Tot to Train/Test

s = typ .== "train"
Xtrain = Xp[s, :] 
Ytrain = Y[s, :]
Xtest = rmrow(Xp, s)
Ytest = rmrow(Y, s)
ntrain = nro(Xtrain)
ntest = nro(Xtest)
ntot = ntrain + ntest
(ntot = ntot, ntrain, ntest)

(ntot = 178, ntrain = 115, ntest = 63)

Working response y

namy = names(Y)[1:3]
j = 2  
nam = namy[j]    # work on the j-th y-variable
ytrain = Ytrain[:, nam]
ytest = Ytest[:, nam]

63-element Vector{Float64}:
 29.8
  1.4
  4.6
 11.0
 17.0
 22.4
 27.9
 46.5
  6.1
  2.0
  ⋮
 18.1
 19.4
 24.8
 27.2
 28.4
 31.3
 33.8
 35.5
 42.5

Model fitting

gamma = 100
nlv = 15
model = kplsr(; gamma, nlv)   # For direct KPLSR (Bennett & Embrechts, 2003), use as below:
#model = dkplsr(; gamma, nlv)  
fit!(model, Xtrain, ytrain)
@names model
fitm = model.fitm
@names fitm

(:X, :Kt, :T, :C, :U, :R, :DKt, :vtot, :xscales, :ymeans, :yscales, :weights, :iter, :kwargs, :par)

Predictions

@head pred = predict(model, Xtest).pred

3×1 Matrix{Float64}:
 30.27263234256255
  1.2491123455562665
  4.4710390881702065
... (63, 1)

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

3×1 Matrix{Float64}:
 29.32135166200425
  2.162005442116719
  5.685860978803078
... (63, 1)

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

3-element Vector{Matrix{Float64}}:
 [18.358260869565218; 18.358260869565218; … ; 18.358260869565218; 18.358260869565218;;]
 [26.864333587051988; 4.1217670557698955; … ; 36.33040026195786; 45.001339268847815;;]
 [29.32135166200425; 2.162005442116719; … ; 39.208838218241105; 46.42134811761902;;]

rmsep(pred, ytest)

1×1 Matrix{Float64}:
 0.5511112753976203

bias(pred, ytest)

1×1 Matrix{Float64}:
 -0.15413101976451102

mse(pred, ytest)

@head r = residreg(pred, ytest) # residuals

3×1 Matrix{Float64}:
 -0.47263234256254805
  0.15088765444373342
  0.1289609118297932
... (63, 1)

Plotting predictions vs. observed data

plotxy(pred, ytest; size = (500, 400), color = (:red, .5), bisect = true, title = string("Test set - variable ", nam), 
    xlabel = "Prediction", ylabel = "Observed").f

Residuals

plotxy(ytest, r; size = (500, 400), color = (:red, .5), zeros = true, title = string("Test set - variable ", nam), 
    xlabel = "Observed (Test)", ylabel = "Residuals").f

Adding a smoothing

zpred = vec(pred)
zr = vec(r)

63-element Vector{Float64}:
 -0.47263234256254805
  0.15088765444373342
  0.1289609118297932
 -0.010875159198107553
  1.4489775124476338
  0.290389476158623
  0.2844942468389391
 -0.06310162831640298
  0.22063704374842352
  0.24732308729945274
  ⋮
  0.12741438205125277
 -0.1781680431607029
  0.9736514052289138
 -0.3613955588708784
 -0.02422859343892725
  1.34464900017079
  0.09247884206364887
 -0.5730256984462514
  0.16607155085765157

Predictions

model_lo = loessr(span = 1/2) 
fit!(model_lo, zpred, ytest)
pred_lo = predict(model_lo, sort(zpred)).pred
f, ax = plotxy(zpred, ytest; size = (500, 400), bisect = true, title = string("Test set - variable ", nam), 
    xlabel = "Predicted", ylabel = "Observed")
lines!(ax, sort(zpred), vec(pred_lo); color = :red)
f

Residuals

model_lo = loessr(span = 1/2) 
fit!(model_lo, zpred, zr)
pred_lo = predict(model_lo, sort(zpred)).pred
f, ax = plotxy(zpred, zr; size = (500, 400), title = string("Test set - variable ", nam), 
    xlabel = "Predictions", ylabel = "Residuals")
lines!(ax, sort(zpred), vec(pred_lo); color = :red)
hlines!(ax, 0; color = :grey)
f

---