Predicted probabilities from a multinomial regression model using zelig and R

I have a multinomial model estimated with the zelig package in R. Whenever I try to use the setx() command, I get an error message saying there is more than one mode. So in stead of using Zelig, I thought I would do it the hard way. I used the instructions here, but I am not sure if I trust these instructions as they give me some weird results.

• Could anyone explain how to get predicted probabilities from a multinomial logit in R?

PS: here are the coefficients from a simple model (there are 6 choice catagories):

 structure(c(3.68021133487111, -0.903496528862169, -1.56339830041814,
-1.13238307296064, -1.67706243532044, -0.177585202845615, 0.0611115470557421,
-0.0458373863009504, 0.0881133593132653, -0.0686190052488972,
0.0163917121907627, 0.0165232098847022, 0.0373815294869855, -0.0353209839724262,
-0.00698911507852077), .Names = c("(Intercept):1", "(Intercept):2",
"(Intercept):3", "(Intercept):4", "(Intercept):5", "Deviance:1",


and here are some values for the independent variables:

structure(c(0.71847390030784, 1.01838748408701, 1.01838748408701,
1.20499277373001, 0.71847390030784, 1.20499277373001, 0.56393315893118,
1.20499277373001, 0.71847390030784, 0.56393315893118, 2, 5, 4,
7, 10, 27, 5, 9, 17, 16), .Dim = c(10L, 2L))

• What "link function" did you use to fit the model? May 29, 2011 at 7:50
• Hi, I used the following command: zelig(choice ~ Deviance + Votes, data = data, model = "mlogit") May 29, 2011 at 7:57

The first thing to do is to construct the "linear predictors" or "logits" for each category for each prediction. So you have your model equation:

$$\eta_{ir}=\sum_{j=1}^{p}X_{ij}\hat{\beta}_{jr}\;\; (i=1,\dots,m\;\; r=1,\dots,R)$$

Where for notational convenience, the above is to be understood to have $\hat{\beta}_{jR}=\eta_{iR}=0$ (as this is the reference category), and $\hat{\beta}_{jr}=0$ if variable $j$ was not included in the linear predictor for class $r$. So you will have an $m\times R$ matrix of logits. You then exponentiate to form predicted odds ratios and renormalise to form predicted probabilities. Note that the predicted odds ratios can be calculated by a simple matrix operation if your data is sufficiently organised: $$\bf{O}=\exp(\boldsymbol{\eta})=\exp(\bf{X}\boldsymbol{\beta})$$ of the $m\times p$ prediction matrix $\bf{X}\equiv\it{\{X_{ij}\}}$ with the $p\times R$ estimated co-efficient matrix $\boldsymbol{\beta}\equiv\it{\{\hat{\beta}_{jr}\}}$, and $\exp(.)$ is defined component wise (i.e. not the matrix exponential). The matrix $\bf{O}$ is an "odds ratio" matrix, with the last column should be all ones. If we take the $m\times 1$ vector $\bf{T}=O1_{R}$ This gives the normalisation constant for each prediction "row" of odds ratios. Now create the $(m\times m)$ diagonal matrix defined by $W_{kk}\equiv T_{k}^{-1}$, and the predicted probability matrix is given by:

$$\bf{P}=\bf{W}\bf{O}$$

So in the example you give the matrix $\boldsymbol{\beta}$ would look like this:

$$\begin{array}{c|c} Int:1 & Int:2 & Int:3 & Int:4 & Int:5 & 0 \\ \hline Dev:1 & Dev:2 & Dev:3 & Dev:4 & Dev:5 & 0 \\ \hline Vote:1 & Vote:2 & Vote:3 & Vote:4 & Vote:5 & 0 \\ \hline \end{array}$$

With the (roughly rounded) values plugged in we get: $$\boldsymbol{\beta}=\begin{array}{c|c} 3.68 & -0.90 & -1.56 & -1.13 & -1.68 & 0 \\ \hline -0.18 & 0.06 & -0.04 & 0.08 & -0.07 & 0 \\ \hline 0.02 & 0.02 & 0.04 & -0.04 & -0.01 & 0 \\ \hline \hline \end{array}$$

And the $\bf{X}$ matrix would look like this:

$$\bf{X}=\begin{array}{c|c} 1 & 0.72 & 2\\ \hline 1 & 1.02 & 5\\ \hline 1 & 1.02 & 4\\ \hline 1 & 1.20 & 7\\ \hline 1 & 0.72 & 10\\ \hline 1 & 1.20 & 27\\ \hline 1 & 0.56 & 5\\ \hline 1 & 1.20 & 9\\ \hline 1 & 0.72 & 17\\ \hline 1 & 0.56 & 16\\ \hline \end{array}$$

So some R-code to do this would simply be (with the matrices $\bf{X}$ and $\boldsymbol{\beta}$ defined as above). The main parts are reading in the data, and padding it with 1s and 0s for $\bf{X}$ and $\boldsymbol{\beta}$:

 beta<-cbind(as.matrix( structure(c(3.68021133487111, -0.903496528862169, -1.56339830041814, -1.13238307296064, -1.67706243532044, -0.177585202845615, 0.0611115470557421, -0.0458373863009504, 0.0881133593132653, -0.0686190052488972, 0.0163917121907627, 0.0165232098847022, 0.0373815294869855, -0.0353209839724262, -0.00698911507852077), .Names = c("(Intercept):1", "(Intercept):2", "(Intercept):3", "(Intercept):4", "(Intercept):5", "Deviance:1", "Deviance:2", "Deviance:3", "Deviance:4", "Deviance:5", "Votes:1", "Votes:2", "Votes:3", "Votes:4", "Votes:5") , .Dim=c(3L,5L)) ),0) X<-cbind(1,as.matrix( structure(c(0.71847390030784, 1.01838748408701, 1.01838748408701, 1.20499277373001, 0.71847390030784, 1.20499277373001, 0.56393315893118, 1.20499277373001, 0.71847390030784, 0.56393315893118, 2, 5, 4, 7, 10, 27, 5, 9, 17, 16), .Dim = c(10L, 2L)) )) P<-diag(as.vector(exp(X %*% beta) %*% as.matrix(rep(1,ncol(beta))))^-1) %*% exp(X %*% beta)