Introduction

\[ \text{y} = {\color{red} f}(\text{DRS spectrum}) + \epsilon \]

What are we trying to do?

\[\text{y} = {\color{red} f}({\color{dodgerblue} g}(\text{DRS spectrum})) + \epsilon\]

What are we trying to do?

A simplified view of the chemometric modelling pipeline:

What are we trying to do?

A simplified view of the prediction pipeline:

Why open-source chemometrics? [1/4]

Proprietary chemometrics software…

Cost and access

• Licence fees

• Institutional dependency

• Risk of discontinuation

• Limited user base → limited community support

Why open-source chemometrics? [2/4]

Proprietary chemometrics software…

Transparency and trust

• Closed source

• Bugs hard to spot or report

• Reproducibility limited

• No audit trail

• Less amenable to peer scrutiny

Why open-source chemometrics? [3/4]

Proprietary chemometrics software…

Flexibility and integration

• Poor interoperability

• No version control

• Cannot script end-to-end pipelines

• Disconnected from modern AI/ML tools

• Model deployment requires additional licences

• Slow response time to new developments

Why open-source chemometrics? [4/4]

Chemometrics in R

https://cran.r-project.org/web/views/ChemPhys.html

Chemometrics in R: prospectr

Released in 2013…

https://CRAN.R-project.org/package=prospectr

Chemometrics in R: resemble

Released in 2014…

https://CRAN.R-project.org/package=resemble

Chemometrics in R: proximetricsR

[First implemented: 2020] [Released today🙂] [Soon on CRAN]

https://github.com/l-ramirez-lopez/proximetricsR

Chemometrics with Python

https://scikit-learn.org/

Comprehensive machine learning framework:

Chemometrics with Julia

https://julialang.org/

High-performance numerical computing for chemometrics:

What are we trying to do?

\[\text{y} = {\color{red} f}({\color{dodgerblue} g}(\text{DRS spectrum})) + \epsilon\]

The data [1/1]

A dataset of soil samples from Spain

• Target properties: organic carbon (OC) and organic matter

• Spectral data/predictors: NIR (1350 - 2550 nm or 7407.4 - 3921.5 cm\(^{-1}\)) collected with a ProxiScout (BUCHI) sensor

• 364 samples for calibration of models

• 150 samples for validation of models

The data [1/2]

In the folder /data/local-samples-spain/ you will find:

• local-spectra-cal.xlsx: raw spectra for the calibration set

• local-properties-cal.xlsx: reference properties for the calibration set

• local-spectra-val.xlsx: raw spectra for the validation set

• local-properties-val.xlsx: reference properties for the validation set

Before we start: organise your workspace

Create a folder manually on your computer and organise it like this:

my_folder/
├── data/
│   └── local-samples-spain/
│       ├── local-spectra-cal.xlsx
│       ├── local-properties-cal.xlsx
│       ├── local-spectra-val.xlsx
│       └── local-properties-val.xlsx
├── R/
└── outputs/

Set the working directory in R

Tell R to use your folder as the working directory:

my_folder <- "C:/path/to/my_folder"
setwd(my_folder)

Verify (you should get the path to your folder):

getwd()

For R we use RStudio

Install and load the proximetricsR package in R

Install proximetricsR from GitHub:

remotes::install_github("l-ramirez-lopez/proximetricsR")

Install and load the proximetricsR package in R

Install proximetricsR from GitHub:

remotes::install_github("l-ramirez-lopez/proximetricsR")

Once installed, load it:

library(proximetricsR)

Load the calibration data

Read the calibration data (samples with spectra and their corresponding reference values):

mlocal_cal <- proxiscout_read_data(
  file = "data/local-samples-spain/local-spectra-cal.xlsx", 
  references_file = "data/local-samples-spain/local-properties-cal.xlsx"
)

Inspect the loaded caliration data

Check how many rows and variables are in the loaded data:

dim(mlocal_cal)

[1] 364   6

Check the name of the variables in the loaded data:

colnames(mlocal_cal)

[1] "sampleName" "deviceId"   "ref"        "OM"         "OC"        
[6] "spc"       

View the table:

View(mlocal_cal)

View the spectra only:

View(mlocal_cal$spc)

Structure of the loaded data (summary):

Now load the validation data

Read the validation data (samples with spectra and their corresponding reference values):

mlocal_val <- proxiscout_read_data(
  file = "data/local-samples-spain/local-spectra-val.xlsx", 
  references_file = "data/local-samples-spain/local-properties-val.xlsx"
)

Check the dimensions of the validation data:

dim(mlocal_val)

[1] 150   6

colnames(mlocal_val)

[1] "sampleName" "deviceId"   "ref"        "OM"         "OC"        
[6] "spc"       

Obtain the vector of wavenumbers

The wavenumbers are stored as the column names of the spectra matrix in the calibration data. We can extract them and convert them to numeric values:

wavs <- as.numeric(colnames(mlocal_cal$spc))
wavs

The spectra

matplot(
  wavs,
  t(mlocal_cal$spc), 
  type = "l", 
  lty = 1,
  col = rgb(1, 0, 0, 0.5),
  xlab = expression("Wavenumber ("*cm^{-1}*")"),
  ylab = "Reflectance",
  main = "Calibration data",
  xlim = c(7500, 4000)
)

matplot(
  wavs,
  t(mlocal_val$spc), 
  type = "l", 
  lty = 1,
  col = rgb(1, 0, 0, 0.5),
  xlab = expression("Wavenumber ("*cm^{-1}*")"),
  ylab = "Reflectance",
  main = "Validation data",
  xlim = c(7500, 4000)
)

The spectra

What are we trying to do?

\[\text{y} = {\color{red} f}({\color{dodgerblue} g}(\text{DRS spectrum})) + \epsilon\]

What are we trying to do?

1. Define what to model (e.g. OC, OM, etc.)

2. Define a set of preprocessing methods \({\color{dodgerblue} g}\) to apply to the spectra

3. Define a set of modelling methods \({\color{red} f}\) to apply to the preprocessed spectra

4. Define additional aspects of the calibtration pipeline (e.g. cross-validation, hyperparameter tuning, etc.)

a) Define what to model

Here you can define in a list a set of formulas to model, e.g. OC, OM, etc. For example, if you want to model OC, you can define the following formula:

my_formulas <- list(OC ~ spc)

Note you can also include more than one formula in the list, e.g. if you want to model both OC and OM:

my_formulas_all <- list(OC ~ spc, OM ~ spc)

b) Define a set of preprocessing methods [1/4]

In proximetricsR, you can define a set of preprocessing methods/recipes.

We use the concept of “constructors” to define preprocessing methods. A constructor is a function that takes some parameters and returns a preprocessing method. The preprocessing constructors have a prep_* prefix. For example:

prep_snv()

- prep_snv 

b) Define a set of preprocessing methods [2/4]

Other examples of construtors:

prep_resample(grid = "proxiscout")
prep_smooth(w = 5, p = 1, algorithm = "savitzky-golay")
prep_derivative(m = 1, w = 11, algorithm = "savitzky-golay")
prep_transform(to = "absorbance")

Check more:

help(proximetricsR)

b) Define a set of preprocessing methods [3/4]

Let’s create a recipe of preprocessing methods to apply to the spectra:

recipe_1 <- preprocess_recipe(
  prep_resample(grid = "proxiscout"),
  prep_snv(),
  prep_derivative(m = 1, w = 11, p = 1, algorithm = "savitzky-golay"),
  device = "proxiscout" # !!!
)
recipe_1

Spectral preprocessing recipe (device: "proxiscout"): 
 - Step 1: prep_resample
 - Step 2: prep_snv
 - Step 3: prep_derivative
    m: 1; w: 11; p: 1; algorithm: 'savitzky-golay'

b) Define a set of preprocessing methods [4/4]

Let’s create a list of various preprocessing recipes (to test later):

recipe_2 <- preprocess_recipe(
  prep_resample(grid = "proxiscout"),
  prep_snv(),
  prep_derivative(m = 2, w = 11, p = 2, algorithm = "savitzky-golay"),
  device = "proxiscout" 
)

recipe_3 <- preprocess_recipe(
  prep_resample(grid = "proxiscout"),
  prep_detrend(p = 2),
  prep_derivative(m = 2, w = 11, p = 2, algorithm = "savitzky-golay"),
  device = "proxiscout" 
)

recipe_4 <- preprocess_recipe(
  prep_resample(grid = "proxiscout"),
  prep_transform(to = "absorbance"),
  prep_derivative(m = 2, w = 11, p = 2, algorithm = "savitzky-golay"),
  device = "proxiscout" 
)

recipe_5 <- preprocess_recipe(
  prep_resample(grid = "proxiscout"),
  prep_transform(to = "absorbance"),
  prep_derivative(m = 1, w = 11, p = 2, algorithm = "savitzky-golay"),
  prep_snv(),
  device = "proxiscout" 
)

recipe_6 <- preprocess_recipe(
  prep_resample(grid = "proxiscout"),
  prep_transform(to = "absorbance"),
  prep_derivative(m = 2, w = 11, p = 2, algorithm = "savitzky-golay"),
  prep_snv(),
  device = "proxiscout" 
)

my_preprocessings <- list(recipe_1, recipe_2, recipe_3, recipe_4, recipe_5)

c) Define a set of modelling methods

Here we also use the concept of constructors:

help(fit_constructors)

Create a list of modeling methods (to test later):

fit_1 <- fit_plsr(ncomp = 10, type = "standard")
fit_2 <- fit_plsr(ncomp = 9, type = "modified") 

my_fittings <- list(fit_1, fit_2)

NOTE: check the vignette made for regression methods

d) Define additional aspects of the calibtration pipeline

help(calibration_control)

my_control <- calibration_control(
  validation_type = "kfold", # kfold cross validation 
  number = 5, # 5 folds
  folds = "random", # folds made at random
  remove_outliers = 1 # fit, remove outliers and fit again (one refit iteration)
)

Calibrate!!!

Now, the calibration of a model becomes a simple integration of instructions:

my_model <- calibrate_models(
  formulas = my_formulas,
  data = mlocal_cal,
  preprocess_recipes = my_preprocessings,
  methods = my_fittings,
  control = my_control
)

Plot detailed modelling results:

help(plot.spectral_model)

Inspect the results in plots:

plot(
  my_model$final_models$`OC ~ spc`,  
  spectral = "all", 
  cv = "all", 
  validation = "all"
)

Validate (externally) [1/3]

one_model <- my_model$final_models$`OC ~ spc`
y_hat <- predict(one_model, newdata = mlocal_val)

axes_range <- range(c(y_hat$predictions, mlocal_val$OC), na.rm = TRUE)

plot(
  y_hat$predictions, 
  mlocal_val$OC, 
  xlab = "Predicted OC (%)", 
  ylab = "Reference OC (%)", 
  main = "External validation",
  xlim = axes_range,
  ylim = axes_range, 
  pch = 16, 
  cex = 2, 
  col = "#F0AD00"
)
abline(0, 1, col = "red")

Validate (externally) [2/3]

Validate (externally) [3/3]

cat("R2 =", round(cor(y_hat$predictions, mlocal_val$OC)^2, 3), "\n")

R2 = 0.734 

cat("RMSE =", round(mean((y_hat$predictions - mlocal_val$OC)^2)^0.5, 3), "\n")

RMSE = 1.24 

Or

val_results <- validate_prediction(
  prediction = y_hat, 
  reference = mlocal_val$OC
)
val_results$validation$ncomp_4$val_stats

    rsq    rmse max_res 
  0.734   1.246   5.444 

What are we trying to do?

Deploy the model

dir.create("outputs", showWarnings = FALSE)

proxiscout_write_model(
  one_model, 
  file = "outputs/local_oc_model.json"
)

• Collect the file and upload it…

What are we trying to do now?

Imagine we only have few samples with reference values (as if we were just starting with a new project)…

• Local samples: Only 30 samples with reference values (OC) and their corresponding spectra. However, many more with spectra but without reference values.

• Spectral library: We have a spectral library with a considerable amount of calibration samples.

GOAL: We want to select the most relevant samples from the spectral library to build a calibration model for our local samples.

Two new files

These files come from a soil spectral library from USA:

my_folder/
├── data/
│   ├── local-samples-spain/
│   │   ├── local-spectra-cal.xlsx
│   │   ├── local-properties-cal.xlsx
│   │   ├── local-spectra-val.xlsx
│   │   └── local-properties-val.xlsx
│   └── spectral-library-usa/
│       ├── spectral-library-spectra.xlsx
│       └── spectral-library-properties.xlsx
├── R/
└── outputs/

https://docs.soilspectroscopy.org/neospectra.html

Background

Ramirez-Lopez et al. (2026)

A trivia example of gesearch

• Topic: beers of the world

• Groups: several groups of 5 people at the bar

• Sessions: every Friday

• Gens: persons

• Individuals: groups

• Evaluation period: every year

Searching for the best training samples

The gesearch algorithm is based on a novel evolutionary search method that allow us to identify the most relevant training samples from a large spectral library to be used to build linear models for a given target domain.

In the resemble package, the gesearch() function implements this search method.

help(gesearch)

Read the new files

mlibrary <- proxiscout_read_data(
  "data/spectral-library-usa/spectral-library-spectra.xlsx",
  references_file = "data/spectral-library-usa/spectral-library-properties.xlsx"
)

dim(mlibrary)

[1] 2106   27

Inspect the spectral library data

colnames(mlibrary)

 [1] "sampleName"                 "deviceId"                  
 [3] "OC"                         "oc_usda.c729_w.pct"        
 [5] "c.tot_usda.a622_w.pct"      "n.tot_usda.a623_w.pct"     
 [7] "s.tot_usda.a624_w.pct"      "ph.h2o_usda.a268_index"    
 [9] "bd_usda.a4_g.cm3"           "clay.tot_usda.a334_w.pct"  
[11] "silt.tot_usda.c62_w.pct"    "sand.tot_usda.c60_w.pct"   
[13] "caco3_usda.a54_w.pct"       "efferv_usda.a479_class"    
[15] "cec_usda.a723_cmolc.kg"     "ca.ext_usda.a722_cmolc.kg" 
[17] "mg.ext_usda.a724_cmolc.kg"  "k.ext_usda.a725_cmolc.kg"  
[19] "na.ext_usda.a726_cmolc.kg"  "wr.33kPa_usda.a415_w.pct"  
[21] "wr.1500kPa_usda.a417_w.pct" "al.dith_usda.a65_w.pct"    
[23] "p.ext_usda.a652_mg.kg"      "p.ext_usda.a1070_mg.kg"    
[25] "k.ext_usda.a1065_mg.kg"     "ec_usda.a364_ds.m"         
[27] "spc"                       

View(mlibrary)

Raw spectral data

Second derivative

second_derivative_recipe <- preprocess_recipe(
  prep_transform(to = "absorbance"),
  prep_derivative(m = 2, w = 9, p = 2, algorithm = "savitzky-golay"), 
  device = "proxiscout"
)
my_second_derivative <- process(mlibrary$spc, second_derivative_recipe)

Simulation

mlocal_cal$ref
mlocal_cal$OC

mlocal_cal <- proxiscout_read_data(
  file = "data/local-samples-spain/local-spectra-cal.xlsx", 
  references_file = "data/local-samples-spain/local-properties-cal.xlsx"
)

mlocal_cal$OC[mlocal_cal$ref == 0] <- NA

mlocal_val <- proxiscout_read_data(
  file = "data/local-samples-spain/local-spectra-val.xlsx", 
  references_file = "data/local-samples-spain/local-properties-val.xlsx"
)

Prepare to execute the serach in parallel in R

The search is computationally intensive. The gesearch() function supports parallel computing.

Load the doParallel package

# Parallel processing
library(doParallel)

# !!!!!
parallel::detectCores() 

n_cores <- 4
cl <- makeCluster(n_cores)
registerDoParallel(cl)

Wavenumber trimming

Here, we define the wavenumber from which we will be trimming our spectral data:

wav_lim <- 6000

and extract (again, just in case) the vector of wavenumbers:

wavs <- as.numeric(colnames(mlocal_cal$spc))

Get a vector of TRUE for the wavenumbers to keep and FALSE for the wavenumbes to ignore:

wav_sel <- wavs < wav_lim
wav_sel

Search!!!

search_recipe <- preprocess_recipe(
  prep_resample(grid = "proxiscout"),
  prep_transform(to = "absorbance"),
  prep_derivative(m = 2, w = 11, p = 2, algorithm = "savitzky-golay"),
  prep_snv(),
  device = "proxiscout" 
)

search_results <- gesearch(
  Xr = process(mlibrary$spc[, wav_sel], search_recipe), 
  Yr = mlibrary$OC,
  Xu = process(mlocal_cal$spc[, wav_sel], search_recipe), 
  Yu = mlocal_cal$OC,
  Yu_lims = c(0.02, 15), # New!!
  k = 30, b = 200, retain = 0.95, # the key gesearch parameters
  target_size = 80,
  fit_method = fit_pls(ncomp = 10, method = "mpls"),
  optimization = c("similarity", "reconstruction", "response", "range"),
  control = gesearch_control(retain_by = "probability"),
  seed = 1410, 
  pchunks = 5
)

Use the gesearch model

The gesearch() function builds a model with the samples found in combination with the samples from the local domain used to guide the search. SO we can use that model to predict the validation samples:

my_gesearch_predictions <- predict(
  search_results, 
  newdata = process(mlocal_val$spc[, wavs < wav_lim], search_recipe)
)[[1]]

my_gesearch_predictions

We could go on and do the validation directly here, but this will be using resemble… Since we want to deploy our model, this is way easier if we switch back to proximetricsR.

Build the hybrid model with proximetricsR

search_results$indices contains the indices of the samples found in the spectral library, so, create a hybrid modelling dataset (the samples found + the few samples from the local domain):

mlibrary_selected <- mlibrary[search_results$indices, ]

# add variables that are in the local that are not in the library
mlibrary_selected$ref <- NA
mlibrary_selected$OM <- NA

mlibrary_selected <- mlibrary_selected[, colnames(mlocal_cal)]

mhybrid_cal <- rbind(
  mlocal_cal, 
  mlibrary_selected
)