Classification from scratch, bagging and forests 10/8

Tenth post of our series on classification from scratch. Today, we’ll see the heuristics of the algorithm inside bagging techniques.

Often, bagging is associated with trees, to generate forests. But actually, it is possible using bagging for any kind of model. Recall that bagging means “boostrap aggregation”. So, consider a model $m:\mathcal{X}\rightarrow \mathcal{Y}$. Let $\widehat{m}_{S}$ denote the estimator of $m$ obtained from sample $S=\{y_i,\mathbf{x}_i\}$ with $i=\{1,\cdots,n\}$.

Consider now some boostrap sample, $S_b=\{y_i,\mathbf{x}_i\}$ with $i$ is randomly drawn from $\{1,\cdots,n\}$ (with replacement). Based on that sample, estimate $\widehat{m}_{S_b}$. Then draw many samples, and consider the agregation of the estimators obtained, using either a majority rule, or using the average of probabilities (if a probabilist model was considered). Hence$$\widehat{m}^{bag}(\mathbf{x})=\frac{1}{B}\sum_{b=1}^B \widehat{m}_{S_b}(\mathbf{x})$$

Bagging logistic regression #1

Consider the case of the logistic regression. To generate a bootstrap sample, it is natural to use the technique describe above. I.e. draw pairs $(y_i,\mathbf{x}_i)$ randomly, uniformly (with probability $1/n$) with replacement. Consider here the small dataset, just to visualize. For the b part of bagging, use the following code

L_logit = list() n = nrow(df) for(s in 1:1000){ df_s = df[sample(1:n,size=n,replace=TRUE),] L_logit[[s]] = glm(y~., df_s, family=binomial)}

Then we should aggregate over the 1000 models, to get the agg part of bagging,

p = function(x){ nd=data.frame(x1=x[1], x2=x[2]) unlist(lapply(1:1000,function(z) predict(L_logit[[z]],newdata=nd,type="response")))}

We now have a prediction for any new observation

vu = seq(0,1,length=101) vv = outer(vu,vu,Vectorize(function(x,y) mean(p(c(x,y))))) image(vu,vu,vv,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+(df$y=="1")],cex=1.5) contour(vu,vu,vv,levels = .5,add=TRUE) Bagging logistic regression #2 Another technique that can be used to generate a bootstrap sample is to keep all $\mathbf{x}_i$‘s, but for each of them, to draw (randomly) a value for $y$, with$$Y_{i,b}\sim\mathcal{B}(\widehat{m}_{S}(\mathbf{x}_i))$$since$$\widehat{m}(\mathbf{x})=\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}].$$Thus, the code for the b part of bagging algorithm is now L_logit = list() n = nrow(df) reg = glm(y~x1+x2, df, family=binomial) for(s in 1:100){ df_s = df df_s$y = factor(rbinom(n,size=1,prob=predict(reg,type="response")),labels=0:1) L_logit[[s]] = glm(y~., df_s, family=binomial) }

The agg part of bagging algorithm remains unchanged. Here we obtain

vu = seq(0,1,length=101) vv = outer(vu,vu,Vectorize(function(x,y) mean(p(c(x,y))))) image(vu,vu,vv,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+(df$y=="1")],cex=1.5) contour(vu,vu,vv,levels = .5,add=TRUE) Of course, we can use that code we check the prediction obtain on the observations we have in our sample. Just to change, consider here the myocarde data. The entiere code is here L_logit = list() reg = glm(as.factor(PRONO)~., myocarde, family=binomial) for(s in 1:1000){ myocarde_s = myocarde myocarde_s$PRONO = 1*rbinom(n,size=1,prob=predict(reg,type="response")) L_logit[[s]] = glm(as.factor(PRONO)~., myocarde_s, family=binomial) } p = function(x){ nd=data.frame(FRCAR=x[1], INCAR=x[2], INSYS=x[3], PRDIA=x[4], PAPUL=x[4], PVENT=x[5], REPUL=x[6]) unlist(lapply(1:1000,function(z) predict(L_logit[[z]],newdata=nd,type="response")))}

For the first observation, with our 1000 simulated datasets, and our 1000 models, we obtained the following estimation for the probability to die.

Fronm bags to forest

Here, we grew a lot of trees, but it is not stricto sensus a random forest algorithm, as introduced in 1995, in Random decision forests. Actually, the difference is in the creation of decision trees. To understand what happens, get back to the previous post on classification trees. As we’ve seen, when we have a node, we look at possible splits : we consider all possible variable, and all possible threshold. The startegy here will be to draw randomly $k$ variables out of $p$ (with of course $k, for instance $k=\sqrt{p}$). That's interesting in high dimension, because at each split, we should look for all variables, all cutoffs, and that can take quite some time (especially with the bootstrap procedure, where the goal will be to grow 1000 trees).

To be continued…

Classification from scratch, trees 9/8

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 Tree

Decision trees are easy to read. So easy to read that they are everywhere

We 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 dataset

library(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 $x$-axis and the $y$-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. $\{y_i,\mathbf{x})i\})$ for some $i\in\mathcal{I}\subset\{1,\cdots,n\}$. Hence, a leaf a caracterized by $\mathcal{I}$. For instance, the first node in the tree is $\mathcal{I}=\{1,\cdots,n\}$. A (binary) split is based on one specific variable – say $x_j$ – and a cutoff, say $s$. Then, there are two options:

• either $x_{i,j}\leq s$, then observation $i$ goes on the left, in $\mathcal{I}_L$
• or $x_{i,j}> s$, then observation $i$ goes on the right, in $\mathcal{I}_R$

Thus, $\mathcal{I}=\mathcal{I}_L\cup\mathcal{I}_R$.

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 $\mathcal{I}$ is defined as $$G(\mathcal{I})=-\sum_{y\in\{0,1\}}p_y(1-p_y)$$where $p_y$ is the proportion of individuals in the leaf of type $y$. 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 $p_y(1-p_y)$? 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 $p_0\approx1$ or $p_1\approx1$, we have strong homogenity. If we want an index to maximize, $-p_y(1-p_y)$ might be an interesting candidate. Further more, the worst case would be a leaf with $p_0\approx1/2$, which is exactly what we have here. Note that we can also write$$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)$$where $n_{y,\mathcal{I}}$ is the number of individuals of type $y$ in the leaf $\mathcal{I}$, and $n_{\mathcal{I}}$ is the number of individuals in the leaf $\mathcal{I}$.

If we do not split, we have index$$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)$$while if we split, define index$$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)$$The code to compute is would be

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)}

Actually, one can consider other indices, like the entropic measure$$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)$$while if we split, $$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)$$

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 variable

CLASSE = myocarde[,1] &lt;=100 table(CLASSE) CLASSE FALSE 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.4640415 Initially, without any split, it was -2*mean(myocarde$PRONO)*(1-mean(myocarde$PRONO)) [1] -0.4832375 which can actually also be obtained with CLASSE = myocarde[,1] gini(y=myocarde$PRONO,classe=CLASSE) [1] -0.4832375

There is a net gain in spliting of

gini(y=myocarde$PRONO,classe=(myocarde[,1]&lt;=100))- gini(y=myocarde$PRONO,classe=(myocarde[,1]&lt;=Inf)) [1] 0.01919591

Now, 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, using

sort(unique(myocarde[,1]))

But in massive datasets, it can be very long. Here, I prefer

seq(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&lt;=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 get 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, 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 is gini(y=myocarde$PRONO,classe=(myocarde[,3]&lt;19))- gini(y=myocarde$PRONO,classe=(myocarde[,3]&lt;=Inf)) [1] 0.2832801 which 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]&lt;19))- gini(y=myocarde$PRONO,classe=(myocarde[,3]&lt;=Inf)))/ gini(y=myocarde$PRONO,classe=(myocarde[,3]&lt;=Inf)) [1] 0.5862131

If 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&lt;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&lt;=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 code

idx = which(myocarde$INSYS&gt;=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&lt;=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 code

Note that the case of categorical explanatory variables has been discussed in a previous post, a few years ago.

Application on our small dataset

On 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 cut

With less observations in the leaves, we can easily get a perfect model here

tree = 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…

Pruning

I 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 Importance

An 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 function

cart = rpart(PRONO~., myocarde) split = summary(cart)$splits If we look at the first part of that object, we get split count ncat improve index adj INSYS 71 -1 0.58621312 18.850 0.0000000 REPUL 71 1 0.55440034 1094.500 0.0000000 INCAR 71 -1 0.54257020 1.690 0.0000000 PRDIA 71 1 0.27284114 17.000 0.0000000 PAPUL 71 1 0.20466714 23.250 0.0000000 So 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.0000000 PVENT 27 -1 0.10803571 14.500 0.0000000 PRDIA 27 1 0.10803571 18.500 0.0000000 PAPUL 27 1 0.10803571 22.500 0.0000000 INCAR 27 1 0.04705882 1.195 0.0000000 On 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 have cart$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 using

barplot(t(cart$variable.importance),horiz=TRUE) To be continued with more trees… Data Science for Actuaries, Regression Models with R After an introduction to Advanced R, we will discuss for the last part of our crash course visualization and graphs (from the previous set of slides), and I just uploaded additional slides on regression models (including some pdf version) Classification on the German Credit Database In our data science course, this morning, we’ve use random forrest to improve prediction on the German Credit Dataset. The dataset is > url="http://freakonometrics.free.fr/german_credit.csv" > credit=read.csv(url, header = TRUE, sep = ",") Almost all variables are treated a numeric, but actually, most of them are factors, > str(credit) 'data.frame': 1000 obs. of 21 variables:$ Creditability   : int  1 1 1 1 1 1 1 1 1 1 ...
$Account.Balance : int 1 1 2 1 1 1 1 1 4 2 ...$ Duration        : int  18 9 12 12 12 10 8  ...
$Purpose : int 2 0 9 0 0 0 0 0 3 3 ... (etc). Let us convert categorical variables as factors, > F=c(1,2,4,5,7,8,9,10,11,12,13,15,16,17,18,19,20) > for(i in F) credit[,i]=as.factor(credit[,i]) Let us now create our training/calibration and validation/testing datasets, with proportion 1/3-2/3 > i_test=sample(1:nrow(credit),size=333) > i_calibration=(1:nrow(credit))[-i_test] The first model we can fit is a logistic regression, on selected covariates > LogisticModel <- glm(Creditability ~ Account.Balance + Payment.Status.of.Previous.Credit + Purpose + Length.of.current.employment + Sex...Marital.Status, family=binomial, data = credit[i_calibration,]) Based on that model, it is possible to draw the ROC curve, and to compute the AUC (on ne validation dataset) > fitLog <- predict(LogisticModel,type="response", + newdata=credit[i_test,]) > library(ROCR) > pred = prediction( fitLog, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCLog1=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCLog1,"\n")
AUC:  0.7340997

An alternative is to consider a logistic regression on all explanatory variables

> LogisticModel <- glm(Creditability ~ .,
+  family=binomial,
+  data = credit[i_calibration,])

We might overfit, here, and we should observe that on the ROC curve

> fitLog <- predict(LogisticModel,type="response",
+                   newdata=credit[i_test,])
> pred = prediction( fitLog, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCLog2=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCLog2,"\n") AUC: 0.7609792 There is a slight improvement here, compared with the previous model, where only five explanatory variables were considered. Consider now some regression tree (on all covariates) > library(rpart) > ArbreModel <- rpart(Creditability ~ ., + data = credit[i_calibration,]) We can visualize the tree using > library(rpart.plot) > prp(ArbreModel,type=2,extra=1) The ROC curve for that model is > fitArbre <- predict(ArbreModel, + newdata=credit[i_test,], + type="prob")[,2] > pred = prediction( fitArbre, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCArbre=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCArbre,"\n")
AUC:  0.7100323

As expected, a single has a lower performance, compared with a logistic regression. And a natural idea is to grow several trees using some boostrap procedure, and then to agregate those predictions.

> library(randomForest)
> RF <- randomForest(Creditability ~ .,
+ data = credit[i_calibration,])
> fitForet <- predict(RF,
+                     newdata=credit[i_test,],
+                     type="prob")[,2]
> pred = prediction( fitForet, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCRF=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCRF,"\n") AUC: 0.7682367 Here this model is (slightly) better than the logistic regression. Actually, if we create many training/validation samples, and compare the AUC, we can observe that – on average – random forests perform better than logistic regressions, > AUC=function(i){ + set.seed(i) + i_test=sample(1:nrow(credit),size=333) + i_calibration=(1:nrow(credit))[-i_test] + LogisticModel <- glm(Creditability ~ ., + family=binomial, + data = credit[i_calibration,]) + summary(LogisticModel) + fitLog <- predict(LogisticModel,type="response", + newdata=credit[i_test,]) + library(ROCR) + pred = prediction( fitLog, credit$Creditability[i_test])
+   AUCLog2=performance(pred, measure = "auc")@y.values[[1]]
+   RF <- randomForest(Creditability ~ .,
+   data = credit[i_calibration,])
+   fitForet <- predict(RF,
+                       newdata=credit[i_test,],
+                       type="prob")[,2]
+   pred = prediction( fitForet, credit$Creditability[i_test]) + AUCRF=performance(pred, measure = "auc")@y.values[[1]] + return(c(AUCLog2,AUCRF)) + } > A=Vectorize(AUC)(1:200) > plot(t(A)) Econometrics vs. Machine Learning with Temporal Patterns A few months ago, I did publish a (long) post entitled ‘some thoughts on economics, mathematics, econometrics, machine learning, etc‘. In that post, I was discussing possible differences between foundations of econometrics, and machine learning. I wanted to get back today on an important point, related to training/sampling datasets, when we have temporal data. I was discussing this morning, with a student of the Data Science for Actuaries program, an interesting point related to claim frequency models, for insurance ratemaking. Since the goal is to predict claims frequency (to assess the level of the insurance premium), he suggested to use old data to train the model, and more recent one to test it. The problem is that the model did not incorporate any temporal pattern, and we got surprising results. Consider here a simple dataset, > set.seed(1) > n=50000 > X1=runif(n) > T=sample(2000:2015,size=n,replace=TRUE) > L=exp(-3+X1-(T-2000)/20) > E=rbeta(n,5,1) > Y=rpois(n,L*E) > B=data.frame(Y,X1,L,T,E) Claims frequency is driven by a Poisson process, with one covariate, X1, and we assume that the intensity decreases (with an exponential rate). Consider here a standard linear regression, without any time effect > reg=glm(Y~X1+offset(log(E)),data=B, + family=poisson) We can also compute the empirical annualized claims frequency > u=seq(0,1,by=.01) > v=predict(reg,newdata=data.frame(X1=u,E=1)) > p=function(x){ + B=B[abs(B$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp=Vectorize(p)(seq(.05,.95,by=.1))

and plot the two curves on the same graph,

> plot(seq(.05,.95,by=.1),vp,type="b")
> lines(u,exp(v),lty=2,col="red")

This is what we usually do in econometrics. In machine learning, and more specifically to assess the quality of the model, and for model selection, it is common to split the dataset in two parts. A training sample, and a validation sample. Consider some randomized training/validation samples, then fit a model on the training sample, and finally use it to get a prediction,

> idx=sample(1:nrow(B),size=nrow(B)*7/8)
> B_a=B[idx,]
> B_t=B[-idx,]
> reg=glm(Y~X1+offset(log(E)),data=B_a,
+ family=poisson)
> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B_a[abs(B_a$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_a=Vectorize(p)(seq(.05,.95,by=.1)) > plot(seq(.05,.95,by=.1),vp_a,col="blue") > lines(u,exp(v),lty=2) > p=function(x){ + B=B_t[abs(B_t$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_t=Vectorize(p)(seq(.05,.95,by=.1))
> lines(seq(.05,.95,by=.1),vp_t,col="red")

The blue curve is the prediction on the training sample (as we usually do in econometrics), but then the red curve is the prediction on the testing sample. Here, volatility probably comes from the small size of the testing sample (1 observation out of 8).

Now, what if we use the year as a splitting criteria : we fit a model on old years to fit a model, and we test it on recent years,

> B_a=subset(B,T<2014)
> B_t=subset(B,T>=2014)
> reg=glm(Y~X1+offset(log(E)),data=B_a,family=poisson)
> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B_a[abs(B_a$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_a=Vectorize(p)(seq(.05,.95,by=.1)) > plot(seq(.05,.95,by=.1),vp_a,col="blue") > lines(u,exp(v),lty=2) > p=function(x){ + B=B_t[abs(B_t$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_t=Vectorize(p)(seq(.05,.95,by=.1))
> lines(seq(.05,.95,by=.1),vp_t,col="red")

Clearly, we miss something here…

We were looking at such a graph this morning, and it took me some time to understand how training and validation samples were designed, and that there was a possible temporal effect (actually, this morning, it was based on a 3 year training sample, and a 1 year validation sample).

Since there is a temporal pattern, let us capture it. As an econometrician, let me use a regression model

> reg=glm(Y~X1+T+offset(log(E)),data=B,
+ family=poisson)
> C=coefficients(reg)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C[1]+C[3]*(2000:2015)))
> lines(u,v,lty=2,col="red")

(I focus only on the evolution of the temporal variate on that graph).

Here, we use a linear model, but there are usually no reason to assume linearity. So we might consider splines

> library(splines)
> reg=glm(Y~X1+bs(T)+offset(log(E)),
+ data=B,family=poisson)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> v2=predict(reg,newdata=data.frame(X1=0,
+ T=2000:2015,E=1))
> plot(2000:2015,exp(v2),type="b")
> lines(u,v,lty=2,col="red")

But here again, why should we assume that there is an underlying smooth function? There might be some ruptures… So let us consider a regression on factors

> reg=glm(Y~0+X1+as.factor(T)+offset(log(E)),
+ data=B,family=poisson)
> C=coefficients(reg)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C[2:17]),type="b")
> lines(u,v,lty=2,col="red")

An alternative might be to consider some more general model, like a regression tree

> library(rpart)
> reg=rpart(Y~X1+T+offset(log(E)),data=B,
+ method="poisson",cp=1e-4)
> p=function(t){
+   B=B[B$T==t,] + B$E=1
+   mean(predict(reg,newdata=B))
+ }
> y_m=Vectorize(function(t) p(t))(2000:2015)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3+.5)
> plot(2000:2015,y_m,ylim=c(.02,.085),type="b")
> lines(u,v,lty=2,col="red")

Here, it seems that something went wrong. I guess it’s coming from the exposure. So consider a simplier model, on the annualized frequency, and with weights that are related to the exposure

> reg=rpart(Y/E~X1+T,data=B,weights=B$E,cp=1e-4) > p=function(t){ + B=B[B$T==t,]
+   B$E=1 + mean(predict(reg,newdata=B)) + } > y_m=Vectorize(function(t) p(t))(2000:2015) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3+.5) > plot(2000:2015,y_m,ylim=c(.02,.085),type="b") > lines(u,v,lty=2,col="red") That was for the econometrician perspective. With a machine learning perspective, consider a training sample (here based on old data) and a validation sample (based on more recent ones) > B_a=subset(B,T<2014) > B_t=subset(B,T>=2014) If we consider a model, it is easy to get a prediction on recent years, even if the model was designed to model older ones, > reg_a=glm(Y~X1+T+offset(log(E)), + data=B_a,family=poisson) > C=coefficients(reg_a) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > plot(2000:2015,exp(C[1]+C[3]*c(2000:2013, + NA,NA)),type="b") > lines(u,v,lty=2,col="red") > points(2014:2015,exp(C[1]+C[3]*2014:2015), + pch=19,col="blue") But if we use years as factors, things are more complicated. > reg_a=glm(Y~0+X1+as.factor(T)+offset(log(E)), + data=B_a,family=poisson) > C=coefficients(reg_a) > RMSE=function(A){ + L=exp(C[1]*B_t$X1+ A[1]*(B_t$T==2014) + A[2]*(B_t$T==2015))
+   Y_t=L*B_t$E + sum( (Y_t - B_t$Y )^2)}
> i=optim(c(.4,.4),RMSE)$par > plot(2000:2015,c(exp(C[2:15]),NA,NA),) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > lines(u,v,lty=2,col="red") > points(2014:2015,exp(i),pch=19,col="blue") becase we need to get a prediction on levels that were not in our training sample. Here, we minimize the RMSE to quantify factor levels for recent years. And the output is not that bad. So yes, it is possible to get a training dataset on older data, and test it on recent years. But one should be careful, and take into account, properly, temporal patterns. How Could Classification Trees Be So Fast on Categorical Variables? I think that over the past months, I have been saying non-correct things about classification with categorical covariates. Because I never took time to look at it carefuly. Consider some simulated dataset, with a logistic regression, > n=1e3 > set.seed(1) > X1=runif(n) > q=quantile(X1,(0:26)/26) > q[1]=0 > X2=cut(X1,q,labels=LETTERS[1:26]) > p=exp(-.1+qnorm(2*(abs(.5-X1))))/(1+exp(-.1+qnorm(2*(abs(.5-X1))))) > Y=rbinom(n,size=1,p) > df=data.frame(X1=X1,X2=X2,p=p,Y=Y) Here, we use some continuous covariate$X_1$, except that$X_1$ is considered as not-observed. Instead, we have a categorical covariate with 26 categories. The (theoretical) relationship between the covariate and the probability is given below, > vx1=seq(0,1,by=.001) > vp=exp(-.1+qnorm(2*(abs(.5-vx1))))/(1+exp(-.1+qnorm(2*(abs(.5-vx1))))) > plot(vx1,vp,type="l") and the empirical probability, for each modality is If we run a classification tree, we get > library(rpart) > tree=rpart(Y~X2,data=df) > library(rpart.plot) > prp(tree, type=2, extra=1) To be more specific, the output is here > tree 1) root 1000 249.90000 0.4900000 2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.3 0.302 4) X2=J,K,L,M,N,O,P,Q,R 346 65.12 0.25144 * 5) X2=F,G,H,I 153 37.22876 0.4183007 * 3) X2=A,B,C,D,E,S,T,U,V,W,X,Y,Z 501 109.61 0.67 6) X2=B,C,D,E,S,T,U,V,W,X 385 90.38 0.623 * 7) X2=A,Y,Z 116 14.50862 0.8534483 * Note that it takes less than a second to get that output. So clearly, we did not look for all combinations between modalities. For the first node, there are like $2^{26}$ possible groups, i.e. > 67108864 It is big… not huge, but too big to try all combinations, since that’s only the first node, and we have to do it again on the two leaves, etc. Antoine (aka @ly_antoine) told me – while we were having a coffee after lunch today – the trick to get a fast algorithm, on categories. And as usual, the idea is very clever… First, we need a function to compute 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)} For the first node, the idea is very simple: • Compute empirical averages $p_i=\mathbb{E}[Y \vert X=\{x_i\}]$ > cond_prob=aggregate(df$Y,by=list(df$X2),mean) • Then sort those values, $p_{1:k}\leq p_{2:k} \leq \cdots \leq p_{k:k}$, • Based on that ordering, consider $\{x_{1:k},x_{2:k},\cdots,x_{k:k}\}$ > Group_Letters=cond_prob[order(cond_prob$x),2]

• Then consider (only) $k$ possible partitions,

$\{x_{1:k},x_{2:k},\cdots,x_{i:k}\}$ against $\{x_{i+1:k},x_{i+2:k},\cdots,x_{k:k}\}$

> v_gini=rep(NA,26)
> for(v in 1:26){
+   CLASSE=df$X2 %in% Group_Letters[1:v] + v_gini[v]=gini(y=df$Y,classe=CLASSE)
+ }

If we plot them, we get

> plot(1:26,v_gini,type="b)

As for continuous variables, we seek for the maximum value, and then, we have our two groups,

> sort(Group_Letters[1:which.max(v_gini)])
[1] F G H I J K L M N O P Q R

That’s exactly what we got with the tree function in R,

1) root 1000 249.90000 0.4900000
2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.30 0.30

Now, consider the leaf on the left (for instance)

> sub_df=df[df$X2 %in% sort(Group_Letters[1:which.max(v_gini)]),] Then use the same algorithm as before: sort the conditional means, > cond_prob=aggregate(sub_df$Y,by=
+ list(sub_df$X2),mean) > s_Group_Letters=cond_prob[order(cond_prob$x),2]

Then compute Gini indices based on groups obtained from that ordering,

> v_gini=rep(NA,length(sub_Group_Letters))
> for(v in 1:length(sub_Group_Letters)){
+   CLASSE=sub_df$X2 %in% s_Group_Letters[1:v] + v_gini[v]=gini(y=sub_df$Y,classe=CLASSE)
+ }

If we plot it, we get our two groups,

> plot(1:length(s_Group_Letters),v_gini,type="b")


And the first group is here

> sort(sub_Group_Letters[1:which.max(v_gini)])
[1] J K L M N O P Q R

Again, that’s exactly what we got with the R function

1) root 1000 249.90000 0.4900000
2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.30 0.30
4) X2=J,K,L,M,N,O,P,Q,R 346  65.12 0.25144  *

Clever, isn’t?

Actuariat de l’Assurance Non-Vie #2

Pour le second cours d’actuariat de l’assurance non-vie à l’ENSAE, qui aura lieu lundi après midi, les slides présentant les modèles classiques pour prédire des variables factorielles (classification) sont en ligne,

Computational Time of Predictive Models

Tuesday, at the end of my 5-hour crash course on machine learning for actuaries, Pierre asked me an interesting question about computational time of different techniques. I’ve been presenting the philosophy of various algorithm, but I forgot to mention computational time. I wanted to try several classification algorithms on the dataset used to illustrate the techniques

> rm(list=ls())
"http://freakonometrics.free.fr/myocarde.csv",
> levels(myocarde$PRONO)=c("Death","Survival") But the dataset is rather small, with 71 observations and 7 explanatory variables. So I decided to replicate the observations, and to add some covariates, > levels(myocarde$PRONO)=c("Death","Survival")
> idx=rep(1:nrow(myocarde),each=100)
> TPS=matrix(NA,30,10)
> myocarde_large=myocarde[idx,]
> k=23
> M=data.frame(matrix(rnorm(k*
+ nrow(myocarde_large)),nrow(myocarde_large),k))
> names(M)=paste("X",1:k,sep="")
> myocarde_large=cbind(myocarde_large,M)
> dim(myocarde_large)
[1] 7100   31
> object.size(myocarde_large)
2049.064 kbytes

The dataset is not big… but at least, it does not take 0.0001 sec. to run a regression.  Actually, to run a logistic regression, it takes 0.1 second

> system.time(fit< glm(PRONO~.,
+ data=myocarde_large, family="binomial"))
user      system     elapsed
0.114       0.016       0.134
> object.size(fit)
9,313.600 kbytes

And I was surprised that the regression object was 9Mo, which is more than four times the size of the dataset. With a large dataset, 100 times larger,

> dim(myocarde_large_2)
[1] 710000     31

it takes 20 sec.

> system.time(fit<-glm(PRONO~.,
+ data=myocarde_large_2, family="binomial"))
utilisateur     système      écoulé
16.394       2.576      19.819
> object.size(fit)
90,9025.600 kbytes

and the object is ‘only’ ten times bigger.

Choosing a Classifier

In order to illustrate the problem of chosing a classification model consider some simulated data,

> n = 500
> set.seed(1)
> X = rnorm(n)
> ma = 10-(X+1.5)^2*2
> mb = -10+(X-1.5)^2*2
> M = cbind(ma,mb)
> set.seed(1)
> Z = sample(1:2,size=n,replace=TRUE)
> Y = ma*(Z==1)+mb*(Z==2)+rnorm(n)*5
> df = data.frame(Z=as.factor(Z),X,Y)

A first strategy is to split the dataset in two parts, a training dataset, and a testing dataset.

> df1 = training = df[1:300,]
> df2 = testing  = df[301:500,]
• The Holdout Method: Training and Testing Datasets

The two datasets can be visualised below, with the training dataset on top, and the testing dataset below

> plot(df1$X,df1$Y,pch=19,col=c(rgb(1,0,0,.4),
+ rgb(0,0,1,.4))[df1$Z]) Variable Selection using Cross-Validation (and Other Techniques) A natural technique to select variables in the context of generalized linear models is to use a stepŵise procedure. It is natural, but contreversial, as discussed by Frank Harrell in a great post, clearly worth reading. Frank mentioned about 10 points against a stepwise procedure. • It yields R-squared values that are badly biased to be high. • The F and chi-squared tests quoted next to each variable on the printout do not have the claimed distribution. • The method yields confidence intervals for effects and predicted values that are falsely narrow (see Altman and Andersen (1989)). • It yields p-values that do not have the proper meaning, and the proper correction for them is a difficult problem. • It gives biased regression coefficients that need shrinkage (the coefficients for remaining variables are too large (see Tibshirani (1996)). • It has severe problems in the presence of collinearity. • It is based on methods (e.g., F tests for nested models) that were intended to be used to test prespecified hypotheses. • Increasing the sample size does not help very much (see Derksen and Keselman (1992)). • It allows us to not think about the problem. • It uses a lot of paper. ‘Variable Importance Plot’ and Variable Selection 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. > set.seed(1) > n=500 > library(clusterGeneration) > library(mnormt) > S=genPositiveDefMat("eigen",dim=15) > S=genPositiveDefMat("unifcorrmat",dim=15) > X=rmnorm(n,varcov=S$Sigma)
> library(corrplot)
> corrplot(cor(X), order = "hclust")

See Gosh & Hendersen (2003) for more details on the methodology.

Spliting a Node in a Tree

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",
> library(rpart)
> cart<-rpart(PRONO~.,data=MYOCARDE)

we get

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

Regression Models, It’s Not Only About Interpretation

Yesterday, I did upload a post where I tried to show that “standard” regression models where not performing bad. At least if you include splines (multivariate splines) to take into accound joint effects, and nonlinearities. So far, I do not discuss the possible high number of features (but with boostrap procedures, it is possible to assess something related to variable importance, that people from machine learning like).

But my post was not complete: I was simply plotting the prediction obtained by some model. And it “looked like” the regression was nice, but so were the random forrest, the $k$-nearest neighbour and boosting algorithm. What if we compare those models on new data?

On Some Alternatives to Regression Models

When you start discussing with people in machine learning, you quickly hear something like “forget your econometric models, your GLMs, I can easily find a machine learning ‘model’ that can beat yours”. I am usually very sceptical, especially when I hear “easily” or “always“. I have no problem about the fact that I use old econometric models, but I had the feeling that things aren’t that easy. I can understand that we might have problems when we do have a lot of features (I am still working on that, I’ll get back to this point soon), but I have the feeling that I can still capture interactions, and non-linearities with standard econometric models as well as any machine learning algorithm.

Just to illustrate, consider the following ‘model

$\mathbb{E}[Y\vert\boldsymbol{X}=\boldsymbol{x}]=m(\boldsymbol{x})$

where $m(\cdot)$ is (just to illustrate)

> n <- 5000
> rtf <- function(x1, x2) { sin(x1+x2)/(x1+x2) }
> xgrid <- seq(1,6,length=31)
> ygrid <- seq(1,6,length=31)
> zgrid <- outer(xgrid,ygrid,rtf)
> persp(xgrid,ygrid,zgrid,theta=30, phi=30,
+ col="green", ticktype="detailed",shade=TRUE)

Growing one Tree

Consider the following toy dataset, with some spam/ham information, and two words, “viagra” and “lottery”.

> load(spam.RData)
Y viagra lottery
27 spam      0       1
37  ham      0       1
57 spam      0       0
89  ham      0       0
20 spam      1       0
86  ham      0       0

For the first node, compute Gini index for the two variables,

> gini=function(variable){
+ T=table(db\$Y,db[,variable])
+ nx=apply(T,2,sum)
+ ProbCond=T/matrix(rep(nx,each=2),2,2)
+ ProbCond
+ Gini=-ProbCond*(1-ProbCond)
+ sum(matrix(rep(nx,each=2),2,2)/sum(nx)*Gini)}
> gini("viagra")
[1] -0.44
> gini("lottery")
[1] -0.487

Here Gini index is maximal for “viagra”, so that will be the first node.