Random Forests: Properties and applications

FINAL DEGREE PROJECT Random Forests: Properties and applications Presented by: Fátima del Pilar Villalba Pizarro Supervised by: DR. EMILIO CARRIZOSA PRIEGO FACULTY OF MATHEMATICS Statistics and Operational Research Department Seville, September 2020 Índice general Resumen 5 Abstract 7 Introducción 9 Introduction 13 1. From CART to random forests 15 1.1. CART (Classification and Regression Trees) . . . . . . . . . . . . . . 15 1.1.1. Structure of decision trees . . . . . . . . . . . . . . . . . . . 16 1.1.2. Construction of decision trees . . . . . . . . . . . . . . . . . 17 1.1.2.1. Topology and type of splitting . . . . . . . . . . . . 18 1.1.2.2. Spliting rules . . . . . . . . . . . . . . . . . . . . . 18 1.1.2.3. Stop-splitting rule . . . . . . . . . . . . . . . . . . 22 1.1.2.4. Designation of class label at terminal nodes . . . . 25 1.1.3. Confusion matrix and accuracy . . . . . . . . . . . . . . . . 25 1.1.4. Regression Trees . . . . . . . . . . . . . . . . . . . . . . . . 26 1.2. Randomforests............................. 26 1.2.1. Construction of random forests . . . . . . . . . . . . . . . . . 27 1.2.1.1. Bagging(Bootstrap-aggregating) . . . . . . . . . . 27 1.2.1.2. Randomness in splitting . . . . . . . . . . . . . . . 28 1.2.1.3. Unprunning . . . . . . . . . . . . . . . . . . . . . 29 1.2.2. Proximity matrix . . . . . . . . . . . . . . . . . . . . . . . . 29 1.2.2.1. Dissimilarity matrix [13][7] . . . . . . . . . . . . . 30 1.2.2.2. Outliers ....................... 31 2. Importance of variables 33 3. Applications of random forests 35 3.1. Random forests in R........................... 35 3.2. ExperimentI .............................. 37 3.3. ExperimentII.............................. 45 3 4 Índice general 3.3.1. Tecator ............................. 45 3.3.1.1. Classification of tecator . . . . . . . . . . . . . . . 45 3.3.1.2. Regression of tecator . . . . . . . . . . . . . . . . 53 3.3.2. Growth ............................. 54 4. Conclusions 61 Glossary 63 Bibliography 64 Resumen Los bosques aleatorios son una herramienta fundamental en el ámbito del aprendizaje supervisado, por lo que son empleados en multitud de disciplinas, demostrando grandes resultados y múltiples ventajas en problemas de clasificación y regresión. En este trabajo se exponen los fundamentos y bases de los bosques aleatorios, destacando sus ventajosas características, para posteriormente hacer diferentes aplicaciones con ellos, con las que contrastar dichas ventajas, observar distintos aspectos que se harán patentes e incluso analizar su comportamiento cuando extendemos su uso al caso de los datos funcionales. 5 Abstract Random forests are considered a fundamental tool in supervised learning. Consequently, random forests are used in a wide range of disciplines, yielding great results and demonstrating many advantages in classification and regression problems. Along this work, the bases and foundations of random forests are exposed, emphasizing their advantageous properties, to subsequently make different applications with them, demonstrating their advantages, analyzing distinct aspects that become noticeable, and even studying the way they behave when extending their use to functional data. 7 Introducción En un mundo en constante avance y desarrollo, es importante usar técnicas que modelen la realidad de manera eficiente. Muchos problemas de regresión y clasificación se pueden resolver mediante aprendizaje supervisado, es decir, mediante una técnica en la que, dada una muestra con distintos datos en forma de pares, uno de entrada (normalmente un vector) y otro de salida, se les atribuye un valor numérico o clase a cada par, con el fin de poder predecir la salida de futuros datos. Este método ha ido cobrando importancia en los últimos años tomando un papel relevante como puede verse en la Figura 1, que nos muestra el interés suscitado por este tema en relación con el número de búsquedas del término “aprendizaje supervisado" (en ingles “supervised learning") en Google. Figura 1: Evolución del número de búsquedas del término “aprendizaje supervisado" (supervised learning) en Google desde 2004 Una herramienta particular de aprendizaje supervisado es la conocida como bosques aleatorios, o del inglés, random forests, desarrollados en el año 2001 por Leo Breiman [6]. Esta técnica se usa en distintas áreas como son la econometría [20], la quimioinformática [18], la bioinformática [10] o la medicina. En cuanto a esta última, se tiene que los bosques aleatorios son útiles en diferentes campos como son la selección de marcadores genéticos responsables de enfermedades, la microbiología o la epidemiología genética, e incluso para predecir la replicación de un virus como el HIV-1 [17]. La relevancia de los bosques aleatorios se debe a que, a diferencia de otras técnicas, los bosques aleatorios son eficaces incluso cuando se trabaja con variables con bastante 9 16 1.1. CART (Classification and Regression Trees) Decision trees are considered leaders in their area because of their easy interpretability. The reason therefore resides in the scheme that follows such a tree: an if-then rule. [9] 1.1.1. Structure of decision trees Let us explain now how decision trees work. We will begin with classification the trees and then go on with an in depth review of the regression trees. Decision trees can handle any type of variables, categorical (i.e., they take values in a finite set not having any natural order) or numerical/continuous (i.e. it is a real number). Let us see first which is the structure of such a tree. A classification tree is a collection of nodes and branches that display a partition of the original set. Let t0be such original set. This set can be split into a compound of sets that form a partition of t0. For simplicity, let us think we are working with the simplest type of tree, the binary tree, which means at each split two new sets are obtained. If we name the new sets obtained from t0,t1and t2it is fulfilled t0=t1tt2. This process is repeated with such new sets. All these sets are called nodes. Three types of nodes are found: The original set t0, called root node. The subsets that are not split, called terminal nodes. These form a partition of t0. Each terminal node is assigned a class label, and it is possible to find more than just one terminal node with the same class label. The remaining nodes, called nonterminal nodes. Such structure of a decision tree can be seen in Figure 1.1. Figure 1.1: Strcture of a decision tree Chapter 1. From CART to random forests 17 Marital status Age Working Married 28 Yes Single 37 No Divorced 42 No Married 41 Yes Table 1.1: Example of learning sample L 1.1.2. Construction of decision trees For the construction of the classification tree, a learning sample must be given. Let L={(Xi, Yi)i=1,..,N }denote the learning sample, where Nis the number of observations, Xi∈Xa vector of Mvariables, named measurement vector, in the measurement space X, and Yithe corresponding observation from among Kpossible classes C=C1, ..., CK. Example 1.1.1. An example of what is explained above can be seen in Table 1.1, where it can be identified: N= 4,M= 2,K= 2 X1is the variable Maritalstatus that takes values in Married, Single, Divorced X2is the variable Age that takes values in 18,19,20, . . .  L=(X1 1=Married, X2 1= 28, Y1=Y es),(X1 2=Single, X2 2= 37, Y2=No),(X1 3=Divorced, X2 3= 42, Y3=No),(X1 4=Married, X2 4= 41, Y4=Y es) According to the given learning sample the classification tree can be built. The construction of a tree is based on four points: 1. Topology and type of splitting 2. Selection of splits 3. Declaration of terminal nodes/Stop-splitting rule 18 1.1. CART (Classification and Regression Trees) 4. Designation of class label at terminal nodes The main goal is constructing the classification tree from the training set L. We now outline some details of the above mentioned elements. 1.1.2.1. Topology and type of splitting According to the number of new nodes resulting from the splitting, two types of splitting can be distinguished: binary splitting, i.e., there are obtained two new nodes, or multi-splitting, that induces more than two new nodes. For continuous variables the second one is not really useful. It can also be set a difference between two types of splits according to the number of variables involved in the splitting. If only one variable takes part in the split, it is called univariate split, otherwise, it is named multivariate split. In what follows, just binary trees, with univariate split will be considered. 1.1.2.2. Spliting rules In this section we will assume we are working with a standard structure, i.e., all measurement vectors Xiare of fixed dimensionality. Many different criteria of splitting have been adopted. Out of then, Leo Breimann highlights two, namely: Gini criterion Twoing criterion The Gini criterion is the one that is going to be explained, because it is also the one that uses the program Rthat is going to be used in the last chapter. Before going down to the explanation of how this criterion operates, we must first introduce a few fundamental concepts in two initial steps, that will converge in the third and last step. 1. Explanation of concepts related with probability 2. Introduction of sets of all possible splitting 3. Explanation of selecting the best split: Gini criterion Probability concepts Beginning with the first step: Remember we said we had the training sample L= (Xi, Yi)i=1,..,N , with Xi∈Xand dimension M, and Yi∈C=C1, ..., CK. Chapter 1. From CART to random forests 19 The prior class probabilities π(k),1≤k≤Kis defined as: π(k) = P(Yi=Ck). Let Nkbe the number of cases of Lin class Ck. Often the prior probabilities can be assessed like: π(k) = Nk N,1≤k≤K. Suppose now we are at node t. Maintaining the same idea as above, let N(t)be the total number of observations in Lthat have reached the node t, and let Nk(t)denote the number of observations of class Ckin t. So, the quotient Nk(t) Nk is interpreted as the proportion of observations of class Ckfalling into t. Using both definitions, we obtain the probability of an observation in Lfulfilling two aspects simultaneously: falling into tand being of class Ck, given by p(k, t) = π(k)Nk(t) Nk . If we sum in k, we obtain the probability of an observation of any class Ckreaching the node t p(t) = X k p(k, t), and we can define now the probability the observation is of class Cksubject to having reached node t p(k|t) = p(k, t) p(t) satisfying: X k p(k|t) = 1 According to the approach done for π(k)we can also assume that p(k|t) = Nk(t) N(t).(1.1) 20 1.1. CART (Classification and Regression Trees) Figure 1.2: Parent node and descendant nodes Splitting sets Let us now continue with the second step and talk about two necessary sets Qand S. Supposed we want to do univariate binary splitting, just one single variable Xm, that matches with a coordinate of the measurement vector X, is responsible for the split. In case Xmis categorical taking values in the set B={b1, . . . , bH}, we can define a set Qlike Q=Questions |Xm i∈A?, A ⊂B and in case Xmis a continuous variable: Q=Questions |Xm i< c?. This way, every question of Qsuggests a possible split, as depicted in Figure 1.2. With regard to the observations that have reached the node t, the question of Qhas to be formulated, and there are just two possible answers: YES or NO. Depending on the answer, the observation reaches now the node to the left, let us call it tleft =tY ES, or the one to the right, namely tright =tNO. Each question generates one possible split. The set of splits is the one called S. If the variable Xm iis categorical, then there are 2H−1−1possible splits, because the number of possible splits is given by the number of possible questions, i.e., by |Q|, taking into account that tleft =tY ES and tright =tNO is symmetrical to tleft =tNO and tright =tY ES, in other words, Xm i∈Ais symmetrical to Xm i/∈A. It is known that Xm itakes values in B={b1, . . . , bH}so |B|=H, and thus the number of possible subsets of Bis 2H, but it must be considered the symmetry above explained, and it must also be considered that ∅(and by symmetry also the total) have to be factored out. So it follows that the number of questions in Q is |Q|= 2H−1−1, and therefore the number of possible splits. In case the variable Xm iis continuous, contrary to what may seem, the number of Chapter 1. From CART to random forests 21 distinct splits is also finite. There are at most as many splits as number of observations N, because this type of splits: 1. Consider all the distinct values of Xm ithat appear in L 2. Sort them 3. Take cihalfway between two of the ordered values. Gini criterion Next action is already the third step, and it consists in selecting the best split, among all the possibilities. Hence, we need to introduce two new concepts: impurity function and impurity measure. Let us begin by defining the impurity function. Definition 1.1.1. A function φ:P→R, where P is defined as P=(p1, . . . , pK)| Pkpk= 1, pk≥0, k = 1, .., K, is said to be an impurity function if it satisfies the following properties: 1. φachieves its maximum only at the point (1/K, . . . , 1/K) 2. φachieves its minimum at the points of the form (1,0,...,0),(0,1,0,...,0), ...,(0,...,0,1) 3. φis a symmetric function of p1, . . . , pK. Let us continue now explaining the impurity measure of a node t,i(t). It is defined as (remember we defined p(k|t)in (1.1) as the probability of the observation being of class Cksubject to have reached node t): i(t) = φ(p(1 |t), . . . , p(K|t)) From the above, the decrease in impurity of a split sin a node tis defined as: ∆i(s, t) = i(t)−prighti(tright)−plefti(tleft)(1.2) where pright and pleft are the proportion of items of the node tsent by the split sto the node tright and to the node tleft, respectively. Now, once the impurity concepts have been explained, according to the probability concepts exposed above, we can finally talk about the criterion we mentioned before: Gini Criterion. 22 1.1. CART (Classification and Regression Trees) This criterion uses the impurity measure of a node tdescribed by the following formula, called Gini Index: X k6=l p(k|t)p(l|t) or what is the same: i(t) = (X k p(k|t))2−X k p2(k|t)=1−X k p2(k|t).(1.3) Remember: p(k|t)refers to the probability of assigning an item selected at random from node t to class Ck, and p(l|t)expresses the probability of the object belonging to the class Cl. Using this impurity measure, it must be calculated now the decrease in impurity ∆i(s, t)using (1.2) and (1.3). So, ∆i(s, t) = 1−X k p2(k|t)−pright "1−X k p2(k|tright)#−pleft "1−X k p2(k|tleft)# ∆i(s, t) = −X k p2(k|t) + pright X k p2(k|tright) + pleft X k p2(k|tleft) The aim is to maximize ∆i(s, t). The better split is the one that achieves it. Any split sverifies ∆i(s, t)≥0. But, it can be demonstrated that ∆i(s, t)is a concave function. Hence, ∆i(s, t) = 0 if and only if p(k|tleft) = p(k|tright) = p(k|t), k = 1, . . . , K. 1.1.2.3. Stop-splitting rule Once the splitting rule has been selected, the next action is to decide the stopsplitting rule, i.e, a rule that indicates when to declare a node a terminal node. One of the initial stopping rules was based on fixing a threshold, call it β > 0, for the maximum decrease in tree impurity. Definition 1.1.2. The tree impurity I:T−→ R, where Tis a set of trees, is a function defined from the impurity measure ias: I(T) = X t∈T0 I(t) Chapter 1. From CART to random forests 23 where T0is a set of terminal nodes achieved after some splitting. Let I(t) = i(t)p(t) then the tree impurity is: I(T) = X t∈T0 I(t) = X t∈T0 i(t)p(t) So, the decrease in tree impurity is given by: ∆I(s, t) = I(t)−I(tright)−I(tleft) Hence, the early stop-splitting is: max s∈S∆I(s, t)< β Nevertheless, this rule was shown not to be fully satisfactory and instead, the tree should be pruned, once it has grown much too large. Therefore is also important to explain what "grown much too large" means. Hence, let us explain, before going on, three new concepts: misclassification cost C(k|l), resubstitution estimate of the expected missclassification cost of a node t,r(t), and the resubstitution estimate of the missclassification cost of a tree T,R(T). We are going to tackle these question from a general perspective and particularize for the specific case of classification trees. Definition 1.1.3. Let dbe a function d:X→C, called classifier. The true misclassification rate of the function d, denoted R∗(d), constructed from the learning sample L, indicates how accurate a classifier is, i.e., the probability of dmisclassifying a new sample drawn from the same distribution as L, by testing the classifier on subsequent cases whose correct classification has been observed. This function is given by: R∗(d) = P(d(X)6=Y) where X∈X, Y ∈C. This means that a good classificator has a low value of R∗(d). It follows that the prunning criterion must minimize R∗(d). The question is how can R∗(d)be estimated. One of the techniques is using the resubstitution estimate. The resubstitution estimate is the proportion of cases misclassified, and is given by: R(d) = 1 N N X i=1 X(d(Xi)6=Yi) 24 1.1. CART (Classification and Regression Trees) where Xis the indicator function, which has value equal to 1 if the argument is true and 0 if it is false. It is used Lto construct dand also in order to calculate R(d). So, if we use the value of R(d)to estimate R∗(d), the result will be unrealistic good, even being possible that the value of R(d) = 0 and so R∗(d) = 0, too. One way of enhancing the estimation, in case the learning sample Lis large enough, is to divide Linto two subsets L1and L2, so that with L1we construct dand with L2calculate R(d). The second subset L2is called test sample.L2can be considered independent of L1and from the same distribution. However for samples Lof small size, the cross-validation method can be used [16]. Now we can use what is explained above in the case of classification trees. Definition 1.1.4. C(k|l)is defined as the cost of misclassifying, as class Ck, an item that indeed corresponds to class Cl, satisfying: C(k|l)≥0if k 6=l = 0 if k =l Definition 1.1.5. The resubstitution estimate of the expected misclassification cost of a node t,r(t)is defined as: r(t) = min kX l C(k|l)p(l|t), where PlC(k|l)p(l|t)is the estimated expected misclassification cost of an unknown class item that reaches node tand is classified as class Ck. Definition 1.1.6. The resubstitution estimate of the misclassification cost of a tree T, R(T)is given by: R(T) = X t∈T0 r(t)p(t) = X t∈T0 R(t) where R(t) = r(t)p(t) The objective is to minimize the value of R(T). Therefore an equilibrium between two opposing properties must be found: 1. If the number of splits increases, then the value of R(T)decreases. In other words, if the number of terminal nodes increases, the value of R(T)decreases. 2. A tree that is grown much too large, so that just one item of the learning sample reaches that terminal node, is overfitted, i.e. it classifies perfectly the given learning sample, but it is likely it does not classify correctly a new sample. If the number of terminal nodes increases too much,R(T)increases. Chapter 1. From CART to random forests 25 Once it has been understood how to measure "how good" a tree is, we can now explain the procedure that must be follow in order to prune the tree properly. First of all, a tree must grow until all terminal nodes are pure: In this first step, the tree must grow until, for every terminal node, it is satisfied that only one object has reached the terminal node in question. Let us call this tree Tmax. Next step is pruning upward, i.e., cut off nodes successively. The final step consists of choosing from the collection of subtrees formed in the previous step, the "optimum-sized" tree. Observation 1.1.1. In fact, a value Nmin (in general it takes the values 1 or 5) can be set, and the tree grows just until N(t)⩽Nmin is fulfilled. (Remember: N(t)is the total number of observations in Lthat have reached the node t.) 1.1.2.4. Designation of class label at terminal nodes Now it is time to assign a class label at all the terminal nodes of the tree. Therefore the class assignment rule C∗ k(t)is used, where C∗ k(t)denotes that Ckis the class given to the terminal node t. The rule is: C∗ k(t) = Ckfor the k satisfying p(k|t) = max lp(l|t) If the maximum is attained at more than one value of k, any of such kcan be taken. 1.1.3. Confusion matrix and accuracy Until this moment it has been explained how to construct a classifier type called classification tree. The goodness of fit of the method is defined through a confusion matrix and the accuracy. The accuracy measures the proportion of well-classified items, so a high accuracy indicates that the classifier is good. The higher the accuracy, the better the classifier. Assume there are two possible classes: positive and negative. When predicting the class of an observation, there are four possible endings: 1. Predict class positive, being the truth class positive 2. Predict class positive, being the truth class negative 3. Predict class negative, being the truth class negative 4. Predict class negative, being the truth class positive The four cases receive following names: true positive (TP), false positive (FP), true negative (TN) and false negative (FN), respectively. Chapter 2 Importance of variables As already pointed out, random forests lose interpretability compared to decision trees, because the last ones show a direct influence of the predictor variable at its position at the tree. While the decision trees show the influence that the predictor variables have, directly on their position at the tree [17], random forests are constructed in a more complex way, and they do not reveal the dependency directly. In order to interpret random forest, one of the most important tasks in random forests are to be done: to measure the importance of the variables[3]. There are different ways of addressing this issue. 1. One of them, and probably the simplest one, consists in just counting how many times each variable is selected by the trees that conform the forest[17]. 2. Another, more complex, measure of the importance of variables is the Mean Decrease Impurity(MDI). In case the selected splitting criterion is the Gini criterion, than the Mean Decrease Impurity is called Gini importance[17]. This measure of variable importance is based on the following assumption: the Gini Index for a parent node has a higher value than the one of its descendants [13]. As we explained previously, when doing a split, the "best" variable for the split must be selected according to the decrease in impurity measure. So averaging the decrease in impurity over all the nodes of the trees in the forest, where the variable is the one selected for the splitting, the MDI is obtained. So, the higher the value of the MDI, the more important the variable [15]. 3. A further way is the permutation of the values of the variable in consideration, called Mean Decrease Accuracy (MDA) [13]. It consists in permuting for i= 1, ..., N the predictor variable Xm i(m−th coordinate of the measurement vector Xi), dissociating this variable from its original observation Yi. 33 34 Marital status Working Married Yes Single No Divorced No Married Yes Table 2.1: Zoom on variable X1 iMarital status of Table 1.1 Marital status Working Single Yes Married No Married No Divorced Yes Table 2.2: Permutation of the variable Marital status up to Table 2.1 For a better understanding of the permutation, let us go back to the example of Table 1.1. The variable X1 iwas for the different i0s X1 1=Married, X1 2= Single, X1 3=Divorced, X1 4=Married, and the correspondent observations Y1=Y es, Y2=Y es, Y3=No, Y4=Y es. Once the permutation of the variable is done, the appearance of Table 2.1 becomes, for instance, the one in Table 2.2. If now the permuted variable, together with the non permuted variables, is used to predict the response, there are two possible scenarios. The first consists in anything changing with respect to the original scenario, where no permutation had be done. This would mean that the permutation of the variable has not influenced on the results, so this variable is not really important for the prediction. The second possible scenario is that after the permutation the number of observations misclassified increases, i.e. the accuracy decreases. The only justification for this outcome is that the variable has a significant importance for the prediction [17]. (It is worth noting that the accuracy is obtained by classifying the OOB sample, so the permutation of the variable takes place just in the OOB) [19]. If the decrease in accuracy, by permuting the variable, is averaged over all the trees in the forest, the MDA is obtained. So, the higher the MDA, the more important is the variable, because more decreases the accuracy if we permute the variable [15]. Chapter 3 Applications of random forests Chapters 1 and 2 have been focused on the theory of random forests. Now, in this chapter, some applications of the above explained theory is shown, running different experiments with the software R. Therefore, in Section 3.1 the main ideas of computing random forests in Rare introduced, in order to subsequently use this software for different applications with the data set Wisconsin breast cancer diagnosis in Section 3.2, and with two different data sets with functional data, namely, tecator and growth, in Section 3.3. 3.1. Random forests in R Before presenting the different experiments, how random forests are implemented in Ris outlined in this section. First step is to install the package related with this item and then call such library: install.packages("randomForest") library(randomForest) Once this has been done, it is possible to use the functions of such package with the data available. Hence, it is important to understand correctly how these functions work. Let us start focusing on the function randomForest. The help-environment of Rprovides following: ## S3 method for class ’formula’ randomForest(formula, data=NULL, ..., subset, na.action=na.fail) ## Default S3 method: 35 36 3.1. Random forests in R randomForest(x, y=NULL, xtest=NULL, ytest=NULL, ntree=500, mtry=if (!is.null(y) && !is.factor(y)) max(floor(ncol(x)/3), 1) else floor(sqrt(ncol(x))), replace=TRUE, classwt=NULL, cutoff, strata, sampsize = if (replace) nrow(x) else ceiling(.632*nrow(x)), nodesize = if (!is.null(y) && !is.factor(y)) 5 else 1, maxnodes = NULL, importance=FALSE, localImp=FALSE, nPerm=1, proximity, oob.prox=proximity, norm.votes=TRUE, do.trace=FALSE, keep.forest=!is.null(y) && is.null(xtest), corr.bias=FALSE, keep.inbag=FALSE, ...) As most of the arguments of the function randomForest are explained by themselves, the most relevant ones are presented in Table 3.1. Argument Meaning and working formula data frame or matrix of predictors, or formula describing the model to be fitted data data frame containing the variables ntree Number of trees to grow mtry Number of variables randomly sampled as candidates at each split. By default for classification it takes p(p), for regression p 3, with pthe number of variables importance if it is set TRUE it calculates the importance of the predictors proximity if it is set TRUE it calculates the proximity matrix Table 3.1: Meaning and working of the most important arguments of randomForest Once the main points of the computing have been presented, two different experiments are going to be exposed. For the first one, a data set conformed by continuous predictor variables, and a categorical response variable is going to be analyzed, focusing on the variable importance. For the second experiment the use of random forests is going to be extended to functional data. Chapter 3. Applications of random forests 37 3.2. Experiment I The data that is used along this application of random forests, has been obtained from the UCI Machine Learning Repository [11], and corresponds to the Wisconsin breast cancer diagnosis. This data is conformed by 31 predictor variables, and a categorical response variable. From among the predictors variables we distinguish the ID number (set in column one) and the ten principal ones (set from column 3 to 12), that represent the means: 1)ID number 3) radius (mean of distances from center to points on the perimeter) 4) texture (standard deviation of gray-scale values) 5) perimeter 6) area 7) smoothness (local variation in radius lengths) 8) compactness (perimeter^2 / area - 1.0) 9) concavity (severity of concave portions of the contour) 10) concave points (number of concave portions of the contour) 11) symmetry 12) fractal dimension ("coastline approximation" - 1) The remaining variables (13 to 30) are the standard error and the worst case of the 10 already exposed predictor variables. The response variable is: 2)Diagnosis (M = malignant, B = benign) The first step consists in reading, and setting a headline to the data. headline=c("Id","Classification","Radius","Texture", "Perimeter","Area","Smoothness","Compactness", "Concavity","Concave Points","Symmetry", "Fractal dimension","Radius_se","Texture_se", "Perimeter_se","Area_se","Smoothness_se", "Compactness_se","Concavity_se","Concave Points_se", "Symmetry_se","Fractal dimension_se","Radius_worst","Texture_worst", "Perimeter_worst","Area_worst", 38 3.2. Experiment I "Smoothness_worst","Compactness_worst", "Concavity_worst","Concave Points_worst", "Symmetry_worst","Fractal dimension_worst") wisconsin=read.table("wdbc.data",sep = ",", col.names = headline) and in order to be able to call the variables by their names: attach(wisconsin) Now the data is correctly implemented and can be used for the experiment. In order to allow the replication of this experiment, it is important to set a seed, and "fix" the randomness. Next, the randomForest is executed like it can be seen, using five variables (√31) at each split and growing 1000 trees, and the results are printed: set.seed(1081) w.rf=randomForest(Classification~.,data=wisconsin, ntree=1000,mtry=5,importance=TRUE) print(w.rf) So, it is obtained: Call: randomForest(formula = Classification ~ ., data = wisconsin,ntree = 1000,mtry = 5, importance = TRUE) Type of random forest: classification Number of trees: 1000 No. of variables tried at each split: 5 OOB estimate of error rate: 3.34% Confusion matrix: B M class.error B 350 7 0.01960784 M 12 200 0.05660377 that can be interpreted as follows. The OOB error is 3.34 %, i.e. this is the percentage of items of the OOB that have been wrongly classified, or, in other words, the random forest classifies 96.66%correctly. Also the confusion matrix is represented, so it can Chapter 3. Applications of random forests 39 be seen how many items have been correctly classified or not, i.e., the false positives (5.66%) and the false negatives (1.96%). Above can be seen, that in order to obtain information about the importance of the variables, the entry importance=TRUE has also been included as an argument of the function randomForest. Just using now importance(w.rf) a list with the MDA and the MDI for each variable of the model is obtained: MeanDecreaseAccuracy MeanDecreaseGini Id 4.918394 1.3460295 Radius 13.210533 10.2696696 Texture 17.081849 4.0069453 Perimeter 14.565036 14.0129864 Area 14.726296 10.8019264 Smoothness 10.008769 1.7488165 Compactness 8.583551 2.3981941 Concavity 16.030364 10.3178601 Concave.Points 21.449959 29.4106724 Symmetry 5.431248 1.1330892 Fractal.dimension 5.452440 0.9137801 Radius_se 14.810328 4.0634558 Texture_se 5.709398 1.2645055 Perimeter_se 13.998440 4.2729203 Area_se 19.755545 8.5036794 Smoothness_se 4.534937 1.1640510 Compactness_se 8.130864 1.3049383 Concavity_se 8.041339 1.6496968 Concave.Points_se 7.312380 1.1211280 Symmetry_se 5.618074 1.0870049 Fractal.dimension_se 4.237092 1.3857821 Radius_worst 24.541105 31.6980086 Texture_worst 19.268786 5.0469406 Perimeter_worst 23.814279 32.1185216 Area_worst 25.123753 29.3425631 Smoothness_worst 16.132225 3.2753597 Compactness_worst 11.639484 4.2820763 Concavity_worst 19.114623 8.8452845 Concave.Points_worst 25.047172 34.6431429 40 3.2. Experiment I Symmetry_worst 9.118319 2.2282686 Fractal.dimension_worst 8.311042 1.8237781 To interpret this more easily, let us show it graphically in Figure 3.1 from: varImpPlot(w.rf,n.var = 31) Figure 3.1: Importance of variables As explained in Chapter 2, the higher the MDA (or the MDI), the more important the variable (according to the respective criteria). Although the results of the two meth- Chapter 3. Applications of random forests 41 ods are, in general, not exactly the same, both will identify very important or very superfluous variables. In this case, according to the MDA the most important variables, in descendant order of importance are: Concave.Points_worst, Area_worst, Radius_worst and Perimeter _worst. For MDI criterion these four are also the most important variables, permuting the order of Concave.Points_worstand Perimeter_worst. When considering variables that are not important for the random forest, it stands out that the MDA considers as completely superfluous variables less variables than MDI. For MDA we have, if we consider for example MDA under 10 as unimportant, in increasing order of importance: Smoothness_se, Id, Symmetry, Fractal dimension, Symmetry_se, Texture_se and a series of not that important variables as can be seen in the graphic. Nevertheless, according to MDI, there are a lot of superfluous variables, it suffices to take a look at the graphic and see how many variables have a value MDI=5 (not even under 10 as was done with MDA), outstanding Fractal.dimension,Symmetry_se, Concave.Points_se and Symmetry among others. As might be expected, in both cases, Id is considered an unimportant variable, because a number of identification assigned to a patient cannot be a reason to detect cancer. Figure 3.2: Ranking of the variables according to MDA and MDI From the above information, it can be deduced that there exist some similarities between the ranking of the variables, when listing them according to the MDA or the MDI. So if we take a look at Figure 3.2 the relations between the position in variable importance that the different variables take, in concordance to both criteria, can be observed. Every variable is given the position they take in the ranking of variable 48 3.3. Experiment II Vari1 Vari2 Vari3 Vari4 Vari5 Vari6 Vari7 Vari8 Vari9 Vari10 Vari11 Vari12 Vari13 Vari14 Vari15 Vari16 Vari17 Vari18 Vari19 Vari20 Vari21 Vari22 Vari23 Vari24 Vari25 Vari26 Vari27 Vari28 Vari29 Vari30 Vari31 Vari32 Vari33 Vari34 Vari35 Vari36 Vari37 Vari38 Vari39 Vari40 Vari41 Vari42 Vari43 Vari44 Vari45 Vari46 Vari47 Vari48 Vari49 Vari50 Vari51 Vari52 Vari53 Vari54 Vari55 Vari56 Vari57 Vari58 Vari59 Vari60 Vari61 Vari62 Vari63 Vari64 Vari65 Vari66 Vari67 Vari68 Vari69 Vari70 Vari71 Vari72 Vari73 Vari74 Vari75 Vari76 Vari77 Vari78 Vari79 Vari80 Vari81 Vari82 Vari83 Vari84 Vari85 Vari86 Vari87 Vari88 Vari89 Vari90 Vari91 Vari92 Vari93 Vari94 Vari95 Vari96 Vari97 Vari98 Vari99 Vari100 0 5 10 15 value variables type MDA MDI Figure 3.5: Variable Importance Chapter 3. Applications of random forests 49 (a) Curves of first derivative as functions of the wavelength (b) Curves of second derivative as functions of the wavelength Figure 3.6: First two derivatives of the absorbances with color as class-distinction The first derivative is built as follows: Deriv=matrix(,nrow=nrow(Explicativ), ncol=ncol(Explicativ)-1) for (i in 1:nrow(Explicativ)){ for (j in 1:ncol(Explicativ)-1){ Deriv[i,j]=Explicativ[i,j+1]-Explicativ[i,j] }} So now it can be computed the random forest with the first derivative, the same way as it was done with the original function just replacing Explicativ for Deriv and naming also the new variables. DataDeriv1=data.frame(Deriv,Response) set.seed(1081) rfD1=randomForest(Res~.,data=DataDeriv1,ntree=1000, importance=TRUE,proximity=TRUE) print(rfD1) In this case a OOB error of 3,26%is obtained, i.e. it has been reduced considerably. When introducing the second derivative, i.e., using: Deriv2=matrix(,nrow=nrow(Deriv), ncol=ncol(Deriv)-1) for (i in 1:nrow(Deriv)){ for (j in 1:ncol(Deriv)-1){ Deriv2[i,j]=Deriv[i,j+1]-Deriv[i,j] }} 50 3.3. Experiment II NamesDeri2=c("Deri1",...,"Res") DataDeriv2=data.frame(Deriv2,Response) colnames(DataDeriv2)<-NamesDeri2 attach(DataDeriv2) set.seed(1081) rfD2=randomForest(Res~.,data=DataDeriv2,ntree=1000, importance=TRUE,proximity=TRUE) print(rfD2) the OOB error is reduced until 0,93%. If the results of the other combinations of function and derivatives before mentioned are to be calculated, the following must be computed: Namestogether1=c("Vari1",...,"Deri1",...,"Res") VariablesTogether1=cbind(Explicativ,Deriv) DatosTogether1=data.frame(VariablesTogether1,Response) colnames(DatosTogether1)<-Namestogether1 set.seed(1081) rfT1=randomForest(Res~.,data=DatosTogether1,ntree=1000, importance=TRUE,proximity=TRUE) print(rfT1) IT1=importance(rfT1) varImpPlot(rfT1) ################# Namestogether=c("Vari1",...,"2Deri1",..., "Res") VariablesTogether2=cbind(Explicativ,Deriv2) DatosTogether2=data.frame(VariablesTogether2,Response) colnames(DatosTogether2)<-Namestogether2 attach(DatosTogether2) set.seed(1081) rfT2=randomForest(Res~.,data=DatosTogether2,ntree=1000, importance=TRUE,proximity=TRUE) print(rfT2) IT2=importance(rfT2) varImpPlot(rfT2) ################## Chapter 3. Applications of random forests 51 Namestogether=c("Vari1",...,"Deri1",...,"2Deri1",..., "Res") VariablesTogether=cbind(Explicativ,Deriv,Deriv2) DatosTogether=data.frame(VariablesTogether,Response) colnames(DatosTogether)<-Namestogether attach(DatosTogether) set.seed(1081) rfT=randomForest(Res~.,data=DatosTogether,ntree=1000, importance=TRUE,proximity=TRUE) print(rfT) IT=importance(rfT) varImpPlot(rfT) ################## NamestogetherD=c("Deri1",...,"2Deri1",...,"Res") VariablesTogetherD=cbind(Deriv,Deriv2) DatosTogetherD=data.frame(VariablesTogetherD,Response) colnames(DatosTogetherD)<-NamestogetherD attach(DatosTogetherD) set.seed(1081) rfTD=randomForest(Res~.,data=DatosTogetherD,ntree=1000, importance=TRUE,proximity=TRUE) print(rfTD) ITD=importance(rfTD) varImpPlot(rfTD) Let us now see the OOB errors for the different combinations in Table 3.2. Functions used for the explicative variables OOB error(%) f 16,28 D1 3,26 D2 0,93 f+D1+D2 0,93 f+D1 2,79 f+D2 0,93 Table 3.2: OOB errors for the different combinations in the classification problem of tecator 52 3.3. Experiment II The result obtained just with the second derivative (D2) is exactly the same result that is obtained when using the second derivative and the function with (f+D1+D2) or without (f+D2) the first derivative, i.e. the variables of the second derivative are the ones that better the results. When using the original function and the first derivative (f+D1) the results are much better than using only the original function (f) and a few better than using only the first derivative (D1), but not as good as when using the second derivative (D2), i.e. the first derivative betters the results obtained with the original function but not as much as the second derivative does. So, answering the question posed above, for this data the use of the derivatives is a clear way of improving the OOB. (a) MDA for original function, first derivative and second derivative (b) MDI for original function, first derivative and second derivative Figure 3.7: Variable importance different derivatives Another question that can be tried to be answered is if the variables that are impor- Chapter 3. Applications of random forests 53 tant for the original function, the first derivative and the second derivative correspond to the same wavelength. When representing the MDA and MDI of the random forests obtained up to the original function, the first derivative (D1) and the second derivative (D2), approximately the variables corresponding to the same wavelengths are important, as can be seen in Figures3.7(a) and 3.7(b). Nevertheless, it is important to mention that, when speaking about a range, what has been just said is correctly, but when speaking about the values inside that range it follows that not exactly the same wavelengths are important for the distinct functions. Taking a look at the range between 900nm and 950nm approximately, it stands out that this range contains the most important wavelengths, independent of speaking of the original function, the first or the second derivative. Nevertheless, inside the range it must be emphasized that the original function and the second derivative have as most important wavelengths more or less the same values around 930nm, but at this point, the first derivative achieves a minimum. It follows that the first derivative is complementary to the original function and the second derivative. 3.3.1.2. Regression of tecator In the beginning it was said that the data set tecator originally consists of a compound of observations with the variable response fat index, a continuous variable, that we decided to transform into a qualitative variable grouping the fat indexes in two classes: high and low fat index, in order to address this problem as a classification problem. Nevertheless, it is also possible to use random forest techniques when using fat index as a continuous variable, considering this time the problem as a regression problem. Proceding this way, the only difference with respect to the classification problem (when computing it with R) is to keep the original response variable like we can see for the case with the original function (For the rest of the programming just keep on (with Response coded like here) the same way as was done with the classification problem): Explicativ=tecator$absorp.fdata$data Response=tecator$y$Fat Data=data.frame(Explicativ,Response) Let us see now what happens when using the different combinations of original function and the first two derivatives as was done with the classification problem. Therefore, see Table 3.3 (The results are obtained with seed 1081). In relation with what can been observed in the Table 3.3 it can be said that when using the original function and the two derivatives (f+D1+D2) the best result is obtained. Nevertheless, the results are very similar to the ones of just using the second derivative or he second derivative and the original function (D2 or f+D2). When just 54 3.3. Experiment II Functions used for the explicative variables Variable explained(%) Mean squared residuals f 70,41 47,80386 D1 98,68 2,135725 D2 99,05 1,526871 f+D1+D2 99,32 1,093536 f+D1 98,68 2,134657 f+D2 99,01 1,595362 Table 3.3: Different combinations in the regression problem of tecator using the first derivative (D1) or the first derivative and the original function (f+D1) the same result is obtained for the percentage of the explained variability of the response variable, probably because the best results when splitting are obtained for the variables of the first derivative (in contrast to the variables of the original function), so using for the random forest just the first derivative (D1) or its combination with the original function (f+D1) does not present any difference. The worst result is obtained when using just the original function (f). 3.3.2. Growth The dataset called growth, available with the library(fda) of Rhas also been addressed as a classification problem. This data set represents the function of the height of 39 boys and 54 girls between the ages of 1 and 18, discretized into 31 points, not equally spaced (1, 1.25, 1.5, 1.75, 2, 3, 4, 5, 6, 7, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 14, 14.5, 15, 15.5, 16, 16.5, 17, 17.5, 18). So this way for each individual ithere are 31 variables Xm i, with m= 1, ..., 31 that can be used in order to classify with random forests the response variable, the sex of the individual. The different curves of height for the different ages are drawn in Figure 3.8, in different colors depending if it corresponds to a boy or a girl. As done with the tecator data set, the two classes of the response variable of the growth data are going to be coded: • "1" if it is a boy • "-1" if it is a girl So following must be computed: Chapter 3. Applications of random forests 55 Figure 3.8: Height at the ages between 1 and 18 library(fda.usc) library(randomForest) #################################################### data("growth") names=c("Vari1","Vari2","Vari3","Vari4","Vari5","Vari6", "Vari7","Vari8","Vari9","Vari10","Vari11", "Vari12","Vari13","Vari14","Vari15","Vari16", "Vari17","Vari18","Vari19","Vari20","Vari21", "Vari22","Vari23","Vari24","Vari25","Vari26", "Vari27","Vari28","Vari29","Vari30","Vari31", "Res") boys=rep(1,39) girls=rep(-1,54) Explicativ=rbind(t(growth$hgtm),t(growth$hgtf)) Response1=c(boys,girls) Response=factor(Response1) Data=data.frame(Explicativ,Response) colnames(Data)<-names attach(Data) So, when computing the random forest with the original function 56 3.3. Experiment II set.seed(33) rf=randomForest(Res~.,data=Data,ntree=1000, importance=TRUE,proximity=TRUE) print(rf) that will have as output: Type of random forest: classification Number of trees: 1000 No. of variables tried at each split: 5 OOB estimate of error rate: 8.6% Confusion matrix: -1 1 class.error -1 48 6 0.11111111 1 2 37 0.05128205 So, the OOB error obtained is of 8,6%that we will try to improve as we did with the tecator data set. The confusion matrix is also shown in the output, with 48 "true positive", i.e. in this case, 48 girls classified as girls, and 37 "true negative", i.e. 37 boys classified by the random forest as boys. Figure 3.9: Variable Importance Let us now see the importance of variables. According to Figure 3.8, it might be expected that the most important variables are those of ages round 15, because until Chapter 3. Applications of random forests 57 that moment, the curves of boys and girls are very similar, and only up to the mentioned age, the curves are different for most of the individuals of both classes. So, taking a look at Figure 3.9, it can been observed that, the most important variables, according to MDI and MDA are Vari27 until Vari31, which correspond to the ages 16 to 18, i.e. the result obtained is coherent with what was expected. Next step is to see what happens when introducing the first and the second derivatives. Therefore, let us see first the curves of these functions in Figures 3.10(a) and 3.10(b). From Figure 3.10(a) can it can be deduced that the most important ages in order to distinguish if the person is a boy or a girl, are approximately the ones over 13 years, because up to that moment, the curves of boys and girls are different. When focusing on the second derivative in Figure 3.10(b) it cannot be done such a deduction. (a) Curves of first derivative of the height as functions of the age (b) Curves of second derivative of the height as functions of the age Figure 3.10: First two derivatives of the height with color as class-distinction Taking a look at Figures3.11(a) and 3.11(b) it can be seen now which are the most important variables for the random forests when using the first and the second derivatives. Figure 3.11(a) shows that the variables Deri19 until Deri30 are the most important ones, i.e. the ages 12.25 until 17.75 as was predicted by Figure 3.10(a). Nevertheless, when looking at the important ages according to the second derivative different "valleys" and "mountains" can be identified, standing out the mountain between Deri17 and Deri20, and the one between Deri22 and Deri27, corresponding to the ages between 11.5and 13, and between 14 and 16.5, respectively. It can also be seen that in the corresponding Figure for the first as also for the second derivative, some MDA values are negative. This means that when permuting the variable values over the different observations the accuracy of the random forests for new data is improved. Nevertheless, in this case we are working with a small database and additionally these values of MDA are very low in contrast to the variables that are considered as important, so they can be considered as zero. 64 a node t TTree T0Set of terminal nodes I(T)Tree Impurity R∗(d)True misclassification rate of the function d R(d)Resubstitution estimate C(k|l)Cost of misclassifying r(t)Resubstitution estimate of the expected misclassification cost of the node t R(T)Resubstitution estimate of the tree T C∗ k(t)Class assignment rule FNumber of trees of the forest Prox Proximity matrix Prox(i, j)Cell (i,j) of the proximity matrix ui,k Measure of outlyingness ˜ui,k Normalized measure of outlyingness mkMean of measure of outlyingness over all elements of class Ck JkNumber of observations of class Ck Bibliography [1] Nelson Lee Afanador, Agnieszka Smolinska, Thanh N. Tran, and Lionel Blanchet. Unsupervised random forest: a tutorial with case studies. Chemometrics, 30:232–241, 2016. [2] Philippe Saint-Pierre Baptiste Gregorutti, Bertrand Michel. Correlation and variable importance in random forests. Stat Comput, 27:659–678, 2017. [3] Gérard Biau and Erwan Scornet. A random forest guided tour. Test, 25(2):197– 227, 2016. [4] R. Blanquero, E. Carrizosa, A. Jiménez-Cordero, and B. Martín-Barragán. Functional-bandwidth kernel for support vector machine with functional data: An alternating optimization algorithm. European Journal of Operational Research, 275(1):195 – 207, 2019. [5] Leo Breiman. Random forests. Machine Learning, 45. [6] Leo Breiman. Statistical modeling: The two cultures. Statistical Science, 16(3):199–215, 2001. [7] Leo Breiman. Manual on setting up, using, and understanding random forests v4.0. Machine Learning. Journal, 2003. [8] Leo Breiman, Jerome H. Friedman, Richard A. Olshen, and Charles J. Stone. Classification And Regression Trees. CHAPMAN HALL/CRC, 1984. [9] Emilio Carrizosa and Dolores Romero Morales. Supervised classification and mathematical optimization. Computers Operations Research,40(1),150-165, 2013. [10] Ramón Díaz Uriarte and Sara Alvarez de Andrés. Gene selection and classification of microarray data using random forest. BMC Bioinformatics, 7:3, 2006. [11] Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. [12] F. Hillier and G. Lieberman. Introduccion a la investigacion de operaciones. Mexico, D.F: McGraw-Hill/Interamericana, 2010. 65 66 Bibliography [13] Alan Julian Izenman. Modern Multivariate Statistical Techniques. Regression,Classification, and Manifold Learning. Springer, 2013. [14] Robert C. Messenger James N. Morgan. THAID,a Sequential Analysis Program for the Analysis of Nominal Scale Dependent Variables. Survey Research Center, Institute for Social Research, University of Michigan, 1973. [15] M.R. Ortiz-Posadas. Pattern Recognition Techniques Applied to Biomedical Problems. Springer. [16] Martha E. Stone. Cross-validation: A review. Math. Operationforsch. Statist. Ser Statist, 9:127–139, 1977. [17] Carolin Strobl, Anne-Laure Boulesteix, Achim Zeileis, and Torsten Hothorn. Bias in random forest variable importance measures: Illustrations, sources and a solution.”. BMC bioinformatics, 8(1):25, 2007. [18] Vladimir Svetnik, Andy Liaw, Christopher Tong, J. Christopher Culberson, Robert P. Sheridan, and Bradley P. Feuston. Random forest: a classification and regression tool for compound classification and qsar modeling. Journal of Chemical Information and Computer Sciences, 43(6):1947–1958, 2003. PMID: 14632445. [19] Jerome Friedman Trevor Hastie, Robert Tibshirani. The Elements of Statistical Learning: Data Mining, Inference and Prediction. Springer, 2009. [20] Hal R. Varian. Big data: New tricks for econometrics. Journal of Economic Perspectives, 28(2):3–28, 2014.