Classification from scratch, neural nets 6/8

Sixth post of our series on classification from scratch. The latest one was on the lasso regression, which was still based on a logistic regression model, assuming that the variable of interest $Y$ has a Bernoulli distribution. From now on, we will discuss technique that did not originate from those probabilistic models, even if they might still have a probabilistic interpretation. Somehow. Today, we will start with neural nets.

Maybe I should start with a disclaimer. The goal is not to replicate well designed R functions, used for predictive modeling. It is simply to get a basic understanding of what’s going on.

Networs, nodes and edges

First of all, neurals nets are nets, or networks. I will skip the parallel with “neural” stuff because it does not help me understanding what is happening (all apologies for my poor knowledge on biology, and cells)

So, it’s about some network. Networks have nodes, and edges (possibly connected) that connect nodes,

or maybe, to more specific (at least it helped me understanding what’s going on), some sort of flow network,

In such a network, we usually have sources (here multiple) sources (here $\color{red}\{s_1,s_2,s_3\}$), on the left, on a sink (here $\{\color{blue}t\}$), on the right. To continue with this metaphorical introduction, information from the sources should reach the sink. An usually, sources are explanatory variables, $\{\mathbf{x}_1,\cdots,\mathbf{x}_p\}$, and the sink is our variable of interest $\mathbf{y}$. And we want to create a graph, from the sources to the sink. We will have directed edges, with only one (unique) direction, where we will put weights. It is not a flow, the parallel with flow will stop here. For instance, the most simple network will be the following one, with no layer (i.e no node between the source and the sink)

The output here is a binary variable $y\in\{0,1\}$ (it can also be $y\in\{-1,+1\}$ but here, it’s not a big deal). In our network, our output will be $y\in(0,1)$, because it is more easy to handly. For instance, consider $y=f($something$)$, for some function $f$ taking values in $(0,1)$. One can consider the sigmoid function$$f(x)=\frac{1}{1+e^{-x}}=\frac{e^{x}}{e^{x}+1}$$which is actually the logistic function (so we should not be surprised to have results somehow close the logistic regression…). This function $f$ is called the activation function, and there are thousands of such functions. If $y\in\{-1,+1\}$, people consider the hyperbolic tangent$$f(x)=\tanh(x)={\frac {(e^{x}-e^{-x})}{(e^{x}+e^{-x})}}$$or the inverse tangent function
$$f(x)=\tan ^{-1}(x)$$And as input for such function, we consider a weighted sum of incoming nodes. So here$$y_i=f\left(\sum_{j=1}^p\omega_j x_{j,i}\right)$$We can also add a constant actually$$y_i=f\left(\omega_0+\sum_{j=1}^p\omega_j x_{j,i}\right)$$So far, we are not far away from the logistic regression. Except that our starting point was a probabilistic model, in the sense that the later was interpreted as a probability (the probability that $Y=1$) and we wanted the model with the highest likelihood. But we’ll talk about selection of weights later one. First, let us construct our first (very simple) neural network. First, we have the sigmoid function

 sigmoid = function(x) 1 / (1 + exp(-x))

The consider some weights. In our model with seven explanatory variables, with need 7 weights. Or 8 if we include the constant term. Let us consider $\mathbf{\omega}=\mathbf{1}$,

 weights_0 = rep(1,8) X = as.matrix(cbind(1,myocarde[,1:7])) y_5_1 = sigmoid(X %*% weights_0)

that’s kind of stupid because all our predictions are 1, here. Let us try something else. Like $\mathbf{\omega}=\widehat{\mathbf{\beta}}^{ols}$. It is optimized, somehow, but we needed something to visualize what’s going on

 weights_0 = lm(PRONO~.,data=myocarde)$coefficients then use  y_5_1 = sigmoid(X %*% weights_0) In order to see if we get a “good” prediction, let use plot the ROC curve, and compare it with the one we got with a (simple) logistic regression  library(ROCR) pred = ROCR::prediction(y_5_1,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) reg = glm(PRONO~.,data=myocarde,family=binomial(link = "logit")) y_0 = predict(reg,type="response") pred0 = ROCR::prediction(y_0,myocarde$PRONO) perf0 = ROCR::performance(pred0,"tpr", "fpr") plot(perf0,add=TRUE,col="red") That’s not bad for a very first attempt. Except that we’ve been cheating here, since we did use $\mathbf{\omega}=\widehat{\mathbf{\beta}}^{ols}$. How, for real, should we choose those weights? Using a loss function Well, if we want an “optimal” set of weights, we need to “optimize” an objective function. So we need to quantify the loss of a mistake, between the prediction, and the observation. Consider here a quadratic loss function  loss = function(weights){ mean( (myocarde$PRONO-sigmoid(X %*% weights))^2) }

It might be stupid to use a quadratic loss function for a classification, but here, it’s not the point. We just want to understand what is the algorithm we use, and the loss function $\ell$ is just one parameter. Then we want to solve$$\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n\ell\left(y_i,f(\omega_0+\mathbf{x}_i^T\mathbf{\omega})\right)\right\rbrace$$Thus, consider

 weights_1 = optim(weights_0,loss)$par (where the starting point is the OLS estimate). Again, to see what’s going on, let us visualize the ROC curve  y_5_2 = sigmoid(X %*% weights_1) pred = ROCR::prediction(y_5_2,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) plot(perf0,add=TRUE,col="red")

That’s not amazing, but again, that’s only a first step.

A single layer

Let us add a single layer in our network.

Those nodes are connected to the sources (incoming from sources) from the left, and then connected to the sink, on the right. Those nodes are not inter-connected. And again, for that network, we need edges (i.e series of weights). For instance, on the network above, we did add one single layer, with (only) three nodes.

For such a network, the prediction formula is $$\mathbf{y}=f\left( \omega_0+ \sum_{h=1}^3\omega_h f_h\left(\omega_{h,0}+ \sum_{j=1}^p \omega_{h,j} x_j\right)\right)$$or more synthetically$$\mathbf{y}=f\left( \omega_0+ \sum_{h=1}^3 \omega_hf_h\left(\omega_{h,0}+ \mathbf{x}^T\mathbf{\omega}_h\right)\right)$$Usually, we consider the same activation function everywhere. Don’t ask me why, I find that weird.

Now, we have a lot of weights to choose. Let us use again OLS estimates

 weights_1 <- lm(PRONO~1+FRCAR+INCAR+INSYS+PAPUL+PVENT,data=myocarde)$coefficients X1 = as.matrix(cbind(1,myocarde[,c("FRCAR","INCAR","INSYS","PAPUL","PVENT")])) weights_2 <- lm(PRONO~1+INSYS+PRDIA,data=myocarde)$coefficients X2=as.matrix(cbind(1,myocarde[,c("INSYS","PRDIA")])) weights_3 <- lm(PRONO~1+PAPUL+PVENT+REPUL,data=myocarde)$coefficients X3=as.matrix(cbind(1,myocarde[,c("PAPUL","PVENT","REPUL")])) In that case, we did specify edges, and which sources (explanatory variables) should be used for each additional node. Actually, here, other techniques could be have been used, like using a PCA. Each node will then be one of the components. But we’ll use that idea later one…  X = cbind(sigmoid(X1 %*% weights_1), sigmoid(X2 %*% weights_2), sigmoid(X3 %*% weights_3)) But we’re not done here. Those were weights from the source to the know nodes, in the layer. We still need the weights from the nodes to the sink. Here, let use use a simple average  weights = c(1/3,1/3,1/3) y_5_3 <- sigmoid(X %*% weights) Again, we can plot the ROC curve to see what we’ve done…  pred = ROCR::prediction(y_5_3,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) plot(perf0,add=TRUE,col="red")

On back propagation

Now, we need some optimal selection of those weights. Observe that with only 3 nodes, there are already $(7+1)\times3+3=27$ parameters in that model! Clearly, parcimony is not the major issue when you start using neural nets! If $$p(\mathbf{x})=f\left( \omega_0+ \sum_{h=1}^3 \omega_hf_h\left(\omega_{h,0}+ \mathbf{x}^T\mathbf{\omega}_h\right)\right)$$we want to solve$$\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n\ell\left(y_i,p(\mathbf{x}_i)\right)\right\rbrace$$for some loss function, which is$$\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n (y_i-p(\mathbf{x}_i))^2 \right\rbrace$$for the quadratic norm, or$$\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n (y_i\log p(\mathbf{x}_i)+[1-y_i]\log [1-p(\mathbf{x}_i)]) \right\rbrace$$if we want to use cross-entropy.

For convenience, let us center all the variable we create, otherwise, we get numerical problems.

 center = function(z) (z-mean(z))/sd(z) loss = function(weights){ weights_1 = weights[0+(1:7)] weights_2 = weights[7+(1:7)] weights_3 = weights[14+(1:7)] weights_ = weights[21+1:4] X1=X2=X3=as.matrix(myocarde[,1:7]) Z1 = center(X1 %*% weights_1) Z2 = center(X2 %*% weights_2) Z3 = center(X3 %*% weights_3) X = cbind(1,sigmoid(Z1), sigmoid(Z2), sigmoid(Z3)) mean( (myocarde$PRONO-sigmoid(X %*% weights_))^2)} Now that we have our objective function, consider some starting points. We can consider weights from a PCA, and then use a gradient descent algorithm,  pca = princomp(myocarde[,1:7]) W = get_pca_var(pca)$contrib weights_0 = c(W[,1],W[,2],W[,3],c(-1,rep(1,3)/3)) weights_opt = optim(weights_0,loss)$par The prediction is then obtained using  weights_1 = weights_opt[0+(1:7)] weights_2 = weights_opt[7+(1:7)] weights_3 = weights_opt[14+(1:7)] weights_ = weights_opt[21+1:4] X1=X2=X3=as.matrix(myocarde[,1:7]) Z1 = center(X1 %*% weights_1) Z2 = center(X2 %*% weights_2) Z3 = center(X3 %*% weights_3) X = cbind(1,sigmoid(Z1), sigmoid(Z2), sigmoid(Z3)) y_5_4 = sigmoid(X %*% weights_) And as previously, why not plot the ROC curve of that model  pred = ROCR::prediction(y_5_4,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) plot(perf,add=TRUE,col="red")

That’s not too bad. But with 27 coefficients, that’s what we would expect, no?

Using nnet() function

That’s more or less what is done in neural nets functions. Let us now have a look at some dedicated R functions.

 library(nnet) myocarde_minmax = myocarde minmax = function(z) (z-min(z))/(max(z)-min(z)) for(j in 1:7) myocarde_minmax[,j] = minmax(myocarde_minmax[,j])

Here, variables are linearly transformed, to take values in $(0,1)$. Then we can construct a neural network with one single layer, and three nodes,

 model_nnet = nnet(PRONO~.,data=myocarde_minmax,size=3) summary(model_nnet) a 7-3-1 network with 28 weights options were - b->h1 i1->h1 i2->h1 i3->h1 i4->h1 i5->h1 i6->h1 i7->h1 -9.60 -1.79 21.00 14.72 -20.45 -5.05 14.37 -17.37 b->h2 i1->h2 i2->h2 i3->h2 i4->h2 i5->h2 i6->h2 i7->h2 4.72 2.83 -3.37 -1.64 1.49 2.12 2.31 4.00 b->h3 i1->h3 i2->h3 i3->h3 i4->h3 i5->h3 i6->h3 i7->h3 -0.58 -6.03 25.14 18.03 -1.19 7.52 -19.47 -12.95 b->o h1->o h2->o h3->o -1.32 29.00 -10.32 26.27

Here, it is the complete full network. And actually, there are (online) some functions that can he used to visualize that network

 library(devtools) source_url('https://gist.githubusercontent.com/fawda123/7471137/raw/466c1474d0a505ff044412703516c34f1a4684a5/nnet_plot_update.r') plot.nnet(model_nnet)

Nice, isn’t it? We clearly see the intermediary layer, with three nodes, and on top the constants. Edges are the plain lines, the darker, the heavier (in terms of weights).

Using neuralnet()

Other R functions can actually be considered.

 library(neuralnet) model_nnet = neuralnet(formula(glm(PRONO~.,data=myocarde_minmax)), myocarde_minmax,hidden=3, act.fct = sigmoid) plot(model_nnet)

Again, for the same network structure, with one (hidden) layer, and three nodes in it.

Network with multiple layers

The good thing is that it’s not possible to add more layers. Like two layers. Nodes from the first layer are no longuer connected with the sink, but with nodes in the second layer. And those nodes will then be connected to the sink. We now have something like
$$p(\mathbf{x})=f\left( \omega_0+ \sum_{h=1}^3 \omega_h f_h\left(\omega_{h,0}+ \mathbf{z}_h^T\mathbf{\omega}_h\right)\right)$$where$$\mathbf{z}_h=f\left( \omega_{h,0}+ \sum_{j=1}^{k_h} \omega_{h,j} f_{h,j}\left(\omega_{h,j,0}+ \mathbf{x}^T\mathbf{\omega}_{h,j}\right)\right)$$I may be rambling here (a little bit) but that’s a lot of parameters. Here is the visualization of such a network,

 library(neuralnet) model_nnet = neuralnet(formula(glm(PRONO~.,data=myocarde_minmax)), myocarde_minmax,hidden=3, act.fct = sigmoid) plot(model_nnet)

Application

Let us get back on our simple dataset, with only two covariates.

 library(neuralnet) df_minmax =df df_minmax$y=(df_minmax$y=="1")*1 minmax = function(z) (z-min(z))/(max(z)-min(z)) for(j in 1:2) df_minmax[,j] = minmax(df[,j]) X = as.matrix(cbind(1,df_minmax[,1:2]))

Consider only one layer, with two nodes

 model_nnet = neuralnet(formula(lm(y~.,data=df_minmax)), df_minmax,hidden=c(2)) plot(model_nnet)

Here, we did not specify it, but the activation function is the sigmoid (actually, it is called logistic here)

 model_nnet$act.fct function (x) { 1/(1 + exp(-x)) } attr(,"type") [1] "logistic" f=model_nnet$act.fct

The weights (on the figure) can be obtained using

 w0 = model_nnet$weights[[1]][[2]][,1] w1 = model_nnet$weights[[1]][[1]][,1] w2 = model_nnet$weights[[1]][[1]][,2] Now, to get our prediction, we should use$$p(\mathbf{x})=f\left( \omega_0+ \omega_1 f(\omega_{1,0}+ \mathbf{x}_h^T\mathbf{\omega}_{1,1:2})+\omega_1 f(\omega_{2,0}+ \mathbf{x}_h^T\mathbf{\omega}_{2,1:2})\right)$$which can be obtained using  f(cbind(1,f(X%*%w1),f(X%*%w2))%*%w0) [,1] [1,] 0.7336477343 [2,] 0.7317999050 [3,] 0.7185803540 [4,] 0.7404005280 [5,] 0.7518482779 [6,] 0.4939774149 [7,] 0.4965876378 [8,] 0.7101714888 [9,] 0.5050760026 [10,] 0.5049877644 Unfortunately, it is not the output of the model here,  neuralnet::prediction(model_nnet) Data Error: 0;$rep1 x1 x2 y 1 0.1250 0.0000000000 0.02030470787 2 0.0625 0.1176470588 0.89621706711 3 0.9375 0.2352941176 0.01995171956 4 0.0000 0.4705882353 1.10849420363 5 0.5000 0.4705882353 -0.01364966058 6 0.3125 0.5294117647 -0.02409150561 7 0.6875 0.8235294118 0.93743057765 8 0.3750 0.8823529412 1.01320924782 9 1.0000 0.9058823529 1.04805134309 10 0.5625 1.0000000000 1.00377379767

If anyone has a clue, I’d be glad to know what went wrong here… I find that odd to have outputs outside the $(0,1)$ interval, but the output is neither$$p(\mathbf{x})=\omega_{0,0}+ \omega_{0,1} f(\omega_{1,0}+ \mathbf{x}_h^T\mathbf{\omega}_{1,1:2})+\omega_{0,2} f(\omega_{2,0}+ \mathbf{x}_h^T\mathbf{\omega}_{2,1:2})$$

 cbind(1,f(X%*%w1),f(X%*%w2))%*%w0 [,1] [1,] 1.01320924782 [2,] 1.00377379767 [3,] 0.93743057765 [4,] 1.04805134309 [5,] 1.10849420363 [6,] -0.02409150561 [7,] -0.01364966058 [8,] 0.89621706711 [9,] 0.02030470787 [10,] 0.01995171956

Supervised Classification, beyond the logistic

In our data-science class, after discussing limitations of the logistic regression, e.g. the fact that the decision boundary line was a straight line, we’ve mentioned possible natural extensions. Let us consider our (now) standard dataset

clr1 <- c(rgb(1,0,0,1),rgb(0,0,1,1))
clr2 <- c(rgb(1,0,0,.2),rgb(0,0,1,.2))
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(x,y,z)
plot(x,y,pch=19,cex=2,col=clr1[z+1])

One can consider a quadratic function of the covariates (instead of a linear one)

reg=glm(z~x+y+I(x^2)+I(y^2)+I(x*y),
data=df,family=binomial)
summary(reg)

pred_1 <- function(x,y){
predict(reg,newdata=data.frame(x=x,
y=y),type="response")>.5 }

x_grid<-seq(0,1,length=101)
y_grid<-seq(0,1,length=101)
z_grid <- outer(x_grid,y_grid,pred_1)
image(x_grid,y_grid,z_grid,col=clr2)
points(x,y,pch=19,cex=2,col=clr1[z+1])

Supervised Classification, Logistic and Multinomial

We will start, in our Data Science course,  to discuss classification techniques (in the context of supervised models). Consider the following case, with 10 points, and two classes (red and blue)

> clr1 <- c(rgb(1,0,0,1),rgb(0,0,1,1))
> clr2 <- c(rgb(1,0,0,.2),rgb(0,0,1,.2))
> 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(x,y,z)
> plot(x,y,pch=19,cex=2,col=clr1[z+1])

To get a prediction, i.e. a partition of the space in two parts, consider some logistic regression

> reg=glm(z~x+y,data=df,family=binomial)
> summary(reg)

Call:
glm(formula = z ~ x + y, family = binomial, data = df)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.6593  -0.4400   0.2564   0.5830   1.5374

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.706      1.999  -0.854    0.393
x             -5.489      5.360  -1.024    0.306
y              8.568      5.515   1.554    0.120

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 13.4602  on 9  degrees of freedom
Residual deviance:  8.1445  on 7  degrees of freedom
AIC: 14.144

Number of Fisher Scoring iterations: 5

Given some point, the predicted class is obtained using

> pred_1 <- function(x,y){
+ predict(reg,newdata=data.frame(x=x,
+ y=y),type="response")>.5
+ }

(here, the predicted class is simply the one that is the most likely). To visualize it use

> x_grid<-seq(0,1,length=101)
> y_grid<-seq(0,1,length=101)
> z_grid <- outer(x_grid,y_grid,pred_1)
> image(x_grid,y_grid,z_grid,col=clr2)
> points(x,y,pch=19,cex=2,col=clr1[z+1])

Since the logistic regression is a (generalized) linear model, the line that separate the two regions is a straight line.

Modeling individual losses with mixtures

Usually, the sentence that I keep saying in my regression classes is “please, look at your data“. In our previous post, we’ve been playing like most econometricians: we did not look at the data. Actually, if we look at the distribution of individual losses, in the dataset, we see the following,

> n=nrow(couts)
> plot(sort(couts$cout),(1:n)/(n+1),xlim=c(0,10000),type="s",lwd=2,col="green") It looks like there are fixed costs claims in our database. How do we deal with it in the standard case (e.g. in Loss Models textbook) ? We can use a mixture of – at least – three distributions here, $f(y\boldsymbol) = p_1 {\color{Blue} f_1(}y{\color{Blue} )} + p_2 {\color{Magenta} \delta_{\kappa}(}y{\color{Magenta} )} + p_3 {\color{Red} f_3(}y{\color{Red} )}$ with • a distribution for small claims, ${\color{Blue} f_1(}\cdot{\color{Blue} )}$, e.g. an exponential distribution • a Dirac mass in ${\color{Magenta} \kappa}$, i.e. ${\color{Magenta} \delta_{\kappa}(}\cdot{\color{Magenta} )}$ • a distribution for larger claims, ${\color{Red} f_3(}\cdot{\color{Red} )}$, e.g. a Gamma, or a lognormal, distribution > I1=which(couts$cout<1120)
>  I2=which((couts$cout>=1120)&(couts$cout<1220))
>  I3=which(couts$cout>=1220) > (p1=length(I1)/nrow(couts)) [1] 0.3284823 > (p2=length(I2)/nrow(couts)) [1] 0.4152807 > (p3=length(I3)/nrow(couts)) [1] 0.256237 > X=couts$cout
>  (kappa=mean(X[I2]))
[1] 1171.998
>  X0=X[I3]-kappa
>  u=seq(0,10000,by=20)
>  F1=pexp(u,1/mean(X[I1]))
>  F2= (u>kappa)
>  F3=plnorm(u-kappa,mean(log(X0)),sd(log(X0))) * (u>kappa)
>  F=F1*p1+F2*p2+F3*p3
>  lines(u,F)

In our previous post, we’ve discussed the idea that all parameters might be related to some covariates, i.e.

$f(y|\boldsymbol{X}) = p_1(\boldsymbol{X}) {\color{Blue} f_1(}y|\boldsymbol{X}{\color{Blue} )} + p_2(\boldsymbol{X}) {\color{Magenta} \delta_{\kappa}(}y{\color{Magenta} )} + p_3(\boldsymbol{X}) {\color{Red} f_3(}y|\boldsymbol{X}{\color{Red} )}$

which yield the following premium model,

$\mathbb{E}(Y|\boldsymbol{X}) = {\color{Blue} {\underbrace{\mathbb{E}(Y|\boldsymbol{X},Y\leq s_1)}_{A} \cdot {\underbrace{\mathbb{P}(Y\leq s_1|\boldsymbol{X})}_{D}}}}\\+{\color{Purple} {{\underbrace{\mathbb{E}(Y|Y\in( s_1,s_2], \boldsymbol{X}) }_{B}}\cdot {\underbrace{\mathbb{P}(Y\in( s_1,s_2]| \boldsymbol{X})}_{D}}}}\\+{\color{Red} {{\underbrace{\mathbb{E}(Y|Y> s_2, \boldsymbol{X}) }_{C}}\cdot {\underbrace{\mathbb{P}(Y> s_2| \boldsymbol{X})}_{D}}}}$

For the ${\color{Blue} A}$${\color{Magenta} B}$ and ${\color{Red} C}$ terms, that’s easy, we can use standard models we’ve seen in the course. For the probability, we should use a multinomial model. Recall that for the logistic regression model, if $(\pi,1-\pi)=(\pi_1,\pi_2)$, then

$\log \frac{\pi}{1-\pi}=\log \frac{\pi_1}{\pi_2} =\boldsymbol{X}'\boldsymbol{\beta}$

i.e.

$\pi_1 = \frac{\exp(\boldsymbol{X}'\boldsymbol{\beta})}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta})}$

and

$\pi_2 = \frac{1}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta})}$

To derive a multivariate extension, write

$\pi_1 = \frac{\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)+\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)}$

$\pi_2 = \frac{\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)+\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)}$

and

$\pi_3 = \frac{1}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)+\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)}$

Again, maximum likelihood techniques can be used, since

$\mathcal{L}(\boldsymbol{\pi},\boldsymbol{y})\propto \prod_{i=1}^n \prod_{j=1}^3 \pi_{i,j}^{Y_{i,j}}$

where here, variable $Y_{i}$  – which take three levels – is splitted in three indicators (like any categorical explanatory variables in standard regression model). Thus,

$\log \mathcal{L}(\boldsymbol{\beta},\boldsymbol{y})\propto \sum_{i=1}^n \sum_{j=1}^2 \left(Y_{i,j} \boldsymbol{X}_i'\boldsymbol{\beta}_j\right) - n_i\log\left[1+1+\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)+\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)\right]$

and, as for the logistic regression, then use Newton Raphson’ algorithm to compute numerically the maximum likelihood. In R, first we have to define the levels, e.g.

> seuils=c(0,1120,1220,1e+12)
> couts$tranches=cut(couts$cout,breaks=seuils,
+ labels=c("small","fixed","large"))
nocontrat    no garantie    cout exposition zone puissance agevehicule
1      1870 17219      1RC 1692.29       0.11    C         5           0
2      1963 16336      1RC  422.05       0.10    E         9           0
3      4263 17089      1RC  549.21       0.65    C        10           7
4      5181 17801      1RC  191.15       0.57    D         5           2
5      6375 17485      1RC 2031.77       0.47    B         7           4
ageconducteur bonus marque carburant densite region tranches
1            52    50     12         E      73     13    large
2            78    50     12         E      72     13    small
3            27    76     12         D      52      5    small
4            26   100     12         D      83      0    small
5            46    50      6         E      11     13    large

Then, we can run a multinomial regression, from

> library(nnet)

using some selected covariates

> reg=multinom(tranches~ageconducteur+agevehicule+zone+carburant,data=couts)
# weights:  30 (18 variable)
initial  value 2113.730043
iter  10 value 2063.326526
iter  20 value 2059.206691
final  value 2059.134802
converged

The output is here

> summary(reg)
Call:
multinom(formula = tranches ~ ageconducteur + agevehicule + zone +
carburant, data = couts)

Coefficients:
(Intercept) ageconducteur agevehicule      zoneB      zoneC
fixed  -0.2779176   0.012071029  0.01768260 0.05567183 -0.2126045
large  -0.7029836   0.008581459 -0.01426202 0.07608382  0.1007513
zoneD      zoneE      zoneF   carburantE
fixed -0.1548064 -0.2000597 -0.8441011 -0.009224715
large  0.3434686  0.1803350 -0.1969320  0.039414682

Std. Errors:
(Intercept) ageconducteur agevehicule     zoneB     zoneC     zoneD
fixed   0.2371936   0.003738456  0.01013892 0.2259144 0.1776762 0.1838344
large   0.2753840   0.004203217  0.01189342 0.2746457 0.2122819 0.2151504
zoneE     zoneF carburantE
fixed 0.1830139 0.3377169  0.1106009
large 0.2160268 0.3624900  0.1243560

To visualize the impact of a covariate (one, only), one can use also spline functions

> library(splines)
> reg=multinom(tranches~agevehicule,data=couts)
# weights:  9 (4 variable)
initial  value 2113.730043
final  value 2072.462863
converged
> reg=multinom(tranches~bs(agevehicule),data=couts)
# weights:  15 (8 variable)
initial  value 2113.730043
iter  10 value 2070.496939
iter  20 value 2069.787720
iter  30 value 2069.659958
final  value 2069.479535
converged

For instance, if the covariate is the age of the car, we do have the following probabilities

> predict(reg,newdata=data.frame(agevehicule=5),type="probs")
small     fixed     large
0.3388947 0.3869228 0.2741825

and for all ages from 0 to 20,

For instance, for new cars, the proportion of fixed costs is rather small (here in purple), and keeps increasing with the age of the car. If the covariate is the density of population in the area the driver lives, we do obtain the following probabilities

> reg=multinom(tranches~bs(densite),data=couts)
# weights:  15 (8 variable)
initial  value 2113.730043
iter  10 value 2068.469825
final  value 2068.466349
converged
> predict(reg,newdata=data.frame(densite=90),type="probs")
small     fixed     large
0.3484422 0.3473315 0.3042263

Based on those probabilities, it is then possible to derive the expected cost of a claims, given some covariates (e.g. the density). But first, define subsets of the whole dataset

> sbaseA=couts[couts$tranches=="small",] > sbaseB=couts[couts$tranches=="fixed",]
> sbaseC=couts[couts$tranches=="large",] with a threshold given by > (k=mean(sousbaseB$cout))
[1] 1171.998

Then, let us run our four models,

> reg=multinom(tranches~bs(densite),data=couts)

We can now compute predictions based on those models,

> nouveau=data.frame(densite=seq(10,100))
> proba=predict(reg,newdata=nouveau,type="probs")
> predA=predict(regA,newdata=nouveau,type="response")
> predB=predict(regB,newdata=nouveau,type="response")
> predC=predict(regC,newdata=nouveau,type="response")+k
> pred=cbind(predA,predB,predC)

To visualize the impact of each component on the premium, we can compute probabilities, are well as expected costs (given a cost in each subset),

> cbind(proba,pred)[seq(10,90,by=10),]
small     fixed     large    predA    predB    predC
10 0.3344014 0.4241790 0.2414196 423.3746 1171.998 7135.904
20 0.3181240 0.4471869 0.2346892 428.2537 1171.998 6451.890
30 0.3076710 0.4626572 0.2296718 438.5509 1171.998 5499.030
40 0.3032872 0.4683247 0.2283881 451.4457 1171.998 4615.051
50 0.3052378 0.4620219 0.2327404 463.8545 1171.998 3961.994
60 0.3136136 0.4417057 0.2446807 472.3596 1171.998 3586.833
70 0.3279413 0.4056971 0.2663616 473.3719 1171.998 3513.601
80 0.3464842 0.3534126 0.3001032 463.5483 1171.998 3840.078
90 0.3652932 0.2868006 0.3479061 440.4925 1171.998 4912.379

Now, it is possible to plot those figures in a graph,

> barplot(t(proba*pred))
> abline(h=mean(couts\$cout),lty=2)

(the dotted horizontal line is the average cost of a claim, in our dataset).