If we grow a tree with standard functions in R, on the same dataset used to introduce classification tree in some previous post,

> MYOCARDE=read.table( + "http://freakonometrics.free.fr/saporta.csv", + head=TRUE,sep=";") > library(rpart) > cart<-rpart(PRONO~.,data=MYOCARDE)

we get

> library(rpart.plot) > library(rattle) > prp(cart,type=2,extra=1)

The first step is to split the first node (based on the whole dataset). To split it, we can use either Gini index

> gini=function(y,classe){ + T=table(y,classe) + nx=apply(T,2,sum) + n=sum(T) + pxy=T/matrix(rep(nx,each=2),nrow=2) + omega=matrix(rep(nx,each=2),nrow=2)/n + g=-sum(omega*pxy*(1-pxy)) + return(g)}

or the entropy

> entropie=function(y,classe){ + T=table(y,classe) + nx=apply(T,2,sum) + n=sum(T) + pxy=T/matrix(rep(nx,each=2),nrow=2) + omega=matrix(rep(nx,each=2),nrow=2)/n + g=sum(omega*pxy*log(pxy)) + return(g)}

For instance, if we choose to split according to the first variable, with threshold 2.5, Gini index would be

> CLASSE=MYOCARDE[,1]<=2.5 > gini(y=MYOCARDE$PRONO,classe=CLASSE) [1] -0.4832375

To get the “optimal” split, we consider all variable, and all threshold (according to some constraint, perhaps, e.g. at leat 5 observations per node)

> mat_gini=mat_v=matrix(NA,7,101) > for(v in 1:7){ + variable=MYOCARDE[,v] + v_seuil=seq(quantile(MYOCARDE[,v], + 6/length(MYOCARDE[,v])), + quantile(MYOCARDE[,v],1-6/length( + MYOCARDE[,v])),length=101) + mat_v[v,]=v_seuil + for(i in 1:101){ + CLASSE=variable<=v_seuil[i] + mat_gini[v,i]= + gini(y=MYOCARDE$PRONO,classe=CLASSE)}} > par(mfrow=c(2,3)) > for(v in 2:7){ + plot(mat_v[v,],mat_gini[v,],type="l", + ylim=range(mat_gini), + main=names(MYOCARDE)[v]) + abline(h=max(mat_gini),col="blue") + }

It looks like we should split according to the second variable (INSYS), as seen on the tree graph, above. Of course, we could be using the entropy

(here we have the same spliting criteria).

Now that we’ve how to split the first node. We keep it.

> idx=which(MYOCARDE$INSYS>=19)

Let us get to the node on the right, when the second variable was above the threshold (here 19). The idea is to run the same code as before, but on that subset,

> mat_gini=mat_v=matrix(NA,7,101) > for(v in 1:7){ + variable=MYOCARDE[idx,v] + v_seuil=seq(quantile(MYOCARDE[idx,v], + 6/length(MYOCARDE[idx,v])), + quantile(MYOCARDE[idx,v],1-6/length( + MYOCARDE[idx,v])), length=101) + mat_v[v,]=v_seuil + for(i in 1:101){ + CLASSE=variable<=v_seuil[i] + mat_gini[v,i]= + gini(y=MYOCARDE$PRONO[idx], + classe=CLASSE)}} > par(mfrow=c(2,3)) > for(v in 2:7){ + plot(mat_v[v,],mat_gini[v,],type="l", + ylim=range(mat_gini), + main=names(MYOCARDE)[v]) + abline(h=max(mat_gini),col="blue")

Here we see that we should split according to the last one (REPUL), exactly as we got on the tree graph.

Hi

5 stars for explaining difficult concepts – the code worked a dream – took a day to understand but I have got it.

You really tackled the key questions of how to generate the first split and proved that it worked by reconciling against R

So many books leave the beginner lost here pointing to Java C++ or to SAS with endless lines of code.

You managed to do it in 6 even though

Really appreciated your work on decision trees -thank you

Please keep up the posts Awesome.

When I run the code without any modification and plot the graphs, I got the following message:

“Error in plot.window(…) : need finite ‘ylim’ values”.

When I check the code without indexing ( ind = pxy > 0), I found that mat_ent contain NaN values for certain variables with certain thresholds. The NaN values were caused by zero elements in the square matrix pxy when we took the log of it. For instance, if we limit our code to variables INCAR and INSYS and we run the entropie function step by step. the log(pxy) ,

> log(pxy)

classe

y FALSE TRUE

DECES -Inf -0.8070914

SURVIE 0.0000000 -0.5908683

which cause the entropie function to generate NaN values.

But, when I remove the zero values from matrix pxy by using index ind (ind = pxy > 0) and plot the graphs, I get different plots. The range of mat_ent are also smaller : between -0.65 and 0.35. My question is the following: if the code need no modification, then how I can get rid off the R error message?.

Thanks for your help and time!.

Salut Arthur,

Thanks for this great blog!.

For the entropy method I used an index to avoid taking the log of zero values;

ind 0)

g = sum(omega[ind]*pxy[ind]*log(pxy[ind])).

But I am still unable to reproduce your graphs???.

Merci!.

how come you get zeros ? in my post, I start with values so that I always have at least 5 people in each node…

Thanks for the nice illustration. Although you point out that this just pertains to the “rpart” package (i.e., CART algorithm), it may be worth mentioning that there are other split variable and split point selection criteria. Specifically, it has been shown that the exhaustive search strategy of CART can lead to a bias towards variables with more splitting points. In this data set this does not matter much because all variables are numeric with similarly many split points. But if there are categorical variables with fewer possible splits, the CART tree might artificially prefer the numeric variables.

Hence, various algorithms have been suggested in the literature that separate variable and split point selection by using statistical tests (e.g., QUEST and GUIDE by Loh and co-workers or CTree and MOB in our “partykit” R package).

As a final remark: While the “pretty rpart plot” is certainly prettier than the default in the “rpart” package, there are other display options that some users find even easier to read. One is to use our generic “partykit” framework to display the “cart” object we have fitted:

library("partykit")

plot(as.party(cart))