{"word_count": 1594, "character_count": 13659}
Nineth post of our series on classification from scratch. Today, we'll see the heuristics of the algorithm inside classification trees. And yes, I promised eight posts in that series, but clearly, that was not sufficient... sorry for the poor prediction.Decision TreeDecision trees are easy to read. So easy to read that they are everywhereWe start from the top, and we go down, with a binary choice, at each stop, each node. Let us see how it works on our datasetlibrary(rpart)cart = rpart(PRONO~.,data=myocarde)library(rpart.plot)prp(cart,type=2,extra=1)We start here with one single leaf. If we have two explanatory variable (the [latex]x[/latex]-axis and the [latex]y[/latex]-axis if we want to plot it), we will check what happens if we cut the leaf accoring to the value of the first variable (and there will be two subgroups, the one on the left and the one on the right)or if we cut according to the second one (and there will be two subgroups, the one on top and the one below).Why and where do we cut? Let us formalize a little bit. A node (a leaf) constains observations, i.e. [latex]\{y_i,\mathbf{x})i\})[/latex] for some [latex]i\in\mathcal{I}\subset\{1,\cdots,n\}[/latex]. Hence, a leaf a caracterized by [latex]\mathcal{I}[/latex]. For instance, the first node in the tree is [latex]\mathcal{I}=\{1,\cdots,n\}[/latex]. A (binary) split is based on one specific variable - say [latex]x_j[/latex] - and a cutoff, say [latex]s[/latex]. Then, there are two options: either [latex]x_{i,j}\leq s[/latex], then observation [latex]i[/latex] goes on the left, in [latex]\mathcal{I}_L[/latex] or [latex]x_{i,j}&#62; s[/latex], then observation [latex]i[/latex] goes on the right, in [latex]\mathcal{I}_R[/latex]Thus, [latex]\mathcal{I}=\mathcal{I}_L\cup\mathcal{I}_R[/latex].Now, define some impurity index, in some node. In the context of a classification tree, the most popular index used (the so-called impurity index) is Gini for node [latex]\mathcal{I}[/latex] is defined as [latex display="true"]G(\mathcal{I})=-\sum_{y\in\{0,1\}}p_y(1-p_y)[/latex]where [latex]p_y[/latex] is the proportion of individuals in the leaf of type [latex]y[/latex]. I use this notation here because it can be extended to the case of more than one class. Here, we consider only binary classification. Now, why [latex]p_y(1-p_y)[/latex]? Because we want leaves that are extremely homogeneous. In our dataset, out of 71 individuals, 42 died, 29 survived. A perfect classification would be obtained if we can split in two, with the 29 survivors on the left, and the 42 dead on the right. In that case, leaves would be perfectly homogneous. So, when [latex]p_0\approx1[/latex] or [latex]p_1\approx1[/latex], we have strong homogenity. If we want an index to maximize, [latex]-p_y(1-p_y)[/latex] might be an interesting candidate. Further more, the worst case would be a leaf with [latex]p_0\approx1/2[/latex], which is exactly what we have here. Note that we can also write[latex display="true"]G(\mathcal{I})=-\sum_{y\in\{0,1\}}\frac{n_{y,\mathcal{I}}}{n_{\mathcal{I}}}\left(1-\frac{n_{y,\mathcal{I}}}{n_{\mathcal{I}}}\right)[/latex]where [latex]n_{y,\mathcal{I}}[/latex] is the number of individuals of type [latex]y[/latex] in the leaf [latex]\mathcal{I}[/latex], and [latex]n_{\mathcal{I}}[/latex] is the number of individuals in the leaf [latex]\mathcal{I}[/latex].If we do not split, we have index[latex display="true"]G(\mathcal{I})=-\sum_{y\in\{0,1\}}\frac{n_{y,\mathcal{I}}}{n_{\mathcal{I}}}\left(1-\frac{n_{y,\mathcal{I}}}{n_{\mathcal{I}}}\right)[/latex]while if we split, define index[latex display="true"]G(\mathcal{I}_L,\mathcal{I}_R)=-\sum_{x\in\{L,R\}}\frac{n_x}{n_{\mathcal{I}_x}}{n_{\mathcal{I}}}\sum_{y\in\{0,1\}}\frac{n_{y,\mathcal{I}_x}}{n_{\mathcal{I}_x}}\left(1-\frac{n_{y,\mathcal{I}_x}}{n_{\mathcal{I}_x}}\right)[/latex]The code to compute is would begini = 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)/ng. = -sum(omega*pxy*(1-pxy))return(g)}Actually, one can consider other indices, like the entropic measure[latex display="true"]E(\mathcal{I})=-\sum_{y\in\{0,1\}}\frac{n_{y,\mathcal{I}}}{n_{\mathcal{I}}}\log\left(\frac{n_{y,\mathcal{I}}}{n_{\mathcal{I}}}\right)[/latex]while if we split, [latex display="true"]E(\mathcal{I}_L,\mathcal{I}_R)=-\sum_{x\in\{L,R\}}\frac{n_x}{n_{\mathcal{I}_x}}{n_{\mathcal{I}}}\sum_{y\in\{0,1\}}\frac{n_{y,\mathcal{I}_x}}{n_{\mathcal{I}_x}}\log\left(\frac{n_{y,\mathcal{I}_x}}{n_{\mathcal{I}_x}}\right)[/latex]entropy = 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)}This index was used originally in C4.5 algorithm.Dividing a leaf (or not)For instance, consider the very first split. Assume that we want to split according to the very first variableCLASSE = myocarde[,1] &#60;=100table(CLASSE)CLASSEFALSE  TRUE    13    58 In that case, there will be 13 invididuals on one side (the left, say), and 58 on the other side (the right).gini(y=myocarde$PRONO,classe=CLASSE)[1] -0.4640415Initially, without any split, it was-2*mean(myocarde$PRONO)*(1-mean(myocarde$PRONO))[1] -0.4832375which can actually also be obtained withCLASSE = myocarde[,1] gini(y=myocarde$PRONO,classe=CLASSE)[1] -0.4832375There is a net gain in spliting ofgini(y=myocarde$PRONO,classe=(myocarde[,1]&#60;=100))-gini(y=myocarde$PRONO,classe=(myocarde[,1]&#60;=Inf))[1] 0.01919591Now, how do we split? Which variable and which cutoff? Well... let's try all possible splits... Here, we have 7 variables. We can consider all possible values, usingsort(unique(myocarde[,1]))But in massive datasets, it can be very long. Here, I preferseq(min(myocarde[,1]),max(myocarde[,1]),length=101) so that we try 101 values of possible cutoff. Overall, the number of computations is rather low, with 707 Gini indices to compute. Again, I won't get back here on the motivations for such a technique to create partitions, I will keep that for the course in Barcelona, but it is fast.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&#60;=v_seuil[i]mat_gini[v,i]=  gini(y=myocarde$PRONO,classe=CLASSE)}}Actually, the range of possible values is slightly different: I do not want cutoff too much on the left or on the right... having a leaf with one or two observations is not the idea, here. Not, if we plot all the functions, we getpar(mfrow=c(3,2))for(v in 2:7){  plot(mat_v[v,],mat_gini[v,],type="l",  ylim=range(mat_gini),xlab="",ylab="",  main=names(myocarde)[v])   abline(h=max(mat_gini),col="blue")}Here, the most homogenous leaves obtained using a cut in two parts is when we use variable 'INSYS'. And the optimal cutoff variable is close to 19. So far, that's the only information we use. Well, actually no. If the gain is sufficiently large, we go for a split. Here, the gain isgini(y=myocarde$PRONO,classe=(myocarde[,3]&#60;19))-gini(y=myocarde$PRONO,classe=(myocarde[,3]&#60;=Inf))[1] 0.2832801which is large. Sufficiently large to go for it, and to split in two. Actually, we look at the relative gain-(gini(y=myocarde$PRONO,classe=(myocarde[,3]&#60;19))-gini(y=myocarde$PRONO,classe=(myocarde[,3]&#60;=Inf)))/gini(y=myocarde$PRONO,classe=(myocarde[,3]&#60;=Inf))[1] 0.5862131If that gain exceed 1% (the default value in R), we split in two.Then, we do it again. Twice. First, on go on the leaf on the left, with 27 observations. And we try to see if we can split it.idx = which(myocarde$INSYS&#60;19)mat_gini = mat_v = matrix(NA,7,101)for(v in 1:7){  variable = myocarde[idx,v]  v_seuil = seq(quantile(myocarde[idx,v],7/length(myocarde[idx,v])),quantile(myocarde[idx,v],1-7/length(myocarde[idx,v])), length=101)  mat_v[v,] = v_seuil  for(i in 1:101){    CLASSE = variable&#60;=v_seuil[i]    mat_gini[v,i]=      gini(y=myocarde$PRONO[idx],classe=CLASSE)}}par(mfrow=c(3,2))for(v in 2:7){  plot(mat_v[v,],mat_gini[v,],type="l",       ylim=range(mat_gini),xlab="",ylab="",       main=names(myocarde)[v])   abline(h=max(mat_gini),col="blue")}The graph is here the following,and observe that the best split is obtained using 'REPUL', with a cutoff around 1585. We check that the (relative) gain is sufficiently large, and then we go for it.And then, we consider the other leaf, and we run the same codeidx = which(myocarde$INSYS&#62;=19)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&#60;=v_seuil[i]    mat_gini[v,i]=      gini(y=myocarde$PRONO[idx],           classe=CLASSE)}}par(mfrow=c(3,2))for(v in 2:7){  plot(mat_v[v,],mat_gini[v,],type="l",       ylim=range(mat_gini),xlab="",ylab="",       main=names(myocarde)[v])   abline(h=max(mat_gini),col="blue")}Here, we should split according to 'REPUL', and the cutoff is about 1094. Here again, we have to make sure that the split is worth it. And we cut.Now we have four leaves. And we should run the same code, again. Actually, not on the very first one, which is homogenous. But we should do the same for the other three. If we do it, we can see that we cannot split them any further. Gains will not be sufficiently interesting.Now guess what... that's exactly what we have obtained with our initial codeNote that the case of categorical explanatory variables has been discussed in a previous post, a few years ago.Application on our small datasetOn our small dataset, we obtain (after changing the default values since in R, we should not have leaves with less than 10 observations... and here, the dataset is too small).tree = rpart(y ~ x1+x2,data=df, control = rpart.control(cp = 0.25,minsplit = 7))prp(tree,type=2,extra=1) u = seq(0,1,length=101)p = function(x,y){predict(tree,newdata=data.frame(x1=x,x2=y),type="prob")[,2]}v = outer(u,u,p)image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10)points(df$x1,df$x2,pch=19,cex=1.5,col="white")points(df$x1,df$x2,pch=c(1,19)[1+z],cex=1.5)contour(u,u,v,levels = .5,add=TRUE) We have a nice and simple cutWith less observations in the leaves, we can easily get a perfect model heretree = rpart(y ~ x1+x2,data=df, control = rpart.control(cp = 0.25,minsplit = 2))prp(tree,type=2,extra=1)u = seq(0,1,length=101)p = function(x,y){predict(tree,newdata=data.frame(x1=x,x2=y),type="prob")[,2]}v = outer(u,u,p)image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10)points(df$x1,df$x2,pch=19,cex=1.5,col="white")points(df$x1,df$x2,pch=c(1,19)[1+z],cex=1.5)contour(u,u,v,levels = .5,add=TRUE)Nice, isn't it? Now, just two little additional comments before growing some more trees...PruningI did not mention pruning here. Because there are two possible strategies when growing trees. Either we keep spliting, until we obtain only homogeneous leaves. Once we have a big, deep tree, we go for pruning. Or we use the stategy mentionned here : at each step, we check if the split is worth it. If not, we stop.Variable ImportanceAn interesting tool is the variable importance function. The heuristic idea is that if we use variable 'INSYS' to split, it is an important variable. And its importance is related to the gain in Gini index. If we get back to the visualization of the tree, it seems that two variables are interesting here: 'INSYS' and 'REPUL'. And we should get back to previous computation to quantify how important both are. This will be used in our next post, on random forests. But actually it is not the case here, with one single tree. Let us get back to the graph on the initial node.Indeed, 'INSYS' is important, since we decided to use it. But what about 'INCAR' or 'REPUL'? They were very close...  And actually, in R, those surrogate splits are considered in the computation, as briefly explained in the vignette. Let us look more carefully at the output of the R functioncart = rpart(PRONO~., myocarde)split = summary(cart)$splitsIf we look at the first part of that object, we getsplit      count ncat    improve    index       adjINSYS    71   -1 0.58621312   18.850 0.0000000REPUL    71    1 0.55440034 1094.500 0.0000000INCAR    71   -1 0.54257020    1.690 0.0000000PRDIA    71    1 0.27284114   17.000 0.0000000PAPUL    71    1 0.20466714   23.250 0.0000000So indeed, 'INSYS' was the most important variable, but surrogate splits can also be considered, and 'INCAR' and 'REPUL' are indeed very important. The gain was 58% (as we obtained) using 'INSYS' but there were gains of 55% (nothing to be ashamed of). So it would be unfair to claim that they have no importance, at all. And it is the same for the other leaves that we split,REPUL    27    1 0.18181818 1585.000 0.0000000PVENT    27   -1 0.10803571   14.500 0.0000000PRDIA    27    1 0.10803571   18.500 0.0000000PAPUL    27    1 0.10803571   22.500 0.0000000INCAR    27    1 0.04705882    1.195 0.0000000On the left, we did use 'REPUL' (with 18% gain), but 'PVENT', 'PRDIA' and 'PAPUL' were not that bad, with (almost) 11% gain... We can obtain variable importance by summing all those values, and we havecart$variable.importance     INSYS      REPUL      INCAR      PAPUL      PRDIA      FRCAR      PVENT 10.3649847 10.0510872  8.2121267  3.2441501  2.8276121  1.8623046  0.3373771 that we can visualize usingbarplot(t(cart$variable.importance),horiz=TRUE)To be continued with more trees...