Load spectra

sim2 contains 100 simulated 1D NMR spectra split evenly into groups A and B. Five of every 25 peaks per spectrum differ between groups by 10 % in area. The group labels are attached as an attribute; see ?sim2 for details.

library(metabodeconplus)
x <- sim2
y <- attr(sim2, "group")
n <- length(x)

We use one half of the data for training and one half for testing.

set.seed(1)
tr <- sort(sample(n, round(0.5 * n)))
te <- setdiff(seq_len(n), tr)
true_x0 <- attr(sim2, "true_x0")   # ppm of the discriminating peaks

The pipeline, step by step

We pick a small group-balanced subset (4 + 4 spectra) just for the plots; the model itself is trained on all of x[tr].

abtr <- c(which(y[tr] == "A")[1:4], which(y[tr] == "B")[1:4])
yab <- y[tr][abtr]

decons <- deconvolute(x[tr], nfit=10, smit=2, smws=5, delta=10, npmax=0, verbose=FALSE)
plot_spectra(decons[abtr])
heat_spectra(decons[abtr], y=yab)

aligns <- clupa(decons, maxShift=50, verbose=FALSE)
ref <- attr(aligns, "ref")
plot_spectra(aligns[abtr])

snapped <- snap_to_ref(aligns, maxCombine=5)

X <- peak_mat(snapped)
peakPos <- attr(X, "peakPos")
dim(X)

## [1] 50 21

heat_spectra(X, y=y[tr], true_x0=true_x0)
heat_spectra(X, y=y[tr], true_x0=true_x0, scale_cols=TRUE)

rf <- ranger::ranger(x=X, y=y[tr], probability=TRUE, num.trees=500, seed=1)
cat(sprintf("OOB error: %.1f%%\n", 100 * rf$prediction.error))

## OOB error: 11.8%

The one-shot call

fit_mdm() performs all five steps – deconvolute, align, snap, featurize, fit – and, when any of npmax / maxShift / maxCombine is a vector, searches the grid and returns the best model. Choose the backend with model = "ranger".

md <- fit_mdm(x[tr], y[tr], model="ranger",
              npmax=0L, maxShift=50L, maxCombine=5L,
              verbosity=0, nworkers=1)
print(md)

## metabodeconplus model (mdm)
##   model:         ranger
##   npmax:         0
##   maxShift:      50
##   maxCombine:    5
##   acc:           88.0%
##   auc:           92.0%

Predict held-out spectra

predict() on an mdm object mirrors the training pipeline for new data, reusing the stored reference and feature columns.

# Small rank-based AUC helper (positive class = second factor level).
auc <- function(y, prob) {
    pos <- y == levels(y)[2]; r <- rank(prob)
    n1 <- sum(pos); n0 <- sum(!pos)
    if (n1 == 0 || n0 == 0) NA_real_ else (sum(r[pos]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}
preds <- predict(md, x[te], type="all", verbosity=0)
acc <- mean(preds$class == y[te])
au  <- auc(y[te], preds$prob)
cat(sprintf("Test accuracy: %.1f%%\n", 100 * acc))

## Test accuracy: 84.0%

cat(sprintf("Test AUC:      %.3f\n", au))

## Test AUC:      0.825

Tune the preprocessing

Passing a vector for npmax, maxShift or maxCombine makes fit_mdm() evaluate the cartesian product and keep the best-scoring cell (accuracy, ties broken by AUC). The augmented grid is returned in md$mog.

mt <- fit_mdm(x[tr], y[tr], model="ranger",
              npmax=c(0L, 30L), maxShift=c(20L, 50L), maxCombine=c(2L, 5L),
              verbosity=0, nworkers=1)
knitr::kable(head(mt$mog[order(-mt$mog$auc), ], 5), row.names=FALSE,
             caption="Top parameter combinations by AUC.")

For an honest generalization estimate on small datasets, wrap the whole search in outer cross-validation with benchmark() (not run here because it repeats the grid search for every fold):

bm <- benchmark(x, y, model="ranger", npmax=0L, maxShift=50L, maxCombine=5L, k=5)
mean(bm$predictions$true == bm$predictions$pred)