# From multinomial regression to binary classification on some Siamese data

There are two kinds of people in the world: people who think there are two kinds of people in the world and people who don’t

(borrowed from Menand (2018)). Because things are always simpler when we face only binary choice, aren’t they? But consider here the case were multiple options are possible, and let us see if we cannot get back to simpler binary choices. Consider a collection of observations $(y_i,\boldsymbol{x}_i)$ where $y_i$ is some categorical variable, $y_i\in\mathcal{A}$ where $\mathcal{A}=\lbrace A_1,\cdots,A_\kappa \rbrace$, with $\kappa$ possible categories. Let $\mathcal{I}_k=\lbrace i:y_i\in A_k \rbrace$.

In a classical multinomial logistic regression, suppose that $A_1$ is the reference, then $$\mathbb{P}[Y=A_j|\boldsymbol{X}=\boldsymbol{x}]=\frac{\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_j]}{1+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_2]+\cdots+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k]}$$With a lot a categories, and a small number of observations, inference can be complicated, and non-robust.

• the Siamese dataset

The name Siamese I use, here, comes from Siamese Networks. Or sort of… As we say in French, it is an « histoire de l’homme qui a vu l’homme qui a vu l’ours » (story of the man who saw the man who saw the bear). A few years ago, a student tried to explain to me the idea of Siamese Networks and this is what I understood. I might be completely wrong, but the idea I got from it did make sense, in my mind at least. That is the story of that blog post…

The idea of the siamese algorithm will be to consider all pairs of observations, $(y_i,\boldsymbol{x}_i)$ and $(y_j,\boldsymbol{x}_j)$ :

1. $\tilde y_{i,j}=\boldsymbol{1}(y_i=y_j)$ indicating if individuals $i$ and $j$ are in the same category
2. $\tilde{\boldsymbol{x}}_{i,j}$ is a collection of $p-1$ variables,
• $\tilde {x}_{k:i,j}={x}_{k:i}-{x}_{k:j}$ if $x_k$ is continuous, or $\tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}|$ (we can use another metric, e.g. $\tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}|^2$, and this is why I decided to use some GAM model in the logistic regression on the Siamese dataset)
• $\tilde {x}_{k:i,j}=({x}_{k:i},{x}_{k:j})\in\mathcal{X}_k\times\mathcal{X}_k$ if $x_k$ is a categorical variables (taking values in the set $\mathcal{X}_k$), or $\tilde {x}_{k:i,j}=\boldsymbol{1}({x}_{k:i}\neq{x}_{k:j})\in\{0,1\}$

The original dataset was a $n\times p$ matrix, and (if there are no categorical variable), it becomes a $n(n-1)/2\times p$ matrix. The key point is that if the original variable $y_i$ was multinomial, $y_{i,j}$ is now binomial. For instance, if our initial dataset was the following, with two covariates, one continuous and one categorical

its siamese counterpart is the following

• Classification step

On the dataset $(\tilde y_{i,j},\tilde {x}_{k:i,j})_{i,j}$, fit a logistic regression, $$\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}]=\frac{\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}$$(or any classification model – CART, random forest, etc). But that is the easy part (unless $n$ is large, because the siamese dataset has (roughly) $n^2/2$ rows). The difficult task is the prediction

• Prediction step

Consider a new input variable $\boldsymbol{x}_{\cdot}$, and define its siamese version, $\tilde{\boldsymbol{x}}_{\cdot}=(\tilde{\boldsymbol{x}}_{\cdot,j})_j$, i.e. a database with $n$ rows. Then compute
$$p_{\cdot,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{\cdot,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]}$$ where $p_{\cdot,j}$ is the probability that $(y_j,\boldsymbol{x}_{j})$ and $(y_{\cdot},\boldsymbol{x}_{\cdot})$ are in the same category, as well as $$p_{i,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{i,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}$$
Let $\boldsymbol{p}_{\cdot}=(p_{\cdot,j})$, and similarly $\boldsymbol{p}_{i}=(p_{i,j})$. Then several techniques can be used to predict $y_{\cdot}$.

1. $\widehat{y}_{\cdot}=y_{j^\star}$ where ${j^\star}=\underset{j=1,\cdots,n}{\text{argmax}}\{p_{\cdot,j}\}$: the predicted class is the one of the observation the most likely to be other same class
2. $\widehat{y}_{\cdot}=y_{j^\star}$ where ${j^\star}=\underset{\ell=1,\cdots,k}{\text{argmax}}\{\overline{p}_{\ell}\}$, where$$\overline{p}_{\ell} = \frac{1}{n_{\ell}}\sum_{j\in\mathcal{I}_s} \boldsymbol{1}(y_j =y_{\ell}),\text{ where }\mathcal{I}_s=\lbrace i:p_{\cdot,i}>s\rbrace$$consider only probabilities sufficiently high to be considered, and predict the most important class (majority rule)
3. $\widehat{y}_{\cdot}=y_{j^\star}$ where $${j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{\theta_i\}$$ where $$\theta_i=\cos(\boldsymbol{p}_{\cdot},\boldsymbol{p}_{i})=\displaystyle{\frac{\boldsymbol{p}_{\cdot}\cdot\boldsymbol{p}_{i}}{\|\boldsymbol{p}_{\cdot}\|\|\boldsymbol{p}_{i}\|}}$$
4. $\widehat{y}_{\cdot}=y_{j^\star}$ where $${j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{KL_{\cdot|i}\}$$ and $KL_{\cdot|i}=\displaystyle{\sum_{j=1}^n p_{\cdot,j}\log\frac{p_{\cdot,j}}{p_{i,j}}}$ (but one can select another metric)
5. $\widehat{y}_{\cdot}=y_{j^\star}$ where ${j^\star}=\underset{j\in\mathcal{J}}{\text{argmax}}\{p_{\cdot,j}\}$ and $\mathcal{J}$ is a sample of $k$ observations, chosen randomly, one in each group (one-shot procedure): the predicted class is the one of the observation the most likely to be other same class

Heuristically, it can be related to some $k$ nearest neighbors strategy: we give the attribute that most neighbors have. The total distance is a weighted sum of the componentwise distances (for the logistic regression).

• Simulation study

In order to test that technique, let us generate some multinomial model where $y$ has 10 possible labels, with 6 (independent) covariates $x_1,\cdots,x_6$, and $\mathbb{P}[Y=A_k|\boldsymbol{X}=\boldsymbol{x}]\propto \exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k]$ (where coefficients $\boldsymbol{\beta}_k$ where generated randomly) for $k\in\{1,2,\cdots,10\}$ (there were 10 categories) and with $n=700$ observations.

n=700 X1=rnorm(n) X2=rnorm(n) X3=rnorm(n) X4=rnorm(n) X5=rnorm(n) X6=rnorm(n) X=cbind(1,X1,X2,X3,sqrt(abs(X4)),X5*X1,X6) k = 10 PARAM = matrix(rnorm(k*6),k,6) PARAM[,1]=PARAM[,1]-1 PARAM=cbind(PARAM,0) P=matrix(NA,n,k-1) for(j in 1:(k-1)) P[,j] = X %*% (PARAM[j,])+rnorm(n) P=cbind(P,0) S=apply(exp(P),1,sum) Pb = exp(P)/S tirage = function(i){ sample(1:10,size=1,prob = Pb[i,]) } Y = LETTERS[Vectorize(tirage)(1:n)] dbase = data.frame(Y=as.factor(Y),X1,X2,X3,X4,X5)

In the paragraph previously, I suggested to take the most likely one. Being wrong means that it was not the first choice. But perhaps being the second or the third choice is not that bad, actually. So in my simulations, I look at the proportion of predictions were our prediction is the good one (top 1), if the true one is either the most likely or the second most likely (top 2), or in the top 3. That will be on my $x$-axis. I draw some line, but we simply have three points (top 1, top 2 and top 3). I compute the proportion of good prediction, using cross-validation techniques (10-fold). The black lines are one of the methods described above. The red one is the standard multinomial model (with a logistic link function). For the Siamese model, I tried several models. I tried a logistic regression, and some smooth version (GAM) on top

and a classification tree, on the left, as well as some random forest on the right, below.

It looks like the multinomial approach performs always better than any Siamese one… and to be honest, I am disappointed.

Here is the code I did use when I considered a logistic regression on the Siamese dataset,

set.seed(1) kfold = sample(rep(1:10,n/10))   KFOLDglm = function(i){ i_test=which(kfold==i) i_calibration=which(kfold!=i) y=credit[i_calibration,"Y"] tirage = function(){ v=c(sample(i_calibration[y==levels(y)[1]],size=1), sample(i_calibration[y==levels(y)[2]],size=1), sample(i_calibration[y==levels(y)[3]],size=1), sample(i_calibration[y==levels(y)[4]],size=1), sample(i_calibration[y==levels(y)[5]],size=1), sample(i_calibration[y==levels(y)[6]],size=1), sample(i_calibration[y==levels(y)[7]],size=1), sample(i_calibration[y==levels(y)[8]],size=1), sample(i_calibration[y==levels(y)[9]],size=1), sample(i_calibration[y==levels(y)[10]],size=1)) names(v)=levels(y) return(v) }   LogisticModel &lt;- multinom(Y ~ ., data = credit[i_calibration,], trace=FALSE)   comparaisonx = function(base,x=base[1,]){ mix_base = base for(j in 1:ncol(base)){ xj = as.numeric(x[j])-base[,j] mix_base[,j] = (xj) } mix_base } comparaisony = function(base,y=base[1]){ as.factor(base == y) } creditx = credit[,-which(names(credit) == "Y")] nc=length(i_calibration) B=comparaisonx(base = creditx[i_calibration[2:nc],],x=creditx[i_calibration[1],]) B$Y=comparaisony(base = credit[i_calibration[2:nc],"Y"],y=credit[i_calibration[1],"Y"]) for(i in 2:(nc-1)){ B0=comparaisonx(base = creditx[i_calibration[(i+1):nc],],x=creditx[i_calibration[i],]) B0$Y=comparaisony(base = credit[i_calibration[(i+1):nc],"Y"],y=credit[i_calibration[i],"Y"]) B=rbind(B,B0) } credit_mix = B   OneShotLogisticModel &lt;- glm(Y ~ ., data = credit_mix, family=binomial) A_ref = table(credit[i_calibration,"Y"])/length(i_calibration)   vect_oneshot = function(i){ B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,]) predict(OneShotLogisticModel,type="response",newdata=B2) }   prediction_oneshot = function(i,type=1){ B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,]) p=predict(OneShotLogisticModel,type="response",newdata=B2) y=credit[i_calibration,"Y"] base = data.frame(p,y) base = base[rev(order(base$p)),] if(type==1){T = table(base$y[1:11]) return(names(which.max(T)))} if(type==2){return(base$y[1])} if(type==3){A=table(base$y[1:10])/10 T=A/A_ref return(names(which.max(T)))} if(type==4){ costheta = rep(NA,length(i_calibration)) for(j in 1:length(i_calibration)){ vecteur_proba = vect_oneshot(i_calibration[j]) costheta[j] = sum(vecteur_proba*p)/(sqrt(sum(vecteur_proba^2))*sqrt(sum(p^2))) } return(y[which.max(costheta)])} if(type==5){ kl = rep(NA,length(i_calibration)) for(j in 1:length(i_calibration)){ vecteur_proba = vect_oneshot(i_calibration[j]) kl[j] = sum(p*log(vecteur_proba/p)) } return(y[which.max(as.vector(kl))])} if(type==6){ ## one shot : tirer au hasard un de chaque, et dire lequel est plus credible !   y=credit[i_calibration,"Y"] tirage = function(){ v=c(sample(i_calibration[y==levels(y)[1]],size=1), sample(i_calibration[y==levels(y)[2]],size=1), sample(i_calibration[y==levels(y)[3]],size=1), sample(i_calibration[y==levels(y)[4]],size=1), sample(i_calibration[y==levels(y)[5]],size=1), sample(i_calibration[y==levels(y)[6]],size=1), sample(i_calibration[y==levels(y)[7]],size=1), sample(i_calibration[y==levels(y)[8]],size=1), sample(i_calibration[y==levels(y)[9]],size=1), sample(i_calibration[y==levels(y)[10]],size=1)) names(v)=levels(y) return(v) } pd=rep(NA,101) for(ix in 1:101){ ids = tirage() B2=comparaisonx(base = creditx[ids,],x=creditx[i,]) p=predict(OneShotLogisticModel,type="response",newdata=B2) pd[ix]=levels(y)[which.max(p)] } levels(y)[which.max(table(pd))] } } PRED0=as.character(credit[i_test,"Y"]) PRED1=as.character(predict(LogisticModel,type = "class", newdata=credit[i_test,])) PRED21=as.character(Vectorize(function(i) prediction_oneshot(i,type=1)) (i_test)) PRED22=as.character(Vectorize(function(i) prediction_oneshot(i,type=2))(i_test)) PRED23=as.character(Vectorize(function(i) prediction_oneshot(i,type=3))(i_test)) PRED24=as.character(Vectorize(function(i) prediction_oneshot(i,type=4))(i_test)) PRED25=as.character(Vectorize(function(i) prediction_oneshot(i,type=5))(i_test)) PRED26=as.character(Vectorize(function(i) prediction_oneshot(i,type=6))(i_test)) B=data.frame(PRED0,PRED1,PRED21,PRED22,PRED23,PRED24,PRED25,PRED26) B }   s=1/100;setTxtProgressBar(pb, s*2) for(i in 2:10){ PREDICTION = rbind(PREDICTION,KFOLDglm(i)) s=s+1/100;setTxtProgressBar(pb, s*2)} for(j in 1:5) PREDICTION[,j]=as.character(PREDICTION[,j]) L=list() v=mean(PREDICTION[,1]!=PREDICTION[,2]) names(v)="logistic" L[["logistic"]]=v v=c(mean(PREDICTION[,1]!=PREDICTION[,3]), mean(PREDICTION[,1]!=PREDICTION[,4]), mean(PREDICTION[,1]!=PREDICTION[,5]), mean(PREDICTION[,1]!=PREDICTION[,6]), mean(PREDICTION[,1]!=PREDICTION[,7]), mean(PREDICTION[,1]!=PREDICTION[,8])) names(v)=c("top10","max","10norm","cos","KL","OS") L[["glm"]]=v

# Multinomial Logit as an Iterated Logit Regression

For the second section of the course at ENSAE, yesterday, we’ve seen how to run a multinomial logistic regression model. It is simply an extension of the binomial logistic regression. But actually, it is also possible to consider iterative binomial regressions.

Consider here a response variable $Y$ with a multinomial distribution (3 factors to have something more general than the binomial), taking values $\{A,B,C\}$, with respective probabilities $\mathbf{p}=(p_A,p_B,p_C)$. Here is a code to generate some multinomial variables

msample=function(A,B,C){ Y=rep(NA,B) for(i in 1:B){Y[i]=sample(A,size=1,prob=C[i,])} return(Y) }

and here is a code to generate a dataset with $n$ rows,

generate3=function(n,x,pb=c(-2,0)){ set.seed(x) X1=runif(n) X2=runif(n) X3=runif(n) s1=pb[1]+X1+X2 s2=pb[2]-X1+X2 P1=exp(s1)/(1+exp(s1)+exp(s2)) P2=exp(s2)/(1+exp(s1)+exp(s2)) Y=msample(0:2,n,cbind(1-P1-P2,P1,P2)) df=data.frame(Y=Y,X1=X1,X2=X2,X3=X3) return(df) }

Let us generate a training dataset and a validation one

 pb=c(.31,.42) DF1=generate3(1000,1,pb=pb) DF2=generate3(500,2,pb=pb)

With a multivariate logistic regression
$$\mathbb{P}[Y=A|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}$$
$$\mathbb{P}[Y=B|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}$$
$$\mathbb{P}[Y=B|\mathbf{x}]=\frac{1}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}$$

For convenience, consider the most popular factor in our training dataset

modalite=names(sort(table(DF1$Y),decreasing = TRUE)) Consider a regression model on the simulated dataset (with several covariates), let us estimate it, and let us get predictions.  library(nnet) reg=multinom(as.factor(Y) ~ ., data = DF1) mp1=predict (reg, DF1, "probs") mp2=predict (reg, DF2, "probs") An alternative can be the following. consider a first regression model on the Bernoulli variable $Y_A=\mathbf{1}(Y=A)$. Actually, we will consider the most important factor, but for convenience, assume that it is $A$. $$\mathbb{P}[Y_A=A|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{a}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{a}]}$$ On our dataset, estimate that model, and get predictions. In the case where $Y\neq A$, define another Bernoulli variable $Y_B=\mathbf{1}(Y=B|Y\neq A)$. We can estimate that model and derive two probabilities, $\mathbb{P}(Y=B|Y\neq A)$ and $\mathbb{P}(Y=C|Y\neq A)$ (the sum of the two being equal to 1). Based on those two models, it is possible to compute the three probabilities we are looking for. $\mathbb{P}[Y=A]$ is obtained from the first model, and we can derive the other two from $\mathbb{P}[Y=B|Y\neq A]\cdot\mathbb{P}[Y\neq A]$ and $\mathbb{P}[Y=C|Y\neq A]\cdot\mathbb{P}[Y\neq A]$.  reg1=glm((Y==modalite[1])~.,data=DF1,family=binomial) reg2=glm((Y==modalite[2])~.,data=DF1[-which(DF1$Y==modalite[1]),],family=binomial) p11=predict (reg1, newdata=DF1, type="response") p12=predict (reg2, newdata=DF1, type="response") p21=predict (reg1, newdata=DF2, type="response") p22=predict (reg2, newdata=DF2, type="response") mmp1=cbind(p11,(1-p11)*p12,(1-p11)*(1-p12)) mmp2=cbind(p21,(1-p21)*p22,(1-p21)*(1-p22)) colnames(mmp1)=colnames(mmp2)=modalite

Let us compare the predicted probabilites, on the same dataset (here the training dataset)

> mmp1[1:9,c("0","1","2")] 0 1 2 1 0.19728737 0.4991805 0.3035321 2 0.17244580 0.5648537 0.2627005 3 0.19291753 0.5971058 0.2099767 4 0.09087176 0.7787304 0.1303978 5 0.23400225 0.4083022 0.3576955 6 0.18063647 0.6637352 0.1556283 7 0.13188881 0.7402710 0.1278401 8 0.13776970 0.6524959 0.2097344 9 0.12325864 0.6790336 0.1977078 > mp1[1:9,c("0","1","2")] 0 1 2 1 0.19691036 0.5022692 0.3008205 2 0.17123189 0.5680647 0.2607034 3 0.19293066 0.5984402 0.2086291 4 0.08821851 0.7813318 0.1304497 5 0.23470739 0.4109990 0.3542936 6 0.18249687 0.6602168 0.1572863 7 0.13128711 0.7400898 0.1286231 8 0.13525341 0.6553618 0.2093848 9 0.12090016 0.6815915 0.1975084 

The two are very close. So yes, it is possible to see the multinomial regression as some sequential binomial regressions.

# Applications of Chi-Square Tests

This morning, in our mathematical statistical class, we’ve seen the use of the chi-square test. The first one was related to some goodness of fit of a multinomial distribution. Assume that $\boldsymbol{N}=(N_1,\cdots,N_k)\sim\mathcal{M}(n,\boldsymbol{p})$. In order to test $H_0:\boldsymbol{p}=\boldsymbol{p}_0$ against $H_1:\boldsymbol{p}\neq\boldsymbol{p}_0$, use the statistic

$Q=\sum_{j=1}^k \frac{[N_j-np_{0,j}]^2}{np_{0,j}}$Under $H_0$$Q\sim\chi^2(k-1)$. For instance, we have the number of weddings, in a large city, per season,

> n=c(301,356,413,262)

We want to test if weddings are celebrated uniformely over the year, i.e. $H_0:\boldsymbol{p}=\boldsymbol{1}/4$.

> np=rep(sum(n)/4,4)
> cbind(n,np)
n  np
[1,] 301 333
[2,] 356 333
[3,] 413 333
[4,] 262 333
> Q=sum( (n-np)^2/np  )
> Q
[1] 39.02102

This quantity should be compared with the quantile of the chi-square distribution

> qchisq(.95,df=4-1)
[1] 7.814728

but it is also possible to compute the p-value,

> 1-pchisq(Q,df=4-1)
[1] 1.717959e-08

Here, we reject the assumption that weddings are celebrated uniformly over the year.

# Inference for the Multinomial Distribution

This morning, in our mathematical statistical class, we’ve seen briefly the multinomial distribution, and statistical inference. $\boldsymbol{N}=(N_1,\cdots,N_k)$ has a $\mathcal{M}(n,\boldsymbol{p})$ distribution if its probability function is

$\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})=\frac{n!}{n_1!\cdots n_k!}p_1^{n_1}\cdots p_k^{n_k}$

with $n_1+\cdots+n_k=n$ and $p_1+\cdots+p_k=1$.

The maximum likelihood estimator is then the optimum of

$\widehat{\boldsymbol{p}}=\text{argmax}\{\log\mathcal{L}(\boldsymbol{p}); \boldsymbol{p}^\text{\sffamily T}\boldsymbol{1}=1\}$

We use Lagrange multiplier to solve this constrained optimization problem,

$\widehat{\boldsymbol{p}}=\text{argmax}\{\log\mathcal{L}(\boldsymbol{p})-\lambda[ \boldsymbol{p}^\text{\sffamily T}\boldsymbol{1}-1]\}$

First order conditions are here

$\frac{\partial \log\mathcal{L}(\boldsymbol{p})-\lambda[ \boldsymbol{p}^\text{\sffamily T}\boldsymbol{1}-1]}{\partial p_j}=\frac{n_j}{p_j}-\lambda=0$

and

$\frac{\partial \log\mathcal{L}(\boldsymbol{p})-\lambda[ \boldsymbol{p}^\text{\sffamily T}\boldsymbol{1}-1]}{\partial\lambda}=\boldsymbol{p}^\text{\sffamily T}\boldsymbol{1}-1=0$

Thus,

$\frac{n_j}{\widehat{p}_j}=\widehat{\lambda}$

From

$\sum_{j=1}^k \widehat{p}_j=1$

we can easily get that Lagrande multiplier is $\widehat{\lambda}=n$. And then

$\widehat{p}_j=\frac{n_j}{n}$

One can easily get that this maximum likelihood estimator is unbiased, since $\mathbb{E}[\widehat{\boldsymbol{p}}]=\boldsymbol{p}$. Actually, we can easily prove that

$\mathbb{E}[N_j]=np_j$

and that $\text{Var}[N_j]=np_j(1-p_j)$, while $\text{Cov}[N_j,N_k]=-p_j p_k$. The trick to get the later is simple,

$\text{Var}[N_j+N_k]=\text{Var}[N_j]+\text{Var}[N_k]+2\text{Cov}[N_j,N_k]$

and $\text{Var}[N_j+N_k]=n(p_j+p_k)(1-p_j-p_k)$. Thus, we can easily get the covariance. From that term, we can write that

$\text{Var}[\widehat{\boldsymbol{p}}]=\Sigma$

with

$\Sigma_{j,j}=\frac{p_j(1-p_j)}{n}$

while

$\Sigma_{j,k}=-\frac{p_jp_k}{n}$

# 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)
> regC=glm((cout-k)~bs(densite),data=sousbaseC,family=Gamma(link="log"))

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