% \VignetteIndexEntry{Addiction - Multinomial Model}
% \VignetteDepends{nnet}
%\VignetteEngine{knitr::knitr}
%\VignetteEncoding{UTF-8}
\documentclass[a4paper]{article}

\title{Addiction - Multinomial Model}

\begin{document}

\maketitle


<<echo=FALSE,eval=FALSE>>=
options(width=85)
@

First the "addiction" data are loaded and attached.

<<results='hide',eval=FALSE>>=
library(catdata)
data(addiction)
attach(addiction)
@

For the multinomial logit model the function "multinom" from the "nnet"--package is used.

<<eval=FALSE>>=
library(nnet)
@

The response "ill" has to be used as factor.

<<eval=FALSE>>=
ill <- as.factor(ill)
addiction$ill<-as.factor(addiction$ill)
@

The first model is a model with the covariates "gender", "university" and a linear effect of "age"

<<eval=FALSE>>=
multinom0 <- multinom(ill ~ gender + age + university, data=addiction)
summary(multinom0)
@

Another possibility to fit multinomial response models is given by the function "vglm" from the package "VGAM". 
<<eval=FALSE>>=
library(VGAM)
multivgam0<-vglm(ill ~ gender + age + university, multinomial(refLevel=1), 
                 data=addiction)
summary(multivgam0)
@
Both models yield the same parameter estimates. 
\\ \\
The second model includes an additional quadratic effect of "age".

<<eval=FALSE>>=
addiction$age2 <- addiction$age^2
multinom1 <- update(multinom0, . ~ . + age2)
summary(multinom1)

multivgam1<-vglm(ill ~ gender + age + university + age2, multinomial(refLevel=1), 
                 data=addiction)
summary(multivgam1)
@
It should be noted that the standard errors for the models generated by "nnet" and "VGAM" differ when age is included quadratically. The parameter estimates are equal again. 
\\ \\
Now the necessity of the quadratic term is tested by using the function "anova".

<<eval=FALSE>>=
anova(multinom0,multinom1)
multinom1$dev - multinom0$dev
@

Now we  plot the probabilities for the  responses against age. First a sequence within the range of age has to be created.

<<eval=FALSE>>=
minage <- min(na.omit(age))
maxage <- max(na.omit(age))

ageindex <- seq(minage, maxage, 0.1)
n <- length(ageindex)
@

Now the   vectors for the other covariates  and the  data sets for men and women 
are  built. 

<<eval=FALSE>>=
ageindex2 <- ageindex^2

gender1 <- rep(1, n)
gender0 <- rep(0, n)
university1 <- rep(1, n)

datamale <- as.data.frame(cbind(gender=gender0,age=ageindex,university=
  university1,age2=ageindex2)) 
datafemale <- as.data.frame(cbind(gender=gender1,age=ageindex,university=
  university1,age2=ageindex2))
@

Now for the built data sets the probabilities based on model "multinom1" are computed.

<<eval=FALSE>>=
probsmale <- predict(multinom1, datamale, type="probs")
probsfemale <- predict(multinom1, datafemale, type="probs")
@

Now the probabilities can be plotted.

<<echo=TRUE,eval=FALSE>>=
par(cex=1.4, lwd=2)

plot(ageindex, probsmale[,1], type="l", lty=1, ylim=c(0,1), main=
"men with university degree", ylab="probabilities")
lines(ageindex, probsmale[,2], lty="dotted") 
lines(ageindex, probsmale[,3], lty="dashed") 
legend("topright", legend=c("Weak-willed", "diseased", "both"), lty=c("solid", 
"dotted", "dashed"))
@

<<echo=TRUE,eval=FALSE>>=
par(cex=1.4, lwd=2)

plot(ageindex, probsfemale[,1], type="l", lty=1, ylim=c(0,1), main=
  "women with university degree", ylab="probabilities")
lines(ageindex, probsfemale[,2], lty="dotted")                                                  
lines(ageindex, probsfemale[,3], lty="dashed") 
legend("topright", legend=c("Weak-willed", "diseased", "both"), 
       lty=c("solid", "dotted", "dashed"))
@
\end{document}