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.

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