### Lecture 11 - Multicategory Logit Models
### Travel Example

### Data: Greene (1995)
## Travel Choices:
# Measure travelers' decisions regarding travel mode
# between Sydney and Melbourne, Australia.
# Modes of travel (choices):
# mode = 0 --> air
# mode = 1 --> train
# mode = 2 --> bus
# mode = 3 --> car
# (Note that this is a nominal response variable -
# the numbers have no quantitative meaning)

## Primary concern is to quantify the association between explanatory variables of interest and travel choice:
# hinc = household income ($1000)
# psize = number of individuals traveling in a party
# Can you think of other unmeasured confounders?

## Source in Stat 539 R functions (glmCI and summ.mfit):
source("http://www.math.montana.edu/shancock/courses/stat539/r/GillenRFunctions.R")


## Read in data and attach:
travel = read.table("http://www.math.montana.edu/shancock/courses/stat539/data/TravelChoices.txt",header=T)

attach(travel)  
# Allows us to refer to variable names directly (without preceding by travel$)
# Make sure to detach at the end of your R session!

head(travel)
dim(travel)
# Sample size = 210

### Start with descriptive statistics:

## Univariate statistics on response:
## How many and what proportion are in each travel category?
cbind(table(mode), table(mode)/sum(table(mode)))
# --> travel by train most common (30%), then
# car (28%)
# air (27.6%)
# bus (14.2%)

## Bivariate statistics of response vs. explanatory variables:
# Travel mode vs. party size?
table(mode,psize)
# As psize increases, probability of traveling by air decreases.
# Note sparse data for large party sizes (only one observation with psize = 6)

# Travel mode vs. hinc?
mode.fac = factor(mode,levels=c(0:3), labels=c("air","train","bus","car"))
boxplot(hinc~mode.fac,ylab="Household Income",)
lapply( split(hinc,mode), FUN=summary)
# Median household income among travelers choosing air or car appears
# to be higher than those traveling by train or bus.

### For demonstration purposes only (learning experience),
### let's fit separate logistic regressions to each choice
### using air as reference category:
##
#####	Model 1:  Air vs. Train
##
fit1 <- glm( (mode==1) ~ hinc + psize, data=travel, family=binomial, subset=(mode==0 | mode==1) )
glmCI( fit1 )
##
#####	Model 2:  Air vs. Bus
##
fit2 <- glm( (mode==2) ~ hinc + psize, data=travel, family=binomial, subset=(mode==0 | mode==2) )
glmCI( fit2 )
##
#####	Model 3:  Air vs. Car
##
fit3 <- glm( (mode==3) ~ hinc + psize, data=travel, family=binomial, subset=(mode==0 | mode==3) )
glmCI( fit3 )
##

### Interpretations?


### Fit multinomial regression model.
## Need nnet library for "multinom" function:
library(nnet)

mfit <- multinom( mode.fac ~ hinc + psize, data=travel)
summary(mfit) # Not too helpful

# Gillen function sourced at beginning of session:
summ.mfit(mfit)
# rrr = "relative risk ratio" = est. odds of level ___ to level ___
# Interpretations of rrr for psize in each equation?

## Visualize fitted model:
# Plot fitted probabilities vs. hinc for given party size:
co <- coef(mfit) # Note is a 3x3 matrix
is.matrix(co)
dim(co)

# For party size = 2
ps <- 2
# Probability of air:
curve(1/(1+exp(co[1,1]+co[1,2]*x+co[1,3]*ps) +
				exp(co[2,1]+co[2,2]*x+co[2,3]*ps) +
				exp(co[3,1]+co[3,2]*x+co[3,3]*ps)), from=0, to=80,
	xlab="Household Income ($k)", ylab="Fitted Probabilities",
	main="Fitted Probabilities of Travel Type for Parties of Size 2",
	col="red", lty=1, ylim=c(0,1), lwd=2)
# Probability of train:
curve(exp(co[1,1]+co[1,2]*x+co[1,3]*ps)/
			(1+exp(co[1,1]+co[1,2]*x+co[1,3]*ps) +
				exp(co[2,1]+co[2,2]*x+co[2,3]*ps) +
				exp(co[3,1]+co[3,2]*x+co[3,3]*ps)), add=T,
				col="blue",lty=2, lwd=2)
# Probability of bus:
curve(exp(co[2,1]+co[2,2]*x+co[2,3]*ps)/
			(1+exp(co[1,1]+co[1,2]*x+co[1,3]*ps) +
				exp(co[2,1]+co[2,2]*x+co[2,3]*ps) +
				exp(co[3,1]+co[3,2]*x+co[3,3]*ps)), add=T,
				col="darkgreen",lty=3, lwd=2)
# Probability of car:
curve(exp(co[3,1]+co[3,2]*x+co[3,3]*ps)/
			(1+exp(co[1,1]+co[1,2]*x+co[1,3]*ps) +
				exp(co[2,1]+co[2,2]*x+co[2,3]*ps) +
				exp(co[3,1]+co[3,2]*x+co[3,3]*ps)), add=T,
				col="purple",lty=4, lwd=2)
legend(60,.8,c("Air","Train","Bus","Car"),
	col=c("red","blue","darkgreen","purple"),
	lty=c(1,2,3,4))


### What if we want to compare train to car?
V <- vcov(mfit)
row.names(V)

est <- co[1,2]-co[3,2]
se <- sqrt(V[2,2]+V[8,8]-2*V[2,8])


# CI for beta_T1 - beta_C1
est + c(-1,1)*1.96*se

# CI for RRR of train to car for a one unit increase in hinc:
exp(est + c(-1,1)*1.96*se)
# Interpret?


# Any interaction between psize and hinc on mode?
mfit.int = multinom( mode.fac ~ hinc * psize, data=travel)
summary(mfit.int)
summ.mfit(mfit.int)

### Interpretations using interaction model:
# Use Gillen function LinContr.mfit.

## What is the effect of party size on the RRR comparing car to air among those making $20k per year? among those making $50k per year?
LinContr.mfit(contr.names=c("car:psize","car:hinc:psize"), contr.coef=c(1,20), mfit.int)

LinContr.mfit(contr.names=c("car:psize","car:hinc:psize"), contr.coef=c(1,50), mfit.int)


## What is the effect of income on the RRR comparing car to air among those with parties of size 1? among those with parties of size 3?
LinContr.mfit(contr.names=c("car:hinc","car:hinc:psize"), contr.coef=c(1,1), mfit.int)

LinContr.mfit(contr.names=c("car:hinc","car:hinc:psize"), contr.coef=c(1,3), mfit.int)



## LRT comparing H0: mfit to Ha: mfit.int
anova(mfit,mfit.int)
# Conclusion?



detach(travel)