6
$\begingroup$

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", 
"Deviance:2", "Deviance:3", "Deviance:4", "Deviance:5", "Votes:1", 
"Votes:2", "Votes:3", "Votes:4", "Votes:5"))

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))
$\endgroup$
2
  • $\begingroup$ What "link function" did you use to fit the model? $\endgroup$ May 29, 2011 at 7:50
  • $\begingroup$ Hi, I used the following command: zelig(choice ~ Deviance + Votes, data = data, model = "mlogit") $\endgroup$ May 29, 2011 at 7:57

1 Answer 1

3
$\begingroup$

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)

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.