# Insurance data science : Pictures

At the Summer School of the Swiss Association of Actuaries, in Lausanne, following the part of Jean-Philippe Boucher (UQAM) on telematic data, I will start talking about pictures this Wednesday. Slides are available online

Ewen Gallic (AMSE) will present a tutorial on satellite pictures, and a simple classification problem, related to Alzeimher detection.

We will try to identify what is on the following pictures, starting with the car

(we will see that the car is indeed identified)

a skier,

and a fire,

We will also discuss previous pictures from the summer school

# Classification from scratch, boosting 11/8

Eleventh post of our series on classification from scratch. Today, that should be the last one… unless I forgot something important. So today, we discuss boosting.

## An econometrician perspective

I might start with a non-conventional introduction. But that’s actually how I understood what boosting was about. And I am quite sure it has to do with my background in econometrics.

The goal here is to solve something which looks like$$m^\star=\underset{m\in\mathcal{M}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \ell(y_i,m(\mathbf{x}_i))\right\rbrace$$for some loss function $\ell$, and for some set of predictors $\mathcal{M}$. This is an optimization problem. Well, optimization is here in a function space, but still, that’s simply an optimization problem. And from a numerical perspective, optimization is solve using gradient descent (this is why this technique is also called gradient boosting). And the gradient descent can be visualized like below

Again, the optimum is not some some real value $x^\star$, but some function $m^\star$. Thus, here we will have something like$$m^{(k)}=m^{(k-1)}+\underset{h\in\mathcal{H}}{\text{argmin}}\left\lbrace \sum_{i=1}^n \ell(y_i,m^{(k-1)}(\mathbf{x}_i)+h(\mathbf{x}_i))\right\rbrace$$(as they write it is serious articles) where the term on the right can also be written$$m^{(k)}=m^{(k-1)}+\underset{h\in\mathcal{H}}{\text{argmin}}\left\lbrace \sum_{i=1}^n \ell(\underbrace{y_i-m^{(k-1)}(\mathbf{x}_i)}_{\varepsilon_{k,i}},h(\mathbf{x}_i))\right\rbrace$$I prefer the later, because we see clearly that $f$ is some model we fit on the remaining residuals.

We can rewrite it like that: define$$r_{i,k}=-\left.\frac{\partial \ell(y_i,m(\mathbf{x}_i))}{\partial m(\mathbf{x}_i)}\right\vert_{m(\mathbf{x}_i)=m^{(k-1)}(\mathbf{x}_i)}$$for all $i=1,\cdots,n$. The goal is to fit a model so that $r_{i,k}=h^\star(\mathbf{x}_i)$, and when we have that optimal function, set $m_k(\mathbf{x})=m_{k-1}(\mathbf{x})+\gamma_k h^\star(\mathbf{x})$ (yes, we can include some shrinkage here).

Two important comments here. First of all, the idea should be weird to any econometrician. First, we fit a model to explain $y$ by some covariates $\mathbf{x}$. Then consider the residuals $\widehat{\varepsilon}$, and to explain them with the same covariate $\mathbf{x}$. If you try that with a linear regression, you’d done at the end of step 1, since residuals $\widehat{\varepsilon}$ are orthogonal to covariates $\mathbf{x}$: no way that we can learn from them. Here it works because we consider simple non linear model. And actually, something that can be used is to add a shrinkage parameter. Do not consider $\widehat{\varepsilon}=y-\widehat{m}(\mathbf{x})$ but $\widehat{\varepsilon}=y-\gamma\widehat{m}(\mathbf{x})$. The idea of weak learners is extremely important here. The more we shrink, the longer it will take, but that’s not (too) important.

I should also mention that it’s nice to keep learning from our mistakes. But somehow, we should stop, someday. I said that I will not mention this part in this series of posts, maybe later on. But heuristically, we should stop when we start to overfit. And this can be observed either using a split training/validation of the initial dataset or to use cross validation. I will get back on that issue later one in this post, but again, those ideas should probably be dedicated to another series of posts.

## Learning with splines

Just to make sure we get it, let’s try to learn with splines. Because standard splines have fixed knots, actually, we do not really “learn” here (and after a few iterations we get to what we would have with a standard spline regression). So here, we will (somehow) optimize knots locations. There is a package to do so. And just to illustrate, use a Gaussian regression here, not a classification (we will do that later on). Consider the following dataset (with only one covariate)

n=300 set.seed(1) u=sort(runif(n)*2*pi) y=sin(u)+rnorm(n)/4 df=data.frame(x=u,y=y)

For an optimal choice of knot locations, we can use

library(freeknotsplines) xy.freekt=freelsgen(df$x, df$y, degree = 1, numknot = 2, 555)

With 5% shrinkage, the code it simply the following

v=.05 library(splines) xy.freekt=freelsgen(df$x, df$y, degree = 1, numknot = 2, 555) fit=lm(y~bs(x,degree=1,knots=xy.freekt@optknot),data=df) yp=predict(fit,newdata=df) df$yr=df$y - v*yp YP=v*yp for(t in 1:200){ xy.freekt=freelsgen(df$x, df$yr, degree = 1, numknot = 2, 555) fit=lm(yr~bs(x,degree=1,knots=xy.freekt@optknot),data=df) yp=predict(fit,newdata=df) df$yr=df$yr - v*yp YP=cbind(YP,v*yp)} nd=data.frame(x=seq(0,2*pi,by=.01)) viz=function(M){ if(M==1) y=YP[,1] if(M&gt;1) y=apply(YP[,1:M],1,sum) plot(df$x,df$y,ylab="",xlab="") lines(df$x,y,type="l",col="red",lwd=3) fit=lm(y~bs(x,degree=1,df=3),data=df) yp=predict(fit,newdata=nd) lines(nd$x,yp,type="l",col="blue",lwd=3) lines(nd$x,sin(nd$x),lty=2)}

To visualize the ouput after 100 iterations, use

viz(100)

Clearly, we see that we learn from the data here… Cool, isn’t it?

## Learning with stumps (and trees)

Let us try something else. What if we consider at each step a regression tree, instead of a linear-by-parts regression (that was considered with linear splines).

library(rpart) v=.1 fit=rpart(y~x,data=df) yp=predict(fit) df$yr=df$y - v*yp YP=v*yp for(t in 1:100){ fit=rpart(yr~x,data=df) yp=predict(fit,newdata=df) df$yr=df$yr - v*yp YP=cbind(YP,v*yp)}

Again, to visualise the learning process, use

viz=function(M){ y=apply(YP[,1:M],1,sum) plot(df$x,df$y,ylab="",xlab="") lines(df$x,y,type="s",col="red",lwd=3) fit=rpart(y~x,data=df) yp=predict(fit,newdata=nd) lines(nd$x,yp,type="s",col="blue",lwd=3) lines(nd$x,sin(nd$x),lty=2)}

This time, with those trees, it looks like not only we have a good model, but also a different model from the one we can get using a single regression tree.

What if we change the shrinkage parameter?

viz=function(v=0.05){ fit=rpart(y~x,data=df) yp=predict(fit) df$yr=df$y - v*yp YP=v*yp for(t in 1:100){ fit=rpart(yr~x,data=df) yp=predict(fit,newdata=df) df$yr=df$yr - v*yp YP=cbind(YP,v*yp)} y=apply(YP,1,sum) plot(df$x,df$y,xlab="",ylab="") lines(df$x,y,type="s",col="red",lwd=3) fit=rpart(y~x,data=df) yp=predict(fit,newdata=nd) lines(nd$x,yp,type="s",col="blue",lwd=3) lines(nd$x,sin(nd$x),lty=2)}

There is clearly an impact of that shrinkage parameter. It has to be small to get a good model. This is the idea of using weak learners to get a good prediction.

## Classification and Adaboost

Now that we understand how bootsting works, let’s try to adapt it to classification. It will be more complicated because residuals are usually not very informative in a classification. And it will be hard to shrink. So let’s try something slightly different, to introduce the adaboost algorithm.

In our initial discussion, the goal was to minimize a convex loss function. Here, if we express classes as $\{-1,+1\}$, the loss function we consider is $e^{-y\cdot m(\mathbf{x})}$ (this product $y\cdot m(\mathbf{x})$) was already discussed when we’ve seen the SVM algorithm. Note that the loss function related to the logistic model would be $\log(1+e^{-y\cdot m(\mathbf{x})})$.

What we do here is related to gradient descent (or Newton algorithm). Previously, we were learning from our errors. At each iteration, the residuals are computed and a (weak) model is fitted to these residuals. The the contribution of this weak model is used in a gradient descent optimization process. Here things will be different, because (from my understanding) it is more difficult to play with residuals, because null residuals never exist in classifications. So we will add weights. Initially, all the observations will have the same weights. But iteratively, we ill change them. We will increase the weights of the wrongly predicted individuals and decrease the ones of the correctly predicted individuals. Somehow, we want to focus more on the difficult predictions. That’s the trick. And I guess that’s why it performs so well. This algorithm is well described in wikipedia, so we will use it.

We start with $\mathbf{\omega}_0=\mathbf{1}/n$, then at each step fit a model (a classification tree) with weights $\mathbf{\omega}_k$(we did not discuss weights in the algorithms of trees, but it is straigtforward in the formula actually). Let $\widehat{h}_{\mathbf{\omega}_k}$ denote that model (i.e. the probability in each leaves). Then consider the classifier $2~\mathbf{1}[\widehat{h}_{\mathbf{\omega}_k}(\cdot)>0.5]-1$ which returns a value in $\{-1,+1\}$. Then set $$\varepsilon_k=\sum_{i\in\mathcal{I}_k}\omega_i$$where $\mathcal{I}_k$ is the set of misclassified individuals,$$\mathcal{I}_k=\big\lbrace i:2~\mathbf{1}[\widehat{h}_{\mathbf{\omega}_k}(\mathbf{x}_i)>0.5]-1\neq y_i\big\rbrace$$Then set $$\alpha_k = \frac{1}{2} \ln \left(\frac{1-\epsilon_k}{\epsilon_k}\right)$$and update finally the model using$$m_{k=1}=m_k+\alpha_k\widehat{h}_{\mathbf{\omega}_k}$$as well as the weights$$\mathbf{\omega}_{k+1}=\mathbf{\omega}_k e^{-\mathbf{y} \alpha_k \widehat{h}_{\mathbf{\omega}_k}(\mathbf{x}_i)}$$(of course, devide by the sum to insure that the total sum is then 1). And as previously, one can include some shrinkage. To visualize the convergence of the process, we will plot the total error on our dataset.

n_iter = 100 y = (myocarde[,"PRONO"]==1)*2-1 x = myocarde[,1:7] error = rep(0,n_iter) f = rep(0,length(y)) w = rep(1,length(y)) # alpha = 1 library(rpart) for(i in 1:n_iter){ w = exp(-alpha*y*f) *w w = w/sum(w) rfit = rpart(y~., x, w, method="class") g = -1 + 2*(predict(rfit,x)[,2]&gt;.5) e = sum(w*(y*g&lt;0)) alpha = .5*log ( (1-e) / e ) alpha = 0.1*alpha f = f + alpha*g error[i] = mean(1*f*y&lt;0) } plot(seq(1,n_iter),error,type=&quot;l&quot;, ylim=c(0,.25),col=&quot;blue&quot;, ylab=&quot;Error Rate&quot;,xlab=&quot;Iterations&quot;,lwd=2)

Here we face a classical problem in machine learning: we have a perfect model. With zero error. That is nice, but not interesting. It is also possible in econometrics, with polynomial fits: with 10 observations, and a polynomial of degree 9, we have a perfect fit. But a poor model. Here it is the same. So the trick is to split our dataset in two, a training dataset, and a validation one

set.seed(123) id_train = sample(1:nrow(myocarde), size=45, replace=FALSE) train_myocarde = myocarde[id_train,] test_myocarde = myocarde[-id_train,]

We construct the model on the first one, and we check on the second one that it’s not that bad…

y_train = (train_myocarde[,"PRONO"]==1)*2-1 x_train = train_myocarde[,1:7] y_test = (test_myocarde[,"PRONO"]==1)*2-1 x_test = test_myocarde[,1:7] train_error = rep(0,n_iter) test_error = rep(0,n_iter) f_train = rep(0,length(y_train)) f_test = rep(0,length(y_test)) w_train = rep(1,length(y_train)) alpha = 1 for(i in 1:n_iter){ w_train = w_train*exp(-alpha*y_train*f_train) w_train = w_train/sum(w_train) rfit = rpart(y_train~., x_train, w_train, method="class") g_train = -1 + 2*(predict(rfit,x_train)[,2]&gt;.5) g_test = -1 + 2*(predict(rfit,x_test)[,2]&gt;.5) e_train = sum(w_train*(y_train*g_train&lt;0)) alpha = .5*log ( (1-e_train) / e_train ) alpha = 0.1*alpha f_train = f_train + alpha*g_train f_test = f_test + alpha*g_test train_error[i] = mean(1*f_train*y_train&lt;0) test_error[i] = mean(1*f_test*y_test&lt;0)} plot(seq(1,n_iter),test_error,col='red') lines(train_error,lwd=2,col='blue')

Here, as previously, after 80 iterations, we have a perfect model on the training dataset, but it behaves badly on the validation dataset. But with 20 iterations, it seems to be ok…

## R function

Of course, it’s possible to use R functions,

library(gbm) gbmWithCrossValidation = gbm(PRONO ~ .,distribution = "bernoulli", data = myocarde,n.trees = 2000,shrinkage = .01,cv.folds = 5,n.cores = 1) bestTreeForPrediction = gbm.perf(gbmWithCrossValidation)

Here cross-validation is considered, and not training/validation, as well as forests instead of single trees, but overall, the idea is the same… Off course, the output is much nicer (here the shrinkage is a very small parameter, and learning is extremely slow)

# 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… # Classification from scratch, linear discrimination 8/8 Eighth post of our series on classification from scratch. The latest one was on the SVM, and today, I want to get back on very old stuff, with here also a linear separation of the space, using Fisher’s linear discriminent analysis. ## Bayes (naive) classifier Consider the follwing naive classification rule$$m^\star(\mathbf{x})=\text{argmin}_y\{\mathbb{P}[Y=y\vert\mathbf{X}=\mathbf{x}]\}$$or$$m^\star(\mathbf{x})=\text{argmin}_y\left\{\frac{\mathbb{P}[\mathbf{X}=\mathbf{x}\vert Y=y]}{\mathbb{P}[\mathbf{X}=\mathbf{x}]}\right\}$$(where $\mathbb{P}[\mathbf{X}=\mathbf{x}]$ is the density in the continuous case). In the case where $y$ takes two values, that will be standard $\{0,1\}$ here, one can rewrite the later as$$m^\star(\mathbf{x})=\begin{cases}1\text{ if }\mathbb{E}(Y\vert \mathbf{X}=\mathbf{x})>\displaystyle{\frac{1}{2}}\\0\text{ otherwise}\end{cases}$$and the set$$\mathcal{D}_S =\left\{\mathbf{x},\mathbb{E}(Y\vert \mathbf{X}=\mathbf{x})=\frac{1}{2}\right\}$$is called the decision boundary. Assume that$$\mathbf{X}\vert Y=0\sim\mathcal{N}(\mathbf{\mu}_0,\mathbf{\Sigma})$$and$$\mathbf{X}\vert Y=1\sim\mathcal{N}(\mathbf{\mu}_1,\mathbf{\Sigma})$$then explicit expressions can be derived.$$m^\star(\mathbf{x})=\begin{cases}1\text{ if }r_1^2< r_0^2+2\displaystyle{\log\frac{\mathbb{P}(Y=1)}{\mathbb{P}(Y=0)}+\log\frac{\vert\mathbf{\Sigma}_0\vert}{\vert\mathbf{\Sigma}_1\vert}}\\0\text{ otherwise}\end{cases}$$where $r_y^2$ is the Manalahobis distance, $$r_y^2 = [\mathbf{X}-\mathbf{\mu}_y]^{\text{{T}}}\mathbf{\Sigma}_y^{-1}[\mathbf{X}-\mathbf{\mu}_y]$$ Let $\delta_y$be defined as$$\delta_y(\mathbf{x})=-\frac{1}{2}\log\vert\mathbf{\Sigma}_y\vert-\frac{1}{2}[{\color{blue}{\mathbf{x}}}-\mathbf{\mu}_y]^{\text{{T}}}\mathbf{\Sigma}_y^{-1}[{\color{blue}{\mathbf{x}}}-\mathbf{\mu}_y]+\log\mathbb{P}(Y=y)$$the decision boundary of this classifier is $$\{\mathbf{x}\text{ such that }\delta_0(\mathbf{x})=\delta_1(\mathbf{x})\}$$which is quadratic in ${\color{blue}{\mathbf{x}}}$. This is the quadratic discriminant analysis. This can be visualized bellow. The decision boundary is here But that can’t be the linear discriminant analysis, right? I mean, the frontier is not linear… Actually, in Fisher’s seminal paper, it was assumed that $\mathbf{\Sigma}_0=\mathbf{\Sigma}_1$. In that case, actually, $$\delta_y(\mathbf{x})={\color{blue}{\mathbf{x}}}^{\text{T}}\mathbf{\Sigma}^{-1}\mathbf{\mu}_y-\frac{1}{2}\mathbf{\mu}_y^{\text{T}}\mathbf{\Sigma}^{-1}\mathbf{\mu}_y+\log\mathbb{P}(Y=y)$$ and the decision frontier is now linear in ${\color{blue}{\mathbf{x}}}$. This is the linear discriminant analysis. This can be visualized bellow Here the two samples have the same variance matrix and the frontier is ## Link with the logistic regression Assume as previously that$$\mathbf{X}\vert Y=0\sim\mathcal{N}(\mathbf{\mu}_0,\mathbf{\Sigma})$$and$$\mathbf{X}\vert Y=1\sim\mathcal{N}(\mathbf{\mu}_1,\mathbf{\Sigma})$$then$$\log\frac{\mathbb{P}(Y=1\vert \mathbf{X}=\mathbf{x})}{\mathbb{P}(Y=0\vert \mathbf{X}=\mathbf{x})}$$is equal to $$\mathbf{x}^{\text{{T}}}\mathbf{\Sigma}^{-1}[\mathbf{\mu}_y]-\frac{1}{2}[\mathbf{\mu}_1-\mathbf{\mu}_0]^{\text{{T}}}\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]+\log\frac{\mathbb{P}(Y=1)}{\mathbb{P}(Y=0)}$$which is linear in $\mathbf{x}$$$\log\frac{\mathbb{P}(Y=1\vert \mathbf{X}=\mathbf{x})}{\mathbb{P}(Y=0\vert \mathbf{X}=\mathbf{x})}=\mathbf{x}^{\text{{T}}}\mathbf{\beta}$$Hence, when each groups have Gaussian distributions with identical variance matrix, then LDA and the logistic regression lead to the same classification rule. Observe furthermore that the slope is proportional to $\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]$, as stated in Fisher’s article. But to obtain such a relationship, he observe that the ratio of between and within variances (in the two groups) was$$\frac{\text{variance between}}{\text{variance within}}=\frac{[\mathbf{\omega}\mathbf{\mu}_1-\mathbf{\omega}\mathbf{\mu}_0]^2}{\mathbf{\omega}^{\text{T}}\mathbf{\Sigma}_1\mathbf{\omega}+\mathbf{\omega}^{\text{T}}\mathbf{\Sigma}_0\mathbf{\omega}}$$which is maximal when $\mathbf{\omega}$ is proportional to $\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]$, when $\mathbf{\Sigma}_0=\mathbf{\Sigma}_1$. ## Homebrew linear discriminant analysis To compute vector $\mathbf{\omega}$ m0 = apply(myocarde[myocarde$PRONO=="0",1:7],2,mean) m1 = apply(myocarde[myocarde$PRONO=="1",1:7],2,mean) Sigma = var(myocarde[,1:7]) omega = solve(Sigma)%*%(m1-m0) omega [,1] FRCAR -0.012909708542 INCAR 1.088582058796 INSYS -0.019390084344 PRDIA -0.025817110020 PAPUL 0.020441287970 PVENT -0.038298291091 REPUL -0.001371677757 For the constant – in the equation $\omega^T\mathbf{x}+b=0$ – if we have equiprobable probabilities, use b = (t(m1)%*%solve(Sigma)%*%m1-t(m0)%*%solve(Sigma)%*%m0)/2 ## Application (on the small dataset) In order to visualize what’s going on, consider the small dataset, with only two covariates, x = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) y = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) z = c(1,1,1,1,1,0,0,1,0,0) df = data.frame(x1=x,x2=y,y=as.factor(z)) m0 = apply(df[df$y=="0",1:2],2,mean) m1 = apply(df[df$y=="1",1:2],2,mean) Sigma = var(df[,1:2]) omega = solve(Sigma)%*%(m1-m0) omega [,1] x1 -2.640613174 x2 4.858705676 Using R regular function, we get library(MASS) fit_lda = lda(y ~x1+x2 , data=df) fit_lda Coefficients of linear discriminants: LD1 x1 -2.588389554 x2 4.762614663 which is the same coefficient as the one we got with our own code. For the constant, use b = (t(m1)%*%solve(Sigma)%*%m1-t(m0)%*%solve(Sigma)%*%m0)/2 If we plot it, we get the red straight line plot(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")]) abline(a=b/omega[2],b=-omega[1]/omega[2],col="red")

As we can see (with the blue points), our red line intersects the middle of the segment of the two barycenters

points(m0["x1"],m0["x2"],pch=4) points(m1["x1"],m1["x2"],pch=4) segments(m0["x1"],m0["x2"],m1["x1"],m1["x2"],col="blue") points(.5*m0["x1"]+.5*m1["x1"],.5*m0["x2"]+.5*m1["x2"],col="blue",pch=19)

Of course, we can also use R function

predlda = function(x,y) predict(fit_lda, data.frame(x1=x,x2=y))$class==1 vv=outer(vu,vu,predlda) contour(vu,vu,vv,add=TRUE,lwd=2,levels = .5) One can also consider the quadratic discriminent analysis since it might be difficult to argue that $\mathbf{\Sigma}_0=\mathbf{\Sigma}_1$ fit_qda = qda(y ~x1+x2 , data=df) The separation curve is here plot(df$x1,df$x2,pch=19, col=c("blue","red")[1+(df$y=="1")]) predqda=function(x,y) predict(fit_qda, data.frame(x1=x,x2=y))$class==1 vv=outer(vu,vu,predlda) contour(vu,vu,vv,add=TRUE,lwd=2,levels = .5) # Computational Actuarial Science, with R, in Barcelona This Wednesday, I will give a graduate crash course on computational actuarial science, with R, which will be the second part of the lecture of Tuesday. Slides are now available, # Graduate Course on Classification Techniques, with R This Monday, I will give a crash course on Classification Techniques with R (an overview). Slides are below (there are some typos in the pdf slides, and many graphs missing) and available in pdf format (incluing animated graphs) # Visualising a Classification in High Dimension So far, when discussing classification, we’ve been playing on my toy-dataset (actually, I should no claim it’s mine, it is inspired by the one used in the introduction of Boosting, by Robert Schapire and Yoav Freund). But in ral life, there are more observations, and more explanatory variables.With more than two explanatory variables, it starts to be more complicated to visualise. For instance, consider MYOCARDE=read.table( "http://freakonometrics.free.fr/saporta.csv", head=TRUE,sep=";") where we have observations from people in E.R., for infarctus, and we want to understand who did survive, to get a predictive model. But before running some classifier, let us visualise our data. Since we have seven explanatory variables and our class (survival or death), we can go for a PCA. library(FactoMineR) # ACP (sur les var continues) X=MYOCARDE[,1:7] acp=PCA(X) To add the death/survival variable, treat it as numerical 0/1 variable (at least to get a direction) MYOCARDE2=MYOCARDE MYOCARDE2$PRONO=(MYOCARDE2\$PRONO=="SURVIE")*1
acp=PCA(MYOCARDE2,quanti.sup=8,graph=TRUE)

The nice thing is that we see here where variables are colinear with that one. It is also possible to visualise individuals, and classes, too

acp=PCA(MYOCARDE,quali.sup=8,graph=TRUE)
plot(acp, habillage = 8,col.hab=c("red","blue"))

# Les Arbres de Classification

J’animerai une formation lundi 28 de 14:00 à 16:00 au local N-6320 de l’UQAM sur le thème introduction aux arbres de classification. Cette formation est organisée dans le cadre des séminaires en méthodes d’analyses quantitatives et qualitatives qui se tiennent régulièrement depuis un peu plus d’un mois. animé par le collectif pour le développement et les applications en mesure et évaluation (Cdame). Les slides sont disponibles en pdf (il y a quelques animations, qui ne passent qu’avec Acrobat)

La base utilisée tout au long des exposés est la suivante
> MYOCARDE=read.table("http://freakonometrics.free.fr/saporta.csv",head=TRUE,sep=";")

# Introduction aux arbres de classification

Dans quelques semaines, je ferais une introduction aux arbres de classification dans le cadre d’un séminaire de deux heures organisé par le collectif pour le développement et les applications en mesure et évaluation (Cdame). Je mettrais du matériel en ligne très bientôt. Je ne peux m’empêcher de mentionner les autres séminaires des semaines à venir, car le thème de cette session sera très statistique,

• Lundi 10 mars de 14:00 à 16:00: Karim Oualkacha (UQAM), Introduction aux modèles linéaires généralisés et leurs applications [Résumé]
• Lundi 17 mars de 14:00 à 16:00: André Achim (UQAM), La détermination du nombre de dimensions en analyse factorielle exploratoire: bien mieux que l’analyse parallèle [Résumé]
• Lundi 7 avril de 14:00 à 16:00: Gérald Boutin (UQAM), Les entretiens de recherche qualitatifs : des savoirs théoriques à la pratique [Résumé]
• Lundi 14 avril de 14:00 à 16:00: Jean-François Angers (Université de Montréal), Bayes 101. [Résumé]
• Lundi 28 avril de 14:00 à 16:00: Arthur Charpentier (UQAM), Les arbres de classification. [Résumé]