\documentclass[nojss]{jss}

%\VignetteIndexEntry{party with the mob}

%\VignetteDepends{mlbench,colorspace}

%\VignetteKeywords{parametric models, object-orientation, recursive partitioning}

%\VignettePackage{party}

%% packages

\usepackage{amsmath}

\SweaveOpts{engine=R}

%% neet no \usepackage{Sweave}

<<setup, echo = FALSE, results = hide>>=

require("party")

options(useFancyQuotes = FALSE)

@

%% commands

\newcommand{\ui}{\underline{i}}

\newcommand{\oi}{\overline{\imath}}

\newcommand{\argmin}{\operatorname{argmin}\displaylimits}

%% author/title

\author{Achim Zeileis\\Universit\"at Innsbruck \And

        Torsten Hothorn\\Ludwig-Maximilians-\\Universit\"at M\"unchen \And

    Kurt Hornik\\WU Wirtschafts-\\universit\"at Wien}

\Plainauthor{Achim Zeileis, Torsten Hothorn, Kurt Hornik}

\title{\pkg{party} with the \code{mob}: Model-Based Recursive Partitioning in \proglang{R}}

\Plaintitle{party with the mob: Model-Based Recursive Partitioning in R}

\Shorttitle{\pkg{party} with the \texttt{mob}}

\Keywords{parametric models, object-orientation, recursive partitioning}

%% abstract

\Abstract{

  The \pkg{party} package \citep{Hothorn+Hornik+Zeileis:2006a} provides

  the function \code{mob()} implementing a recently suggested algorithm

  for \underline{mo}del-\underline{b}ased recursive partitioning

  \citep{Zeileis+Hothorn+Hornik:2005}. The basic steps are: (1)~fit a

  parametric model to a data set, (2)~test for parameter instability over

  a set of partitioning variables, (3)~if there is some overall parameter

  instability, split the model with respect to the variable associated with

  the highest instability, (4)~repeat the procedure in each of the child nodes.

  It is discussed how these steps of the conceptual algorithm are translated

  into computational tools in an object-oriented manner, allowing the user to

  plug in various types of parametric models. The outcome is a tree where

  each node is associated with a fitted parametric model that can be

  effectively visualized and summarized.

}

\Address{

  Achim Zeileis\\

  Department of Statistics\\

  Faculty of Economics and Statistics\\

  Universit\"at Innsbruck\\

  Universit\"atsstr.~15\\

  6020 Innsbruck, Austria\\

  E-mail: \email{Achim.Zeileis@R-project.org}\\

  URL: \url{http://eeecon.uibk.ac.at/~zeileis/}\\

  Torsten Hothorn\\

  Institut f\"ur Statistik\\

  Ludwig-Maximilians-Universit\"at M\"unchen\\

  Ludwigstra{\ss}e 33\\

  80539 M\"unchen, Germany\\

  E-mail: \email{Torsten.Hothorn@R-project.org}\\

  URL: \url{http://www.stat.uni-muenchen.de/~hothorn/}\\

  Kurt Hornik\\

  Institute for Statistics and Mathematics\\

  WU Wirtschaftsuniversit\"at Wien\\

  Augasse 2--6\\

  1090 Wien, Austria\\

  E-mail: \email{Kurt.Hornik@R-project.org}\\

  URL: \url{http://statmath.wu.ac.at/~hornik/}\\

}

\begin{document}

\section{Motivation}

\label{sec:motivation}

Consider a parametric model $\mathcal{M}(Y, \theta)$

with (possibly vector-valued) observations $Y$ and a

$k$-dimensional vector of parameters $\theta$. This model 

could be a (possibly multivariate) normal distribution for $Y$, or some

kind of regression model when $Y = (y, x)$ can be split up into a dependent variable

$y$ and regressors $x$. An example for the latter could be a linear regression

model $y = x^\top \theta$ or a generalized linear model (GLM) or a survival 

regression.

Given $n$ observations $Y_i$ ($i = 1, \dots, n$) the model can be fitted

by minimizing some objective function $\Psi(Y, \theta)$, e.g., a residual sum of squares

or a negative log-likelihood leading to ordinary least squares (OLS) or maximum

likelihood (ML) estimation.

If a global model for all $n$ observations does not fit well and further

covariates $Z_1, \dots, Z_\ell$ are available, it might be possible to partition

the $n$ observations with respect to these variables and find a fitting model

in each cell of the partition. The algorithm described here tries to find

such a partition adaptively using a greedy forward search. This procedure is

implemented in the function \code{mob()} and described in more detail in the following

section. However, we we will state some goals and design principles in advance.

To translate the model-based partitioning problem into \proglang{R}, we start with

a formula description of the variables involved. This formula should be of type

\verb:y ~ x1 + ... + xk | z1 + ... + zl: where the variables on the

left of the \code{|} specify the data $Y$ and the variables on the right specify the

partitioning variables $Z_j$. Classical regression trees usually have a univariate

response $Y$ and various partitioning variables, i.e., could be specified as

\verb:y ~ 1 | z1 + ... + zl:. Structural change models, on the other hand, are usually

regression models that are segmented with respect to a single partitioning variable,

typically time: \verb:y ~ x1 + ... + xk | z:.

The type of models $\mathcal{M}$ to be used with \code{mob()} should not be

confined (by the implementation), hence we have written an object-oriented 

implementation. The idea is that $\mathcal{M}$ is translated into software

by a model of class ``\code{StatModel}'' as provided by the \pkg{modeltools}

package. The algorithm the relies on various methods being available for these

models. The ``\code{StatModel}'' objects \code{linearModel} and \code{glinearModel},

implementing (generalized) linear regression models, are readily available in

\pkg{modeltools}, others can easily be user-defined.

\section{The model-based recursive partitioning algorithm}

\label{sec:algorithm}

The basic idea is to grow a tee in which every node is associated with a model

of type $\mathcal{M}$. To assess whether splitting of the node is necessary a

fluctuation test for parameter instability is performed. If there is significant instability

with respect to any of the partitioning variables $Z_j$, the node is splitted

into $B$ locally optimal segments (currently only $B = 2$ is implemented)

and then the procedure is repeated in each of the $B$ children.

If no more significant instabilities can be found, the recursion stops.

More precisely, the steps of the algorithm are

\begin{enumerate}

\item Fit the model once to all observations in the current node.

\item Assess whether the parameter estimates are stable with respect to

  every partitioning variable $Z_1, \dots, Z_\ell$. If there is some overall instability,

  select the variable $Z_j$ associated with the highest parameter instability, otherwise

  stop.

\item Compute the split point(s) that locally optimize the objective function $\Psi$.

\item Split the node into child nodes and repeat the procedure.

\end{enumerate}

The details for steps 1--3 are specified in the following.

\subsection{Parameter estimation}

This step of the algorithm is common practice, the only additional 

requirement is (as previously noted) that model has to be of the class

``\code{StatModel}'' as provided by \pkg{modeltools}. Looking at the source

code for the \code{linearModel} provided by this package illustrates how

a simple wrapper to existing \proglang{R} functionality can be written.

In particular, a method to the generic function \code{reweight()} has to

be available. The reason is that it is inefficient to fit a brand-new model

\code{modelobj} (including formula-parsing) in every node---much computation time

is saved if simply \code{reweight(modelobj, weights)} is called in each of the

child nodes. The \code{weights} argument controls which observations go into which

of the child nodes.

\subsection{Testing for parameter instability}

The task in this step of the algorithm is to find out whether the parameters

of the fitted model are stable over each particular ordering implied by

the partitioning variables $Z_j$ or whether splitting the sample with respect

to one of the $Z_j$ might capture instabilities in the parameters and thus improve the fit.

The tests used in this step belong to the class of generalized M-fluctuation

tests \citep{ZeileisHornik2003,Zeileis2005}. For numerical partitioning variables

$Z_j$ the $\sup LM$~statistic is used which is the maximum over all single split

$LM$ statistics. For categorical partitioning variables, a $\chi^2$~statistic is

employed which captures the fluctuation within each of the categories of $Z_j$.

For computing the test statistics and corresponding $p$~values $p_j$ for 

each of the partitioning variables $Z_j$ in \proglang{R}, the only requirement

is that there are methods for the extractor functions \code{estfun()} and 

\code{weights()}. The \code{estfun()} method extracts the empirical estimating

functions (model scores) from the fitted \code{modelobj}, these are the 

main ingredient for M-fluctuation tests. The \code{weights()} method is used

to determine which observations are in the current node (i.e., have a weight

greater than zero) and which are not (i.e., have zero weight).

To determine whether there is some overall instability, it is checked whether

the minial $p$~value $p_{j^*} = \min_{j = 1, \dots, \ell} p_j$ falls below a

pre-specified significance level $\alpha$ (by default $\alpha = 0.05$) or not.

To adjust for multiple testing, the $p$ values can be Bonferroni adjusted

(which is the default). If there is significant instability, the variable $Z_{j^*}$

associated with the minimal $p$~value is used for splitting the node.

\subsection{Splitting}

In this step, the observation in the current node are split with respect

to the chosen partitioning variable $Z_{j^*}$ into $B$ child nodes. Currently,

the infrastructure in \pkg{party} only supports binary splits, i.e., $B = 2$.

For deterimining the split point, an exhaustive search procedure is adopted:

For each conceivable split point, the $B$ child node models are fit and the split

associated with the minimal value of the objective function $\Psi$ is chosen.

Computationally, this means that the fitted model \code{modelobj} is \code{reweight()}ed

for each of the child nodes. The observations entering a child node keep their

current weight while those observations that go into different child nodes receive

zero weight. To compare the objective function $\Psi$, an extractor function

is required to compute it from the fitted \code{modelobj}. This extractor function

can be user-specified and set in \verb:mob_control():, it defaults to \code{deviance()}.

This concludes one iteration of the recursive partitioning algorithm and steps~1--3

are carried out again in each of the $B$ daughter nodes until no significant 

instability is detected in step~2.

\section{Illustrations}

\label{sec:illustration}

\subsection{Boston housing data}

Since the analysis by \cite{BreimanFriedman1985}, the Boston housing data are 

a popular and well-investigated empirical basis for illustrating non-linear 

regression methods both in machine learning and statistics

\citep[see][for two recent examples]{Gama2004,Samarovetal2005} and we follow 

these examples by segmenting a bivariate linear regression model for the house

values.

Thus, the model $\mathcal{M}$ used is \code{linearModel} from the \pkg{modeltools}

package which is automatically loaded together with \pkg{party}.

<<>>=

library("party")

@

The data set is available in package \pkg{mlbench} via

<<>>=

data("BostonHousing", package = "mlbench")

@

and provides $n = \Sexpr{NROW(BostonHousing)}$ observations of the median value of owner-occupied

homes in Boston (in USD~1000) along with $\Sexpr{NCOL(BostonHousing)}$ covariates including in particular

the number of rooms per dwelling (\code{rm}) and the percentage

of lower status of the population (\code{lstat}). A segment-wise linear relationship between

the value and these two variables is very intuitive, whereas the shape of the influence

of the remaining covariates is rather unclear and hence should be learned from the data.

Therefore, a linear regression model for median value explained by \verb:rm^2:

and \verb:log(lstat): with $k = 3$ regression coefficients is employed and

partitioned with respect to all $\ell = 11$ remaining variables. To facilitate subsequent

commands, the transformations are explicitely stored in \code{BostonHousing}:

<<>>=

BostonHousing$lstat <- log(BostonHousing$lstat)

BostonHousing$rm <- BostonHousing$rm^2

@

Choosing appropriate

transformations of the dependent variable and the regressors that enter the linear

regression model is important to obtain a well-fitting model in each segment and

we follow in our choice the recommendations of \cite{BreimanFriedman1985}. Monotonous

transformations of the partitioning variables do not affect the recursive partitioning

algorithm and hence do not have to be performed. However, it is important to distinguish

between numerical and categorical variables for choosing an appropriate parameter 

stability test. The variable \code{chas} is a dummy indicator variable (for tract bounds

with Charles river) and should thus be turned into a factor. Furthermore, the variable

\code{rad} is an index of accessibility to radial highways and takes only 9 distinct

values. Thus it is most appropriately treated as an ordered factor.

<<>>=

BostonHousing$chas <- factor(BostonHousing$chas, levels = 0:1, labels = c("no", "yes"))

BostonHousing$rad <- factor(BostonHousing$rad, ordered = TRUE)

@

Both transformations only affect the parameter stability test chosen (step~2), not the splitting

procedure (step~3).

The model is estimated

by OLS, the instability is assessed using a Bonferroni-corrected

significance level of $\alpha = 0.05$ and the nodes are split with a required minimal

segment size of $40$ observations. The control parameters are thus set to

<<>>=

ctrl <- mob_control(alpha = 0.05, bonferroni = TRUE, minsplit = 40,

  objfun = deviance, verbose = TRUE)

@

Actually, all of these settings are the defaults except \code{minsplit = 40} and

\code{verbose = TRUE} which causes some information about the fitting process

being written to the screen. The objective function \code{deviance()} extracts in

this case the residual sum of squares from a fitted \code{linearModel} object.

Having collected all building blocks, we can now call the function \code{mob()}

that takes the model specification of the linear regression model \verb:medv ~ lstat + rm:

plus all partitioning variables, along with the \code{data} set, the \code{control}

settings and the \code{model} to be used.

<<>>=

fmBH <- mob(medv ~ lstat + rm | zn + indus + chas + nox + age + dis + rad + tax + crim + b + ptratio,

  data = BostonHousing, control = ctrl, model = linearModel)

@

The result is the fitted model \code{fmBH} of class ``\code{mob}'' that contains

the tree with a fitted linear regression associated with every node. Printing

this object will show the splits, their $p$ values and call the \code{print()} method

for the model in each terminal node (i.e., this simply relies on a \code{print()}

method being available for the fitted model and re-uses it).

<<>>=

fmBH

@

Looking at the printed output is typically rather tedious, a visualization via

the \code{plot()} method

<<eval=FALSE>>=

plot(fmBH)

@

is much easier to interpret. By default, this produces partial scatter plots of the

variable $y$ against each of the regressors $x_i$ in the terminal nodes. Each scatter

plot also shows the fitted values, i.e., a project of the fitted hyperplane.

\setkeys{Gin}{width=\textwidth}

\begin{figure}[p]

\begin{center}

<<BostonHousing-plot,echo=FALSE,results=hide,fig=TRUE,height=8,width=12,include=FALSE>>=

plot(fmBH)

@

\includegraphics[width=18cm,keepaspectratio,angle=90]{MOB-BostonHousing-plot}

\caption{\label{fig:BostonHousing} Linear-regression-based tree for the Boston housing data.

The plots in the leaves give partial scatter plots for \code{rm} (upper panel) and 

\code{lstat} (lower panel).}

\end{center}

\end{figure}

From this visualization, it can be seen that in the nodes~4, 6, 7 and 8 the increase of

value with the number of rooms dominates the picture (upper panel) whereas in node~9 the

decrease with the lower status population percentage (lower panel) is more pronounced. 

Splits are performed in the variables \code{tax} (poperty-tax rate) and

\code{ptratio} (pupil-teacher ratio).

Various quantities of interest can be computed, provided that the \code{model} used

provides the corresponding methods, e.g., \code{predict()}, \code{residuals()}, \code{logLik()},

\code{coef()} and \code{summary()}. The latter two by default try to extract information

for the terminal nodes, but a \code{node} argument can be set to the node IDs of interest.

As an example, the regression coefficients for the terminal node models can be easily 

extracted by

<<>>=

coef(fmBH)

@

reflecting the differences of the models that can also be seen in the the associated

\code{plot()}. Even more information is available in a \code{summary()}, e.g., for node 7:

<<>>=

summary(fmBH, node = 7)

@

The test statistics and $p$~values computed in each node, can be extracted analogously

by using the method for the function \code{sctest()} (for performing \underline{s}tructural

\underline{c}hange \underline{test}s).

<<>>=

sctest(fmBH, node = 7)

@

For summarizing the quality of the fit, we could compute the mean squared error, log-likelihood

or AIC:

<<>>=

mean(residuals(fmBH)^2)

logLik(fmBH)

AIC(fmBH)

@

<<echo=FALSE>>=

nt <- NROW(coef(fmBH))

nk <- NCOL(coef(fmBH))

@

As the \code{logLik()} method simply re-uses the method for \code{linearModel} objects,

this does not only report $\Sexpr{(nk+1)*nt-1}$ estimated parameters ($\Sexpr{nk}$ parameters in

each of the $\Sexpr{nt}$ terminal nodes plus $\Sexpr{nt} - 1$ split points) but

$\Sexpr{(nk+2)*nt-1}$ parameters because each terminal node is additionally associated with a 

variance estimate. However, for the fitting process, the variance was treated as a nuisance 

parameter as we employed OLS estimation (rather than fully-specified ML estimation). 

\subsection{Pima Indians diabetes data}

Another popular benchmark data set for binary classifications is the Pima Indians

diabetes database which is also available from \pkg{mlbench}:

<<>>=

data("PimaIndiansDiabetes2", package = "mlbench")

PimaIndiansDiabetes <- na.omit(PimaIndiansDiabetes2[,-c(4, 5)])

@

After omitting missing values (and the variables \verb:triceps: and \verb:insulin: which

are missing for most women), the data set provides diabetes test results for

$n = \Sexpr{NROW(PimaIndiansDiabetes)}$ women

along with $\Sexpr{NCOL(PimaIndiansDiabetes)}$ covariates including in particular

the plasma glucose concentration \code{glucose} as an important predictor for diabetes.

Fitting a logistic regression model \verb:diabetes ~ glucose: seems to be straightforward,

whereas the influence of the remaining variables should again be learned by recursive

partitioning. This will yield a model tree with $k = 2$ regression coefficients in each terminal

node, partitioned with respect to the remaining $\ell = 5$ remaining variables.

The model is estimated by ML employing the \code{glinearModel}, 

the instability is assessed using a Bonferroni-corrected

significance level of $\alpha = 0.05$ and the nodes are split with a required minimal

segment size of $20$ observations. Hence, all control parameters correspond to the

default values in \verb:mob_control(): and do not have to be set explicitely in 

the \code{mob()} call:

<<>>=

fmPID <- mob(diabetes ~ glucose | pregnant + pressure + mass + pedigree + age,

  data = PimaIndiansDiabetes, model = glinearModel, family = binomial())

@

To visualize this, we simply call again:

<<eval=FALSE>>=

plot(fmPID)

@

\setkeys{Gin}{width=\textwidth}

\begin{figure}[bth]

\begin{center}

<<PimaIndiansDiabetes-plot,echo=FALSE,results=hide,fig=TRUE,height=6,width=9>>=

plot(fmPID)

@

\caption{\label{fig:PimaIndiansDiabetes} Logistic-regression-based tree for the Pima Indians

diabetes data. The plots in the leaves give spinograms for \code{diabetes} versus 

\code{glucose}.}

\end{center}

\end{figure}

which produces again a plot of the dependent variable $y$ against the only regressors

$x$ in the terminal nodes. As $y$ is \code{diabetes}, a binary variable, and $x$ is 

\code{glucose}, a numeric variable, a spinogram is chosen for visualization. The breaks

in the spinogram are the five-point summary of \code{glucose} on the full data set. The

fitted lines are the mean predicted probabilities in each group.

The model tree distinguishes three different groups:

\begin{itemize}

  \item[\#2] Women with low body mass index that have on average a low risk of

    diabetes, however this increases clearly with glucose level.

  \item[\#4] Women with average and high body mass index, younger than 30 years,

    that have a higher avarage risk that also increases with glucose level.

  \item[\#5] Women with average and high body mass index, older than 30 years,

    that have a high avarage risk that increases only slowly with glucose level.

\end{itemize}

The same interpretation can also be drawn from the coefficient estimates

and the corresponding odds ratios (with respect to glucose):

<<>>=

coef(fmPID)

exp(coef(fmPID)[,2])

@  

<<echo=FALSE>>=

risk <- round(100 * (exp(coef(fmPID)[,2])-1), digits = 1)

@

i.e., the odds increase by \Sexpr{risk[1]}\%, \Sexpr{risk[2]}\% and \Sexpr{risk[3]}\%

with respect to glucose in the three groups.

\section{Conclusion}

\label{sec:conclusion}

The function \code{mob()} in the \pkg{party} package provides a flexible and object-oriented

implementation of the general algorithm for model-based recursive partitioning.

Models of class ``\code{StatModel}''---that employ a formula interface and are equipped with

methods for the generic functions \code{reweight()}, \code{weights()}, \code{estfun()} plus

some function for extracting the value of the objective function---can be easily partitioned.

The resulting ``\code{mob}'' tree can be flexibly summarized, both numerically and graphically,

and used for predictions on new data.

\bibliography{partyrefs}

\end{document}