mental <- read.table(file="http://sites.williams.edu/bklingen/files/2013/12/mental.dat", header=TRUE)
mental$impair = factor(mental$impair, labels=c("well","mild","moderate","impaired"), ordered=TRUE)
head(mental)
## Ordinal Package
require(ordinal)
fit <- clm(impair ~ life + ses, data=mental)
summary(fit)
fit0 <- clm(impair ~ ses, data=mental)
anova(fit,fit0)

fit1 <- clm(impair ~ ses*life, data=mental)
summary(fit1)
anova(fit1,fit)


## VGAM Package
require(VGAM)
fit <- vglm(impair ~ life + ses, family=cumulative(parallel=TRUE), data=mental)
summary(fit)
#coefficients(fit,matrix=TRUE)
maxl=logLik(fit)
maxl
fit0 <- vglm(impair ~ ses, family=cumulative(parallel=TRUE), data=mental)
maxl0 = logLik(fit0)
maxl0
LR.stat = -2*(maxl0 - maxl)
LR.stat
1 - pchisq(LR.stat,df=1)

fit1 <- vglm(impair ~ life + ses + life*ses, family=cumulative(parallel=TRUE), data=mental)
LR.stat = -2*(maxl - logLik(fit1))
LR.stat
1 - pchisq(LR.stat,df=1)

### fitted logits
life1 <- seq(0,9,1)
fit.logit.ses0 <- predict(fit,newdata=data.frame(life=life1,ses=0))
fit.logit.ses1 <- predict(fit,newdata=data.frame(life=life1,ses=1))
name <- colnames(fit.logit.ses1)
plot.data <- data.frame(life=rep(life1,3),ses=rep(c("SES = low","SES = high"),each=3*10), type=rep(name,each=10), logit=c(fit.logit.ses0,fit.logit.ses1))
head(plot.data)
require("lattice")
xyplot(logit~life|ses, group=type, data=plot.data, type="l", auto.key=list(points=FALSE, lines=TRUE, columns=3))

### fitted cumulative probabilities
fit.cumprob.ses0 <- predict(fit,newdata=data.frame(life=life1,ses=0), untransform=TRUE)
fit.cumprob.ses1 <- predict(fit,newdata=data.frame(life=life1,ses=1), untransform=TRUE)
name <- colnames(fit.cumprob.ses1)
plot.data <- data.frame(life=rep(life1,3),ses=rep(c("SES = low","SES = high"),each=3*10), type=rep(name,each=10), cumProb=c(fit.cumprob.ses0,fit.cumprob.ses1))
head(plot.data)
xyplot(cumProb~life|ses, group=type, data=plot.data, type="l", auto.key=list(points=FALSE, lines=TRUE, columns=3))

### fitted category probabilities
fit.prob.ses0 <- predict(fit,newdata=data.frame(life=life1,ses=0), type="response")
fit.prob.ses1 <- predict(fit,newdata=data.frame(life=life1,ses=1), type="response")
name <- colnames(fit.prob.ses1)
plot.data <- data.frame(life=rep(life1,4),ses=rep(c("SES = low","SES = high"),each=4*10), type=rep(name,each=10), Prob=c(fit.prob.ses0,fit.prob.ses1))
head(plot.data)
xyplot(Prob~life|ses, group=type, data=plot.data, type="l", auto.key = list(points=FALSE, lines=TRUE, columns=3))