soilVAE Workflow

This vignette predicts soil carbon from NIR spectra with soilVAE and compares it with a PLS baseline.

The data are NIRsoil from the prospectr package. They are loaded from prospectr when the vignette runs; soilVAE itself contains code only and bundles no data.

The model runs in Python. Install the backend once with vae_install(); soilVAE then finds it in every session. The model-fitting chunks below run only when vae_available() is TRUE (it never installs anything), and the whole vignette needs the suggested packages prospectr and pls.

library(soilVAE)
has_pkgs <- requireNamespace("prospectr", quietly = TRUE) &&
  requireNamespace("pls", quietly = TRUE)
has_tf <- has_pkgs && vae_available()
knitr::opts_chunk$set(eval = has_pkgs)
# The tuning example fits several models; it runs only outside CRAN checks.
run_tuning <- has_tf && identical(Sys.getenv("NOT_CRAN"), "true")
library(prospectr)
library(pls)

Data

NIRsoil holds 825 soil samples with NIR absorbance spectra (1100–2498 nm) and laboratory values. The target here is total carbon (Ciso, %). The data come with a calibration/test split (train).

data("NIRsoil", package = "prospectr")
d <- NIRsoil[is.finite(NIRsoil$Ciso), ]
table(set = ifelse(d$train == 1, "calibration", "test"))
#> set
#> calibration        test 
#>         548         184

Preprocessing

Standard normal variate (SNV) removes scatter. No further scaling is needed: vae_fit() centres each band and divides all bands by one overall SD (the default scale_x = "global", which keeps the shape of the spectrum), and stores these statistics in the model, so vae_predict() scales new spectra the same way. (PLS centres the data itself.)

X <- standardNormalVariate(d$spc)
cal <- d$train == 1

y <- d$Ciso
X_cal <- X[cal, ];  y_cal <- y[cal]
X_te  <- X[!cal, ]; y_te  <- y[!cal]

wl <- as.numeric(colnames(d$spc))
matplot(wl, t(standardNormalVariate(d$spc[1:30, ])), type = "l", lty = 1,
        col = "grey40", xlab = "Wavelength (nm)", ylab = "SNV absorbance")

Metrics

metrics <- function(obs, pred) {
  rmse <- sqrt(mean((obs - pred)^2))
  c(R2   = 1 - sum((obs - pred)^2) / sum((obs - mean(obs))^2),
    RMSE = rmse,
    RPIQ = IQR(obs) / rmse)
}

Baseline: PLS

The number of components is chosen by 10-fold cross-validation on the calibration set.

set.seed(1)
pls_fit <- plsr(y_cal ~ X_cal, ncomp = 30, validation = "CV")
ncomp <- which.min(RMSEP(pls_fit)$val[1, 1, -1])
pls_pred <- predict(pls_fit, X_te, ncomp = ncomp)[, 1, 1]
round(metrics(y_te, pls_pred), 3)
#>    R2  RMSE  RPIQ 
#> 0.683 0.855 0.999

soilVAE

A fifth of the calibration samples is held out for early stopping: training stops when the prediction error on these samples stops improving, and the best weights are kept. Early stopping starts after a warm-up of patience epochs. (Passing validation_split = 0.2 instead of X_val/y_val does the same split inside vae_fit().)

set.seed(1)
va <- sample(nrow(X_cal), round(0.2 * nrow(X_cal)))

model <- vae_build(input_dim = ncol(X_cal), seed = 1L)
vae_fit(model, X_cal[-va, ], y_cal[-va], X_cal[va, ], y_cal[va],
        epochs = 300L, patience = 20L)

vae_pred <- vae_predict(model, X_te)
round(metrics(y_te, vae_pred), 3)
#>    R2  RMSE  RPIQ 
#> 0.755 0.752 1.137

The training history shows the validation prediction loss that early stopping watched:

h <- vae_history(model)
plot(h$val_y_loss, type = "l", xlab = "Epoch", ylab = "Validation y loss (standardised)")

Tuning

tune_vae_train_val() fits a grid of configurations on a train/validation split; select_best_from_grid() picks one. Each configuration gets its own seed, so a tuning run is reproducible. (This chunk is skipped during CRAN checks because it fits several models.)

grid <- data.frame(latent_dim = c(16L, 32L), dropout = 0.1, lr = 5e-4,
                   beta_kl = 1e-3, alpha_y = 1, epochs = 300L,
                   batch_size = 64L, patience = 20L)
grid$hidden_enc <- list(c(512L, 256L), c(512L, 256L))
grid$hidden_dec <- list(c(256L, 512L), c(256L, 512L))

tuned <- tune_vae_train_val(X_cal[-va, ], y_cal[-va], X_cal[va, ], y_cal[va],
                            seed = 1, grid_vae = grid)
tuned$tuning_df[, c("latent_dim", "R2_val", "RMSE_val")]
#>   latent_dim    R2_val  RMSE_val
#> 1         16 0.9076586 0.6155152
#> 2         32 0.8963394 0.6521497
select_best_from_grid(tuned$tuning_df, "rmse")$best$latent_dim
#> [1] 16

Comparison on the test set

round(rbind(PLS = metrics(y_te, pls_pred), soilVAE = metrics(y_te, vae_pred)), 3)
#>            R2  RMSE  RPIQ
#> PLS     0.683 0.855 0.999
#> soilVAE 0.755 0.752 1.137
plot(y_te, vae_pred, xlab = "Observed Ciso (%)", ylab = "Predicted Ciso (%)",
     pch = 16, col = "grey30")
abline(0, 1, lty = 2)

Saving a model

A fitted model, including the predictor and target scaling, is saved to one .keras file:

f <- tempfile(fileext = ".keras")
vae_save(model, f)
model2 <- vae_load(f)
all.equal(vae_predict(model2, X_te), vae_pred)
#> [1] TRUE