Multinomial Models for Discrete Outcomes

The purpose of this session is to show you how to use R's procedures for doing Multinomial Logit (MNL). Additionally, we look at Ordered Logit and Probit. Note that both STATA and LIMDEP also have “canned” procedures for conditional and nested logit. There is, to date, no canned procedure in R for doing nested logit, but there is a conditional logit procedure in the survival package. LIMDEP also has a “canned” procedure for multinomial Probit. There is an R package (MNP) which will estimate the multinomial Probit model using Markov Chain Monte Carlo methods. However, we do not demonstrate these packages in this lesson.

# This file uses R to compare procedures for multinomial logit and ordered probit and logit for

# data that are naturally ordered. A common mistake is to estimate naturally ordered data with MNL

# or MNP. First, we estimate a multinomial logit (MNL) for the Spector and Mazzeo data with a new

# dependent variable, LETTERS. LETTERS is coded 0=C, 1=B, 2=A. Then we reestimate the model as an

# ordered logit to compare the results.

# Note of caution: Strictly speaking, choice models assume that n individuals make the choices,

# and that the choices are independent. However, in this example it is the instructor making the choices,

# and the choices are obviously not independent. In order to consider this example as a "choice" model,

# we need to think of the instructor as making n independent decisions based on the students' past

# grades, achievement tests, and whether s/he received a treatment, PSI. This may not be substantively correct,

# since the instructor should merely evaluate performance based on course materials, rather

# than paying attention to the student's attributes. However, the objective is to illustrate how to

# estimate these models; # not specify an appropriate substantive model.

library(car)

library(stats)

# First, let's load the data.

Spector <- read.table("c:/users/B. Dan Wood/My Documents/My Teaching/Maximum Likelihood/Data/Letters.txt", header=TRUE)

attach(Spector)

summary(Spector)

# We will first estimate a multinomial logit using the Zelig library (note that GLM could also be used since

# Zelig is partially just an interface for GLM). Note that the "canned" procedures in R use GLM for

# estimating many MLE oriented models.

#Note: The following is to unmask Zelig's simulate function when another library claims that name.

#rm(sim)

library(Zelig)

mlogit.out <- zelig(as.factor(LETTERS) ~ GPA + TUCE + PSI, model = "mlogit", data = Spector)

summary(mlogit.out)

# Following are some interesting values that can be extracted from the model for later use.

# The value of the log-likelihood at the maximum can be used in a likelihood ratio test for

# testin the model specification.

loglik <- -1/2*deviance(mlogit.out)

loglik

# Note that Zelig does not appear to be standardized in the way to extract certain outputs.

# With mlogit the deviance and other interesting values are extracted with e.g.,

# deviance(mlogit.out). However, in other procedures the deviance is obtained by

# mlogit.out$deviance. Note also that this is not consistent with their web site

# help system. Also, names(mlogit.out) returns a null value.

# Same for the other interesting values.

# Model Coefficients

beta <-coefficients(mlogit.out)

beta

# Fitted Probabilities

fitted.probs <- fitted.values(mlogit.out)

fitted.probs

# Let's get the fitted probabilities the hard way now to demonstrate how it is done.

# First create the necessary matrices.

x <- cbind(1,GPA,TUCE,PSI)

betamat<- matrix(coefficients(mlogit.out),4,2, byrow=T)

beta0 <- betamat[,1]

beta1 <- betamat[,2]

# Now calculate the probabilities

p2 <- 1/(1+exp(x %*% beta0)+exp(x %*% beta1))

p1 <- exp(x %*% beta1)*p2

p0 <- exp(x %*% beta0)*p2

# We can also calculate the probabilities with all variables at their means the hard way.

# First, calculate the means.

xbar <- as.matrix(mean(cbind(1,Spector[2:4])))

# Now calculate the probabilities at the means.

p2bar <- 1/(1+exp(t(xbar) %*% beta0)+exp(t(xbar) %*% beta1))

p2bar

p1bar <- exp(t(xbar) %*% beta1)*p2bar

p1bar

p0bar <- exp(t(xbar) %*% beta0)*p2bar

p0bar

# Note that the preceding can be used to calculate first differences simply by replacing

# elements on xbar with interesting units. For example, we might consider moving PSI=0 to PSI=1

# to see the effect on the three probabilities.

xbar[4,1]=0

p2bar <- 1/(1+exp(t(xbar) %*% beta0)+exp(t(xbar) %*% beta1))

p2bar

p1bar <- exp(t(xbar) %*% beta1)*p2bar

p1bar

p0bar <- exp(t(xbar) %*% beta0)*p2bar

p0bar

xbar[4,1]=1

p2bar <- 1/(1+exp(t(xbar) %*% beta0)+exp(t(xbar) %*% beta1))

p2bar

p1bar <- exp(t(xbar) %*% beta1)*p2bar

p1bar

p0bar <- exp(t(xbar) %*% beta0)*p2bar

p0bar

# We can also allow Zelig to calculate first differences.

#Simulate the predicted probabilities of making an A, PSI=0, rest at mean

x.out <- setx(mlogit.out, GPA=2, PSI = 0)

s.out <- sim(mlogit.out, x = x.out, num=10000)

summary(s.out)

#Simulate the predicted probability of making an A, PSI=1, rest at mean

x.out <- setx(mlogit.out, GPA=3, PSI=0)

s.out <- sim(mlogit.out, x = x.out)

summary(s.out)

# We could just take the difference between these to get a first difference. However, Zelig

# can do this in one step.

# First Difference when Shifting GPA from 2 to 3 with PSI=0 and TUCE at mean

gpa.low <- setx(mlogit.out, GPA=2, PSI=0)

gpa.high <- setx(mlogit.out, GPA = 3, PSI=0)

s.out <- sim(mlogit.out, x = gpa.low, x1 = gpa.high)

summary(s.out)

plot(s.out)

# Now let's use Zelig to do an ordered probit and an ordered logit

oprobit.out <- zelig(as.factor(LETTERS) ~ GPA + TUCE + PSI, model = "oprobit", data = Spector)

summary(oprobit.out)

# Simulate ordered probit first differences using PSI=0 to PSI=1

x.low <- setx(oprobit.out, PSI=0)

x.high <- setx(oprobit.out, PSI= 1)

s.out <- sim(oprobit.out, x = x.low, x1 = x.high)

summary(s.out)

plot(s.out)

ologit.out <- zelig(as.factor(LETTERS) ~ GPA + TUCE + PSI, model = "ologit", data = Spector)

summary(ologit.out)

# Simulate ordered logit first differences using PSI=0 to PSI=1

x.low <- setx(ologit.out, PSI=0)

x.high <- setx(ologit.out, PSI= 1)

s.out <- sim(ologit.out, x = x.low, x1 = x.high)

summary(s.out)

plot(s.out)

# Note that there are no "canned" procedures in R for calculating marginal effects or

# testing assumptions about IIA. See Greene 2003 for details that could be programmed.

# Also, there are no "canned" nested logit procedures, or an MNP procedure such as in

# LIMDEP. There is a package for doing MNP with MCMC techniques. However, I have not

# explored its utility for this lesson. Therefore, unless you are willing to do some

# programming work it is probably better to use STATA or LIMDEP when these are needed.