funcml: A Formula-First Framework for Machine Learning in R

Machine learning in R has long suffered from a fragmentation problem. Fitting a random forest, a regularized linear model, and a neural network each requires learning a different package, a different formula syntax, and a different set of conventions for prediction and evaluation. funcml was built to close that gap.

What is funcml?

funcml provides a unified, formula-first interface to supervised learning in R. Six top-level verbs cover the full ML workflow:

A minimal example

Verb	What it does
`fit()`	Train any learner via a common formula interface
`evaluate()`	Resampling-based performance with any metric
`tune()`	Grid or random hyperparameter search
`compare_learners()`	Benchmark multiple models side by side
`interpret()`	14 model-agnostic interpretability methods
`estimate()`	Plug-in g-computation for ATE/CATE estimation

library(funcml)

# fit a random forest
fit_rf <- fit(Sepal.Length ~ ., data = iris, model = "ranger")

# evaluate with 5-fold CV
evaluate(data = iris, formula = Sepal.Length ~ ., model = "ranger", resampling = cv(5), seed = 1)

<funcml_eval> model: ranger | task: regression
  metric       mean          sd n   std_error conf_level   conf_low  conf_high
1   rmse 0.33905736 0.045133466 5 0.020184300       0.95 0.28301676 0.39509796
2    mae 0.27705969 0.030313120 5 0.013556440       0.95 0.23942098 0.31469840
3    mse 0.11658952 0.032166013 5 0.014385078       0.95 0.07665014 0.15652890
4  medae 0.24278001 0.039072187 5 0.017473613       0.95 0.19426548 0.29129454
5   mape 0.04779262 0.004993545 5 0.002233181       0.95 0.04159232 0.05399293
6    rsq 0.82362138 0.055242283 5 0.024705100       0.95 0.75502903 0.89221373

# compare to a regularized regression
compare_learners(
  data = iris,
  formula = Sepal.Length ~ .,
  models = c("ranger", "glmnet"),
  resampling = cv(5),
  seed = 1
)

<funcml_compare> task: regression | tuned: FALSE
    model metric       mean          sd n   std_error conf_level   conf_low
1  ranger   rmse 0.33905736 0.045133466 5 0.020184300       0.95 0.28301676
2  ranger    mae 0.27705969 0.030313120 5 0.013556440       0.95 0.23942098
3  ranger    mse 0.11658952 0.032166013 5 0.014385078       0.95 0.07665014
4  ranger  medae 0.24278001 0.039072187 5 0.017473613       0.95 0.19426548
5  ranger   mape 0.04779262 0.004993545 5 0.002233181       0.95 0.04159232
6  ranger    rsq 0.82362138 0.055242283 5 0.024705100       0.95 0.75502903
7  glmnet   rmse 0.31746014 0.015929563 5 0.007123917       0.95 0.29768098
8  glmnet    mae 0.25600583 0.014843643 5 0.006638279       0.95 0.23757501
9  glmnet    mse 0.10098394 0.010238700 5 0.004578886       0.95 0.08827092
10 glmnet  medae 0.22566286 0.022954893 5 0.010265740       0.95 0.19716060
11 glmnet   mape 0.04394887 0.003418004 5 0.001528578       0.95 0.03970485
12 glmnet    rsq 0.84775805 0.023406028 5 0.010467494       0.95 0.81869563
    conf_high tuned rank
1  0.39509796 FALSE    2
2  0.31469840 FALSE    2
3  0.15652890 FALSE    2
4  0.29129454 FALSE    2
5  0.05399293 FALSE    2
6  0.89221373 FALSE    2
7  0.33723930 FALSE    1
8  0.27443664 FALSE    1
9  0.11369697 FALSE    1
10 0.25416512 FALSE    1
11 0.04819288 FALSE    1
12 0.87682047 FALSE    1

Interpretability built in

interpret() wraps 14 methods, from permutation importance and partial dependence plots to SHAP values and accumulated local effects, behind a single call:

interpret(fit_rf, data = iris, method = "pdp", features = "Petal.Width")$result$curves |> head()

      feature     value     yhat
1 Petal.Width 0.1000000 5.577112
2 Petal.Width 0.2142857 5.591556
3 Petal.Width 0.3285714 5.595045
4 Petal.Width 0.4428571 5.672101
5 Petal.Width 0.5571429 5.662717
6 Petal.Width 0.6714286 5.672528

Causal estimation

When the research question goes beyond prediction, estimate() wraps plug-in g-computation to estimate average and conditional treatment effects:

set.seed(1)
n <- 200
X1 <- rnorm(n); X2 <- rnorm(n)
A <- rbinom(n, 1, plogis(0.3 * X1))
Y <- 1 + 0.5 * A + 0.4 * X1 - 0.2 * X2 + rnorm(n)
mydata <- data.frame(Y, A, X1, X2)

# ATE of treatment A on outcome Y, adjusted for covariates X
estimate(
  data      = mydata,
  formula   = Y ~ A + X1 + X2,
  model     = "ranger",
  treatment = "A",
  estimand  = "ATE",
  seed      = 1
)

<funcml_estimand> ATE via g-computation
Treatment: A (1 vs 0)
Estimate: 0.4531 | SE: 0.0156 | 95% normal CI [0.4225, 0.4836]

The learner registry

funcml ships with 26 learners spanning linear models, tree ensembles, SVMs, neural networks, additive models, boosting, and stacking. All are accessed by name:

list_learners()[, c("learner", "has_tune")]

        learner has_tune
16     adaboost     TRUE
23         bart     TRUE
10          C50     TRUE
19      cforest     TRUE
18        ctree     TRUE
7     e1071_svm     TRUE
12        earth     TRUE
15          fda     TRUE
13          gam     TRUE
9           gbm     TRUE
1           glm     TRUE
3        glmnet     TRUE
11         kknn     TRUE
20          lda     TRUE
22     lightgbm     TRUE
6           mlp     TRUE
14   naivebayes     TRUE
5          nnet     TRUE
17          pls     TRUE
21          qda     TRUE
8  randomForest     TRUE
4        ranger     TRUE
2         rpart     TRUE
25     stacking     TRUE
26 superlearner     TRUE
24      xgboost     TRUE