| Title: | Fitting Interpretable Neural Additive Models Using Orthogonalization |
|---|---|
| Description: | An algorithm for fitting interpretable additive neural networks for identifiable and visualizable feature effects using post hoc orthogonalization. Fit custom neural networks intuitively using established 'R' 'formula' notation, including interaction effects of arbitrary order while preserving identifiability to enable a functional decomposition of the prediction function. For more details see Koehler et al. (2025) <doi:10.1038/s44387-025-00033-7>. |
| Authors: | David Köhler [aut, cre] (ORCID: <https://orcid.org/0009-0006-0027-4046>) |
| Maintainer: | David Köhler <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.1.0 |
| Built: | 2026-11-03 18:03:30 UTC |
| Source: | https://github.com/koehlibert/onam_r |
Returns a function function(inputs) {...} that builds a feed-forward
stack of layer_dense() calls on top of the inputs tensor it is
later called with, and returns the resulting keras3::keras_model().
build_dnn(units = NULL, layer_df = NULL, default_activation = NULL)build_dnn(units = NULL, layer_df = NULL, default_activation = NULL)
units |
Optional integer vector of layer sizes. Ignored if
|
layer_df |
Optional data frame / tibble with one row per layer.
Must contain a
Any column not supplied falls back to
tibble::tibble(
units = c(64, 32, 1),
activation = c("gelu", "gelu", "linear"),
kernel_regularizer = list(keras3::regularizer_l2(0.001), NULL, NULL)
)
|
default_activation |
Activation function to be used if not otherwise
specified in |
Layers can be specified either as a simple vector of unit counts, or as a data frame / tibble giving full per-layer control over activation, use_bias, initializers, regularizers, layer name, etc.
A function of a single argument inputs (a keras input
tensor) that returns a keras3::keras_model().
mod_units <- build_dnn(units = c(16, 8, 1)) mod_df <- build_dnn(layer_df = data.frame( units = c(32, 16, 8, 1), activation = c("gelu", "gelu", "gelu", "linear"), use_bias = c(FALSE, FALSE, FALSE, TRUE) ))mod_units <- build_dnn(units = c(16, 8, 1)) mod_df <- build_dnn(layer_df = data.frame( units = c(32, 16, 8, 1), activation = c("gelu", "gelu", "gelu", "linear"), use_bias = c(FALSE, FALSE, FALSE, TRUE) ))
Compute Generalized Sobol Indices for dependent features
gen_sobol(object, data = NULL)gen_sobol(object, data = NULL)
object |
Either model of class |
data |
Data for which the model is to be evaluated. If |
For details on generalized sobol indices, see Chastaing et al. (2012) doi:10.1214/12-EJS749.
Returns a named vector of percentage of variance explained by each interaction order.
# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- build_dnn(units = c(16, 8, 1)) list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) gen_sobol(mod)# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- build_dnn(units = c(16, 8, 1)) list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) gen_sobol(mod)
Helper function to install Keras and packages necessary for package
functionality into a conda environment. Use this function if
keras3::install_keras() does not work, esp. on windows machines.
install_conda_env( envname = "r-keras", python_version = "python=3.10", overwrite = FALSE )install_conda_env( envname = "r-keras", python_version = "python=3.10", overwrite = FALSE )
envname |
Name for the conda environment to be created. |
python_version |
Python version to be installed in the conda environment. |
overwrite |
Should an existing conda environment of name |
No return value, called for side effects
Load a fitted onam model from disk
load_onam(dir)load_onam(dir)
dir |
Directory as created by save_onam. |
Returns the restored object of class onam.
Fits an interpretable neural additive model with post hoc orthogonalization for a given network architecture and user-specified feature sets.
onam( formula, list_of_deep_models, data, model = NULL, prediction_function = NULL, model_data = NULL, categorical_features = NULL, target = "continuous", n_ensemble = 10, epochs = 500, learning_rate = 0.001, callback = NULL, seed = NULL, progresstext = FALSE, verbose = 0 )onam( formula, list_of_deep_models, data, model = NULL, prediction_function = NULL, model_data = NULL, categorical_features = NULL, target = "continuous", n_ensemble = 10, epochs = 500, learning_rate = 0.001, callback = NULL, seed = NULL, progresstext = FALSE, verbose = 0 )
formula |
Formula for model fitting. Specify deep parts with the same
name as |
list_of_deep_models |
List of named models used in |
data |
Data to be fitted |
model |
Prediction model that is to be explained. Output of the model as
returned from |
prediction_function |
Prediction function to be used to generate the
outcome. Only used if |
model_data |
Data used for generating predictions of |
categorical_features |
Vector of feature names of categorical features. |
target |
Target of prediction task. Can be either "continuous" or "binary". For "continuous"(default), an additive model for the prediction of a continuous outcome is fitted. For "binary", a binary classification with sigmoid activation in the last layer is fitted. |
n_ensemble |
Number of orthogonal neural additive model ensembles |
epochs |
Number of epochs to train the model. See
|
learning_rate |
Learning rate for model fitting. See
|
callback |
Callback to be called during training. See
|
seed |
Random seed used by R, python, numpy, and backend framework.
See |
progresstext |
Show model fitting progress. If |
verbose |
Verbose argument for internal model fitting. Used for
debugging. See |
For more details see Koehler et al. (2025) https://doi.org/10.1038/s44387-025-00033-7.
Returns a model object of class onam, containing all ensemble
members, ensemble weights, and main and interaction effect outputs.
# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- function(inputs) { # `inputs` is a list of shared feature inputs for interaction terms combined <- if (is.list(inputs)) layer_concatenate(inputs) else inputs outputs <- combined %>% layer_dense(units = 16, activation = "gelu") %>% layer_dense(units = 8, activation = "linear", use_bias = TRUE) %>% layer_dense(units = 1, activation = "linear", use_bias = TRUE) keras_model(inputs, outputs) } list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) summary(mod)# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- function(inputs) { # `inputs` is a list of shared feature inputs for interaction terms combined <- if (is.list(inputs)) layer_concatenate(inputs) else inputs outputs <- combined %>% layer_dense(units = 16, activation = "gelu") %>% layer_dense(units = 8, activation = "linear", use_bias = TRUE) %>% layer_dense(units = 1, activation = "linear", use_bias = TRUE) keras_model(inputs, outputs) } list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) summary(mod)
Plot Interaction Effect
plot_inter_effect( object, feature1, feature2, interpolate = FALSE, labs = NULL, custom_colors = "spectral", n_interpolate = 200, include_main = FALSE )plot_inter_effect( object, feature1, feature2, interpolate = FALSE, labs = NULL, custom_colors = "spectral", n_interpolate = 200, include_main = FALSE )
object |
Either model of class |
feature1, feature2
|
Effects to be plotted. |
interpolate |
If TRUE, values will be interpolated for a smooth plot. If FALSE (default), only observations in the data will be plotted. |
labs |
An optional named vector that can contain axis labels and the name of the effect. Expected vector names are 'xlab', 'ylab' and 'effect'. |
custom_colors |
color palette object for the interaction plot. Default is "spectral", returning a color palette based on the spectral theme. |
n_interpolate |
number of values per coordinate axis to interpolate. Ignored if 'interpolate = FALSE'. |
include_main |
If TRUE, main effects for features feature1 and feature2 will be added to the interaction term to give combined effect of main and interaction effects. Default is FALSE. |
Returns a 'ggplot2' object of the specified effect interaction
# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- build_dnn(units = c(16, 8, 1)) list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) plot_inter_effect(mod, "x1", "x2")# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- build_dnn(units = c(16, 8, 1)) list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) plot_inter_effect(mod, "x1", "x2")
Plot Main Effect
plot_main_effect(object, feature, reference_level = NULL, labs = NULL)plot_main_effect(object, feature, reference_level = NULL, labs = NULL)
object |
Either model of class |
feature |
Feature for which the effect is to be plotted, must be present in the model formula. For interaction terms, use plotInteractionEffect |
reference_level |
Reference level to be used when plotting categorical effects for better interpretability. Effect of the reference level will be set to zero and subtracted from all over factor levels. |
labs |
An optional named vector that can contain axis labels. Expected vector names are 'xlab', 'ylab' and 'effect'. If both 'ylab' and 'effect' are specified, 'effect' will be used as y-axis label. |
Returns a ggplot2 object of the specified effect
# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- build_dnn(units = c(16, 8, 1)) list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) plot_main_effect(mod, "x1")# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- build_dnn(units = c(16, 8, 1)) list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) plot_main_effect(mod, "x1")
Plot generalized Sobol indices
## S3 method for class 'gen_sobol' plot(x, ...)## S3 method for class 'gen_sobol' plot(x, ...)
x |
Object of class |
... |
further arguments, currently unused |
For details on generalized sobol indices, see Chastaing et al. (2012) doi:10.1214/12-EJS749.
Returns a 'ggplot2' object showing, for each fitted effect, the
generalized Sobol index split into the effect's own contribution
(gen_sobol_index_1) and its contribution through interactions with
other effects of the same order (gen_sobol_index_2)
Plot variance decomposition
## S3 method for class 'var_decomp' plot(x, ...)## S3 method for class 'var_decomp' plot(x, ...)
x |
Object of class |
... |
further arguments, currently unused |
Returns a 'ggplot2' object showing the fraction of total variance explained by each interaction order
Evaluate orthogonal neural additive model
## S3 method for class 'onam' predict(object, newdata = NULL, ...)## S3 method for class 'onam' predict(object, newdata = NULL, ...)
object |
model of class |
newdata |
Data for which the model is to be evaluated. If NULL (default),
data with which |
... |
some methods for this generic require additional arguments. None are used in this method. |
Returns a list containing data, model output for each observation in
newdata and main and interaction effects obtained by the model
An onam object holds one or more fitted 'keras' submodels per ensemble
member (in object$ensemble[[i]]$model_list), which cannot be
serialized with saveRDS() alone. This function saves each submodel to
its own .keras file and the remaining orthogonalization metadata
(weights, model_info, data, predictions, ...) to a single .rds file,
all inside dir. Use load_onam to restore the object.
save_onam(object, dir, overwrite = FALSE)save_onam(object, dir, overwrite = FALSE)
object |
Object of class |
dir |
Directory the model is saved to. Created if it does not exist. |
overwrite |
Should an existing directory of the same name be overwritten? |
No return value, called for side effects.
generates a summary of a fitted onam object including
information on ensembling strategy and performance metrics such as
correlation and degree of interpretability
## S3 method for class 'onam' summary(object, ...) ## S3 method for class 'summary.onam' print(x, ...)## S3 method for class 'onam' summary(object, ...) ## S3 method for class 'summary.onam' print(x, ...)
object |
onam object of class |
... |
further arguments passed to or from other methods. |
x |
object of class |
For examples see example(onam)
Gives summary of the onam object, including model inputs, number
of ensembles, correlation of model output and original outcome variable, and
interpretability metrics i_1 and i_2
Get variance decomposition of orthogonal neural additive model
var_decomp(object, data = NULL)var_decomp(object, data = NULL)
object |
Either model of class |
data |
Data for which the model is to be evaluated. If |
Returns a named vector of percentage of variance explained by each interaction order.
# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- build_dnn(units = c(16, 8, 1)) list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) var_decomp(mod)# Basic example for a simple ONAM-model # Create training data n <- 1000 x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) y <- sin(x1) + ifelse(x2 > 0, pweibull(x2, shape = 3), pweibull(-x2, shape = 0.5)) + x1 * x2 data_train <- cbind(x1, x2, y) # Define model model_formula <- y ~ mod1(x1) + mod1(x2) + mod1(x1, x2) mod1 <- build_dnn(units = c(16, 8, 1)) list_of_deep_models <- list(mod1 = mod1) # Fit model mod <- onam(model_formula, list_of_deep_models, data_train, n_ensemble = 1, epochs = 10) var_decomp(mod)