coal <- read.table("http://sites.williams.edu/bklingen/files/2012/02/coalminers.txt",header=TRUE)
coal
require("VGAM")
# Proportional odds model:
fit1 <- vglm(cbind(normal, mild, severe) ~ log(exposure), cumulative(parallel=TRUE), data=coal)
#testing proportional odds assumption via LR test:
fit2 <- vglm(cbind(normal, mild, severe) ~ log(exposure), cumulative(parallel=FALSE), data=coal)
deviance(fit1)-deviance(fit2) #very small, so proportional odds assumption reasonable
coef(fit1, matrix=TRUE)

exposure1=5:100
pred <- predict(fit1,type="link",newdata=data.frame(exposure=exposure1))
plot(pred[,1]~log(exposure1),type="l")
lines(pred[,2]~log(exposure1))
#Compare with sample logits to check fit:
n=rowSums(coal[,2:4])
cp1=coal[,2]/n
cp2=rowSums(coal[,2:3])/n
logitcp1=log(cp1/(1-cp1))
logitcp2=log(cp2/(1-cp2))
points(logitcp1~log(coal$exposure))
points(logitcp2~log(coal$exposure))

#fitted probabilities:
prob <- predict(fit1,type="response",newdata=data.frame(exposure=exposure1))
plot(prob[,1]~exposure1,type="l",ylim=c(0,1),col="green",xlab="exposure (years)", ylab="Prob. of reaction")
lines(prob[,2]~exposure1,col="orange")
lines(prob[,3]~exposure1,col="red")

cumprob <- t(apply(prob,1,cumsum))
plot(cumprob[,1]~exposure1,type="l",ylim=c(0,1),col="green",xlab="exposure (years)", ylab="Cum. prob. of reaction")
lines(cumprob[,2]~exposure1,col="orange")