Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/760/data/gator.txt")
attach(Data)

## Descriptive
table(Food)
table(Lake)
table(Size)
table(Sex)

## Set fish as baseline
Food <- relevel(Food, "fish")
Size <- relevel(Size, "small")

## Some model selection
AIC(multinom(Food~1, trace=FALSE))
AIC(multinom(Food~Size, trace=FALSE))
AIC(multinom(Food~Size+Lake, trace=FALSE))
AIC(multinom(Food~Size+Lake+Sex, trace=FALSE))
AIC(multinom(Food~Size*Lake, trace=FALSE))
fit <- multinom(Food~Size+Lake, trace=FALSE) ## Final model

## Plots of probability
New <- data.frame(Size="large", Lake=levels(Lake))
p <- predict(fit, New, type="prob")
col <- c("#FF4E37FF", "#979B00FF", "#00BE5DFF", "#00AFFFFF", "#FF00FFFF")
matplot(p, xaxt="n", pch=19, col=col, type="o", lty=1, ylab="Pr")
axis(1, at=1:4, labels=levels(Lake))
legend("topright", colnames(p), col=col, pch=19, lwd=1)

New <- data.frame(Lake="Trafford", Size=levels(Size))
p <- predict(fit, New, type="prob")
matplot(p, xaxt="n", pch=19, col=col, type="o", lty=1, ylab="Pr")
axis(1, at=1:2, labels=levels(Size))
legend("topright", colnames(p), col=col, pch=19, lwd=1)

## ORs
Food <- relevel(Food, "invert")
fit <- multinom(Food~Size+Lake, trace=FALSE)
exp(cbind(coef(fit)[,2], t(drop(confint(fit, 2)))))

Food <- relevel(Food, "fish")
fit <- multinom(Food~Size+Lake, trace=FALSE)
exp(cbind(coef(fit)[,5], t(drop(confint(fit, 5)))))

## Wald p-values
summ <- summary(fit)
z <- summ$coefficients/summ$standard.errors
2*pnorm(-abs(z))

## LR p-values
fit0 <- multinom(Food~Size, trace=FALSE)
anova(fit0, fit, test="Chisq")

fit0 <- multinom(Food~Lake, trace=FALSE)
anova(fit0, fit, test="Chisq")