trauma <- read.table("http://sites.williams.edu/bklingen/files/2013/12/trauma.dat", header=TRUE)
trauma
library(VGAM)
fit <- vglm(cbind(y1,y2,y3,y4,y5) ~ dose, family=cumulative(parallel=TRUE), data=trauma)
summary(fit)
fitted(fit)

vglm(cbind(y1,y2,y3,y4,y5) ~ 1, family=cumulative(parallel=TRUE), data=trauma)  #  null model
1 - pchisq(2*(53.67903 - 48.87282) , df=1)

fit2 <- vglm(cbind(y1,y2,y3,y4,y5) ~ dose, family=cumulative, data=trauma)  #non-proportional odds
summary(fit2)

### fitted logits
fit.logit <- predict(fit2,se=TRUE)
require("reshape2")
pred.long = melt(as.data.frame(fit.logit$fitted.values))
LB=fit.logit$fitted.values - qnorm(0.975)*fit.logit$se.fit
UB=fit.logit$fitted.values + qnorm(0.975)*fit.logit$se.fit
LB.long=melt(as.data.frame(LB))
UB.long=melt(as.data.frame(UB))
plotdata=data.frame(dose=c(1,2,3,4),type=LB.long$variable, pred=pred.long$value, LB=LB.long$value,UB=UB.long$value)

# plot for fitted logits + Wald confints
require("ggplot2")
ggplot(data=plotdata, aes(x=dose,color=type)) + geom_line(aes(y=pred),size=1.5) + theme_bw() +
theme(legend.position="top", legend.title=element_blank()) + ylab("Cumulative Logits") +
geom_ribbon(aes(ymin=LB, ymax=UB, fill=type),alpha=0.3) + ggtitle("Cumulative Logit Model \n without proportional odds")+
guides(col=guide_legend(nrow=2))

#ggsave(file = "C:\\Users\\Local PC Account\\Research\\presentations\\Ordinal\\FDADec13\\fig8a_logitstrauma_nonpropodds.eps", device=cairo_ps)

#plot for fitted cum. probs + Wald confints
ggplot(data=plotdata, aes(x=dose,color=type)) + geom_line(aes(y=1/(1+exp(-pred))),size=1.5) + theme_bw() +
theme(legend.position="top", legend.title=element_blank()) + ylab("Cumulative Probabilities") +
geom_ribbon(aes(ymin=1/(1+exp(-LB)), ymax=1/(1+exp(-UB)), fill=type),alpha=0.3) + ggtitle("Cumulative Logit Model \n without proportional odds") +
guides(col=guide_legend(nrow=2))
ggsave(file = "C:\\Users\\Local PC Account\\Research\\presentations\\Ordinal\\FDADec13\\fig8b_cumprobtrauma_nonpropodds.eps", device=cairo_ps)

fit.cump <- predict(fit2,untransform=TRUE)
fit.cump.long=cbind(dose=c(1,2,3,4),melt(as.data.frame(fit.cump)))
ggplot(data=fit.cump.long, aes(x=dose,y=value,color=variable)) + geom_line(size=1.5) + theme_bw() +
theme(legend.position="top", legend.title=element_blank()) + ylab("Cumulative Probabilities")



plot(fit.logit[,1]~dose, type="b", col="red", ylim=c(-2,4.5), ylab="Fitted Cumulative Logit", xlab="Dose", main="Cumulative Logit Model \n without proportional odds")
lines(fit.logit[,2]~dose, type="b", col="blue")
## add sample logits:
n <- rowSums(happy[,2:4]) #total sample size
sample.cumprob1 <- happy[,2]/n
sample.cumprob2 <- rowSums(happy[,2:3])/n
sample.logit1 <- logit(sample.cumprob1)
sample.logit2 <- logit(sample.cumprob2)
points(sample.logit1~income, pch="+", col="red")
points(sample.logit2~income, pch="+", col="blue")
legend("top", legend=name, lty=c(1,1), col=c("red","blue"), ncol=2, bty=n)