\documentclass[nojss]{jss}

%\VignetteIndexEntry{party: A Laboratory for Recursive Partytioning}

%\VignetteDepends{coin, TH.data}

%\VignetteKeywords{conditional inference, non-parametric models, recursive partitioning}

%\VignettePackage{party}

%% packages

\usepackage{amstext}

\usepackage{amsfonts}

\usepackage{amsmath}

\usepackage{thumbpdf}

\usepackage{rotating}

%% need no \usepackage{Sweave}

%% commands

\renewcommand{\Prob}{\mathbb{P} }

\renewcommand{\E}{\mathbb{E}}

\newcommand{\V}{\mathbb{V}}

\newcommand{\Var}{\mathbb{V}}

\newcommand{\R}{\mathbb{R} }

\newcommand{\N}{\mathbb{N} }

\newcommand{\C}{\mathbb{C} }

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

\newcommand{\argmax}{\operatorname{argmax}\displaylimits}

\newcommand{\LS}{\mathcal{L}_n}

\newcommand{\TS}{\mathcal{T}_n}

\newcommand{\LSc}{\mathcal{L}_{\text{comb},n}}

\newcommand{\LSbc}{\mathcal{L}^*_{\text{comb},n}}

\newcommand{\F}{\mathcal{F}}

\newcommand{\A}{\mathcal{A}}

\newcommand{\yn}{y_{\text{new}}}

\newcommand{\z}{\mathbf{z}}

\newcommand{\X}{\mathbf{X}}

\newcommand{\Y}{\mathbf{Y}}

\newcommand{\sX}{\mathcal{X}}

\newcommand{\sY}{\mathcal{Y}}

\newcommand{\T}{\mathbf{T}}

\newcommand{\x}{\mathbf{x}}

\renewcommand{\a}{\mathbf{a}}

\newcommand{\xn}{\mathbf{x}_{\text{new}}}

\newcommand{\y}{\mathbf{y}}

\newcommand{\w}{\mathbf{w}}

\newcommand{\ws}{\mathbf{w}_\cdot}

\renewcommand{\t}{\mathbf{t}}

\newcommand{\M}{\mathbf{M}}

\renewcommand{\vec}{\text{vec}}

\newcommand{\B}{\mathbf{B}}

\newcommand{\K}{\mathbf{K}}

\newcommand{\W}{\mathbf{W}}

\newcommand{\D}{\mathbf{D}}

\newcommand{\I}{\mathbf{I}}

\newcommand{\bS}{\mathbf{S}}

\newcommand{\cellx}{\pi_n[\x]}

\newcommand{\partn}{\pi_n(\mathcal{L}_n)}

\newcommand{\err}{\text{Err}}

\newcommand{\ea}{\widehat{\text{Err}}^{(a)}}

\newcommand{\ecv}{\widehat{\text{Err}}^{(cv1)}}

\newcommand{\ecvten}{\widehat{\text{Err}}^{(cv10)}}

\newcommand{\eone}{\widehat{\text{Err}}^{(1)}}

\newcommand{\eplus}{\widehat{\text{Err}}^{(.632+)}}

\newcommand{\eoob}{\widehat{\text{Err}}^{(oob)}}

\newcommand{\bft}{\mathbf{t}}

\hyphenation{Qua-dra-tic}

\title{\pkg{party}: A Laboratory for Recursive Partytioning}

\Plaintitle{party: A Laboratory for Recursive Partytioning}

\author{Torsten Hothorn\\Ludwig-Maximilians-\\Universit\"at M\"unchen

   \And Kurt Hornik\\Wirtschaftsuniversit\"at Wien

   \And Achim Zeileis\\Wirtschaftsuniversit\"at Wien}

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

\Abstract{

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

  at providing a recursive part(y)itioning laboratory assembling

  various high- and low-level tools for building tree-based regression

  and classification models. This includes conditional inference trees

  (\code{ctree}), conditional inference forests (\code{cforest}) and

  parametric model trees (\code{mob}). At the core of the package is

  \code{ctree}, an implementation of conditional inference trees which

  embed tree-structured regression models into a well defined theory of

  conditional inference procedures. This non-parametric class of regression

  trees is applicable to all kinds of regression problems, including

  nominal, ordinal, numeric, censored as well as multivariate response

  variables and arbitrary measurement scales of the covariates.

  This vignette comprises a practical guide to exploiting the flexible

  and extensible computational tools in \pkg{party} for fitting and

  visualizing conditional inference trees.

}

\Keywords{conditional inference, non-parametric models, recursive partitioning}

\Address{

  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, Achim Zeileis\\

  Institute for Statistics and Mathematics\\

  WU Wirtschaftsuniversit\"at Wien\\

  Augasse 2--6\\

  1090 Wien, Austria\\

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

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

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

}

\begin{document}

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

options(width = 70, SweaveHooks = list(leftpar = 

    function() par(mai = par("mai") * c(1, 1.1, 1, 1))))

require("party")

require("coin")

set.seed(290875)

@

\setkeys{Gin}{width=\textwidth}

\section{Introduction}

The majority of recursive partitioning algorithms 

are special cases of a simple two-stage algorithm: 

First partition the observations by univariate splits in a recursive way and

second fit a constant model in each cell of the resulting partition.

The most popular implementations of such algorithms are `CART' 

\citep{classifica:1984} and `C4.5' \citep{Quinlan1993}. Not unlike AID,   

both perform an exhaustive search over all possible splits maximizing an

information measure of node impurity selecting the 

covariate showing the best split.

This approach has two fundamental problems: overfitting and a selection 

bias towards covariates with many possible splits. 

With respect to the 

overfitting problem \cite{Mingers1987} notes that the algorithm

\begin{quote}

[\ldots] has no concept of statistical significance, and so cannot

distinguish between a significant and an insignificant improvement in the

information measure.

\end{quote}

With the \pkg{party} package \citep[see][for a full description of its

methodological foundations]{Hothorn+Hornik+Zeileis:2006a}

we enter at the point where \cite{WhiteLiu1994} demand for

\begin{quote}

[\ldots] a \textit{statistical} approach [to recursive partitioning] which

takes

into account the \textit{distributional} properties of the measures.

\end{quote} 

We present a unified framework embedding recursive binary partitioning into

the well defined theory of permutation tests developed by

\cite{StrasserWeber1999}.

The conditional distribution of statistics measuring the association between

responses and covariates is the basis for an unbiased selection among

covariates measured at different scales. 

Moreover, multiple test procedures are applied to determine whether no

significant

association between any of the covariates and the response can be stated and 

the recursion needs to stop.

\section{Recursive binary partitioning} \label{algo}

We focus on regression models describing the conditional distribution of a

response variable $\Y$ given the status of $m$ covariates by

means of tree-structured recursive partitioning. The response $\Y$ from some

sample space $\sY$ may be multivariate as well. 

The $m$-dimensional covariate vector $\X = (X_1, \dots, X_m)$ is taken 

from a sample space $\sX = \sX_1 \times \cdots \times \sX_m$.

Both response variable and covariates may be measured

at arbitrary scales.

We assume that the conditional distribution $D(\Y | \X)$ of the response 

$\Y$ given the covariates $\X$ depends on a function $f$ of the covariates

\begin{eqnarray*} 

D(\Y | \X) = D(\Y | X_1, \dots, X_m) = D(\Y | f( X_1, \dots,

X_m)),

\end{eqnarray*}

where we restrict ourselves to partition based regression relationships,

i.e., $r$ disjoint cells $B_1, \dots, B_r$ partitioning the covariate space $\sX

= \bigcup_{k = 1}^r B_k$.

A model of the regression relationship is to be fitted based on a learning 

sample $\LS$, i.e., a random sample of $n$ independent and

identically distributed observations, possibly with some covariates $X_{ji}$

missing,

\begin{eqnarray*}

\LS & = & \{ (\Y_i, X_{1i}, \dots, X_{mi}); i = 1, \dots, n \}.

\end{eqnarray*}

For the sake of simplicity, we use a learning sample

<<party-data, echo = TRUE>>=

ls <- data.frame(y = gl(3, 50, labels = c("A", "B", "C")), x1 = rnorm(150) + rep(c(1, 0, 0), c(50, 50,

50)), x2 = runif(150))

@

in the following illustrations.

A generic algorithm for recursive binary partitioning for a given learning

sample $\LS$ can be formulated using non-negative integer valued case weights $\w

= (w_1, \dots, w_n)$. Each node of a tree is represented by a vector of

case weights having non-zero elements when the corresponding observations

are elements of the node and are zero otherwise. The following algorithm

implements recursive binary partitioning:

\begin{enumerate}

\item For case weights $\w$ test the global null hypothesis of independence between

      any of the $m$ covariates and the response. Stop if this 

      hypothesis cannot be rejected. 

      Otherwise select the covariate $X_{j^*}$ with strongest 

      association to $\Y$.

\item Choose a set $A^* \subset \sX_{j^*}$ in order to split $\sX_{j^*}$ into

      two disjoint sets $A^*$ and $\sX_{j^*} \setminus A^*$. 

      The case weights $\w_\text{left}$ and $\w_\text{right}$ determine the

      two subgroups with $w_{\text{left},i} = w_i I(X_{j^*i} \in A^*)$ and 

      $w_{\text{right},i} = w_i I(X_{j^*i} \not\in A^*)$ for all $i = 1,

      \dots, n$ ($I(\cdot)$ denotes the indicator function).

\item Recursively repeat steps 1 and 2 with modified case weights 

      $\w_\text{left}$ and $\w_\text{right}$, respectively.

\end{enumerate}

The separation of variable

selection and splitting procedure into steps 1 and 2 of the algorithm

is the key for the construction of interpretable tree

structures not suffering a systematic tendency towards covariates with many

possible splits or many missing values. In addition, a statistically

motivated and intuitive stopping criterion can be implemented: We stop 

when the global null hypothesis of independence between the

response and any of the $m$ covariates cannot be rejected at a pre-specified

nominal level~$\alpha$. The algorithm induces a partition $\{B_1, \dots, B_r\}$ of

the covariate space $\sX$, where each cell $B \in \{B_1, \dots, B_r\}$ 

is associated with a vector of case weights. 

In package \pkg{party}, the dependency structure and the variables may be

specified in a traditional formula based way

<<party-formula, echo = TRUE, results = hide>>=

library("party")

ctree(y ~ x1 + x2, data = ls)

@

Case counts $\w$ may be specified using the \texttt{weights} argument.

\section{Recursive partitioning by conditional inference} \label{framework}

In the main part of this section we focus on step 1 of the generic algorithm.

Unified tests for independence are constructed by means of the conditional

distribution of linear statistics in the permutation test framework

developed by \cite{StrasserWeber1999}. The determination of the best binary split

in one selected covariate and the handling of missing values

is performed based on standardized linear statistics within the same

framework as well. 

\subsection{Variable selection and stopping criteria}

At step 1 of the generic algorithm given in Section~\ref{algo} we face an 

independence problem. We need to decide whether there is any information

about the response variable covered by any of the $m$  covariates. In each node

identified by case weights $\w$, the

global hypothesis of independence is formulated in terms of the $m$ partial hypotheses

$H_0^j: D(\Y | X_j) = D(\Y)$ with global null hypothesis $H_0 = \bigcap_{j = 1}^m

H_0^j$.

When we are not able to reject $H_0$ at a pre-specified 

level $\alpha$, we stop the recursion.

If the global hypothesis can be rejected, we measure the association

between $\Y$ and each of the covariates $X_j, j = 1, \dots, m$, by

test statistics or $P$-values indicating the deviation from the partial

hypotheses $H_0^j$.  

For notational convenience and without loss of generality we assume that the

case weights $w_i$ are either zero or one. The symmetric group of all

permutations of  the elements of $(1, \dots, n)$ with corresponding case

weights $w_i = 1$ is denoted by $S(\LS, \w)$. A more general notation is

given in the Appendix. We measure the association between $\Y$ and $X_j, j = 1, \dots, m$, 

by linear statistics of the form

\begin{eqnarray} \label{linstat}

\T_j(\LS, \w) = \vec \left( \sum_{i=1}^n w_i g_j(X_{ji})

h(\Y_i, (\Y_1, \dots, \Y_n))^\top \right) \in \R^{p_jq}

\end{eqnarray}

where $g_j: \sX_j \rightarrow \R^{p_j}$ is a non-random transformation of

the covariate $X_j$. The transformation may be specified using the

\texttt{xtrafo} argument. If, for example, a ranking \textit{both}

\texttt{x1} and \texttt{x2} is required,

<<party-xtrafo, echo = TRUE, results = hide>>=

ctree(y ~ x1 + x2, data = ls, xtrafo = function(data) trafo(data,

numeric_trafo = rank))

@

can be used. The \emph{influence function} 

$h: \sY \times \sY^n \rightarrow

\R^q$ depends on the responses $(\Y_1, \dots, \Y_n)$ in a permutation

symmetric way. 

%%, i.e., $h(\Y_i, (\Y_1, \dots, \Y_n)) = h(\Y_i, (\Y_{\sigma(1)}, \dots,

%%\Y_{\sigma(n)}))$ for all permutations $\sigma \in S(\LS, \w)$. 

Section~\ref{examples} explains how to choose $g_j$ and $h$ in different 

practical settings. A $p_j \times q$ matrix is converted into a 

$p_jq$ column vector by column-wise combination using the `vec' operator. 

The influence function can be specified using the \texttt{ytrafo} argument.

The distribution of $\T_j(\LS, \w)$ under $H_0^j$ depends on the joint

distribution of $\Y$ and $X_j$, which is unknown under almost all practical

circumstances. At least under the null hypothesis one can dispose of this

dependency by fixing the covariates and conditioning on all possible 

permutations of the responses. This principle leads to test procedures known

as \textit{permutation tests}. 

The conditional expectation $\mu_j \in \R^{p_jq}$ and covariance $\Sigma_j

\in \R^{p_jq \times p_jq}$ of $\T_j(\LS, \w)$ under $H_0$ given

all permutations $\sigma \in S(\LS, \w)$ of the responses are derived by

\cite{StrasserWeber1999}: 

\begin{eqnarray}

\mu_j & = & \E(\T_j(\LS, \w) | S(\LS, \w)) = 

\vec \left( \left( \sum_{i = 1}^n w_i g_j(X_{ji}) \right) \E(h | S(\LS,

\w))^\top

\right), \nonumber \\

\Sigma_j & = & \V(\T_j(\LS, \w) | S(\LS, \w)) \nonumber \\

& = & 

    \frac{\ws}{\ws - 1}  \V(h | S(\LS, \w)) \otimes

        \left(\sum_i w_i  g_j(X_{ji}) \otimes w_i  g_j(X_{ji})^\top \right) \label{expectcovar}

\\

& - & \frac{1}{\ws - 1}  \V(h | S(\LS, \w))  \otimes \left(

        \sum_i w_i g_j(X_{ji}) \right)

\otimes \left( \sum_i w_i g_j(X_{ji})\right)^\top 

\nonumber

\end{eqnarray}

where $\ws = \sum_{i = 1}^n w_i$ denotes the sum of the case weights,

$\otimes$ is the Kronecker product and the conditional expectation of the

influence function is 

\begin{eqnarray*}

\E(h | S(\LS, \w)) = \ws^{-1} \sum_i w_i h(\Y_i, (\Y_1, \dots, \Y_n)) \in

\R^q

\end{eqnarray*} 

with corresponding $q \times q$ covariance matrix

\begin{eqnarray*}

\V(h | S(\LS, \w)) = \ws^{-1} \sum_i w_i \left(h(\Y_i, (\Y_1, \dots, \Y_n)) - \E(h | S(\LS, \w))

\right) \\

\left(h(\Y_i, (\Y_1, \dots, \Y_n)) - \E(h | S(\LS, \w))\right)^\top.

\end{eqnarray*}

Having the conditional expectation and covariance at hand we are able to

standardize a linear statistic $\T \in \R^{pq}$ of the form

(\ref{linstat}) for some $p \in \{p_1, \dots, p_m\}$. 

Univariate test statistics~$c$ mapping an observed multivariate 

linear statistic $\bft \in

\R^{pq}$ into the real line can be of arbitrary form.  An obvious choice is

the maximum of the absolute values of the standardized linear statistic

\begin{eqnarray*}

c_\text{max}(\bft, \mu, \Sigma)  = \max_{k = 1, \dots, pq} \left| \frac{(\bft -

\mu)_k}{\sqrt{(\Sigma)_{kk}}} \right|

\end{eqnarray*}

utilizing the conditional expectation $\mu$ and covariance matrix

$\Sigma$. The application of a quadratic form $c_\text{quad}(\bft, \mu, \Sigma)  = 

(\bft - \mu) \Sigma^+ (\bft - \mu)^\top$ is one alternative, although

computationally more expensive because the Moore-Penrose 

inverse $\Sigma^+$ of $\Sigma$ is involved.

The type of test statistic to be used can be specified by means of the

\texttt{ctree\_control} function, for example

\begin{Schunk}

\begin{Sinput}

R> ctree(y ~ x1 + x2, data = ls, 

      control = ctree_control(teststat = "max"))

\end{Sinput}

\end{Schunk}

uses $c_\text{max}$ and 

\begin{Schunk}

\begin{Sinput}

R> ctree(y ~ x1 + x2, data = ls, 

      control = ctree_control(teststat = "quad"))

\end{Sinput}

\end{Schunk}

takes $c_\text{quad}$ (the default).

It is important to note that the test statistics $c(\bft_j, \mu_j, \Sigma_j),

j = 1, \dots, m$, cannot be directly compared in an unbiased  

way unless all of the covariates are measured at the same scale, i.e.,

$p_1 = p_j, j = 2, \dots, m$. 

In order to allow for an unbiased variable selection we need to switch to

the $P$-value scale because $P$-values for the conditional distribution of test

statistics $c(\T_j(\LS, \w), \mu_j, \Sigma_j)$ 

can be directly compared among covariates

measured at different scales. In step 1 of the generic algorithm we 

select the covariate with minimum $P$-value, i.e., the covariate $X_{j^*}$

with $j^*  =  \argmin_{j = 1, \dots, m} P_j$, where

\begin{displaymath}

P_j  =  \Prob_{H_0^j}(c(\T_j(\LS, \w), \mu_j,

            \Sigma_j) \ge c(\bft_j, \mu_j, \Sigma_j) | S(\LS, \w))

\end{displaymath}

denotes the $P$-value of the conditional test for $H_0^j$. 

So far, we have only addressed testing each partial hypothesis $H_0^j$, which

is sufficient for an unbiased variable selection. A global test 

for $H_0$ required in step 1 can be constructed via an aggregation of the 

transformations $g_j, j = 1, \dots, m$,

i.e., using a linear statistic of the form

\begin{eqnarray*}

\T(\LS, \w) = \vec \left( \sum_{i=1}^n w_i

\left(g_1(X_{1i})^\top, \dots,

g_m(X_{mi})^\top\right)^\top h(\Y_i, (\Y_1, \dots, \Y_n))^\top \right).

%%\in \R^{\sum_j p_jq}

\end{eqnarray*}

However, this approach is less attractive for learning samples with

missing values. Universally applicable approaches are multiple test 

procedures based on $P_1, \dots, P_m$. 

Simple Bonferroni-adjusted $P$-values ($mP_j$ prior to version 0.8-4, now $1

- (1 - P_j)^m$ is used), available via 

<<party-Bonf, echo = TRUE, results = hide>>=

ctree_control(testtype = "Bonferroni")

@

or a min-$P$-value resampling approach 

<<party-MC, echo = TRUE, results = hide>>=

ctree_control(testtype = "MonteCarlo")           

@

are just examples and we refer

to the multiple testing literature \citep[e.g.,][]{WestfallYoung1993}

for more advanced methods. We reject $H_0$ when the minimum 

of the adjusted $P$-values is less than a pre-specified nominal level $\alpha$

and otherwise stop the algorithm. In this sense, $\alpha$

may be seen as a unique parameter determining the size of the resulting trees.

\subsection{Splitting criteria}

Once we have selected a covariate in step 1 of the algorithm, the split itself can be

established by any split criterion, including those established by

\cite{classifica:1984} or \cite{Shih1999}. Instead of simple binary splits, 

multiway splits can be implemented as well, for example utilizing

the work of \cite{OBrien2004}. However, most splitting criteria are not

applicable to response variables measured at arbitrary scales and we

therefore utilize the permutation test framework described above  

to find the optimal binary split in one selected 

covariate $X_{j^*}$ in step~2 of the generic algorithm. The goodness of a split is

evaluated by two-sample linear statistics which are special cases of the linear 

statistic (\ref{linstat}). 

For all possible subsets $A$ of the sample space $\sX_{j^*}$ the linear

statistic

\begin{eqnarray*}

\T^A_{j^*}(\LS, \w) = \vec \left( \sum_{i=1}^n w_i I(X_{j^*i} \in A) 

                                  h(\Y_i, (\Y_1, \dots, \Y_n))^\top \right) \in \R^{q}

\end{eqnarray*}

induces a two-sample statistic measuring the discrepancy between the samples 

$\{ \Y_i | w_i > 0 \text{ and } X_{ji} \in A; i = 1, \dots, n\}$ 

and $\{ \Y_i | w_i > 0 \text{ and } X_{ji} \not\in A; i = 1, \dots, n\}$. 

The conditional expectation $\mu_{j^*}^A$ and covariance

$\Sigma_{j^*}^A$

can be computed by (\ref{expectcovar}).

The split $A^*$ with a test statistic maximized over all possible

subsets $A$ is established:

\begin{eqnarray} \label{split}

A^* = \argmax_A c(\bft_{j^*}^A, \mu_{j^*}^A, \Sigma_{j^*}^A).

\end{eqnarray}

The statistics $c(\bft_{j^*}^A, \mu_{j^*}^A, \Sigma_{j^*}^A)$ are available

for each node with

<<party-ss, echo = TRUE, results = hide>>=

ctree_control(savesplitstats = TRUE)

@

and can be used to depict a scatter plot of the covariate

$\sX_{j^*}$ against the statistics.

Note that we do not need to compute the distribution of 

$c(\bft_{j^*}^A, \mu_{j^*}^A, \Sigma_{j^*}^A)$ in step~2. 

In order to anticipate pathological splits one can restrict the number of

possible subsets that are evaluated, for example by introducing restrictions

on the sample size or the sum of the case weights in each of the two groups

of observations induced by a possible split. For example,

<<party-minsplit, echo = TRUE, results = hide>>=

ctree_control(minsplit = 20)

@

requires the sum of the weights in both the left and right daughter 

node to exceed the value of $20$.

\subsection{Missing values and surrogate splits}

If an observation $X_{ji}$ in covariate $X_j$ is missing, we set the

corresponding case weight $w_i$ to zero for the computation of $\T_j(\LS, \w)$

and, if we would like to split in $X_j$, in $\T_{j}^A(\LS, \w)$ as well.

Once a split $A^*$ in $X_j$ has been implemented, surrogate splits can be 

established by

searching for a split leading to roughly the same division of the

observations as the original split. One simply replaces the original

response variable by a binary variable $I(X_{ji} \in A^*)$ coding the

split and proceeds as described in the previous part. The number of

surrogate splits can be controlled using

<<party-maxsurr, echo = TRUE, results = hide>>=

ctree_control(maxsurrogate = 3)

@

\subsection{Inspecting a tree}

Once we have fitted a conditional tree via

<<party-fitted, echo = TRUE>>=

ct <- ctree(y ~ x1 + x2, data = ls)

@

we can inspect the results via a \texttt{print} method

<<party-print, echo = TRUE>>=

ct

@

or by looking at a graphical representation as in Figure~\ref{party-plot1}.

\begin{figure}[h]

\begin{center}

<<party-plot, echo = TRUE, fig = TRUE, height = 5, width = 8>>=

plot(ct)

@

\caption{A graphical representation of a classification tree.

\label{party-plot1}}

\end{center}

\end{figure}

Each node can be extracted by its node number, i.e., the root node is

<<party-nodes, echo = TRUE>>=

nodes(ct, 1)

@

This object is a conventional list with elements

<<party-nodelist, echo = TRUE>>=

names(nodes(ct, 1)[[1]])

@

and we refer to the manual pages for a description of those elements.

The \texttt{Predict} function aims at computing predictions in the space of

the response variable, in our case a factor

<<party-Predict, echo = TRUE>>=

Predict(ct, newdata = ls)

@

When we are interested in properties of the conditional distribution of the

response given the covariates, we use

<<party-treeresponse, echo = TRUE>>=

treeresponse(ct, newdata = ls[c(1, 51, 101),])

@

which, in our case, is a list with conditional class probabilities.

We can determine the node numbers of nodes some new observations are falling

into by

<<party-where, echo = TRUE>>=

where(ct, newdata = ls[c(1,51,101),])

@

\section{Examples} \label{examples}

\subsection{Univariate continuous or discrete regression}

For a univariate numeric response $\Y \in \R$, the most natural influence function

is the identity $h(\Y_i, (\Y_1, \dots, \Y_n)) = \Y_i$. 

In case some observations with extremely large or small values have been

observed, a

ranking of the observations may be appropriate:

$h(\Y_i, (\Y_1, \dots, \Y_n)) = \sum_{k=1}^n w_k I(\Y_k \le \Y_i)$ for $i = 1, \dots,

n$.

Numeric covariates can be handled by the identity transformation 

$g_{ji}(x) = x$ (ranks are possible, too). Nominal covariates at levels $1,

\dots, K$ are

represented by $g_{ji}(k) = e_K(k)$, the unit vector of length $K$ with

$k$th element being equal to one. Due to this flexibility, special test 

procedures like the Spearman test, the 

Wilcoxon-Mann-Whitney test or the Kruskal-Wallis test and permutation tests

based on ANOVA statistics or correlation coefficients

are covered by this framework. Splits obtained from (\ref{split}) maximize the

absolute value of the standardized difference between two means of the

values of the influence functions. 

For prediction, one is usually interested in an estimate of 

the expectation of

the response $\E(\Y | \X = \x)$ in each cell, an estimate can be 

obtained by 

\begin{eqnarray*}

\hat{\E}(\Y | \X = \x) = \left(\sum_{i=1}^n w_i(\x)\right)^{-1} \sum_{i=1}^n

w_i(\x) \Y_i.

\end{eqnarray*}

\subsection{Censored regression}

The influence function $h$ may be chosen as 

Logrank or Savage scores taking censoring into

account and one can proceed as for univariate continuous regression. This is

essentially the approach first published by \cite{regression:1988}. 

An alternative is the weighting scheme suggested by

\cite{MolinaroDudiotvdLaan2003}. A weighted Kaplan-Meier curve for the case

weights $\w(\x)$ can serve as prediction. 

\subsection{$J$-class classification}

The nominal response variable at levels $1, \dots, J$ 

is handled by influence

functions\linebreak $h(\Y_i, (\Y_1, \dots, \Y_n)) = e_J(\Y_i)$. Note that for a

nominal covariate $X_j$ at levels $1, \dots, K$ with  

$g_{ji}(k) = e_K(k)$ the

corresponding linear statistic $\T_j$ is a vectorized contingency table.

The conditional class probabilities can be estimated via 

\begin{eqnarray*}

\hat{\Prob}(\Y = y | \X = \x) = \left(\sum_{i=1}^n w_i(\x)\right)^{-1}

\sum_{i=1}^n

w_i(\x) I(\Y_i = y), \quad y = 1, \dots, J.

\end{eqnarray*}

\subsection{Ordinal regression}

Ordinal response variables measured at $J$ levels, and ordinal covariates

measured at $K$ levels, are associated with score vectors $\xi \in

\R^J$ and $\gamma \in \R^K$, respectively. Those scores reflect the

`distances' between the levels: If the variable is derived from an

underlying continuous variable, the scores can be chosen as the midpoints

of the intervals defining the levels. The linear statistic is now a linear

combination of the linear statistic $\T_j$ of the form

\begin{eqnarray*}

\M \T_j(\LS, \w) & = & \vec \left( \sum_{i=1}^n w_i \gamma^\top g_j(X_{ji})

            \left(\xi^\top h(\Y_i, (\Y_1, \dots, \Y_n)\right)^\top \right)

\end{eqnarray*}

with $g_j(x) = e_K(x)$ and $h(\Y_i, (\Y_1, \dots, \Y_n)) = e_J(\Y_i)$.

If both response and covariate are ordinal, the matrix of coefficients

is given by the Kronecker product of both score vectors $\M = \xi \otimes \gamma \in

\R^{1, KJ}$. In case the response is ordinal only, the matrix of 

coefficients $\M$ is a block matrix 

\begin{eqnarray*}

\M = 

\left( 

    \begin{array}{ccc}

        \xi_1 &          & 0     \\

              &  \ddots  &       \\

          0   &          & \xi_1 

    \end{array} \right| 

%%\left.

%%    \begin{array}{ccc}

%%       \xi_2 &          & 0     \\

%%              &  \ddots  &       \\

%%          0   &          & \xi_2 

%%    \end{array} \right| 

%%\left.

    \begin{array}{c}

         \\

      \hdots \\

          \\

    \end{array} %%\right.

\left|

    \begin{array}{ccc}

        \xi_q &          & 0     \\

              &  \ddots  &       \\

          0   &          & \xi_q 

    \end{array} 

\right) 

%%\in \R^{K, KJ}

%%\end{eqnarray*}

\text{ or } %%and if one covariate is ordered 

%%\begin{eqnarray*}

%%\M = \left( 

%%    \begin{array}{cccc}

%%        \gamma & 0      &        & 0 \\

%%        0      & \gamma &        & \vdots \\

%%       0      & 0      & \hdots & 0 \\

%%        \vdots & \vdots &        & 0 \\

%%        0      & 0      &        & \gamma

%%    \end{array} 

%%    \right) 

\M = \text{diag}(\gamma)

%%\in \R^{J, KJ}

\end{eqnarray*}

when one covariate is ordered but the response is not. For both $\Y$ and $X_j$

being ordinal, the corresponding test is known as linear-by-linear association

test \citep{Agresti2002}. Scores can be supplied to \texttt{ctree} using the

\texttt{scores} argument, see Section~\ref{illustrations} for an example.

\subsection{Multivariate regression}

For multivariate responses, the influence function is a combination of

influence functions appropriate for any of the univariate response variables

discussed in the previous paragraphs, e.g., indicators for multiple binary

responses \citep{Zhang1998,NohSongPark2004}, Logrank or Savage scores

for multiple failure times 

%%\citep[for example tooth loss times, ][]{SuFan2004}

and the original observations or a rank transformation for multivariate regression 

\citep{Death2002}.

\section{Illustrations and applications} \label{illustrations}

In this section, we present regression problems which illustrate the

potential fields of application of the methodology.  

Conditional inference trees based on $c_\text{quad}$-type test statistics 

using the identity influence function for numeric responses 

and asymptotic $\chi^2$ distribution are applied.

For the stopping criterion a simple

Bonferroni correction is used and we follow the usual convention by choosing

the nominal level of the conditional independence tests as $\alpha = 0.05$.

\subsection{Tree pipit abundance}

<<treepipit-ctree, echo = TRUE>>=

data("treepipit", package = "coin")

tptree <- ctree(counts ~ ., data = treepipit)

@

\begin{figure}[t]

\begin{center}

<<treepipit-plot, echo = TRUE, fig = TRUE, height = 5, width = 8>>=

plot(tptree, terminal_panel = node_hist(tptree, breaks = 0:6-0.5, ymax = 65, 

horizontal = FALSE, freq = TRUE))

@

\caption{Conditional regression tree for the tree pipit data.}

\end{center}

\end{figure}

<<treepipit-x, echo = FALSE, results = hide>>=

x <- tptree@tree

@

The impact of certain environmental factors on the population density of the 

tree pipit \textit{Anthus trivialis} 

%%in Frankonian oak forests 

is investigated by \cite{MuellerHothorn2004}.

The occurrence of tree pipits was recorded several times at

$n = 86$ stands which were established on a long environmental gradient. 

Among nine environmental factors, 

the covariate showing the largest association to the number of tree pipits

is the canopy overstorey $(P = $\Sexpr{round(1 - x$criterion[[3]], 3)}$)$. 

Two groups of stands can be distinguished: Sunny stands with less than $40\%$

canopy overstorey $(n = \Sexpr{sum(x$left$weights)})$ show a

significantly higher density of tree pipits compared to darker stands with more than

$40\%$ canopy overstorey $(n = \Sexpr{sum(x$right$weights)})$.

This result is important for management decisions

in forestry enterprises: Cutting the overstorey with release of old

oaks creates a perfect habitat for this indicator species

of near natural forest environments. 

\subsection{Glaucoma and laser scanning images}

<<glaucoma-ctree, echo = TRUE>>=

data("GlaucomaM", package = "TH.data")

gtree <- ctree(Class ~ ., data = GlaucomaM)

@

<<glaucoma-x, echo = FALSE, results = hide>>=

x <- gtree@tree

@

Laser scanning images taken from the eye background are expected to serve as 

the basis of an automated system for glaucoma diagnosis. 

Although prediction is more important in this application

\citep{new-glauco:2003}, a simple

visualization of the regression relationship is useful for comparing the

structures inherent in the learning sample with subject matter knowledge. 

For $98$ patients and $98$ controls, matched by age

and gender, $62$ covariates describing the eye morphology are available. The

data is part of the 

\pkg{TH.data} package, \url{http://CRAN.R-project.org}).

The first split in Figure~\ref{glaucoma} separates eyes with a volume 

above reference less than

$\Sexpr{x$psplit$splitpoint} \text{ mm}^3$ in the inferior part of the

optic nerve head (\texttt{vari}). 

Observations with larger volume are mostly controls, a finding which

corresponds to subject matter knowledge: The volume above reference measures the

thickness of the nerve layer, expected to decrease with a glaucomatous

damage of the optic nerve. Further separation is achieved by the volume

above surface global (\texttt{vasg}) and the volume above reference in the

temporal part of the optic nerve head (\texttt{vart}).

\setkeys{Gin}{width=.9\textwidth}

\begin{figure}[p]

\begin{center}

<<glaucoma-plot, echo = FALSE, fig = TRUE, height = 6, width = 10>>=

plot(gtree)

@

\caption{Conditional inference tree for the glaucoma data. For each inner node, the

Bonferroni-adjusted $P$-values are given, the fraction of glaucomatous eyes

is displayed for each terminal node. (Note: node 7 was not splitted prior to

version 0.8-4 because of using another formulation of the

Bonferroni-adjustment.) \label{glaucoma}}

<<glaucoma-plot-inner, echo = FALSE, fig = TRUE, height = 7, width = 10>>=

plot(gtree, inner_panel = node_barplot, 

     edge_panel = function(...) invisible(), tnex = 1)

@

\caption{Conditional inference tree for the glaucoma data with the 

         fraction of glaucomatous eyes displayed for both inner and terminal

         nodes. \label{glaucoma-inner}}

\end{center}

\end{figure}

The plot in Figure~\ref{glaucoma} is generated by

<<glaucoma-plot2, echo = TRUE, eval = FALSE>>=

plot(gtree)

@

\setkeys{Gin}{width=\textwidth}

and shows the distribution of the classes in 

the terminal nodes. This distribution can be shown for the inner nodes as

well, namely by specifying the appropriate panel generating function

(\texttt{node\_barplot} in our case), see Figure~\ref{glaucoma-inner}.

<<glaucoma-plot-inner, echo = TRUE, eval = FALSE>>=

plot(gtree, inner_panel = node_barplot, 

     edge_panel = function(...) invisible(), tnex = 1)

@

As mentioned in Section~\ref{framework}, it might be interesting to have a

look at the split statistics the split point estimate was derived from.

Those statistics can be extracted from the \texttt{splitstatistic} element

of a split and one can easily produce scatterplots against the selected

covariate. For all three inner nodes of \texttt{gtree}, we produce such a

plot in Figure~\ref{glaucoma-split}. For the root node, the estimated split point

seems very natural, since the process of split statistics seems to have a

clear maximum indicating that the simple split point model is something

reasonable to assume here. This is less obvious for nodes $2$ and,

especially, $3$.

\begin{figure}[t]

\begin{center}

<<glaucoma-split-plot, echo = TRUE, fig = TRUE, height = 4, width = 11, leftpar = TRUE>>=

cex <- 1.6

inner <- nodes(gtree, c(1, 2, 5))

layout(matrix(1:length(inner), ncol = length(inner)))

out <- sapply(inner, function(i) {

    splitstat <- i$psplit$splitstatistic

    x <- GlaucomaM[[i$psplit$variableName]][splitstat > 0]

    plot(x, splitstat[splitstat > 0], main = paste("Node", i$nodeID),

         xlab = i$psplit$variableName, ylab = "Statistic", ylim = c(0, 10),

         cex.axis = cex, cex.lab = cex, cex.main = cex)

    abline(v = i$psplit$splitpoint, lty = 3)

})

@

\caption{Split point estimation in each inner node. The process of 

         the standardized two-sample test statistics for each possible 

         split point in the selected input variable is show.

         The estimated split point is given as vertical dotted line.

         \label{glaucoma-split}}

\end{center}

\end{figure}

The class predictions of the tree for the learning sample (and for new

observations as well) can be computed using the \texttt{Predict} function. A

comparison with the true class memberships is done by

<<glaucoma-prediction, echo = TRUE>>=

table(Predict(gtree), GlaucomaM$Class)

@

When we are interested in conditional class probabilities, the

\texttt{treeresponse} method must be used. A graphical representation is

shown in Figure~\ref{glaucoma-probplot}.

\setkeys{Gin}{width=.5\textwidth}

\begin{figure}[ht!]

\begin{center}

<<glaucoma-classprob, echo = TRUE, fig = TRUE, height = 4, width = 5>>=

prob <- sapply(treeresponse(gtree), function(x) x[1]) +

               runif(nrow(GlaucomaM), min = -0.01, max = 0.01)

splitvar <- nodes(gtree, 1)[[1]]$psplit$variableName

plot(GlaucomaM[[splitvar]], prob, 

     pch = as.numeric(GlaucomaM$Class), ylab = "Conditional Class Prob.",

     xlab = splitvar)

abline(v = nodes(gtree, 1)[[1]]$psplit$splitpoint, lty = 2)

legend(0.15, 0.7, pch = 1:2, legend = levels(GlaucomaM$Class), bty = "n")

@

\caption{Estimated conditional class probabilities (slightly jittered) 

         for the Glaucoma data depending on the first split variable. 

         The vertical line denotes the first split point. \label{glaucoma-probplot}}

\end{center}

\end{figure}

\subsection{Node positive breast cancer}

<<GBSGS-ctree, echo = TRUE>>=

data("GBSG2", package = "TH.data")  

stree <- ctree(Surv(time, cens) ~ ., data = GBSG2)

@

\setkeys{Gin}{width=\textwidth}

\begin{figure}[t]

\begin{center}

<<GBSG2-plot, echo = TRUE, fig = TRUE, height = 6, width = 10>>=

plot(stree)

@

\caption{Tree-structured survival model for the GBSG2 data and the distribution of 

survival times in the terminal nodes. The median survival time is displayed in

each terminal node of the tree. \label{gbsg2}}

\end{center}

\end{figure}  

Recursive partitioning for censored responses has attracted a lot of interest

\citep[e.g.,][]{regression:1988,relative-r:1992}. 

Survival trees using $P$-value adjusted Logrank statistics

are used by \cite{statistics:2001a} 

for the evaluation of prognostic factors for the German Breast

Cancer Study Group (GBSG2) data,  a prospective controlled clinical trial on the

treatment of node positive breast cancer patients. Here, we use

Logrank scores as well. Complete data of seven prognostic factors of $686$

women are used for prognostic modeling, the

dataset is available within the 

\pkg{TH.data} package.

The number of positive lymph nodes

(\texttt{pnodes}) and the progesterone receptor (\texttt{progrec}) have

been identified as prognostic factors in the survival tree analysis

by \cite{statistics:2001a}.  Here, the binary variable coding whether a

hormonal therapy was applied or not (\texttt{horTh}) additionally is

part of the model depicted in Figure~\ref{gbsg2}.

The estimated median survival time for new patients is less informative

compared to the whole Kaplan-Meier curve estimated from the patients in the

learning sample for each terminal node. We can compute those `predictions'

by means of the \texttt{treeresponse} method

<<GBSG2-KM, echo = TRUE>>=

treeresponse(stree, newdata = GBSG2[1:2,])

@

\subsection{Mammography experience}

<<mammo-ctree, echo = TRUE>>=

data("mammoexp", package = "TH.data")

mtree <- ctree(ME ~ ., data = mammoexp)

@

\begin{figure}[t]

\begin{center}

<<mammo-plot, echo = TRUE, fig = TRUE, height = 5, width = 8>>=

plot(mtree)

@

\caption{Ordinal regression for the mammography experience data with the

fractions of (never, within a year, over one year) given in the nodes.

No admissible split was found for node 4 because only $5$ of $91$ women reported 

a family history of breast cancer and the sample size restrictions would 

require more than $5$ observations in each daughter node. \label{mammoexp}}

\end{center}

\end{figure}

Ordinal response variables are common in investigations where the response

is a subjective human interpretation. 

We use an example given by \cite{HosmerLemeshow2000}, p. 264, 

studying the relationship between the mammography experience (never,

within a year, over one year) and opinions about mammography expressed in

questionnaires answered by $n = 412$ women. The resulting partition based on

scores $\xi = (1,2,3)$ is given in Figure~\ref{mammoexp}. 

Women who (strongly) agree with the question `You do not need a mammogram unless

you develop symptoms' seldomly have experienced a mammography. The variable

\texttt{benefit} is a score with low values indicating a strong agreement with the

benefits of the examination. For those women in (strong) disagreement with the first

question above, low values of \texttt{benefit} identify persons being more likely to have 

experienced such an examination at all. 

\bibliography{partyrefs}

\end{document}

\subsection{Hunting spiders}

Finally, we take a closer look at a challenging dataset on animal abundance

first reported by \cite{AartSmeenk1975} and 

re-analyzed by \cite{Death2002} using regression trees dealing with

multivariate responses. The

abundance of $12$ hunting spider species is regressed on six

environmental variables (\texttt{water}, \texttt{sand}, \texttt{moss},

\texttt{reft}, \texttt{twigs} and \texttt{herbs}) for $n = 28$

observations. Because of the small sample size we allow for a split if at

least $5$ observations are element of a node and approximate the

distribution of a $c_\text{max}$ test statistic by $9999$ Monte-Carlo

replications. The prognostic factors \texttt{twigs} and \texttt{water}

found by \cite{Death2002} are confirmed by the model 

shown in Figure~\ref{spider} which additionally identifies \texttt{reft}.

The data are available in package 

\pkg{mvpart} \citep{mvpart2004} at \url{http://CRAN.R-project.org}.

<<spider-ctree, echo = TRUE, eval = FALSE>>=

data("spider", package = "mvpart")   

sptree <- ctree(arct.lute + pard.lugu + zora.spin + pard.nigr + pard.pull +

           aulo.albi + troc.terr + alop.cune + pard.mont + alop.acce + 

           alop.fabr + arct.peri ~ herbs+reft+moss+sand+twigs+water, 

           controls = ctree_control(teststat = "max", testtype = "MonteCarlo", nresample = 9999, 

           minsplit = 5, mincriterion = 0.9), data = spider)

sptree@tree$criterion

library("coin")

data("spider", package = "mvpart")

it <- independence_test(arct.lute + pard.lugu + zora.spin + pard.nigr + pard.pull +

                       aulo.albi + troc.terr + alop.cune + pard.mont + alop.acce +   

                       alop.fabr + arct.peri ~ herbs+reft+moss+sand+twigs+water, 

                       data = spider, distribution = approximate(B = 19999))

statistic(it, "standardized")

pvalue(it)

@

\begin{figure}[t]

\begin{center}

<<spider-plot, eval = FALSE, echo = TRUE, fig = TRUE>>=

plot(sptree, terminal_panel = node_terminal)

@

\caption{Regression tree for hunting spider abundance. The species with

         maximal abundance is given for each terminal node. \label{spider}}

\end{center}

\end{figure}