pneumo.df <- read.table("data/pneumo.dat",header=T)
attach(pneumo.df)
lodds <- function(x,y) { log((x+.5)/(y+.5)) }
# Proportional odds model fit
#postscript("Rcode/pneumo-po.eps",horizontal=FALSE,onefile=FALSE,height=8,width=6)
par(mfrow=c(2,1))
plot(log(time),lodds(normal,mild+severe),xlab="",ylab="1st CUM.LOGIT")
title("Cumulative logit plots for pneumoconiosis data")
plot(log(time),lodds(mild+normal,severe),xlab="LOG EXPOSURE",ylab="2nd CUM.LOGIT")
#dev.off()
logexpos <- log(rep(time,3))
freq <- c(normal,mild,severe)
response <- as.factor(rep(1:3,each=8))
pnewmo.df <- data.frame(response,logexpos,freq)
library(MASS)
pnewmo.po <- polr(response~logexpos,weights=freq,data=pnewmo.df)
# Continuation ratio model fit
#postscript("Rcode/pneumo-cr.eps",horizontal=FALSE,onefile=FALSE,height=8,width=6)
par(mfrow=c(2,1))
plot(log(time),lodds(mild+severe,normal),xlab="",ylab="1st CR.LOGIT")
title("Continuation-ratio logit plots for pneumoconiosis data")
plot(log(time),lodds(severe,mild),xlab="LOG EXPOSURE",ylab="2nd CR.LOGIT")
#dev.off()
m <- c(normal+mild+severe,mild+severe)
y <- c(mild+severe,severe)
x <- as.factor(rep(1:2,each=8))
lt <- log(rep(time,2))
pneumo.cr <- glm(cbind(y,m-y)~x+lt-1,binomial)