Avner Bar-Hen

Servane Gey

Jean-Michel Poggi

Spatial CART Classification Trees

Keywords: CART, Bivariate Marked Point Process, Spatial CART, Ripley's intertype K-function

We propose to extend CART for bivariate marked point processes to provide a segmentation of the space into homogeneous areas for interaction between marks. While usual CART tree considers marginal distribution of the response variable at each node, the proposed algorithm, SpatCART, takes into account the spatial location of the observations in the splitting criterion. We introduce a dissimilarity index based on Ripley's intertype K-function quantifying the interaction between two populations. This index used for the growing step of the CART strategy, leads to a heterogeneity function consistent with the original CART algorithm. Therefore the new variant is a way to explore spatial data as a bivariate marked point process using binary classification trees. The proposed procedure is implemented in an R package, and illustrated on simulated examples. SpatCART is finally applied to a tropical forest example.

Introduction

We propose an extension of CART for bivariate marked spatial point processes, which takes into account the spatial information in the splitting criterion and providing a segmentation of the space into homogeneous areas for interaction between marks.

CART (Classification And Regression Trees

) is a statistical method, introduced by Breiman et al. [START_REF] Breiman | Classification and regression trees[END_REF], and designing tree predictors for both regression and classification. We restrict our attention on the classification case with two populations. Each observation is characterized by some input variables gathered in vector X and a binary label Y which is the output (or response) variable. The general principle of CART is to recursively partition the input space using binary splits and then determine an optimal partition for prediction. The classical representation of the model relating Y to X is a tree representing the underlying process of construction of the model. If the explanatory variables are spatial coordinates, we get a spatial decision tree and this induces a tessellation of the space (of the X variables). A cell of this tessellation corresponds to a leaf of the decision tree. For a leaf of this tree, the response variable Y is constant and corresponds to the majority label of the observations belonging to this leaf.

The aim of this article is to adapt CART decision trees within the framework of binary marked point processes by considering that the points inside the cells defined by the partition are a realization of a spatial point process, where the response variable Y can be seen as a mark of the point of position given by X. Following this idea, the proposed variant of the CART algorithm takes into account the spatial dimension of the bivariate marked point processes, and thus makes it possible to provide the user a partition of the plan defining homogeneous areas for interaction between the two marks. Hence, the new variant is a way to explore spatial data as a bivariate marked point process using binary classification trees.

While a classical assumption for CART is to consider a sample of i.i.d. observations, our variant takes into account spatial dependence using the Ripley's intertype K-function to suitably redefine the splitting criterion in CART. Instead of being based only on the respective marks' proportions, the splits will be selected according to the strength of interaction between the two populations. When exist, the spatial patterns corresponding to homogeneous interaction sub-windows cannot, in general, be discriminated by traditional CART. Indeed, CART segments space according to the marks' proportions, favoring nodes of individuals of the same mark, without any reference to the spatial positions of the individuals. At the contrary, the new impurity function based on spatial interaction between the marks is adapted to capture homogeneous interaction sub-windows.

Let us mention some akin methods, out of the scope of this paper, since they are related to the regression case. Among the various extensions of CART, Bel et al. ( [START_REF] Bel | CART algorithm for spatial data: application to environmental and ecological data[END_REF]) proposed a variant for spatial data by considering the observations as a sample of regionalized variables modeled as a random field with spatial dependence. Various reweighting schemes to equalize the contribution of each surface unit for the choices of splits are proposed. In [START_REF] Loecher | RgoogleMaps and loa: Unleashing R graphics power on map tiles[END_REF], the authors suggested the use of recursive partitioning scheme like CART for space segmentation for the spatio-temporal hotspot detection problem, where the response variable is continuous. We could also cite some methods able to take into account the spatial effect even if they are not especially designed to cope with spatial data. For example a model-based regression tree with a spatial lag is considered in [START_REF] Wagner | Heterogeneity and spatial dependence of regional growth in the EU: A recursive partitioning approach[END_REF]. Complexifying the approach, by merging the possibility of incorporating covariates and taking into account the spatial effect, we could refer to the GAM regression models framework. For example, the BAMLSS family [START_REF] Umlauf | BAMLSS: Bayesian additive models for location, scale, and shape (and beyond)[END_REF] which allows the parameters of the response distribution to be modeled by explanatory variables. Then ensembles of models like the BART (Bayesian Additive Regression Trees) [START_REF] Ha Chipman | BART: Bayesian additive regression trees[END_REF] method implementing a sum-of-trees model and a regularization prior on the parameters of that model. Finally, other extensions of GAM combined with boosted trees may also be mentioned, see for example [START_REF] Hofner | Model-based boosting in R: a hands-on tutorial using the R package mboost[END_REF].

The paper is organized as follows. In Section 2, we recall some basics about CART decision trees in the classification case. In Section 3, we review some conventional ways to quantify interaction between two point processes and we define the heterogeneity function associated with the Ripley's intertype K-function. Then, in Section 4, we propose a spatial variant of CART dealing with point processes and describe how to implement it. Section 5 illustrates its use on simple simulated examples. Section 6 addresses an application to the spatial distribution of two species of tropical forest. Finally, the concluding Section 7 sketches some perspectives.

CART method

Let us briefly recall some general background on classical settings about Classification And Regression Trees (CART). The data are considered as an independent sample of the random variables (X 1 , . . . , X p , Y ), where the X k s are the explanatory variables (supposed to be numerical in this article) and Y is the categorical variable to be explained. In a general framework, Y is a multinomial variable, but we suppose here that Y is a binary variable with two labels. CART is a rule-based method that generates a binary tree through recursive partitioning. Let us consider a classification tree T , where a class label is assigned to each terminal node (or leaf) of T . Hence T can be viewed as a mapping to assign a value Y i = T (X 1 i , . . . , X p i ) to each observation.

The growing step leading to a deep maximal tree is obtained by recursive partitioning of the training sample by selecting the best split at each node according to some heterogeneity index, such that it is equal to 0 when there is only one class represented in the node to be split, and is maximum when all classes are equally frequent. The two most popular heterogeneity criteria are the Shannon entropy and the Gini index. Among all binary partitions of each set of values of the explanatory variables at a node t, the principle of CART is to split t into two sub-nodes t L and t R according to a threshold on one of the variables, such that the reduction of heterogeneity between a node and the two sub-nodes is maximized. Some variables may be used several times while others may not be used at all. The growing procedure is stopped when there is no more admissible splitting. Each leaf is assigned to the most frequent class label of its observations. Of course, such a maximal tree (denoted by T max ) generally overfits the training data and the associated prediction error R(T max ), with

R(T ) = P(T (X 1 , . . . , X p ) = Y ), (1)
is typically large. Since the goal is to build from the available data a tree T whose prediction error is as small as possible, in a second stage the tree T max is pruned to produce a subtree T whose expected performance is close to the minimum of R(T ) over all binary subtrees T of T max . Since the joint distribution P of (X 1 , . . . , X p , Y ) is unknown, the pruning is based on a penalized empirical risk to balance optimistic estimates of empirical risk by adding a complexity term α|T | that penalizes larger subtrees, with |T | the number of leaves of tree T and α a positive penalty constant. The pruning step produces a decreasing sequence of K subtrees pruned each one from another, denoted by T 1 . . . T K = {t 1 }, and associated with an increasing sequence of K complexities, denoted by 0 = α 1 < . . . < α K .

According to the penalized criterion used in the pruning procedure, the sequences of subtrees pruned from the same maximal tree are quite different. In the most commonly used variant of CART classification trees, the aim is to predict labels Y ∈ {m 1 ; m 2 }. Hence the natural penalized criterion Rpen (T ) used in this context is based on the misclassification rate of a tree T :

Rpen (T ) = 1 n n i=1 1 {T (X 1 i ,...,X p i ) =Y i } + α|T | (2)
where 1 is the indicator function, n the total number of observations, α the positive penalty constant, |T | the number of leaves of the tree T , and Y i the ith random realization of Y .

In the variant of CART class probability trees, the aim is to predict probabilities p(j|x) = P(Y = m j | X = x) rather than labels Y ∈ {m 1 ; m 2 }. Hence this variant prunes the maximal tree through a penalized criterion based on estimation error of the (p(j|x)) j=1,2 rather than misclassification rate. As it will be seen in subsection 4.2, this kind of variant is useful when dealing with marks' spatial intensity rather than marks' majority vote. The penalized criterion used for class probability trees is derived from the mean square error risk as follows. Let (Z 1 , Z 2 ) be defined by Z j = 1 {Y =m j } and let p(j|x

) = P(Y = m j | X = x). Then E(Z j | X = x) = p(j|x), j = 1, 2.
In this context, the criterion is a penalized version of the Gini index. An empirical version of penalized Gini index can be written as [7, sec. 4.6]:

Rpen G (T ) = 1 n t∈ T n t   1 - 2 j=1 p(j|t) 2   + α|T |, (3)
where n t is the number of observations falling in node t, and p(j|t) the proportion of observations with mark j falling in node t.

The maximal tree T max is computed via recursive partitioning, and then pruned with the penalized criterion based on misclassification rate defined in (2) to produce classification trees, or based on Gini index defined in (3) to produce class probability trees. For all α ≥ 0, argmin T Tmax Rpen (T ), or argmin T Tmax Rpen G (T ), belongs to the sequence of nested pruned subtrees (T k ) 1 k K . The final step is to choose a convenient tree among the sequence T 1 . . . T K = {t 1 } or equivalently a convenient complexity among the sequence 0 = α 1 < . . . < α K . This choice can be made via cross-validation.

In this article, p = 2 since we only consider the spatial coordinates as available information.

3 Quantifying the link between two point patterns

Basics on point processes

A point process is a random variable that gives the localization of events in a set W ⊂ R d . Another way to define a given point process is to consider, for each B ⊂ W , the number of events φ(B) occurring in B, where φ is the distribution of the occurrences of the point process.

Since characterization of a spatial repartition is strongly dependent on the scale of observation, the repartition has to be characterized for each scale.

There are classical assumptions about point processes. At first we consider that the probability to observe two points at the same place is null. Up to misrecording, this hypothesis is not very restrictive. Two extra common hypotheses are stationarity and isotropy. Intuitively, it means that the distribution of the occurrences of the point process φ is not affected by translation or rotation around the origin, i.e. the characteristics of the point process are the same for the whole area under study. A point process that is stationary and isotropic is said to be homogeneous. A marked point process is a point process such that a random mark is associated with each localization. In this article, we only consider bivariate point processes, i.e. the mark is a qualitative random variable with two possible issues. Equivalently, the bivariate point process can be viewed as the realization of two point processes (one per level of the mark).

There are several ways to consider the relationships between two clouds of points, mainly related to three aspects: independence, association and random labelling (see [START_REF] Bar-Hen | Simulation study of dissimilarity between point process[END_REF] for example). It ends up that relationships between two clouds of points can be described in various ways and therefore many indices can be defined. Each index will give a specific information about these relationships and will greatly depend on the point process that leads to the observed repartition. For bivariate point processes, many tools based on first-order characteristics of the point processes, may be used to quantify departure from independence (see [START_REF] Cressie | Statistics for Spatial Data[END_REF] for example).

Intertype K-function

Definition

Under the assumptions of stationarity and isotropy, the intertype K ij -function is a bivariate extension of Ripley's K-function, proposed in [START_REF] Hw Lotwick | Methods for analysing spatial processes of several types of points[END_REF], and defined as

K ij (r) = λ -1 j E(
number of points of type j within distance r of a randomly chosen point of type i) (4) where the intensity parameters λ i and λ j correspond to the expected numbers of type i and type j points per unit area, respectively.

The idea of intertype K ij -function is to characterize interaction between two point processes, i.e. second order characteristic of the bivariate point process. The assumption of stationarity is quite a strong hypothesis and could be relaxed using local estimation of the intensity [START_REF] Baddeley | Non-and semiparametric estimation of interaction in inhomogeneous point patterns[END_REF], but this leads to instability of the estimator of the intertype K ij -function. While Ripley's K-function characterizes the spatial structure of a univariate pattern at various scales, the intertype K ij -function characterizes the spatial structure of a bivariate pattern, and more precisely the spatial relationship between two types of points located in the same study area. The intertype K ij -function is defined so that λ j K ij (r) is the expected number of type j points in a circle of radius r centered on an arbitrary type i point of the pattern. Symmetrically, we can define an intertype K ji -function so that λ i K ji (r) is the expected number of type i points in a circle of radius r centered on an arbitrary type j point. If the bivariate spatial process is homogeneous, then K ij (r) = K ji (r). Under independence, the intertype K function is K ij (r) = πr 2 , regardless of the individual univariate spatial patterns of the two types of events.

In order to interpret the observed values of K ij (r), we need to compare them to the theoretical values obtained for simple cases of bivariate patterns, and especially to πr 2 that corresponds to a null hypothesis of absence of interaction between the two types of points. Positive values of K ij (r) -πr 2 indicate attraction between the two processes at distance r while negative values indicate repulsion. However, depending on the context of the study, this appropriate null hypothesis can correspond to at least two different statistical hypotheses (see [START_REF] Pj Diggle | Second-order analysis of spatial clustering for inhomogeneous populations[END_REF]): independence or random labelling. Note that these two hypotheses are different. On one hand, independence of the marginal processes states that the point process that gives the locations of the first point process is independent of the point process that gives the locations of the second point process. The two processes are therefore distributed independently. On the other hand, random labelling states that, conditionally on the locations of points, their mark is assigned to them at random.

Estimation

The estimator of the intertype K ij -function can be defined by:

Kij (r) = ( λi λj A) -1 k,l 1 {d i k ,j l <r} (5)
where d i k ,j l is the distance between the kth location of type i point and the lth location of type j point, A is the area of the region of interest and λi and λj are the estimated intensities. This estimator Kij (r) is a function of the distance between points. It can be integrated to provide single valued indices (see [START_REF] Pj Diggle | Bivariate Cox processes: some models for bivariate spatial point patterns[END_REF] for example).

As the theoretical distributions of the estimators are unknown, confidence intervals are commonly estimated through Monte Carlo simulations of a specified null hypothesis. To test independence, a classical method is to keep the patterns of both point processes unchanged, but to randomize their relative position at each Monte Carlo simulation. On the other hand, a classical approach to test random labelling is to simulate realizations of a point process with the same spatial structure as the overall observed pattern (i.e. without type distinction), and to randomly attribute marks.

Remark 1. Various edge corrections have been suggested; one common example is the extension of Ripley's estimator, leading to:

Kij (r) = ( λi λj A) -1 k,l w(i k , j l )1 {d i k ,j l <r} .
The coefficient w(i k , j l ) is the inverse of the proportion of the perimeter of the circle centered at the kth location of type i point with radius d i k ,j l that lies inside the study area. Ripley ( [START_REF] Bd | Modelling spatial patterns[END_REF]) shows that this corrected estimator is unbiased. When edge corrections are used Kij (r) and Kji (r) are positively correlated and no more equal. In the sequel, in order to obtain explicit formulation for variation of heterogeneity, we use a Ripley's intertype function without edge correction as defined in Equation 5.

4 Spatial CART

Impurity loss based on K ij

The key idea is to take into account the spatial dependence of the data in the CART splitting strategy. It is done by modifying the original impurity loss, which is usually the entropy index. We introduce a dissimilarity index based on Ripley's intertype K-function quantifying the interaction between two populations within a rectangular window of finite area A.

Let us focus on the impurity loss associated with Kij as defined in [START_REF] Bar-Hen | Simulation study of dissimilarity between point process[END_REF]. First, we define the impurity of a node t at radius r as Kt ij (r), the estimation of the Ripley's intertype K-function restricted to node t:

Kt ij (r) = 1 λt i λt j A t i k ,j l ∈t 1 {d i k ,j l <r} , (6)
where one has:

A t the area of node t λt * the estimation of the density of mark * in node t( * = i, j)

d i k j l the euclidean distance between i-marked individual i k and j-marked individual j l
Then, for a split s splitting node t into two child nodes t L and t R , the impurity (6) of node t at radius r is decomposed as the sum of the contribution within nodes and between nodes:

Kt ij (r) = 1 λt i λt j A t i k ,j l ∈t 1 {d i k ,j l <r} = 1 λt i λt j A t i k ,j l ∈t R 1 {d i k ,j l <r} + 1 λt i λt j A t i k ,j l ∈t L 1 {d i k ,j l <r} (7)
+ 1 λt i λt j A t   i k ∈t R j l ∈t L 1 {d i k ,j l <r} + i k ∈t L j l ∈t R 1 {d i k ,j l <r}   = A t L A t λt L i λt L j λt i λt j Kt L ij (r) + A t R A t λt R i λt R j λt i λt j Kt R ij (r) (8)
+ 1 λt i λt j A t   i k ∈t R j l ∈t L 1 {d i k ,j l <r} + i k ∈t L j l ∈t R 1 {d i k ,j l <r}  
This leads to the natural definition of the impurity loss ∆I ij (s, t, r) as the variation of heterogeneity, i. e. the variation of interaction between marks i and j, coming from the split of node t using s, at radius r:

∆I ij (s, t, r) := Kt ij (r) -α s λt L i λt L j λt i λt j Kt L ij (r) -(1 -α s ) λt R i λt R j λt i λt j Kt R ij (r), (9)
where α s = A t L A t . The intensity ratios are similar to marks' proportions within each child node in a spatial context, leading to a classical form for ∆I ij (s, t, r) as an impurity loss, up to the area factor α s . Nonetheless, this area factor is natural when dealing with spatial data since it leads to reweight properly the impurity of the two nodes t R and t L .

Finally, this choice of impurity loss is convenient since it can be rewritten as (see [START_REF] Ha Chipman | BART: Bayesian additive regression trees[END_REF] and ( 9)):

∆I ij (s, t, r) = 1 λt i λt j A t   i k ∈t R j l ∈t L 1 {d i k ,j l <r} + i k ∈t L j l ∈t R 1 {d i k ,j l <r}   . ( 10
)
∆I ij (s, t, r) is positive, which is necessary to define it as impurity loss. In addition, ∆I ij (s, t, r) is null if and only if the children nodes t L and t R are pure at distance r along split s, that is

∀i k ∈ t L , j l ∈ t R and ∀i k ∈ t R , j l ∈ t L d i k ,j l r,
highlighting splits that do not discriminate marks at all. Hence maximizing ∆I ij (s, t, r) leads to increasing spatial purity at fixed scale r.

In addition to the positivity of s → ∆I ij (s, t, r), which is mandatory, a desirable property is the strict concavity, which would ensure that the best split is unique, and then avoids ties. This is an ingredient of the original CART algorithm but this is not the case here as in most extensions of CART. It should be noted that the growing part of the algorithm still works, that is splitting always purifies nodes even if from an algorithmic point of view, the choice of the split could be arbitrary, without any statistical drawback. However, a possibility, not implemented here, could be to look for a resolution r avoiding ties.

Description of the algorithm

The usual penalized misclassification error rate (Equation 2) used in CART classification trees is natural to determine areas that are homogeneous in terms of mark. Our focus is to determine areas of the point process that are homogeneous in terms of marks' intensity. To be able to choose a convenient definition for the risk (as defined in Equation 3), we consider the following property of the intensity of the marks: in the case where the marked spatial point process (X, M ) is stationary, then the intensity of points with mark j ∈ {1; 2} inside surface A can be written as

Λ(A, j) = P(X ∈ A, M = j) = P(M = j | X ∈ A)P(X ∈ A) = P(M = j | X ∈ A)λA
where λ > 0 is the intensity of X (without marks), and A is the area of A. Then estimating Λ(A, j) is equivalent to estimate P(M = j | X ∈ A), since λ can be estimated directly from the point process. Thus the penalized criterion of Equation 3 used in CART class probability trees is natural to estimate the intensities of the marks: intensity for each mark will then be locally estimated on a tessellation of the plane.

We propose an algorithm using impurity loss ∆I ij defined by [START_REF] Cressie | Statistics for Spatial Data[END_REF] to develop the maximal tree T max . The estimator Kij of the intertype K ij -function is computed at each node t, the value of r = r 0 is fixed as the one for which the estimated intertype K-function is the farthest from the one of random labelling.

Recall that a maximal tree T max is computed via recursive partitioning, and then pruned with the penalized criterion defined in (2) or ( 3). This produces a decreasing sequence of K subtrees pruned each one from another, denoted by T 1 . . . T K = {t 1 }, and associated with a sequence of complexities 0 = α 1 < . . . < α K .

Rephrasing the algorithm in spatial terms, r 0 can be considered as the characteristic scale from the region for which we compute r 0 . Then before splitting, we consider for r = r 0 the quantity ∆I ij (s, t, r) and we seek to the best split. After splitting we seek to the best r ≤ r 0 maximizing ∆I ij (ŝ, t, r). Then, after splitting, the two child are considered in parallel in the same way, recomputing r 0,L and r 0,R . It turns out that ∆I ij (s, t, r) is decreasing along any branch of the tree. So the reordered sequence of r i can be used to define a sequence of nested subtrees. Let us remark that, in the algorithm, the value of r t maximizing the impurity criterion for the best split is set to 0 if t is a leaf, what corresponds to the smallest possible value of r that is the best resolution.

The final step is to choose a convenient tree among the sequence T 1 . . . T K = {t 1 }. This choice can be made via cross-validation. Nevertheless, to avoid randomness due to the choice of subsamples, we use two variants of the slope heuristic method proposed by Birgé, Massart in the 2000s (see [START_REF] Arlot | Minimal Penalty and the Slope Heuristic: A Survey (with discussion)[END_REF][START_REF] Jp Baudry | Slope heuristics: overview and implementation[END_REF][START_REF] Gey | Using CART to detect Multiple Change Points in the Mean[END_REF] for more details), and based on the behavior of the number of leaves with respect to the complexity. First, the general idea is that, if a tree is a good predictor, then high energy is necessary to prune it. Hence, the penalty to be chosen in the penalized criterion shall increase a lot before the tree is pruned, leading to a large plateau in the α k → |T k | graph. Second, if there exist plateaus negligible with respect to the largest one, the corresponding trees can be considered as artefacts of the pruning procedure. These trees are then removed from the sequence, and the final tree model is selected through the classical maximal jump in the number of leaves for this new collection of trees. In practice, only jumps occurring before the largest plateau are considered. An artificial example of such a behavior is given in Figure 1, where represents the tree selected via the modified largest jump method, while represents the tree selected via the largest plateau method.

The algorithm used to obtain the final tree is detailed in Table 1 for the calculation of the maximal tree, and in Table 2 for the pruning strategy, and the calculation of the nested sequence of pruned subtrees. It should be noted that this pruning step is the same as for original CART, since it is based either on the minimization of penalized misclassification rate or on the minimization of penalized Gini index. Let us emphasize that this sequence of pruned subtrees and the corresponding complexities allow the user to choose the final partition according to the data and the peculiarity of the studied problem. This is of particular interest if the best subtree is the root.

Remark 2. To avoid too strong edge effects, a minimum distance from the boundaries is imposed to candidate splits in the construction of maximal tree.

The initial resolution r 0

The other parameter of our algorithm is the initial resolution r 0 , which is crucial to define the first split (the most important one). One first idea to provide a default value is to take as initial value for r 0 the one corresponding to the maximum of the difference between the estimated K ij and its theoretical value, that is the one for which the dissimilarity with the homogeneous case is maximal. Depending on the criterion configurations, it appears that critical values can be far from this default value and user should be careful about automatic choice of r 0 .

CART and SpatCART in action: illustration by simulations

An R package spatcart, and the R codes to reproduce experiments of this section, are available on https://github.com/Servane-Gey/Spatial-classification-trees. Package spatcart may also be directly installed with R package devtools from the github repository Servane-Gey/spatcart. Package spatcart requires the following R packages to implement the results: Pruning Sequence (T k ) 1 k K of subtrees pruned from Tmax, and complexities (α k ) 1 k K (see Table 2).

Selection Set k1 = max k | k = argmax 0 j K-1 (αj+1 -αj) , k2 the index of the tree selected via the modified largest jump method, Output Trees Tmax, T k1 , T k2 . Sequences T1 . . . TK = {t1} and 0 = α1 < . . . < αK . Collections rt and Kij(rt) for each node t of Tmax. spatstat to deal with point processes, and in particular to compute ∆I ij in the construction of the maximal tree, tree to deal with tree structures.

Two methods are to be considered: Spatial CART (SpatCART, see Table 1), and CART with Gini splitting criterion. Class probability trees and classification trees are constructed using

Pruning Algorithm

Input Maximal tree Tmax, number of observations n, risk ρ from misclassification error or Gini index.

Initialization k = 1 and α1 = 0, T = T1 the smallest subtree pruned from Tmax at complexity α1, n f the number of leaves of T , Compute for each node t of T nt the number of observations in t, ρ(t) the local risk of t:

ρ(t) = n -1 (x i ,m i )∈t 1 mt =m i for misclassification rate, or ρ(t) = (1 -2 j=1 p(j|t) 2
)nt/n for Gini index ρ(Tt) = {f leaf of T t } ρ(f ) the risk of the branch Tt issued from t, nT t the number of leaves of branch Tt.

While While n f > 1, compute α k+1 = min {t internal node of T } ρ(t) -ρ(Tt) nT t -1 .
Prune all branches Tt of T verifying ρ(Tt) + α k+1 nT t = ρ(t) + α k+1 Set T k+1 the pruned subtree obtained in that way. Set T = T k+1 and k = k + 1. either SpatCART or CART. We could expect the following results. First, since the pruning strategy (see Table 2) is the same for the two methods, based either on the minimization of penalized misclassification rate or on the minimization of penalized Gini index, the major differences between our proposal (Spatial CART) and the existing schemes must be concentrated on the growing step. In addition to the final trees, we also have to look at the comparison of the two maximal trees. Second, since the new splitting criterion reduces more or less to the classical one in the case of spatial homogeneity, we have to consider some artificial homogeneous examples to illustrate this and, what is more crucial, to define some artificial inhomogeneous examples to capture the interest of our proposition. This second point must be considered because the spatial nature of the data must be taken into account, and we can refer to [START_REF] Haining | Bivariate correlation with spatial data[END_REF] and [START_REF] Anselin | Spatial statistical analysis and geographic information systems[END_REF].

Output Sequence (T k , α k ) 1 k K .

Illustrative examples

So, the idea here is to check that CART and SpatCART behave as expected on two typical examples: (a) a control example where the Bayes classifier is a classification tree. In this case, the two algorithms could find a partition close to the partition defined by the underlying model. With these two different scenarios in mind, it remains to analyze the results. But before that, we propose to focus on the kernel of our proposal: the new splitting criterion convenient for detecting a certain kind of spatial heterogeneity and the initial spatial resolution to be selected.

Intertype K-function and splitting criterion

The key tool for defining the splitting criterion is the intertype K-function. More precisely the difference between the estimated one and its theoretical value is evaluated in order to select the split. In each cell of the tesselation, the intensity is assumed to be constant, i.e. the process is assumed to be homogeneous. Figure 3 presents the Chess simulated example: on the left we have the estimated intertype K-function (as a function of r) and its theoretical counterpart and, on the right, the difference between these two functions. In this typical situation, this difference is strictly decreasing and the natural choice is to take r 0 as large as possible. The next object of interest is the impurity function directly connected with the splitting criterion. Figure 4 illustrates the behavior of the impurity function with respect to first split for Chess example, which appears to be clear and expected.

Maximal and optimal partitions

With two different situations depicted by Figure 2, we compare the results of the usual CART strategy (for growing and pruning) and those obtained applying our new proposal.

The homogeneous case

Depending on the run of simulations, the obtained partitions can be roughly the same, or quite different due the ties problem for the split choice, which occurs more frequently in such an example exhibiting high symmetry. Selecting an example for which the maximal trees are similar leads to Figure 5, where SpatCART trees have been obtained from initial resolution r 0 = 0.25 equal to the maximal difference between estimated and theoretical intertype K-functions (see Figure 3). Let us remark that the true underlying frontiers are recovered by the CART and SpatCART optimal classification trees, and that the differences appear inside the blocks of the maximal trees.

Examining the corresponding behaviors of the number of leaves as a function of complexity for SpatCART and CART (see Figure 6), this similarity is visible. So to summarize, as expected, when there is nothing to discover from the spatial viewpoint, SpatCART can behave as the usual CART but in most cases would lead to quite different partitions, due to ties.

The inhomogeneous case

Using roughly the same sequence of plots, let us examine now a more interesting aspect, crucial to illustrate what is the specific contribution of SpatCART. The first step is to compare the maximal trees. We start with an initial resolution r 0 smaller than the locally repulsion scale r = 0.05. Indeed, since the domain consists of a left part where the points of the minority mark push back those of the majority mark to at least 0.05 units, and a right part where the two marks are distributed independently, therefore, to detect this difference, it is necessary to adopt a sufficiently small resolution, which corresponds to values of r 0 lower than the repulsion threshold. data set. The partitions are extremely different. More precisely, since the initial r 0 is small enough, SpatCART does not split the left half of the square, and the right half is split iteratively according to the direction of the first split. This is expected since the precision of the estimation of the intensities is directly related to the number of points and the intensity of the rarer is too low to provide a reliable estimate. A solution could be to impose a minimum number of points for each mark to allow splitting. We have not implemented it since pruning will remove these artificial splits. Note that if the gain of impurity is the same for both direction, the x direction is chosen. On the other hand, CART splits according to the empirical distribution of the blue points, leading to a very different partition. So the generated maximal partitions are, as expected, different in the inhomogeneous situation. The corresponding sequences of complexities for class probability trees are given by Figure 8. A second difference appears: SpatCART seems to generate less splits to provide partitions similar Figure 6: Behavior of number of leaves versus complexity for SpatCART and CART on Chess data set. For each case, represents the tree selected via the modified largest jump method, while represents the tree selected via the largest plateau method.

to CART in terms of spatial repartition, generating less false alarms (useless spurious splits). To end this analysis, let us mention that in this case, results using the pruning strategy are convenient: the selection based on misclassification rate leads to root trees. This is logical since the percentage of blue marks is too low to be recovered. Here, the difference is clear: since there is no heterogeneity in the marginal structures of the response distribution, nothing is inferred using CART.

To summarize, SpatCART seems to provide better partitions in terms of marks spatial repartition (see the maximal trees figures) and SpatCART needs less splits to provide partitions similar to CART in terms of spatial repartition, suggesting that it avoids false alarms in a better way.

Example

Data

We applied these methods to a tropical rain-forest located at Paracou, 40 km west from Kourou in French Guiana (5 • 15'N, 52 • 55'W). It is an experimental site that is devoted to studying the effects of logging damage on stock recovery. Twelve plots (6.25 hectares each) of undisturbed forest were settled in 1984. On each plot, the circumference of every tree with a dbh (diameter at breast height) greater than 10 cm was measured with a precision of 0.5 cm. Its spatial coordinates (± 50 cm) and its specie were noted too. Measurements have been carried annually from 1984 to 1995, and once every two years since. A more precise description of the Paracou plots may be found in [START_REF]Ecology and Management of a Neotropical Rainforest: Lessons Drawn from Paracou, a Long-term Experimental Research Site in French Guiana[END_REF] and the data can be freely downloaded from the website https://paracou.cirad.fr/. Experiments on the Paracou data set may be reproduced with We focus on two species:

• Vouacapoua americana is a hermaphroditic shade tolerant tree specie of mature tropical rainforests and whose distribution spans the eastern part of the Guiana shield. Its local density averages around 10 individuals per hectare for trees with greater than 10 cm dbh, but this varies greatly because of spatial clustering. Individuals are clustered in large patches of a few hectares that are mainly located on hill tops and slopes [START_REF] Traissac | Dynamique spatiale de Vouacapoua americana (Aublet), arbre de forêt tropicale humide à répartition agrégée[END_REF].

• Oxandra asbeckii is a specie of the understorey, the largest individuals of which measure less than 15 cm dbh. It is a shade tolerant specie [START_REF] V Favrichon | Classification des espèces arborées en groupes fonctionnels en vue de la réalisation d'un modèle de dynamique de peuplement en forêt guyanaise[END_REF]. It shows animal dispersal: the seeds are dispersed after passage through the vertebrate gut. Oxandra asbeckii is located on hill tops and slopes. data set, with initial resolution r0 smaller than the locally repulsion scale r = 0.05. For each case, represents the tree selected via the modified largest jump method, while represents the tree selected via the largest plateau method.

This explains the spatial distribution for the two species displayed on Figure 11. The two species Vouacapoua americanaa and Oxandra asbeckii were selected at Paracou because their spatial distribution is linked to the relief: they are both located on hill tops and slopes. Elevation is the environmental factor that drives their spatial distribution and this creates a strong interaction between both repartitions. We focus on one plot represented in Figure 11. The figure also gives the contour lines in order to be able to determine hill tops and slopes of the plot. Seventy trees of Vouacapoua americana and eighty trees of Oxandra asbeckii were referenced for this plot. The data consists of seventy lines (one per tree) and four columns: the 3-D coordinates (longitude, latitude and elevation) as well as the specie indication.

Results and discussion

From the representation of the estimated and theoretical intertype K-functions and their difference evaluated on the Paracou data set (see Figure 9), we can see three regimes: one before the scale value r = 6 (blue dashed line on the right of Figure 9), where there is no interaction between species; one between scale values r = 6 and r = 24 (red dashed line on the right of Figure 9) where species begin to interact; and one after the scale value r = 24 where the interaction between species increases rapidly. From an ecological point of view, r = 24 can be interpreted as the limit of seeds dispersal (see [START_REF] Traissac | Dynamique spatiale de Vouacapoua americana (Aublet), arbre de forêt tropicale humide à répartition agrégée[END_REF]). Hence, we choose the initial median scale value r 0 = 15 to proceed with SpatCART, in order to be sure to catch a sufficiently large scale to capture interaction, and not too large to avoid deeper maximal trees. To better interpret the regions highlighted by the algorithm in terms of interaction between marks, we color each region of the space (each leaf t) using a monotonic gray colormap based on the value of Kt ij (r t ) (see equation 6). Thus Figure 10 colouring the pieces of the partition using this interaction-based colormap, highlights the presence of Oxandra asbeckii at the hill of left top of the plot as well as the competition between both species for the hill at the bottom of the plot. This representation magnifies the regions of high interaction, showing that SpatCART recovers the spatial structure and that the interaction-based grey colormap is meaningful. In addition, the colormap based on interaction allows to visually appreciate the difference between two adjacent regions. This opens a very natural post-processing by ignoring useless frontiers and provides more interpretable tessellations.

The class probability tree has more leaves and if both resulting partitions are coherent, even exactly the same from the right side of the window space, an interesting difference must be noted. On the right side of the window space, the class probability tree partition provides a more detailed segmentation (possibly induced by the slope of the hill) but cleaning the two partitions according to the interaction-based colormap leads to highlight the same areas. A contrario, CART results are really poor with only two leaves, for the modified largest jump variant (left panel of Figure 11). Basically it separates the hill at the bottom of the plot from the rest but cannot catch the mixed structure of species along the hill. The spatial structure as well as the ecology of the two species on this plot cannot be inferred from CART results. CART results using the largest plateau variant (see the right panel of Figure 11) are not more informative for the spatial viewpoint, highlighting the regions only according to the most prevalent specie.

From an ecological perspective, few points can be highlighted. At first, elevation is a key variable to describe interactions between these two species. Even if this two species are shade tolerant, competition for light is a major determinant of forest composition. Hills intercept a greater portion of incoming solar radiation and thus, generate contrasting strategies between species to invest a greater portion of their resources in expanding their crowns. The result is not surprising but the method allows to determine more precisely the interaction areas. On the other hand, below-ground competition for water is also an important driver of tree growth. Slope does not store water and it is a strong disadvantage during drought events. Competition for water among neighbouring trees clearly lead an advantage to Vouacapoua americana compared to Oxandra asbeckii.

Conclusion and perspectives

For a bivariate marked spatial point process, we extend CART to segment the space into homogeneous areas for interaction between marks. The original CART constructs homogeneous zones with respect to the distribution of the variable of interest (here, the mark) conditionally to the explanatory variables (here, the position in space). By modifying the splitting criterion in the CART algorithm, using an empirical version of the intertype K-function, we obtain a new procedure, called SpatCART, adapted to the problem at hand. The intertype function itself depends on a parameter r 0 which must be carefully chosen: not to set it once and for all, but rather to start from a rather large value (at the root of the tree) and gradually decrease it as the tree is built by SpatCART. The proposed variant is a way to explore spatial data as a bivariate marked point process using binary classification trees.

This paper is strongly dependent on the binary case, some extensions to handle the multiclass case can be sketched. One possibility is to adopt the strategy used to handle multiclass Support Vector Machines (SVMs) which are intrinsically two-class classifiers (see [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]). A technique widely used in practice is to build one-versus-rest classifiers. Then we could obtain several tessellations and select the final partition by taking into account some criterion maximizing some global measure of heterogeneity between cells, which could be related to the K ij (r * ) for a suitably chosen scale r * .

Another possible extension is to modify, not only the growing step of the CART algorithm, but also the pruning strategy in order to drive it by point processes properties. This requires to define some additive measure of heterogeneity of the partitions. An intermediate solution could be to use a classical pruning step to simplify and avoid spurious useless splits but the final choice could be to select the final tree among the nested sequence of tree by maximizing the same kind of heterogeneity criterion mentioned before.

Instability is one of the main drawbacks of CART and many classical ways to contain it are available. Ensemble methods like bagging, boosting and random forests (see [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF] for example) can then be similarly defined in the spatial case but the aggregation part should be defined. To introduce and experiment a bagging-like scheme could be of interest for a future work. Note that the sensitivity of CART with respect to rotation could also be used as an alternative way to generate several tessellations. Using such an idea, we could imagine to develop a variant selecting randomly, in case of ties, the splits on each direction to choose, the best candidate per direction and if necessary to randomly select the direction. Such a scheme would lead to equivalent but different partitions which could finally be combined.

Another possible extension is to incorporate covariables. In the specific forestry example, it is possible to think about the introduction of the elevation as a third spatial coordinate. Nevertheless, it is not straightforward to estimate intensity of a point process in higher dimension. Alternatively and more generally, we could imagine to incorporate covariables using a threesteps procedure. First by performing a classical CART using only the additional covariables and then, in each leaf, to apply a SpatCART on spatial variables and finally select the best space tessellation or aggregate them.

R package and source codes

An R package spatcart, and the R codes to reproduce experiments of sections 5 and 6, are available on the repository https://github.com/Servane-Gey/Spatial-classification-trees. Package spatcart may also be directly installed with R package devtools from the github repository Servane-Gey/spatcart.

Figure 1 :

 1 Figure 1: Typical behavior of number of leaves with respect to complexity, and represent the two highlighted trees.

  (b) an example where the classification problem is trivial, and where the interaction between marks plays an important role. In this case, SpatCART must produce splits significantly different from those produced by CART. Let us detail the two illustrative simulated examples.The first example is the Chess bivariate Poisson point process on the unit square, with marks simulated from a blue and red chessboard with 9 squares, represented in the left panel of Figure 2. It is an example of spatial homogeneity for which CART and SpatCART could lead to similar results.The second example is the Locally repulsive bivariate Poisson point process on the unit square, with a majority red mark, and a minority blue one. The marked point process is designed by splitting vertically the unit square into two parts (see the right panel of Figure2): on the left part, the minority blue marked points repulse red ones, with a constant repulsion radius equal to r = 0.05 ; on the right part, blue and red marked points are independently distributed.In this example, a spatial inhomogeneity is introduced, then CART and SpatCART should lead to different results: CART should miss the difference between the two parts and SpatCART should highlight it.

Figure 2 :

 2 Figure 2: Left: Chess simulated data set. Right: Locally repulsive simulated data set.

Figure 3 :

 3 Figure 3: Example of intertype K function and its theoretical value in the homogeneous situation for Chess simulated example. Left : Intertype K-functions. Right : Difference between estimated and theoretical intertype K-functions.

Figure 4 :

 4 Figure 4: Behavior of impurity with respect to first split for Chess example.

Figure 7

 7 contains the SpatCART (top) and the CART (bottom) trees on Locally repulsive

Figure 5 :

 5 Figure 5: Maximal (left) and optimal (right) trees for SpatCART (up) and CART (bottom) on Chess data set with initial resolution r0 = 0.25. The pruned subtrees are selected through the largest plateau method (see Figure 6)

Figure 7 :

 7 Figure 7: Maximal (left) and optimal (right) trees for SpatCART (up) and CART (bottom) on Locally repulsive data set, with initial resolution r0 smaller than the locally repulsion scale r = 0.05.

Figure 8 :

 8 Figure 8: Behavior of number of leaves versus complexity for SpatCART and CART on Locally repulsive

Figure 9 :

 9 Figure 9: Intertype K function and its theoretical value in the homogeneous situation for the Paracou data set. Left : Intertype K-functions. Right : Difference between estimated and theoretical intertype K-functions; blue: scale r = 6, red: scale r = 24.

Figure 10 :

 10 Figure 10: SpatCART class probability (left) and classification (right) optimal trees on the Paracou data set with initial resolution r0 = 15. The grey colored Interaction between marks is the estimation of the intertype K-function restricted to each leaf of the tree as defined by (6).

Figure 11 :

 11 Figure 11: CART common optimal class probability and largest plateau optimal classification tree (left), and modified largest jump variant optimal classification tree (right) on the Paracou data set.

Table 1 :

 1 Spatial Classification Trees (SpatCART).

Table 2 :

 2 Pruning Algorithm.