{"word_count": 368, "character_count": 4485}
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 [latex]Y[/latex] with a multinomial distribution (3 factors to have something more general than the binomial), taking values [latex]\{A,B,C\}[/latex], with respective probabilities [latex]\mathbf{p}=(p_A,p_B,p_C)[/latex]. Here is a code to generate some multinomial variablesmsample=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 [latex]n[/latex] 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+X2s2=pb[2]-X1+X2P1=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[latex display="true"]\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}]}[/latex][latex display="true"]\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}]}[/latex][latex display="true"]\mathbb{P}[Y=B|\mathbf{x}]=\frac{1}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}[/latex]For convenience, consider the most popular factor in our training datasetmodalite=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 [latex]Y_A=\mathbf{1}(Y=A)[/latex]. Actually, we will consider the most important factor, but for convenience, assume that it is [latex]A[/latex].[latex display="true"]\mathbb{P}[Y_A=A|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{a}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{a}]}[/latex]On our dataset, estimate that model, and get predictions. In the case where [latex]Y\neq A[/latex], define another Bernoulli variable [latex]Y_B=\mathbf{1}(Y=B|Y\neq A)[/latex]. We can estimate that model and derive two probabilities, [latex]\mathbb{P}(Y=B|Y\neq A)[/latex] and [latex]\mathbb{P}(Y=C|Y\neq A)[/latex] (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. [latex]\mathbb{P}[Y=A][/latex] is obtained from the first model, and we can derive the other two from [latex]\mathbb{P}[Y=B|Y\neq A]\cdot\mathbb{P}[Y\neq A][/latex] and [latex]\mathbb{P}[Y=C|Y\neq A]\cdot\mathbb{P}[Y\neq A][/latex].  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)=modaliteLet us compare the predicted probabilites, on the same dataset (here the training dataset)&#62; mmp1[1:9,c("0","1","2")]0         1         21  0.19728737 0.4991805 0.30353212  0.17244580 0.5648537 0.26270053  0.19291753 0.5971058 0.20997674  0.09087176 0.7787304 0.13039785  0.23400225 0.4083022 0.35769556  0.18063647 0.6637352 0.15562837  0.13188881 0.7402710 0.12784018  0.13776970 0.6524959 0.20973449  0.12325864 0.6790336 0.1977078&#62; mp1[1:9,c("0","1","2")]0         1         21  0.19691036 0.5022692 0.30082052  0.17123189 0.5680647 0.26070343  0.19293066 0.5984402 0.20862914  0.08821851 0.7813318 0.13044975  0.23470739 0.4109990 0.35429366  0.18249687 0.6602168 0.15728637  0.13128711 0.7400898 0.12862318  0.13525341 0.6553618 0.20938489  0.12090016 0.6815915 0.1975084The two are very close. So yes, it is possible to see the multinomial regression as some sequential binomial regressions.