Introduction to densemlp

What densemlp is

densemlp is a compact dense feedforward neural network (multilayer perceptron) for regression, classification and survival analysis on tabular data. Forward propagation, backpropagation, and Adam optimization are implemented natively in C++ via RcppArmadillo - there is no torch/libtorch dependency, which is what makes it fast to install and fast to train.

Hidden layers are Linear -> BatchNorm -> ReLU -> [gate] -> [dropout] -> [residual] (set batch_norm = FALSE to drop the normalization step), with a plain Linear -> task activation output layer (linear for regression, sigmoid for binary classification, softmax for multiclass, a linear risk score for survival). Numeric predictors are centered/scaled and categorical predictors are one-hot encoded internally, so x can be a mixed-type data frame.

Fitting a single model

set.seed(1)
n <- 300
x <- data.frame(a = rnorm(n), b = rnorm(n), c = rnorm(n))
y <- 2 * x$a - x$b + 0.5 * x$a * x$c + rnorm(n, sd = 0.2)

train <- sample.int(n, floor(0.8 * n))
test <- setdiff(seq_len(n), train)

fit <- densemlp(x[train, ], y[train], hidden_units = c(32, 16), epochs = 60)
pred <- predict(fit, x[test, ])
sqrt(mean((pred - y[test])^2))
#> [1] 0.3286649

task = "auto" (the default) infers "regression" from a numeric outcome, or "binary" / "multiclass" from a factor/character outcome with 2 or more than 2 levels respectively.

y_class <- factor(ifelse(y > median(y), "high", "low"))
fit_class <- densemlp(x[train, ], y_class[train], hidden_units = c(32, 16), epochs = 60)
predict(fit_class, x[test, ], type = "prob")[1:5, ]
#>            high        low
#> [1,] 0.07849898 0.92150102
#> [2,] 0.96028934 0.03971066
#> [3,] 0.80822509 0.19177491
#> [4,] 0.26178535 0.73821465
#> [5,] 0.33715449 0.66284551
predict(fit_class, x[test, ], type = "class")[1:5]
#> [1] low  high high low  low 
#> Levels: high low

Architecture options

All accuracy-oriented options are opt-in and default to off, so a plain densemlp(x, y) call keeps the smallest, fastest architecture.

fit_opts <- densemlp(
  x[train, ], y[train], hidden_units = c(32, 16),
  residual = TRUE,          # learned skip connection per hidden block
  gated = TRUE,              # learned sigmoid gate per hidden block
  dropout = 0.05,            # inverted dropout per hidden layer
  batch_norm = TRUE,         # batch normalization inside each hidden block
  input_projection = 8,      # linear map of the encoded inputs before layer 1
  ema_decay = 0.99,          # moving-average weights for eval/final predictions
  lr_schedule = "cosine",
  epochs = 60
)

input_projection = k inserts a bare linear layer (no activation, no batch normalization) that maps the encoded predictors to k dimensions before the first hidden block; it cannot be combined with interaction.

interaction = TRUE (an efficient learned cross-feature layer) is also available, but is currently numerically unstable on small-n/high-p regression data and is not recommended - see NEWS.md for details. It is left out of tune_densemlp()’s search grid for the same reason.

Internal ensembles

Setting ensemble > 1 fits several bootstrap-resampled members internally and averages their predictions transparently - predict() still returns a single vector/matrix. Use ncores to fit members in parallel (via parallel::mclapply on Unix-alikes; falls back to serial on Windows).

fit_ens <- densemlp(
  x[train, ], y[train], hidden_units = c(32, 16), epochs = 60,
  ensemble = 5, ncores = 2
)
sqrt(mean((predict(fit_ens, x[test, ]) - y[test])^2))
#> [1] 0.3758169

Cross-validation

cv_densemlp() fits densemlp() on each of folds splits and scores the held-out fold with densemlp_metrics() (RMSE/NRMSE/R² for regression; accuracy/balanced accuracy/macro AUC/log loss for classification). Extra arguments are forwarded to every fold’s densemlp() call.

cv <- cv_densemlp(x, y, folds = 5, hidden_units = c(32, 16), epochs = 40)
cv
#> densemlp cross-validation (5 folds, task: regression)
#> 
#> Per-fold metrics:
#>  fold      rmse     nrmse       rsq
#>     1 0.4378918 0.1923046 0.9623922
#>     2 0.4579028 0.2072769 0.9563081
#>     3 0.4443015 0.1739609 0.9692247
#>     4 0.4958937 0.2292346 0.9465608
#>     5 0.4950703 0.2426051 0.9401452
#> 
#> Summary:
#>  metric      mean         sd
#>    rmse 0.4662120 0.02768098
#>   nrmse 0.2090764 0.02760508
#>     rsq 0.9549262 0.01173195

Hyperparameter tuning

tune_densemlp() grid-searches over hidden_units, dropout, residual, gated, ema_decay, lr_schedule, epochs, batch_size, and lr. Each candidate is repeated (repeats) with successive seeds and ranked by mean best-epoch validation loss; the best configuration is refit on the full data by default.

tuned <- tune_densemlp(
  x, y, repeats = 2,
  grid = list(
    hidden_units = list(c(16), c(32, 16)),
    lr = c(1e-3, 3e-3)
  )
)
tuned$best_config
#>   hidden_units dropout residual gated ema_decay lr_schedule epochs batch_size
#> 1           16       0    FALSE FALSE         0        none    100         32
#>      lr      score    score_sd repeats
#> 1 0.003 0.01670958 0.003480581       2

Survival outcomes

With task = "survival" the response is a survival::Surv(time, event) object (or a two-column (time, event) matrix). Two losses are available: loss = "cox" (the default) trains a single linear risk score with a batch-wise Breslow-tie Cox partial likelihood, and loss = "brier" trains a discrete-time hazard head against the IPCW integrated Brier score.

library(survival)
data(lung, package = "survival")
#> Warning in data(lung, package = "survival"): data set 'lung' not found
lung <- na.omit(lung[, c("time", "status", "age", "sex", "ph.ecog", "ph.karno", "wt.loss")])
sy <- Surv(lung$time, lung$status == 2)
sx <- lung[, c("age", "sex", "ph.ecog", "ph.karno", "wt.loss")]

sfit <- densemlp(sx, sy, task = "survival", hidden_units = c(16, 8), epochs = 60)
risk <- predict(sfit, sx, type = "response")
densemlp_metrics(sy, risk, task = "survival")
#> $concordance
#> [1] 0.6801176

Variable importance

perm_importance() gives model-agnostic permutation importance: each predictor is shuffled in turn and the drop in densemlp_metrics() performance is recorded.

imp <- perm_importance(fit, x[test, ], y[test])
imp$data
#>   feature importance
#> 1       a  2.7388069
#> 2       b  1.4793011
#> 3       c  0.6227803

Visualizing training

plot(fit)          # or, equivalently, plot_history(fit)

Learn more

See ?densemlp, ?cv_densemlp, ?tune_densemlp, ?densemlp_metrics and ?perm_importance for full argument documentation, and NEWS.md for the changelog.