{"word_count": 712, "character_count": 6230}
Classification trees are nice. They provide an interesting alternative to a logistic regression.  I started to include them in my courses maybe 7 or 8 years ago. The question is nice (how to get an optimal partition), the algorithmic procedure is nice (the trick of splitting according to one variable, and only one, at each node, and then to move forward, never backward), and the visual output is just perfect (with that tree structure). But the prediction can be rather poor. The performance of that algorithme can hardly compete with a (well specified) logistic regression.Then I discovered forests (see Leo Breiman's page for a detailed presentation). Being a huge fan of boostrap procedures I loved the idea. In regression models, I usually mention boostrap to avoid asymptotic approximations: we boostrap the rows (the observations). In the case of random forest, I have to admit that the idea of selecting randomly a set of possible variables at each node is very clever. The performance is much better, but interpretation is usually more difficult. And something that I love when there are a lot of covariance, the variable importance plot. Which is something that we can hardly get with econometric models (please let me know if I'm wrong).In order to illustrate, let us generate a large dataset. Not necessarily huge, but large, so that we really have to select variables.  Since it is more interesting if we have possibly correlated variables, we need a covariance matrix. There is a nice package in R to randomly generate covariance matrices.&#62; set.seed(1)&#62; n=500&#62; library(clusterGeneration)&#62; library(mnormt)&#62; S=genPositiveDefMat("eigen",dim=15)&#62; S=genPositiveDefMat("unifcorrmat",dim=15)&#62; X=rmnorm(n,varcov=S$Sigma)&#62; library(corrplot)&#62; corrplot(cor(X), order = "hclust")See Gosh &#38; Hendersen (2003) for more details on the methodology.Now that we have a  covariance matrix, let us generate a dataset,&#62; Score=-1+X[,1]/sd(X[,1])+X[,2]/sd(X[,2])++ 0.5*X[,3]/sd(X[,3])-0.25*X[,4]/sd(X[,4])&#62; P=exp(Score)/(1+exp(Score))&#62; Y=rbinom(n,size=1,prob=P)&#62; df=data.frame(Y,X)&#62; allX=paste("X",1:ncol(X),sep="")&#62; names(df)=c("Y",allX)The variable importance plot is obtained by growing some trees,&#62; require(randomForest)&#62; fit=randomForest(factor(Y)~., data=df)Then we can use simple functions&#62; (VI_F=importance(fit))    MeanDecreaseGiniX1          31.14309X2          31.78810X3          20.95285X4          13.52398X5          13.54137X6          10.53621X7          10.96553X8          15.79248X9          14.19013X10         10.02330X11         11.46241X12         11.36008X13         10.82368X14         10.17462X15         10.45530which can also be obtained using&#62; library(caret)&#62; varImp(fit)     OverallX1  31.14309X2  31.78810X3  20.95285X4  13.52398X5  13.54137X6  10.53621X7  10.96553X8  15.79248X9  14.19013X10 10.02330X11 11.46241X12 11.36008X13 10.82368X14 10.17462X15 10.45530But the popular plot that we see in all reports is usually&#62; varImpPlot(fit,type=2)which represents the mean decrease in node impurity (and not the mean decrease in accuracy). One can also visualise Partial Response Plots, as suggested in Friedman (2001), in the context of boosting,&#62; importanceOrder=order(-fit$importance)&#62; names=rownames(fit$importance)[importanceOrder][1:15]&#62; par(mfrow=c(5, 3), xpd=NA)&#62; for (name in names)+   partialPlot(fit, df, eval(name), main=name, xlab=name,ylim=c(-.2,.9))Those variable importance functions can be obtained on simple trees, not necessarily forests. As discussed in a previous post, given an impurity function  such as Gini index we split at some node  if the change in the indexis significant, where  is the node on the left, and  the node on the right.It is possible to evalute the importance of some variable when predicting  by adding up the weighted impurity decreases for all nodes  where is used (averaged over all trees in the forest, but actually, we can use it on a single tree),&#160;where the second sum is only on nodes  based on variable . If  is Gini index, then  is called Mean Decrease Gini function.Consider a single tree, just to illustrate, as suggested in some old post on http://stats.stackexchange.com/&#62; library(rpart)&#62; fit=rpart(factor(Y)~., df)&#62; plot(fit)&#62; text(fit)The idea is look at each node which variable was used to split, and to store it, and then to compute some average (see http://stats.stackexchange.com/)&#62; tmp=rownames(fit$splits)&#62; allVars=colnames(attributes(fit$terms)$factors)  &#62; rownames(fit$splits)=1:nrow(fit$splits)&#62; splits=data.frame(fit$splits)&#62; splits$var=tmp&#62; splits$type=""&#62; frame=as.data.frame(fit$frame)&#62; index=0&#62; for(i in 1:nrow(frame)){+   if(frame$var[i] != "&#60;leaf&#62;"){+   index=index + 1+   splits$type[index]="primary"+   if(frame$ncompete[i] &#62; 0){+   for(j in 1:frame$ncompete[i]){+   index=index + 1+   splits$type[index]="competing"}}+   if(frame$nsurrogate[i] &#62; 0){+   for(j in 1:frame$nsurrogate[i]){+      index=index + 1+      splits$type[index]="surrogate"}}}}&#62; splits$var=factor(as.character(splits$var))&#62; splits=subset(splits, type != "surrogate")&#62; out=aggregate(splits$improve,+     list(Variable = splits$var),+     sum, na.rm = TRUE)&#62; allVars=colnames(attributes(fit$terms)$factors)&#62;  if(!all(allVars %in% out$Variable)){+   missingVars=allVars[!(allVars %in% out$Variable)]+   zeros=data.frame(x = rep(0, length(missingVars)), Variable = missingVars)+   out=rbind(out, zeros)}&#62; out2=data.frame(Overall = out$x)&#62; rownames(out2)=out$Variable&#62; out2      OverallX1  51.024692X10  4.328443X11 19.087255X12 10.399549X13 15.248933X15  9.989834X2  68.758329X3  41.986055X4  15.211913X5  18.247668X7  18.857998X8  43.318540X9  30.299429X6   0.000000X14  0.000000This is the variable influence table we got on our original tree&#62; VI_T=out2&#62; barplot(unlist(VI_T/sum(VI_T)),names.arg=1:15)If we compare we the one on the forest, we get something rather similar&#62; barplot(t(VI_F/sum(VI_F)))This graph is a great tool for variable selection, when we have a lot of variables. And we can get it on a single tree, if it is deep enough.