\documentclass[a4paper]{article}
\usepackage{jmlr2e}
\usepackage{color,hyperref}
\usepackage{amsmath}
\usepackage{amssymb}
\newcommand{\argmin}{\operatorname{argmin}\displaylimits}

%%\VignetteIndexEntry{mboost}
%%\VignetteDepends{mboost, partykit, TH.data}

%% software markup
\let\proglang=\textsf
\let\pkg=\textit
\let\code=\texttt

%% \usepackage{Sweave} is essentially
\RequirePackage[T1]{fontenc}
\RequirePackage{ae,fancyvrb}
\DefineVerbatimEnvironment{Sinput}{Verbatim}{fontshape=sl,fontsize=\small}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{fontsize=\small}
\DefineVerbatimEnvironment{Scode}{Verbatim}{fontshape=sl,fontsize=\small}
\newenvironment{Schunk}{}{}
\setkeys{Gin}{width=.49\textwidth}
\SweaveOpts{engine=R, eps=FALSE, keep.source = TRUE}
<<options, echo=FALSE>>=
options(prompt = "R> ", continue = "+  ", digits = 4)
set.seed(290875)
## use old degrees of freedom for reproducibility
options(mboost_dftraceS = TRUE)
@

%% header declarations
\author{\name Torsten Hothorn \addr LMU M\"unchen \email Torsten.Hothorn@R-project.org
 \AND \name Peter B\"uhlmann \addr ETH Z\"urich \email buhlmann@stat.math.ethz.ch
 \AND \name Thomas Kneib \addr Universit\"at Oldenburg \email thomas.kneib@uni-oldenburg.de
 \AND \name Matthias Schmid \addr FAU Erlangen-N\"urnberg \email Matthias.Schmid@imbe.med.uni-erlangen.de
 \AND \name Benjamin Hofner \addr FAU Erlangen-N\"urnberg \email Benjamin.Hofner@imbe.med.uni-erlangen.de}

\ShortHeadings{Model-based Boosting 2.0}{Hothorn et al.}
\firstpageno{1}

\title{Model-based Boosting 2.0}
\editor{Mikio Braun}
% hyperref setup
\definecolor{Red}{rgb}{0.5,0,0}
\definecolor{Blue}{rgb}{0,0,0.5}
\hypersetup{%
  pdftitle = {Model-Based Boosting 2.0},
  pdfsubject = {submitted to JMLR-MLOSS},
  pdfkeywords = {open-source software, functional gradient descent, statistical learning},
  pdfauthor = {Torsten Hothorn, Peter Buehlmann, Thomas Kneib, Matthias Schmid and Benjamin Hofner},
  %% change colorlinks to false for pretty printing
  colorlinks = {true},
  linkcolor = {Blue},
  citecolor = {Blue},
  urlcolor = {Red},
  hyperindex = {true},
  linktocpage = {true},
}

\begin{document}

\maketitle

\begin{center}
\small
\textbf{This is an extended version of the manuscript}\\
Torsten Hothorn, Peter B\"uhlmann, Thomas Kneib, Mattthias Schmid
and Benjamin Hofner (2010),  Model-based Boosting 2.0. \emph{Journal of Machine Learning Research}, \textbf{11},
2109 -- 2113;\\
URL \url{http://jmlr.csail.mit.edu/papers/v11/hothorn10a.html}.
\end{center}
\bigskip

\begin{abstract}
This is an extended version of the manuscript Torsten Hothorn, Peter B\"uhlmann, Thomas Kneib, Mattthias Schmid
and Benjamin Hofner (2010),  Model-based Boosting 2.0. \emph{Journal of Machine Learning Research}, \textbf{11},
2109 -- 2113; \url{http://jmlr.csail.mit.edu/papers/v11/hothorn10a.html}.

We describe version 2.0 of the \proglang{R} add-on package \pkg{mboost}.
The package implements boosting for optimizing general risk functions utilizing
component-wise (penalized) least squares estimates or regression
trees as base-learners for fitting generalized linear, additive
and interaction models to potentially high-dimensional data.
\end{abstract}

\begin{keywords}
component-wise functional gradient descent, splines, decision trees.
\end{keywords}

\section{Overview}
The \proglang{R} add-on package \pkg{mboost} \citep{PKG:mboost} implements tools for fitting
and evaluating
a variety of regression and classification models that have been suggested in machine learning and statistics.
Optimization within the empirical risk minimization framework is performed
via a component-wise functional gradient descent algorithm.
The algorithm originates from the statistical view on boosting
algorithms \citep{fht00,pbyu03}. The theory and its implementation in
\pkg{mboost} allow for fitting complex prediction models,
taking potentially many interactions of features into account, as well as
for fitting additive and linear models.
The model class the package deals with is best described by so-called structured
additive regression (STAR) models, where some characteristic $\xi$ of the
conditional distribution of a response variable $Y$ given features $X$ is
modeled through a regression function $f$ of the features
$\xi(Y | X = x) = f(x)$.
In order to facilitate parsimonious and interpretable models, the regression function
$f$ is structured, i.e., restricted to additive functions $f(x) = \sum_{j = 1}^p f_j(x)$.
Each model component $f_j(x)$ %may or may not
might take only a subset of the
features into account. Special cases are linear models $f(x) = x^\top \beta$, additive
models $f(x) = \sum_{j = 1}^p f_j(x^{(j)})$, where $f_j$ is a function of the $j$th feature
$x^{(j)}$ only (smooth functions or stumps, for example) or a more complex function where
$f(x)$ is implicitly defined as the sum of multiple decision trees including higher-order
interactions. The latter case corresponds to boosting with trees. Combinations
of these structures are also possible.
The most important advantage of such a decomposition of the regression function is that
each component of a fitted model can be looked at and interpreted separately for gaining
a better understanding of the model at hand.

The characteristic $\xi$ of the distribution depends on the measurement scale
of the response $Y$ and the scientific question to be answered. For binary or numeric
variables, some function of the expectation
may be appropriate, but also quantiles or expectiles may be interesting. The definition
of $\xi$ is determined by defining a loss function $\rho$ whose empirical
risk is to be minimized under some algorithmic constraints (i.e., limited
number of boosting iterations). The model is then fitted using
$$
(\hat{f}_1, \dots, \hat{f}_p) = \argmin_{(f_1, \dots, f_p)}
    \sum_{i = 1}^n w_i \rho\left(y_i, \sum_{j = 1}^p f_j(x)\right).
$$
Here $(y_i, x_i), i = 1, \dots, n$, are $n$ training samples with responses $y_i$ and
potentially high-dimensional feature vectors $x_i$, and $w_i$ are some weights.
The component-wise boosting algorithm starts with some offset for $f$ and iteratively fits residuals defined
by the negative gradient of the loss function evaluated at the current fit by updating
only the best model component in each iteration. The
details have been described by \cite{pbyu03}. Early stopping via resampling approaches or
AIC leads to sparse models in the sense that only a subset of important model components
$f_j$ defines the final model. A more thorough introduction to boosting with applications
in statistics based on version 1.0 of \pkg{mboost} is given by \cite{BuhlmannHothorn06}.

As of version 2.0, the package allows for fitting models to binary, numeric, ordered and
censored responses, i.e., regression of the mean, robust regression, classification
(logistic and exponential loss), ordinal regression\footnote{Model family is new in version
2.0 or was added after the release of \pkg{mboost} 1.0.},
quantile\footnotemark[\value{footnote}] and expectile\footnotemark[\value{footnote}] regression,
censored regression (including Cox, Weibull\footnotemark[\value{footnote}],
log-logistic\footnotemark[\value{footnote}] or lognormal\footnotemark[\value{footnote}] models)
as well as Poisson and negative binomial regression\footnotemark[\value{footnote}] for count
data can be performed. Because the structure of the regression function
$f(x)$ can be chosen independently from the loss function $\rho$, interesting new models can be fitted
\citep[e.g., in geoadditive regression,][]{Kneib+Hothorn+Tutz:2009}.

\section{Design and Implementation}

The package incorporates an infrastructure for representing loss functions
(so-called `families'), base-learners defining the structure of the
regression function and thus the model components $f_j$, and a generic
implementation of component-wise functional gradient descent.  The main
progress in version 2.0 is that only one
implementation of the boosting algorithm is applied to all possible models
(linear, additive, tree-based) and all families.  Earlier versions were
based on three implementations, one for linear models, one for additive
models, and one for tree-based boosting.  In comparison to the 1.0 series,
the reduced code basis is easier to maintain, more robust and regression
tests have been set-up in a more unified way.  Specifically, the new code
basis results in an enhanced and more user-friendly formula interface.  In
addition, convenience functions for hyper-parameter selection, faster
computation of predictions and improved visual model diagnostics are
available.

Currently implemented base-learners include component-wise linear models
(where only one variable is updated in each iteration of the algorithm),
additive models with quadratic penalties \citep[e.g., for fitting smooth
functions via penalized splines, varying coefficients or bi- and trivariate tensor
product splines,][]{Schmid+Hothorn:2008b}, and trees.

As a major improvement over the 1.0 series, computations on larger data sets
(both with respect to the number of observations and the number of variables) are now
facilitated by memory efficient implementations of the base-learners, mostly
by applying sparse matrix techniques \citep[package
\pkg{Matrix},][]{PKG:Matrix} and parallelization for a
cross-validation-based choice of the number of boosting iterations
per default via package \pkg{parallel}. A more
elaborate description of \pkg{mboost} 2.0 features is available
from the \texttt{mboost} vignette\footnote{Accessible via \code{vignette("mboost", package = "mboost")}.}.

%% Regression tests

\section{User Interface by Example}

We illustrate the main components of the user-interface by
a small example on human body fat composition: \cite{garcia2005}
used a linear model for predicting body fat content by means of common anthropometric
measurements that were obtained for $n = 71$ healthy German women.
In addition, the women's body composition was measured by
Dual Energy X-Ray Absorptiometry (DXA). The aim is to describe
the DXA measurements as a function of the anthropometric features.
Here, we extend the linear model by i) an intrinsic variable
selection via early stopping, ii) additional terms allowing for
smooth deviations from linearity where necessary \citep[by means of penalized splines orthogonalized to the linear effect,][]{Kneib+Hothorn+Tutz:2009}, iii) a possible
interaction between two variables with known impact on body fat composition
(hip and waist circumference) and iv) utilizing a robust
median regression approach instead of $L_2$ risk.
For the data (available as data frame \code{bodyfat}), the model structure is specified
via a formula involving the base-learners corresponding to the
different model components (linear terms: \code{bols()}; smooth
terms: \code{bbs()}; interactions: \code{btree()}).
The loss function (here, the check function for the $0.5$ quantile)
along with its negative gradient function are defined by the
\code{QuantReg(0.5)} family \citep{Fenske+Kneib+Hothorn:2009}.
The model structure (specified using the formula \code{fm}), the data and the family
are then passed to function \code{mboost()} for model fitting:
<<setup, echo = FALSE>>=
### make sure package mboost is available
if (!require("mboost"))
    install.packages("mboost", dependencies = TRUE)
library("mboost")
### package partykit is necessary for fitting trees
if (!require("partykit"))
    install.packages("partykit", dependencies = TRUE)
library("partykit")
### speed up things a little bit: if the model is already
### available load it into R
if (file.exists("model.Rda")) {
    load("model.Rda")
    mboost <- function(...) return(model)
    cvrisk <- function(...) return(cvm)
}
@
<<mboost-bodyfat-setup, echo = FALSE>>=
### set-up model formula
### names of features
data("bodyfat", package = "TH.data")
features <- names(bodyfat)[-2]
### set up model structure:
fml <- paste("bols(", features, ")", collapse = " + ") ### linear functions
fms <- paste("bbs(", features, ", center = TRUE, df = 1)",
             collapse = " + ")  ### smooth deviations from linearity
fmt <- "btree(hipcirc, waistcirc, tree_controls = ctree_control(maxdepth = 2, mincriterion = 0))" ### tree-based interaction
fm <- as.formula(paste("DEXfat", paste(fml, fms, fmt, sep = "+"), sep = "~"))
@
<<mboost-fm>>=
library("mboost")               ### attach package `mboost'
print(fm)                       ### model structure
@
<<mboost-bodyfat, results = hide>>=
### fit model for conditional median of DEXfat
model <- mboost(fm,                           ### model structure
                data = bodyfat,               ### 71 observations
                family = QuantReg(tau = 0.5)) ### median regression
@
<<mboost-mstop, echo = FALSE, results = hide>>=
model[1000]                     ### 900 more iterations
@
Once the model has been fitted
%% using $1,000$ initial boosting iterations,
it is important to assess the appropriate number of boosting
iterations via the out-of-sample empirical risk. By default,
$25$ bootstrap samples from the training data are drawn and the
out-of-bag empirical risk is computed (parallel computation if possible):
<<mboost-cvrisk, echo = TRUE, results = hide>>=
### bootstrap for assessing the `optimal' number of boosting iterations
cvm <- cvrisk(model, grid = 1:100 * 10)
model[mstop(cvm)]               ### restrict model to optimal mstop(cvm) iterations
@
Now, the final model is ready for a visual inspection:
<<bodyfat-plot, echo = TRUE, eval = FALSE>>=
plot(cvm); plot(model)          ### depict out-of bag risk & selected components
@
The resulting plots are given in Figure~\ref{mboostplot}. They indicate that
a model based on three components, including a smooth function of
\code{anthro3b} and a bivariate function of hip and waist circumference,
provides the best characterization of the median body fat composition (given
the model specification offered to the boosting algorithm).  A hip
circumference larger than $110$ cm leads to increased body fat but only if
the waist circumference is larger than $90$ cm.

\begin{figure}
\begin{center}
<<mboost-bodyfat-plot-1, echo = FALSE, fig = TRUE, width = 6, height = 6>>=
### plot age and kneebreadth
cex <- 1.3
layout(matrix(c(1, 2, 1, 3), nr = 2))
par(mai = par("mai") * c(0.8, 1.1, 0.8, 0.8))
plot(cvm, cex.lab = cex)
mtext(text = "(A)", side = 3, 1, adj = 1)
plot(model, which = "bols(anthro3b", cex.lab = cex)
mtext(text = "(B)", side = 3, 1, adj = 1)
plot(model, which = "bbs(anthro3b",
     ylim = range(predict(model, which = "bols(anthro3b")),
     cex.lab = cex)
mtext(text = "(C)", side = 3, 1, adj = 1)
@
<<mboost-bodyfat-plot-2, echo = FALSE, results = hide, fig = TRUE>>=
### plot interaction of hip and waist circumference
### first setup grid of hip and waist values
nd <- with(bodyfat,
           expand.grid(hipcirc = h <- seq(from = min(hipcirc),
                                          to = max(hipcirc),
                                          length = 100),
                  waistcirc = w <- seq(from = min(waistcirc),
                                       to = max(waistcirc),
                                       length = 100)))
### define colors for plot
col <-
c("#023FA5", "#1141A4", "#1A44A4", "#2146A4", "#2749A4", "#2C4BA4",
"#304DA4", "#3550A5", "#3852A5", "#3C54A6", "#4056A6", "#4359A7",
"#465BA7", "#495DA8", "#4C5FA9", "#4F61AA", "#5264AA", "#5566AB",
"#5868AC", "#5B6AAD", "#5D6CAE", "#606EAE", "#6270AF", "#6572B0",
"#6775B1", "#6A77B2", "#6C79B3", "#6F7BB4", "#717DB5", "#747FB6",
"#7681B6", "#7883B7", "#7B85B8", "#7D87B9", "#7F89BA", "#828BBB",
"#848DBC", "#868FBD", "#8891BE", "#8A93BE", "#8D94BF", "#8F96C0",
"#9198C1", "#939AC2", "#959CC3", "#979EC4", "#99A0C4", "#9BA1C5",
"#9DA3C6", "#9FA5C7", "#A1A7C8", "#A3A8C9", "#A5AAC9", "#A7ACCA",
"#A9AECB", "#ABAFCC", "#ACB1CC", "#AEB3CD", "#B0B4CE", "#B2B6CF",
"#B4B8CF", "#B5B9D0", "#B7BBD1", "#B9BCD2", "#BABED2", "#BCBFD3",
"#BEC1D4", "#BFC2D4", "#C1C4D5", "#C3C5D6", "#C4C7D6", "#C6C8D7",
"#C7C9D7", "#C9CBD8", "#CACCD9", "#CBCDD9", "#CDCFDA", "#CED0DA",
"#CFD1DB", "#D1D2DB", "#D2D3DC", "#D3D4DC", "#D4D6DD", "#D6D7DD",
"#D7D8DE", "#D8D9DE", "#D9DADF", "#DADBDF", "#DBDCDF", "#DCDCE0",
"#DDDDE0", "#DEDEE0", "#DEDFE1", "#DFDFE1", "#E0E0E1", "#E1E1E2",
"#E1E1E2", "#E2E2E2", "#E2E2E2", "#E2E2E2", "#E2E2E2", "#E2E2E2",
"#E2E2E2", "#E2E1E1", "#E2E0E1", "#E2E0E0", "#E1DFDF", "#E1DEDF",
"#E1DDDE", "#E1DCDD", "#E0DBDC", "#E0DADB", "#E0D9DA", "#DFD8D9",
"#DFD7D8", "#DFD6D7", "#DED5D6", "#DED3D5", "#DDD2D4", "#DDD1D3",
"#DDCFD2", "#DCCED0", "#DCCDCF", "#DBCBCE", "#DBCACD", "#DAC8CB",
"#DAC7CA", "#D9C5C8", "#D9C4C7", "#D8C2C6", "#D8C0C4", "#D7BFC3",
"#D7BDC1", "#D6BBC0", "#D5B9BE", "#D5B8BD", "#D4B6BB", "#D3B4B9",
"#D3B2B8", "#D2B0B6", "#D1AEB4", "#D1ADB3", "#D0ABB1", "#CFA9AF",
"#CEA7AE", "#CEA5AC", "#CDA3AA", "#CCA1A8", "#CB9FA7", "#CB9CA5",
"#CA9AA3", "#C998A1", "#C8969F", "#C7949D", "#C6929C", "#C5909A",
"#C48D98", "#C38B96", "#C38994", "#C28792", "#C18490", "#C0828E",
"#BF808C", "#BE7D8A", "#BD7B88", "#BB7986", "#BA7684", "#B97482",
"#B87180", "#B76F7E", "#B66C7C", "#B56A7A", "#B46777", "#B26575",
"#B16273", "#B05F71", "#AF5D6F", "#AE5A6D", "#AC576B", "#AB5569",
"#AA5266", "#A84F64", "#A74C62", "#A64960", "#A4475E", "#A3445B",
"#A24159", "#A03D57", "#9F3A55", "#9D3752", "#9C3450", "#9A304E",
"#992C4C", "#982949", "#962447", "#942045", "#931B42", "#911640",
"#900F3E", "#8E063B")
### use plot method to draw fitted values of the tree component only
print(plot(model, which = "btree", newdata = nd, col.regions = col,
     at = seq(from = -16, to = 16, length = 100)))
### save model for future use
save(model, cvm, file = "model.Rda")
@
\caption{Out-of-bag empirical risk (A) indicating that
         \Sexpr{mstop(cvm)} iterations are appropriate.
         Fitted model components for variable \code{anthro3b},
         consisting of a linear (B) and smooth term (C).
         The right panel shows the interaction model component
         between hip and waist circumferences. \label{mboostplot}}
\end{center}
\end{figure}

The sources of the \pkg{mboost} package are distributed
at the Comprehensive \proglang{R} Archive Network under GPL-2,
along with binaries for all major platforms as well as documentation and regression tests.
Development versions are available from \url{http://R-forge.R-project.org}.


\section{Overview on 2.0 Series Features}
\label{sec:2.0-series-features}

This additional section gives some further details on new features in the \pkg{mboost}
2.0 series. First of all, new families were included to permit quantile and
expectile regression (\code{QuantReg()}, \code{ExpectReg()}), ordinal regression
(\code{PropOdds()}), censored regression (\code{Weibull()}, \code{Loglog()},
\code{Lognormal()}) and count data regression (\code{NBinomial()}). All new
families, except the quantile and expectile family, have an additional scale
parameter, which is included in the family and is then subsequently estimated in
the boosting algorithm without further need of modifications thereof.

For increased usability of \pkg{mboost}, some changes in the user interface of
the base-learners were necessary. Most prominently, the \code{center} argument
in \code{bols()} was renamed to \code{intercept} to reflect that it is used to
specify whether the base-learner should contain an intercept (\code{intercept =
  TRUE}) or not. The argument \code{by} (formerly denoted by \code{z} in \pkg{mboost} 1.0),
which can be used to specify varying coefficients, is now able to handle factors
with more than two levels in addition to binary and continuous covariates.
Furthermore, the base-learners \code{bss()} (smoothing spline) and \code{bns()}
(penalized natural splines) are deprecated and replaced by \code{bbs()}
(penalized splines), which results in qualitatively the same models but is
computationally much more attractive.

The speed of model estimation and prediction has been notably improved in
version 2.0. Regarding large data sets, a newly added search for duplicated observations of covariates in each
base-learner leads to enormous speed-ups and far less memory consumption in the
base-learner fitting in each boosting iteration. The effective number of
observations in a weighted penalized least-squares problem is the number of
unique covariate observations which is usually considerably smaller than the number of
observations in the learning sample, especially for factors or smooth effects
for millions of observations. The price to pay is an increased pre-processing
time, however, pre-processing is just done once at the beginning of the
algorithm.

In addition, version 2.0 of \pkg{mboost} offers an improved infrastructure:
Firstly, the subset method \code{model[i]} can now be used to
\emph{enhance} or restrict a given boosting model to the specified boosting
iteration \code{i}. Note that in both cases the \emph{original} \code{model}
will be changed to reduce the memory footprint. If the boosting model is
enhanced by specifying an index that is larger than the initial \code{mstop},
only the missing \code{i - mstop} steps are fitted. If the model is restricted,
the spare steps are not dropped, but only hidden from the user, i.e., if we
increase \code{i} again, these boosting steps are immediately available.
The newly introduced function \code{selected()} allows to inspect the path
of selected base-learners and is the basis for stability selection
\citep[in function \code{stabsel()}, cf.][]{mebu10}.
Secondly, an interface for parallel computing (based on \pkg{parallel}) is
automatically used (if available) for cross-validation-based stopping of the
algorithm in \code{cvrisk()}. With some extra work, other parallelization packages
such as \pkg{snow} \citep{PKG:snow} can immediately be integrated.
Moreover, \code{cvrisk()} allows for more flexible definitions of resampling schemes
($k$-fold, cross-validation, bootstrap, subsampling).
Thirdly, the
\code{predict()}, \code{plot()} and \code{coef()} functions were enhanced.
For linear models, the intercept is now adjusted for centered covariates and
non-zero offsets are added to the intercept when using \code{coef(..., off2int = TRUE)}.
Predictions are computed much faster now and for all three functions an
argument \code{which} was included to allow to specify which base-learner(s)
should be used for prediction and plotting or which coefficients should be
extracted. Per default all \emph{selected} base-learners are used. Users can
specify \code{which} as numeric value(s) or as a (vector of) character
string(s):
<<mboost-predict>>=
### new data on a grid on range(anthro3b)
nd <- with(bodyfat, data.frame(anthro3b = seq(min(anthro3b), max(anthro3b),
                                              length = 100)))
### predictions for all base-learners of `anthro3b'
pr <- predict(model, which = "anthro3b", newdata = nd)
pr <- rowSums(pr)    ### aggregate linear and smooth effect
@
The combined effect can now be used for plotting (see Figure~\ref{fig:anthro3b}).

\begin{figure}[ht!]
\begin{center}
<<mboost-anthro3b, echo = TRUE, fig = TRUE, width = 4, height = 4>>=
plot(nd$anthro3b, pr, type = "l", xlab = "anthro3b",
     ylab = "f(anthro3b)")
lines(nd$anthro3b, predict(model, which = "bols(anthro3b", newdata = nd),
      type = "l", lty = "dashed")
@
\caption{Effect of \code{anthro3b}, i.e., the combination of the linear and
  smooth effect (solid line) and for comparison solely the linear effect
  (dashed line). \label{fig:anthro3b}}
\end{center}
\end{figure}

\bibliography{boost}

\end{document}
