RTFBoost: a package for (robust) tree-based functional boosting algorithms ================ Xiaomeng Ju and Matias Salibian Barrera 2026-09-12
This repository contains R code implementing a robust
tree-based boosting algorithm for scalar-on-function regression.
You can install the development version of the package in R using:
devtools::install_github("xmengju/RTFBoost")## DEoptimR (1.1-4 -> 1.2-1) [CRAN]
## mvtnorm (1.3-6 -> 1.4-2) [CRAN]
## Rcpp (1.1.1 -> 1.1.2) [CRAN]
## lme4 (2.0-1 -> 2.0-6) [CRAN]
##
## The downloaded binary packages are in
## /var/folders/wv/nltvdj6j4j5003tn7xxjybxr0000gn/T//RtmpqESttp/downloaded_packages
## ── R CMD build ─────────────────────────────────────────────────────────────────
## checking for file ‘/private/var/folders/wv/nltvdj6j4j5003tn7xxjybxr0000gn/T/RtmpqESttp/remotesa34153b2b8e2/xmengju-RTFBoost-7f4d934/DESCRIPTION’ ... ✔ checking for file ‘/private/var/folders/wv/nltvdj6j4j5003tn7xxjybxr0000gn/T/RtmpqESttp/remotesa34153b2b8e2/xmengju-RTFBoost-7f4d934/DESCRIPTION’
## ─ preparing ‘RTFBoost’:
## checking DESCRIPTION meta-information ... ✔ checking DESCRIPTION meta-information
## ─ checking for LF line-endings in source and make files and shell scripts
## ─ checking for empty or unneeded directories
## NB: this package now depends on R (>= 3.5.0)
## WARNING: Added dependency on R >= 3.5.0 because serialized objects in
## serialize/load version 3 cannot be read in older versions of R.
## File(s) containing such objects:
## ‘RTFBoost/data/fruitfly.RData’
## ─ building ‘RTFBoost_0.1.tar.gz’
##
## Once installed you can load the package with:
library(RTFBoost)Below we illustrate the use of the package with the fruit fly dataset. The original data is provided here. The data set consists of number of eggs laid daily for each of 1000 medflies until time of death. We used a subset of the data including flies that lived at least 30 days, to predict the number of living days. The number of eggs laid in the first 30 days is treated as the functional predictor with values evaluated at each day.
In the data set, fruitfly$eggs denotes the predictor,
and fruitfly$lifetime denotes the response. We plot the
data, highlighting 10 randomly chosen predictor curves in red. We also
plot the distribution of the lifetime variable, measured in days.
data(fruitfly)
matplot(t(fruitfly$eggs), lty = 1, type = "l", ylab = "price", xlab = "days",
main = "number of eggs laid each day (10 days)", col='gray70', ylim=c(0, 150))
set.seed(0)
n <- nrow(fruitfly$eggs)
matplot(t(fruitfly$eggs[sample(n, 10), ]), type='l', col='red', add=TRUE, lty=1)
hist(fruitfly$lifetime, xlab = "days",
main = "Histogram of the lifetime variable")
In order to train our predictor, we split the data set into a
training set (with 60% of the available data), a
validation set and a test set (both with 20%
of the data). We first randomly select the observations for each of
these three sets:
n0 <- floor( 0.2 * n)
set.seed(1)
idx_test <- sample(n, n0)
idx_train <- sample((1:n)[-idx_test], floor( 0.6 * n ) )
idx_val <- (1:n)[ -c(idx_test, idx_train) ] We now create the matrices of explanatory variables (x)
and vectors of responses (y) corresponding to this
partition.
xtrain <- fruitfly$eggs[idx_train, ]
ytrain <- fruitfly$lifetime[idx_train ]
xval <- fruitfly$eggs[idx_val, ]
yval <- fruitfly$lifetime[idx_val ]
xtest <- fruitfly$eggs[idx_test, ]
ytest <- fruitfly$lifetime[idx_test ]The RTFBoost function implements tree-based functional
boosting with four options:
TFBoost(LS): TFBoost with the squared loss;TFBoost(LAD): TFBoost with the L1 loss;RTFBoost(LAD-M): two-stage robust TFBoost that involves
a LAD-stage and an M-stage; andRTFBoost(RR): two-stage robust TFBoost that involves an
S-stage and an M-stage.These options correspond to setting control$type = 'L2',
control$type = 'LAD', control$type = 'LAD-M',
and control$type = ''RR' respectively.
We now explain how to fit a RTFBoost estimator with
different types and compare it with the fgam estimator
proposed in McLean el
al. (2014) and implemented in the refund package.
The RTFBoost functional allows using Type A or Type B
trees as the base learners. To specify the tree type, the user needs to
set tree.type = A or tree.type = B in
control$tree.control. Below, we will fit
RTFBoost with the type B trees, which train much faster
compared to the type A trees.
The following are parameters required for our estimator
tree.type <- "B" # type of the base learner
num.dir <- 200 # number of random directions for type B tree
gg <- 1:30 # specify the grid the functional predictor was evaluated on
tt <- c(0,30) # domain of the functional predictorn
niter <- 1000 # number of boosting iterations
make.prediction <- TRUE # make predictions based on test data
shrinkage <- 0.05 # shrinkage parameter for boosting
nknot <- 3 # the number of interior knots for cubic B-spline basis
init.type <- "median" # intialize with the median of the training responsesThe depth of the base learners in RTFBoost is set with
the argument max.depth in
control$tree.control. We considered max.depth
from 1 to 4, and chose the depth the minimizes the robust MSPE on the
validation set (\(\mathcal{I}_{\text{val}}\)) at early
stopping time. Denote validation residuals as \[
\begin{equation}
r_i = \hat{F}(x_i) - y_i, i \in \mathcal{I}_{\text{val}}
\end{equation}
\], the robust MSPE is defined as \[\hat{\mu}_M^2( \{r_i, i \in
\mathcal{I}_{\text{val}}\}) + \hat{\sigma}_M^2( \{r_i, i \in
\mathcal{I}_{\text{val}}\}), \] where \(\mu_M\) is the M-location estimator and
\(\sigma_M\) is the M-scale estimator.
We use Tukey’s score function and let asymptotic efficiency of \(\mu_M\) and \(\sigma_M\) be 95%. Below is the function to
compute robust MSPE:
cal.rmspe <- function(r){
tmp <- RobStatTM::locScaleM(x=r, psi='bisquare', eff= 0.95)
return( val.errors[dd] <- tmp$mu^2 + tmp$disper^2)
}Below we train RTFBoost with
control$type = 'L2' and select the depth
max.depth. The estimator is initialized at the median of
all training responses.
type <- "L2"
tree.depths <- 1:3
model.list.l2 <- vector('list', length(tree.depths))
val.errors <- rep(NA, length(tree.depths)) # vector of validation robust MSPE
for(dd in 1:length(tree.depths)) {
model.list.l2[[dd]] <-
RTFBoost(x.train = xtrain, y.train = ytrain, x.val = xval, y.val = yval,
x.test = xtest, y.test = ytest, grid = gg, t.range = tt,
control = RTFBoost.control(make.prediction = make.prediction, niter = niter,
tree.control = TREE.control(tree.type = tree.type, num.dir = num.dir, max.depth = dd),
shrinkage = shrinkage, nknot = nknot, type = type, init.type = init.type))
val.errors[dd] <- cal.rmspe(model.list.l2[[dd]]$f.val - yval)
}
model.l2 <- model.list.l2[[which.min(val.errors)]]
rm(model.list.l2)A sanity check
all.equal(mean((model.l2$f.test - ytest)^2), tail(model.l2$err.test,1))## [1] TRUEThen we repeat the same procedure fitting TFBoost(LAD)
and TFBoost(LAD-M).
type <- "LAD"
tree.depths <- 1:3
model.list.lad <- vector('list', length(tree.depths))
val.errors <- rep(NA, length(tree.depths)) # vector of validation robust MSPE
for(dd in 1:length(tree.depths)) {
model.list.lad[[dd]] <-
RTFBoost(x.train = xtrain, y.train = ytrain, x.val = xval, y.val = yval,
x.test = xtest, y.test = ytest, grid = gg, t.range = tt,
control = RTFBoost.control(make.prediction = make.prediction, niter = niter,
tree.control = TREE.control(tree.type = tree.type, num.dir = num.dir, max.depth = dd),
shrinkage = shrinkage, nknot = nknot, type = type, init.type = init.type))
val.errors[dd] <- cal.rmspe(model.list.lad[[dd]]$f.val - yval)
}
model.lad <- model.list.lad[[which.min(val.errors)]]
rm(model.list.lad)type <- "LAD-M"
n.init <- 500 # maximum number of iterations in the LAD-stage
niter <- 1000 # maximum number of iterations in the M-stage
tree.depths <- 1:3
model.list.ladm <- vector('list', length(tree.depths))
val.errors <- rep(NA, length(tree.depths)) # vector of validation robust MSPE
for(dd in 1:length(tree.depths)) {
model.list.ladm[[dd]] <-
RTFBoost(x.train = xtrain, y.train = ytrain, x.val = xval, y.val = yval,
x.test = xtest, y.test = ytest, grid = gg, t.range = tt,
control = RTFBoost.control(make.prediction = make.prediction, niter = niter, n.init = n.init,
tree.control = TREE.control(tree.type = tree.type, num.dir = num.dir, max.depth = dd),
shrinkage = shrinkage, nknot = nknot, type = type, init.type = init.type))
val.errors[dd] <- cal.rmspe(model.list.ladm[[dd]]$f.val - yval)
}
model.ladm <- model.list.ladm[[which.min(val.errors)]]
rm(model.list.ladm)For RTFBoost(RR) we controlled the level of robustness
of the S-stage by setting bb = 0.5 in the S-stage,
corresponding to a M-local estimator achieving 50% breakdown point. To
initialize RTFBoost(RR), we used a single-index Type A LAD
tree. We set the minimum number of observations per node of the initial
tree is to 10 and selected the maximum depth of the tree from 0, 1, 2,
and 3 to minimize the robust MSPE on the validation set. An initial tree
with depth 0 corresponds to initializing RTFBoost(RR) with
the median of all training responses. This step may take several minutes
to run.
n.init <- 200 # maximum number of iterations in the S-stage
niter <- 1000 # maximum number of iterations in the M-stage
type <- "RR"
tree.depths <- 1:3
model.list.rr <- vector('list', length(tree.depths))
val.errors <- rep(NA, length(tree.depths)) # vector of validation robust MSPE
for(dd in 1:length(tree.depths)) {
model.list.rr[[dd]] <-
RTFBoost.validation(x.train = xtrain, y.train = ytrain, x.val = xval, y.val = yval, x.test = xtest, y.test = ytest, grid = gg, t.range = tt,
max.depth.init.set = c(1,2,3), min.leafsize.init.set = c(10), control = RTFBoost.control(make.prediction = make.prediction, niter = niter, n.init = n.init,
tree.control = TREE.control(tree.type = tree.type, num.dir = num.dir, max.depth = dd),
shrinkage = shrinkage, nknot = nknot, type = type, init.type = "LADTree"), save.all = FALSE)
val.errors[dd] <- cal.rmspe(model.list.rr[[dd]]$model.best$f.val - yval)
}## [1] 1 1 1
## [1] 1
## [1] 2
## [1] 3
## [1] 1 1 1
## [1] 1
## [1] 2
## [1] 3
## [1] 1 1 1
## [1] 1
## [1] 2
## [1] 3model.rr <- model.list.rr[[which.min(val.errors)]]$model.best
params <- model.list.rr[[which.min(val.errors)]]$params
print(params)## min_leafs max_depths
## 1 10 1RTFBoost(RR) initialized with the median of training
responses.
n.init <- 200 # maximum number of iterations in the S-stage
niter <- 1000 # maximum number of iterations in the M-stage
type <- "RR"
tree.depths <- 1:3
model.list.rrm <- vector('list', length(tree.depths))
val.errors <- rep(NA, length(tree.depths)) # vector of validation robust MSPE
for(dd in 1:length(tree.depths)) {
model.list.rrm[[dd]] <-
RTFBoost(x.train = xtrain, y.train = ytrain, x.val = xval, y.val = yval, x.test = xtest, y.test = ytest, grid = gg, t.range = tt,
control = RTFBoost.control(type = type, make.prediction = make.prediction, niter = niter, n.init = n.init,
tree.control = TREE.control(tree.type = tree.type, num.dir = num.dir, max.depth = dd),
shrinkage = shrinkage, nknot = nknot, init.type = "median"))
val.errors[dd] <- cal.rmspe(model.list.rrm[[dd]]$f.val - yval)
}
model.rr.m <- model.list.rrm[[which.min(val.errors)]]We fit fgam which estimates the regression function
using tensor-product B-splines with roughness penalties. We set the
marginal basis to be cubic B-splines of dimension 15 and used REML to
select the parameter that penalizes the second order marginal
differences.
library(refund)
nbasis_FGAM <- 15
xx <- xtrain
yy <- ytrain
model.fgam <- fgam(yy ~ af(xx, splinepars = list(k=c(nbasis_FGAM,nbasis_FGAM),
m = list(c(2,2),c(2,2)))),
method ="REML")We compare test errors of different estimators:
err.l2 <-mean((model.l2$f.test- ytest)^2 )
err.lad <-mean( (model.lad$f.test- ytest)^2 )
err.ladm <-mean( (model.ladm$f.test- ytest)^2 )
err.rr <-mean( (model.rr$f.test- ytest)^2 )
err.rr.m <-mean( (model.rr.m$f.test- ytest)^2 )
err.fgam <- mean( (predict(model.fgam,
newdata=list(xx =xtest),
type='response') - ytest)^2 )
print(c(err.l2, err.lad, err.ladm, err.rr, err.rr.m, err.fgam))## [1] 150.3275 145.6992 145.3288 150.9795 146.8934 152.9956err.l2 <-cal.rmspe(model.l2$f.test- ytest)
err.lad <-cal.rmspe(model.lad$f.test- ytest)
err.ladm <-cal.rmspe(model.ladm$f.test- ytest)
err.rr <-cal.rmspe(model.rr$f.test- ytest)
err.rr.m <-cal.rmspe(model.rr.m$f.test- ytest)
err.fgam <- cal.rmspe( (predict(model.fgam,
newdata=list(xx =xtest),
type='response') - ytest))
print(c(err.l2, err.lad, err.ladm, err.rr, err.rr.m, err.fgam))## [1] 116.33195 97.48760 97.96605 103.63262 95.22485 127.61979For RTFBoost, we can also separate the process of
training the predictor and evaluating it on a test set. In this way we
have the flexibility to compute predictions on multiple test sets when
needed. In this case the user needs to use save.tree = TRUE
in the argument control. This option includes the tree
objects at all iterations in the returned object, which are needed for
calculating predictions. The following code illustrates this procedure
running TFBoost(L2) with max.depth = 1:
max.depth <- 1
type <- "L2"
model.l2.separate <- RTFBoost(x.train = xtrain, y.train = ytrain, x.val = xval, y.val = yval,
x.test = xtest, y.test = ytest, grid = gg, t.range = tt,
control = RTFBoost.control(make.prediction = TRUE, niter = niter,
tree.control = TREE.control(tree.type = tree.type, num.dir = num.dir, max.depth = max.depth),
shrinkage = shrinkage, nknot = nknot, type = type, init.type = "median", save.tree = TRUE))
predictions <- RTFBoost.predict(model.l2.separate, newx = xtest, grid = gg, t.range = tt)A sanity check
all.equal(mean((predictions$f.new - ytest)^2), tail(model.l2.separate$err.test, 1))## [1] TRUEall.equal(mean((model.l2.separate$f.test - ytest)^2), tail(model.l2.separate$err.test, 1))## [1] TRUE