MultiNomReg.R

# script originally found here:
# http://www.unt.edu/rss/class/Jon/R_SC/Module9/MultiNomReg.R



#
#
########## Multinomial Logistic Regression example ##########
#
#
# This script assumes you have worked through all the previous notes from 
# the web page and you have downloaded, installed, and updated all available
# R packages. 

# Load the following libraries if you have not already.

library(foreign)
library(Rcmdr)

# The MultiNomReg.sav data file contains one categorical outcome (y) with 3 categories, 
# and three continuous predictors (x1, x2, x3).  
# Load the "MultiNomReg.sav" data file naming is 'mdata1'.

mdata1 <- read.spss("http://www.unt.edu/rss/class/Jon/R_SC/Module9/MultiNomReg.sav", 
  use.value.labels=TRUE, max.value.labels=Inf, to.data.frame=TRUE)
summary(mdata1)

# Next, change the outcome variable (which is currently numeric) to a factor.

mdata2 <- mdata1
mdata2$y <- as.factor(mdata2$y)
summary(mdata2)

## Conduct the multinomial analysis using the mlogit library; which allows you to specify the reference 
## category of the outcome. 

# Load the appropriate library.

library(mlogit)

# Then create a new 'form' of data which can be read by the mlogit function. 
# To do this, we need to expand the outcome variable (y) much like we would 
# for dummy coding a categorical variable for inclusion in standard multiple 
# regression.

mdata3 <- mlogit.data(mdata2, varying=NULL, choice="y", shape="wide")
head(mdata3)

# Run the multinomial analysis with 'mlogit' function, specifying the reference category as "1".

model.1 <- mlogit(y ~ 1 | x1 + x2 + x3, data=mdata3, reflevel="1")
summary(model.1)

# To get the "exp(b)" values (ignore the intercepts which are the first 3 printed and are 
# labeled "alt[T.2]", "alt[T.3]", "alt[T.4]"):

exp(coef(model.1))

# If desired, residuals can be produced.

residuals(model.1)

hist(residuals(model.1))

# Rcmdr can lead to different results because it uses the "nnet" package. 
# If you specify the same model using Rcmdr, you get the following script; which 
# produces slightly different results when compared to what is above. 
 
MLM.1 <- multinom(y ~ x1 +x2 +x3, data=mdata2, trace=FALSE)
summary(MLM.1, cor=FALSE, Wald=TRUE)