###################################################
### chunk number 1: 
###################################################
library(lattice)
options(width=65,show.signif.stars=FALSE)
options(SweaveHooks=list(fig=function() {
ltheme <- canonical.theme(color = FALSE)     ## in-built B&W theme
ltheme$strip.background$col <- "transparent" ## change strip bg
ltheme$plot.symbol$col <- "black"
ltheme$axis.text$cex <- 1.2
ltheme$par.xlab.text$cex <- 1.2
ltheme$par.ylab.text$cex <- 1.2
lattice.options(default.theme = ltheme)      ## set as default
}))


###################################################
### chunk number 2: mice.plot eval=FALSE
###################################################
## mice <- data.frame(dose=c(422,744,948,2069),
## r=c(0,1,3,5),n=c(5,5,5,5))
## mice <- transform(mice,
##                   ldose=log(dose),
##                   p = r/n)
## library(lattice)
## print(xyplot(p~ldose,data=mice,xlab="log dose"))


###################################################
### chunk number 3: 
###################################################
mice <- data.frame(dose=c(422,744,948,2069),
r=c(0,1,3,5),n=c(5,5,5,5))
mice <- transform(mice,
                  ldose=log(dose),
                  p = r/n)
library(lattice)
print(xyplot(p~ldose,data=mice,xlab="log dose"))


###################################################
### chunk number 4: mice.identity
###################################################
try(glm(cbind(r,(n-r))~ldose,
    data=mice,family=binomial(link="identity")))


###################################################
### chunk number 5: 
###################################################
(mice.glm1 <- glm(cbind(r,(n-r))~ldose,data=mice,family=binomial))
#summary(mice.glm1)
gridx <- seq(6,7.75,by=0.01)
fpg <- predict(mice.glm1,new=data.frame(ldose=gridx),type='response')
(mice.glm2 <- update(mice.glm1,family=binomial(link="probit")))
#summary(mice.glm2)
fpp <- predict(mice.glm2,new=data.frame(ldose=gridx),type='response')
(mice.glm3 <- update(mice.glm1,family=binomial(link="cloglog")))
#summary(mice.glm3)
fpc <- predict(mice.glm3,new=data.frame(ldose=gridx),type='response')
#
(mice.glm4 <- update(mice.glm1,family=binomial(link="cauchit")))
fpcauch <- predict(mice.glm4,new=data.frame(ldose=gridx),type='response')
#
mice.new.df <- data.frame(gridx,fpg,fpp,fpc,fpcauch)


###################################################
### chunk number 6: mice.fitted.plots eval=FALSE
###################################################
## library(lattice)
## print(xyplot(fpg+fpp+fpc+fpcauch~gridx,data=mice.new.df,type='l',
## xlab='log dose',ylab='p',lwd=1.5,lty=c(1,2,3,4),
## #key=list(text=c("probit","cloglog","response","cauchit")),
## panel=function(x,y,...){
## panel.xyplot(x,y,...)
## panel.text(mice$ldose,mice$p,"x",cex=1.5,type="p",...)}))
## #legend(x=6.1,y=0.9,legend=c("logit","probit","cloglog"),lty=1:3)}))


###################################################
### chunk number 7: 
###################################################
library(lattice)
print(xyplot(fpg+fpp+fpc+fpcauch~gridx,data=mice.new.df,type='l',
xlab='log dose',ylab='p',lwd=1.5,lty=c(1,2,3,4),
#key=list(text=c("probit","cloglog","response","cauchit")),
panel=function(x,y,...){
panel.xyplot(x,y,...)
panel.text(mice$ldose,mice$p,"x",cex=1.5,type="p",...)}))
#legend(x=6.1,y=0.9,legend=c("logit","probit","cloglog"),lty=1:3)}))


###################################################
### chunk number 8: profile.slope eval=FALSE
###################################################
## beta1 <- seq(0,25,by=0.25)
## pl <- NULL
## #beta0 <- NULL
## for (b1 in beta1){
## z <- b1*mice$ldose
## fit <- glm(cbind(r,(n-r))~offset(z),data=mice,family=binomial)
## logl  <- with(mice,sum(dbinom(r,n,fitted(fit),log=TRUE)))
## #beta0 <- c(beta0,coef(fit))
## pl <- c(pl,logl)
## }
## pl <- pl - max(pl)
## print(xyplot(exp(pl)~beta1,type='l',
## xlab=expression(beta[1]),
## ylab="relative likelihood",lwd=1.5,
## panel=function(...){
## panel.xyplot(...)
## panel.abline(h=exp(3.84/-2),lty=2,lwd=1.5)}))


###################################################
### chunk number 9: 
###################################################
beta1 <- seq(0,25,by=0.25)
pl <- NULL
#beta0 <- NULL
for (b1 in beta1){
z <- b1*mice$ldose
fit <- glm(cbind(r,(n-r))~offset(z),data=mice,family=binomial)
logl  <- with(mice,sum(dbinom(r,n,fitted(fit),log=TRUE)))
#beta0 <- c(beta0,coef(fit))
pl <- c(pl,logl)
}
pl <- pl - max(pl)
print(xyplot(exp(pl)~beta1,type='l',
xlab=expression(beta[1]),
ylab="relative likelihood",lwd=1.5,
panel=function(...){
panel.xyplot(...)
panel.abline(h=exp(3.84/-2),lty=2,lwd=1.5)}))


###################################################
### chunk number 10: 
###################################################
fc2(beta1[which.max(pl)])
#a <- beta1[-2*pl < 3.84]
#print(a[1])
#print(rev(a)[1])


###################################################
### chunk number 11: 
###################################################
confint(mice.glm1)


###################################################
### chunk number 12: profile.prob
###################################################
pseq <- seq(0.1,0.99,by=0.01)
logl <- NULL
xdose <- NULL
x0 <- log(1100)
for (p in pseq){
theta <- rep(log(p/(1-p)),4)
fit <- glm(cbind(r,(n-r))~offset(theta)+I(ldose-x0)-1,data=mice,family=binomial)
logl <- c(logl,with(mice,sum(dbinom(r,n,fitted(fit),log=TRUE))))
}
logl <- logl-max(logl)
#a <- pseq[logl>-1.92]
print(xyplot(exp(logl)~pseq,type="l",
xlab='probability',ylab='relative likelihood',
panel=function(...){
panel.xyplot(...)
panel.abline(h=exp(3.84/-2),lty=2,lwd=1.5)}))


###################################################
### chunk number 13: 
###################################################
a <- pseq[exp(logl)>exp(3.84/-2)]


###################################################
### chunk number 14: 
###################################################
pseq <- seq(0.1,0.99,by=0.01)
logl <- NULL
xdose <- NULL
x0 <- log(1100)
for (p in pseq){
theta <- rep(log(p/(1-p)),4)
fit <- glm(cbind(r,(n-r))~offset(theta)+I(ldose-x0)-1,data=mice,family=binomial)
logl <- c(logl,with(mice,sum(dbinom(r,n,fitted(fit),log=TRUE))))
}
logl <- logl-max(logl)
#a <- pseq[logl>-1.92]
print(xyplot(exp(logl)~pseq,type="l",
xlab='probability',ylab='relative likelihood',
panel=function(...){
panel.xyplot(...)
panel.abline(h=exp(3.84/-2),lty=2,lwd=1.5)}))


###################################################
### chunk number 15: profile.ED50
###################################################
ldseq50 <- log(seq(600,1600,len=1001))
logl <- NULL
fit <- NULL
p <- 0.5
logitp <- rep(log(p/(1-p)),4)
for (ld in ldseq50){
fit <- glm(cbind(r,(n-r))~I(ldose-ld)+offset(logitp)-1,
data=mice,family=binomial)
logl <- c(logl,sum(with(mice,dbinom(r,n,fitted(fit),log=TRUE))))
}
logl <- logl - max(logl)
plot(ldseq50,exp(logl),type='l',xlab="log dose",ylab="relative likelihood")
abline(h=exp(3.84/-2),lty=2,lwd=1.5)


###################################################
### chunk number 16: 
###################################################
a <- ldseq50[exp(logl)>0.15]


###################################################
### chunk number 17: 
###################################################
ldseq50 <- log(seq(600,1600,len=1001))
logl <- NULL
fit <- NULL
p <- 0.5
logitp <- rep(log(p/(1-p)),4)
for (ld in ldseq50){
fit <- glm(cbind(r,(n-r))~I(ldose-ld)+offset(logitp)-1,
data=mice,family=binomial)
logl <- c(logl,sum(with(mice,dbinom(r,n,fitted(fit),log=TRUE))))
}
logl <- logl - max(logl)
plot(ldseq50,exp(logl),type='l',xlab="log dose",ylab="relative likelihood")
abline(h=exp(3.84/-2),lty=2,lwd=1.5)


###################################################
### chunk number 18: profile_ED90
###################################################
ldseq90 <- log(seq(800,8000,len=1001))
logl <- NULL
fit <- NULL
p <- 0.9
logitp <- rep(log(p/(1-p)),4)
for (ld in ldseq90){
fit <- glm(cbind(r,(n-r))~I(ldose-ld)+offset(logitp)-1,
data=mice,family=binomial)
logl <- c(logl,sum(with(mice,dbinom(r,n,fitted(fit),log=TRUE))))
}
logl <- logl - max(logl)
plot(ldseq90,exp(logl),type='l',xlab="log dose",ylab="relative likelihood")
abline(h=exp(3.84/-2),lty=2,lwd=1.5)
aa <- ldseq90[exp(logl)>0.15]


###################################################
### chunk number 19: 
###################################################
vaso <- read.csv("Data/vaso.csv")
vaso$Subject <- factor(vaso$Subject)


###################################################
### chunk number 20: vaso eval=FALSE
###################################################
## data(vaso,package="SMIR")


###################################################
### chunk number 21: vaso.model
###################################################
(vaso.glm <- glm(Y~Rate+Volume,data=vaso,family=binomial))


###################################################
### chunk number 22: 
###################################################
round(resid(vaso.glm,type="pearson"),4)
round(sum(resid(vaso.glm,type="pearson")^2),2)


###################################################
### chunk number 23: vaso.response eval=FALSE
###################################################
## print(xyplot(Rate~Volume,data=vaso,group=Y))


###################################################
### chunk number 24: 
###################################################
print(xyplot(Rate~Volume,data=vaso,group=Y))


###################################################
### chunk number 25: 
###################################################
vaso.glm1 <- update(vaso.glm,.~.+Subject)
coef(vaso.glm1)
anova(vaso.glm,vaso.glm1)
round(resid(vaso.glm1,type="pearson"),4)
round(sum(resid(vaso.glm1,type="pearson")^2),2)


###################################################
### chunk number 26: vaso.linear.model.curves eval=FALSE
###################################################
## print(xyplot(Rate~Volume,data=vaso,group=Y,pch=c(1,3),col="black",
##        panel=function(x,y,...){
##          panel.xyplot(x,y,...)
##          probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
##          logits <- log(probs/(1-probs))
##          line.types <- c(4,3,2,1,2,3,4)
##          for (i in 1:7) 
##            panel.curve((logits[i] - coef(vaso.glm)[1] - coef(vaso.glm)[2] * x)/coef(vaso.glm)[3],
##                        lty=line.types[i])
##          panel.text(c(2,2.6,3.2),c(0.6,0.7,0.8),c("10","50","90"),font=2)
## }))


###################################################
### chunk number 27: 
###################################################
print(xyplot(Rate~Volume,data=vaso,group=Y,pch=c(1,3),col="black",
       panel=function(x,y,...){
         panel.xyplot(x,y,...)
         probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
         logits <- log(probs/(1-probs))
         line.types <- c(4,3,2,1,2,3,4)
         for (i in 1:7) 
           panel.curve((logits[i] - coef(vaso.glm)[1] - coef(vaso.glm)[2] * x)/coef(vaso.glm)[3],
                       lty=line.types[i])
         panel.text(c(2,2.6,3.2),c(0.6,0.7,0.8),c("10","50","90"),font=2)
}))


###################################################
### chunk number 28: vaso.glm2
###################################################
vaso <- transform(vaso,lv=log(Volume),lr=log(Rate))
vaso.glm2 <- glm(Y~lv+lr,data=vaso,family=binomial)
anova(vaso.glm,vaso.glm2)
#summary(vaso.glm2)
sort(round(resid(vaso.glm2,type="pearson"),4))
round(sum(resid(vaso.glm2,type="pearson")^2),4)


###################################################
### chunk number 29: vaso.log.model.curves eval=FALSE
###################################################
## logvolest <- function(x,model,p){
## logitp <- log(p/(1-p))
## exp((logitp-coef(model)[1]-coef(model)[2]*log(x))/coef(model)[3])
## }
## print(xyplot(Rate~Volume,data=vaso,group=Y,pch=c(1,3),col="black",
##        panel=function(x,y,...){
##          panel.xyplot(x,y,...)
##          probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
##          line.types <- c(4,3,2,1,2,3,4)
##          for (i in 1:7)
##          panel.curve(logvolest(x,vaso.glm2,p=probs[i]),lty=line.types[i])
## #         panel.curve(logvolest(x,vaso.glm2,p=0.1),lty=3)
## #         panel.curve(logvolest(x,vaso.glm2,p=0.3),lty=2)
## #         panel.curve(logvolest(x,vaso.glm2,p=0.5),lty=1)
## #         panel.curve(logvolest(x,vaso.glm2,p=0.7),lty=2)
## #         panel.curve(logvolest(x,vaso.glm2,p=0.9),lty=3)
## #         panel.curve(logvolest(x,vaso.glm2,p=0.95),lty=4)
##          panel.text(c(2,2.6,3.2),c(0.55,0.65,0.8),c("10","50","90"),font=2)
##        }))


###################################################
### chunk number 30: 
###################################################
logvolest <- function(x,model,p){
logitp <- log(p/(1-p))
exp((logitp-coef(model)[1]-coef(model)[2]*log(x))/coef(model)[3])
}
print(xyplot(Rate~Volume,data=vaso,group=Y,pch=c(1,3),col="black",
       panel=function(x,y,...){
         panel.xyplot(x,y,...)
         probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
         line.types <- c(4,3,2,1,2,3,4)
         for (i in 1:7)
         panel.curve(logvolest(x,vaso.glm2,p=probs[i]),lty=line.types[i])
#         panel.curve(logvolest(x,vaso.glm2,p=0.1),lty=3)
#         panel.curve(logvolest(x,vaso.glm2,p=0.3),lty=2)
#         panel.curve(logvolest(x,vaso.glm2,p=0.5),lty=1)
#         panel.curve(logvolest(x,vaso.glm2,p=0.7),lty=2)
#         panel.curve(logvolest(x,vaso.glm2,p=0.9),lty=3)
#         panel.curve(logvolest(x,vaso.glm2,p=0.95),lty=4)
         panel.text(c(2,2.6,3.2),c(0.55,0.65,0.8),c("10","50","90"),font=2)
       }))


###################################################
### chunk number 31: vaso.hatvalues eval=FALSE
###################################################
## round(hatvalues(vaso.glm),4)
## print(dotplot(~hatvalues(vaso.glm),
## panel=function(...){
## panel.dotplot(...)
## panel.abline(v=0.154,lty='dotted')}))


###################################################
### chunk number 32: 
###################################################
round(hatvalues(vaso.glm),4)
print(dotplot(~hatvalues(vaso.glm),
panel=function(...){
panel.dotplot(...)
panel.abline(v=0.154,lty='dotted')}))


###################################################
### chunk number 33: vaso.glm3
###################################################
vaso.glm3 <- update(vaso.glm,subset=row.names(vaso)!="32")
summary(vaso.glm3)


###################################################
### chunk number 34: vaso.glm4
###################################################
vaso.glm4 <- update(vaso.glm,subset=!row.names(vaso)%in%c("4","18"))
summary(vaso.glm4)


###################################################
### chunk number 35: vaso.2.model.curves eval=FALSE
###################################################
## volest <- function(x,model,p){
## logitp <- log(p/(1-p))
## (logitp-coef(model)[1]-coef(model)[3]*x)/coef(model)[2]
## }
## print(xyplot(Rate~Volume,data=vaso,group=Y,pch=c(1,3),col="black",
##        panel=function(x,y,...){
##          panel.xyplot(x,y,...)
##          probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
##          line.types <- c(4,3,2,1,2,3,4)
##          for (i in 1:7)
##          panel.curve(volest(x,vaso.glm4,p=probs[i]),lty=line.types[i])
## #         panel.curve(volest(x,vaso.glm4,p=0.05),lty=4)
## #         panel.curve(volest(x,vaso.glm4,p=0.1),lty=3)
## #         panel.curve(volest(x,vaso.glm4,p=0.3),lty=2)
## #         panel.curve(volest(x,vaso.glm4,p=0.5),lty=1)
## #         panel.curve(volest(x,vaso.glm4,p=0.7),lty=2)
## #         panel.curve(volest(x,vaso.glm4,p=0.9),lty=3)
## #         panel.curve(volest(x,vaso.glm4,p=0.95),lty=4)
##          panel.text(c(2,2,2),c(0.4,0.7,1),c("10","50","90"),font=2)
##        }))


###################################################
### chunk number 36: 
###################################################
volest <- function(x,model,p){
logitp <- log(p/(1-p))
(logitp-coef(model)[1]-coef(model)[3]*x)/coef(model)[2]
}
print(xyplot(Rate~Volume,data=vaso,group=Y,pch=c(1,3),col="black",
       panel=function(x,y,...){
         panel.xyplot(x,y,...)
         probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
         line.types <- c(4,3,2,1,2,3,4)
         for (i in 1:7)
         panel.curve(volest(x,vaso.glm4,p=probs[i]),lty=line.types[i])
#         panel.curve(volest(x,vaso.glm4,p=0.05),lty=4)
#         panel.curve(volest(x,vaso.glm4,p=0.1),lty=3)
#         panel.curve(volest(x,vaso.glm4,p=0.3),lty=2)
#         panel.curve(volest(x,vaso.glm4,p=0.5),lty=1)
#         panel.curve(volest(x,vaso.glm4,p=0.7),lty=2)
#         panel.curve(volest(x,vaso.glm4,p=0.9),lty=3)
#         panel.curve(volest(x,vaso.glm4,p=0.95),lty=4)
         panel.text(c(2,2,2),c(0.4,0.7,1),c("10","50","90"),font=2)
       }))


###################################################
### chunk number 37: vaso.glm5
###################################################
vaso.glm5 <- glm(Y~lv+lr,data=vaso,family=binomial)
summary(vaso.glm5)
sort(round(hatvalues(vaso.glm5),4))


###################################################
### chunk number 38: vaso.glm6
###################################################
vaso.glm6 <- update(vaso.glm5,subset=row.names(vaso)!="31")
summary(vaso.glm6)


###################################################
### chunk number 39: vaso.glm7
###################################################
vaso.glm7 <- update(vaso.glm5,subset=!row.names(vaso)%in%c("4","18"))
summary(vaso.glm7)


###################################################
### chunk number 40: vaso.4.model.curves eval=FALSE
###################################################
## print(xyplot(Rate~Volume,data=vaso,group=Y,pch=c(1,3),col="black",
##              subset=!row.names(vaso)%in%c("4","18"),
##        panel=function(x,y,...){
##          panel.xyplot(x,y,...)
##          probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
##          logits <- log(probs/(1-probs))
##          line.types <- c(4,3,2,1,2,3,4)
##          for (i in 1:7)
##          panel.curve(exp((logits[i] + 24.58 - 39.55 * log(x))/31.93),lty=line.types[i])
## #         panel.curve(exp((-2.94 + 24.58 - 39.55 * log(x))/31.93),lty=4)
## #         panel.curve(exp((-2.02 + 24.58 - 39.55 * log(x))/31.93),lty=3)
## #         panel.curve(exp((-0.85 + 24.58 - 39.55 * log(x))/31.93),lty=2)
## #         panel.curve(exp((0 + 24.58 - 39.55 * log(x))/31.93),lty=1)
## #         panel.curve(exp((0.85 + 24.58 - 39.55 * log(x))/31.93),lty=2)
## #         panel.curve(exp((2.02 + 24.58 - 39.55 * log(x))/31.93),lty=3)
## #         panel.curve(exp((2.94 + 24.58 - 39.55 * log(x))/31.93),lty=4)
## }))


###################################################
### chunk number 41: 
###################################################
print(xyplot(Rate~Volume,data=vaso,group=Y,pch=c(1,3),col="black",
             subset=!row.names(vaso)%in%c("4","18"),
       panel=function(x,y,...){
         panel.xyplot(x,y,...)
         probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
         logits <- log(probs/(1-probs))
         line.types <- c(4,3,2,1,2,3,4)
         for (i in 1:7)
         panel.curve(exp((logits[i] + 24.58 - 39.55 * log(x))/31.93),lty=line.types[i])
#         panel.curve(exp((-2.94 + 24.58 - 39.55 * log(x))/31.93),lty=4)
#         panel.curve(exp((-2.02 + 24.58 - 39.55 * log(x))/31.93),lty=3)
#         panel.curve(exp((-0.85 + 24.58 - 39.55 * log(x))/31.93),lty=2)
#         panel.curve(exp((0 + 24.58 - 39.55 * log(x))/31.93),lty=1)
#         panel.curve(exp((0.85 + 24.58 - 39.55 * log(x))/31.93),lty=2)
#         panel.curve(exp((2.02 + 24.58 - 39.55 * log(x))/31.93),lty=3)
#         panel.curve(exp((2.94 + 24.58 - 39.55 * log(x))/31.93),lty=4)
}))


###################################################
### chunk number 42: 
###################################################
drop1(vaso.glm5,test="Chisq")


###################################################
### chunk number 43: 
###################################################
vaso.glm6 <- update(vaso.glm5,.~I(lv+lr))
summary(vaso.glm6)


###################################################
### chunk number 44: 
###################################################
vaso.glm7 <- update(vaso.glm5,.~offset(5*I(lv+lr)))
summary(vaso.glm7)


###################################################
### chunk number 45: vas.glm7
###################################################
#vaso$ofs <- log(0.05)+5*(vaso$lv+vaso$lr)
vaso.glm7 <- update(vaso.glm7,.~offset(log(0.05)+5*I(lv+lr))-1)
summary(vaso.glm7)


###################################################
### chunk number 46: vaso.fitted.lvr eval=FALSE
###################################################
## print(xyplot(Y~I(lv+lr),data=vaso, group=Y,pch=c(1,3), 
##              xlab="log(Volume)+log(Rate)",
## subset=!row.names(vaso)%in%c("32"),
## panel=function(x,y,...){
## panel.xyplot(x,y,...)
## panel.curve(exp(eta <- (log(0.05)+5*x))/(1+exp(eta)))}))


###################################################
### chunk number 47: 
###################################################
print(xyplot(Y~I(lv+lr),data=vaso, group=Y,pch=c(1,3), 
             xlab="log(Volume)+log(Rate)",
subset=!row.names(vaso)%in%c("32"),
panel=function(x,y,...){
panel.xyplot(x,y,...)
panel.curve(exp(eta <- (log(0.05)+5*x))/(1+exp(eta)))}))


###################################################
### chunk number 48:  eval=FALSE
###################################################
## data(bronchit,package="SMIR")


###################################################
### chunk number 49: 
###################################################
bronchit <- read.csv('bronchit.csv')


###################################################
### chunk number 50: bronchitis.response eval=FALSE
###################################################
## print(xyplot(poll~cig,data=bronchit,group=r,pch=c(1,3),col="black"))


###################################################
### chunk number 51: 
###################################################
print(xyplot(poll~cig,data=bronchit,group=r,pch=c(1,3),col="black"))


###################################################
### chunk number 52: bronchit.glm
###################################################
bronchit.glm <- glm(r~cig+poll,data=bronchit,family=binomial)
summary(bronchit.glm)
round(sum(resid(bronchit.glm,type='pearson')^2),1)


###################################################
### chunk number 53: bronchitis_fitted eval=FALSE
###################################################
## print(xyplot(poll~cig,data=bronchit,group=r,pch=c(1,3),col="black",
## panel=function(x,y,...){
## panel.xyplot(x,y,...)
##          probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
##          logits <- log(probs/(1-probs))
##          line.types <- c(4,3,2,1,2,3,4)
## for ( i in 1:7)
## panel.curve((logits[i]-coef(bronchit.glm)[1]-
## coef(bronchit.glm)[2]*x)/coef(bronchit.glm)[3],lty=line.types[i])                
##          panel.text(c(4,14,24),55,c("10","50","90"),font=2)
## }))


###################################################
### chunk number 54: 
###################################################
print(xyplot(poll~cig,data=bronchit,group=r,pch=c(1,3),col="black",
panel=function(x,y,...){
panel.xyplot(x,y,...)
         probs <- c(0.05,0.1,0.3,0.5,0.7,0.9,0.95)
         logits <- log(probs/(1-probs))
         line.types <- c(4,3,2,1,2,3,4)
for ( i in 1:7)
panel.curve((logits[i]-coef(bronchit.glm)[1]-
coef(bronchit.glm)[2]*x)/coef(bronchit.glm)[3],lty=line.types[i])                
         panel.text(c(4,14,24),55,c("10","50","90"),font=2)
}))


###################################################
### chunk number 55: 
###################################################
bronchit$w <- (abs((rbron <- resid(bronchit.glm,type="pearson")))>2)
cbind(bronchit,fitted=round(fitted(bronchit.glm),3),
residual=round(rbron,3))[bronchit$w,]


###################################################
### chunk number 56: poll.cig eval=FALSE
###################################################
## print(xyplot(poll~cig,data=bronchit,group=r,pch=c(1,3),col=1,
## cex=1.3,cex.axis=1.3,cex.lab=1.3,subset=w==TRUE))


###################################################
### chunk number 57: 
###################################################
print(xyplot(poll~cig,data=bronchit,group=r,pch=c(1,3),col=1,
cex=1.3,cex.axis=1.3,cex.lab=1.3,subset=w==TRUE))


###################################################
### chunk number 58: 
###################################################
(glm(r~poll,data=bronchit,family=binomial,subset=cig==0))
(glm(r~cig+poll,data=bronchit,family=binomial,subset=cig!=0))


###################################################
### chunk number 59: bronchit.glm1
###################################################
bronchit.glm1 <- glm(r~cig+poll+I(poll^2)+I(cig^2)+cig*poll,data=
bronchit,family=binomial)
summary(bronchit.glm1)


###################################################
### chunk number 60: 
###################################################
options(digits=4)
bronchit.glm2 <- update(bronchit.glm1,.~.-I(poll^2)-cig:poll)
summary(bronchit.glm2)


###################################################
### chunk number 61: 
###################################################
w <- abs(resid(bronchit.glm2,type="pearson"))>2
cbind(bronchit,fitted=round(fitted(bronchit.glm2),3),
residual=round(resid(bronchit.glm2,type="pearson"),3))[w,]


###################################################
### chunk number 62: 
###################################################
bronchit <- transform(bronchit,
                      fcig=factor(cut(bronchit$cig,c(0,0.1,1,3,5,8,Inf),
                        include.lowest=TRUE),ordered=FALSE),
                      fpol= factor(cut(bronchit$pol,c(50,55,57.5,60,62.5,65,Inf)),
                        ordered=FALSE))


###################################################
### chunk number 63: 
###################################################
(tr <- with(bronchit,tapply(r,list(fcig,fpol),sum,na.rm=TRUE)))
tn <- xtabs(~fcig+fpol,data=bronchit)


###################################################
### chunk number 64: 
###################################################
tn1 <- tn
names(dimnames(tn1)) <- NULL
tn1


###################################################
### chunk number 65: bronchitis.frequencies
###################################################
tbronchit <- data.frame(tn)
tbronchit$tr <- as.vector(tr)


###################################################
### chunk number 66: bronchit.glm3
###################################################
bronchit.glm3 <- glm(cbind(tr,(Freq-tr))~fcig+fpol,
                     data=tbronchit,family=binomial)
summary(bronchit.glm3)
round(resid(bronchit.glm3,type="pearson"),4)


###################################################
### chunk number 67: bronchit.test.main.effects
###################################################
drop1(bronchit.glm3,test="Chisq")


###################################################
### chunk number 68: bronchit.glm4
###################################################
bronchit.glm4 <- update(bronchit.glm3,.~fcig+unclass(fpol))
anova(bronchit.glm4,bronchit.glm3,test="Chisq")
round(coef(summary(bronchit.glm4)),4)
round(resid(bronchit.glm4),4)


###################################################
### chunk number 69: bronchit.glm5
###################################################
tbronchit$mcig <- tbronchit$fcig
levels(tbronchit$mcig) <- levels(tbronchit$mcig)[c(1,2,2,3,3,4)]
bronchit.glm5 <- update(bronchit.glm4,.~.-fcig+mcig)
round(coef(summary(bronchit.glm5)),4)


###################################################
### chunk number 70: bronchit.glm6
###################################################
bronchit$c <- bronchit$fcig
levels(bronchit$c) <- levels(bronchit$c)[c(1,2,2,3,3,4)]
bronchit.glm6 <- update(bronchit.glm1,.~c+unclass(poll))
round(coef(summary(bronchit.glm6)),4)


###################################################
### chunk number 71: 
###################################################
ghq <- read.csv('Data/ghq.csv')


###################################################
### chunk number 72:  eval=FALSE
###################################################
## data(ghq,package="SMIR")


###################################################
### chunk number 73: ghq.plot eval=FALSE
###################################################
## ghq <- transform(ghq,
##                  n=c+nc,p=c/(c+nc),
##                  sex= factor(sex,labels=c("men","women")))
## print(xyplot(p~ghq|sex,data=ghq,cex=1.2))#,
## #panel=function(...){
## #panel.xyplot(...)
## #panel.loess(...)}))


###################################################
### chunk number 74: 
###################################################
ghq <- transform(ghq,
                 n=c+nc,p=c/(c+nc),
                 sex= factor(sex,labels=c("men","women")))
print(xyplot(p~ghq|sex,data=ghq,cex=1.2))#,
#panel=function(...){
#panel.xyplot(...)
#panel.loess(...)}))


###################################################
### chunk number 75: ghq.glm
###################################################
ghq.glm <- glm(cbind(c,nc)~sex+ghq,data=ghq,family=binomial)
summary(ghq.glm)


###################################################
### chunk number 76: ghq.glm1
###################################################
(ghq.glm1 <- update(ghq.glm,.~.-sex))
anova(ghq.glm,ghq.glm1,test="Chisq")


###################################################
### chunk number 77: fig.ghq.fitted eval=FALSE
###################################################
## print(xyplot(p~ghq,data=ghq,col=1,
## panel=function(x,y,...){
## panel.xyplot(x,y,...)
## panel.curve(exp(eta <- (coef(ghq.glm1)[1]+coef(ghq.glm1)[2]*x))/(1+exp(eta)),
## ,lwd=1.5)}))


###################################################
### chunk number 78: 
###################################################
print(xyplot(p~ghq,data=ghq,col=1,
panel=function(x,y,...){
panel.xyplot(x,y,...)
panel.curve(exp(eta <- (coef(ghq.glm1)[1]+coef(ghq.glm1)[2]*x))/(1+exp(eta)),
,lwd=1.5)}))


###################################################
### chunk number 79: ghq.glm2
###################################################
ghq.glm2 <- update(ghq.glm,.~.+sex:ghq)
summary(ghq.glm2)


###################################################
### chunk number 80: ghq.women.glm
###################################################
ghq.women.glm <- glm(cbind(c,nc)~ghq, data=ghq, family=binomial,
subset=sex=="women")#,control=glm.control(maxit=30)) 
summary(ghq.women.glm)


###################################################
### chunk number 81: 
###################################################
ghq.men.glm <- glm(cbind(c,nc)~ghq, data=ghq, family=binomial,
subset=sex=="men",control=glm.control(maxit=30)) 
summary(ghq.men.glm)


###################################################
### chunk number 82: fig.profile.likelihoods eval=FALSE
###################################################
## pseq <- seq(from=0.001,to=0.998,by=0.001)
## rprob <- function(x0,pseq){
## logl <- NULL
## for (p in pseq){
## theta <- rep(log(p/(1-p)),nrow(ghq))
## xx0 <- ghq$ghq-rep(x0,nrow(ghq))
## fit <- glm(cbind(c,nc)~offset(theta)+xx0-1,data=ghq,
## family=binomial,na.action=na.omit)
## logl <- c(logl,with(ghq,sum(dbinom(c,(nc+c),fitted(fit),log=TRUE))))
## }
## exp(logl-max(logl))
## }
## print(xyplot(c(0,1)~c(0,1),type="n",xlab="",ylab="",
## panel=function(x,y,...){
## for (i in 0:6) panel.curve(rprob(i,x),from=0.001,to=0.999)
## panel.text(c(0.05,0.14,0.4,0.75,0.94),rep(1,5),0:6)
## }))


###################################################
### chunk number 83: 
###################################################
pseq <- seq(from=0.001,to=0.998,by=0.001)
rprob <- function(x0,pseq){
logl <- NULL
for (p in pseq){
theta <- rep(log(p/(1-p)),nrow(ghq))
xx0 <- ghq$ghq-rep(x0,nrow(ghq))
fit <- glm(cbind(c,nc)~offset(theta)+xx0-1,data=ghq,
family=binomial,na.action=na.omit)
logl <- c(logl,with(ghq,sum(dbinom(c,(nc+c),fitted(fit),log=TRUE))))
}
exp(logl-max(logl))
}
print(xyplot(c(0,1)~c(0,1),type="n",xlab="",ylab="",
panel=function(x,y,...){
for (i in 0:6) panel.curve(rprob(i,x),from=0.001,to=0.999)
panel.text(c(0.05,0.14,0.4,0.75,0.94),rep(1,5),0:6)
}))


###################################################
### chunk number 84: ghq.profile.likelihood
###################################################
ghqeta <- predict(ghq.glm1,new=data.frame(ghq=c(0,1,2,3,4,5,6)),type="link",se=T)
predp <- 1/(1+exp(-ghqeta$fit))
lleta<- ghqeta$fit-1.96*ghqeta$se.fit
llp <- 1/(1+exp(-lleta))
uleta <- ghqeta$fit+1.96*ghqeta$se.fit
ulp <- 1/(1+exp(-uleta))
predest <- matrix(NA,nrow=7,ncol=3)
for (i in 1:7) {
R <- rprob(i-1,pseq)
predest[i,2] <- pseq[which.max(R)]
a <- pseq[R>0.15]
predest[i,1] <- a[1]
predest[i,3] <- max(a)
}


###################################################
### chunk number 85: conditional_profile eval=FALSE
###################################################
## clog <- data.frame(r=c(11,0),n=c(11,1),group=factor(c(1,2)))
## #clog.glm <- glm(cbind(r,(n-r))~group,data=clog,family=binomial)
## beta1 <- seq(-1,10,by=0.1)
## pl <- NULL
## for (b in beta1){
## z <- b*c(1,0)
## logl <- logLik(glm(cbind(r,(n-r))~offset(z),family=binomial,data=clog))
## pl  <- c(pl,logl)
## }
## pl <- exp(pl - max(pl))
## #
## #b0seq <- seq(-3,3,len=100)
## #b1seq <- seq(-1,10,len=100)
## #L <- matrix(NA,nrow=100,ncol=100)
## #for (i in 1:100){
## #  for (j in 1:100){
## #    eta1 <- b0seq[i]+ b1seq[j]
## #    eta2 <- b0seq[i]
## #    p1 <- exp(eta1)/(1+exp(eta1))
## #    p2 <- exp(eta2)/(1+exp(eta2))
## #    theta <- exp(beta1[j])
## #    L[i,j] <- 11*log((p1/(1-p1))/(p2/(1-p2)))+
## #      11*log(p2/(1-p2))+11*log(1-p1)+log(1-p2)
## #  }
## #}
## #L <- L-max(L)
## #pl <- apply(L,1,max)
## #
## print(xyplot(pl~beta1,type="l",
## panel=function(x,y,...){
## panel.xyplot(x,y,...)
## panel.curve(exp(x)/(11+exp(x)),from=min(beta1),to=max(beta1),lty="dotted")
## }))


###################################################
### chunk number 86: 
###################################################
clog <- data.frame(r=c(11,0),n=c(11,1),group=factor(c(1,2)))
#clog.glm <- glm(cbind(r,(n-r))~group,data=clog,family=binomial)
beta1 <- seq(-1,10,by=0.1)
pl <- NULL
for (b in beta1){
z <- b*c(1,0)
logl <- logLik(glm(cbind(r,(n-r))~offset(z),family=binomial,data=clog))
pl  <- c(pl,logl)
}
pl <- exp(pl - max(pl))
#
#b0seq <- seq(-3,3,len=100)
#b1seq <- seq(-1,10,len=100)
#L <- matrix(NA,nrow=100,ncol=100)
#for (i in 1:100){
#  for (j in 1:100){
#    eta1 <- b0seq[i]+ b1seq[j]
#    eta2 <- b0seq[i]
#    p1 <- exp(eta1)/(1+exp(eta1))
#    p2 <- exp(eta2)/(1+exp(eta2))
#    theta <- exp(beta1[j])
#    L[i,j] <- 11*log((p1/(1-p1))/(p2/(1-p2)))+
#      11*log(p2/(1-p2))+11*log(1-p1)+log(1-p2)
#  }
#}
#L <- L-max(L)
#pl <- apply(L,1,max)
#
print(xyplot(pl~beta1,type="l",
panel=function(x,y,...){
panel.xyplot(x,y,...)
panel.curve(exp(x)/(11+exp(x)),from=min(beta1),to=max(beta1),lty="dotted")
}))


###################################################
### chunk number 87: prenatal
###################################################
prenatal <- data.frame(death=c(3,4,17,2),total=c(179,297,214,25),
clinic=factor(c("A","A","B","B")),
attend=factor(c("<1",">1","<1",">1" )))
prenatal.glm <- glm(cbind(death,(total-death))~attend+clinic,
data=prenatal,family=binomial)
prenatal.glm1 <- update(prenatal.glm,.~clinic+attend)
#anova(prenatal.glm)
#anova(prenatal.glm1)


###################################################
### chunk number 88: 
###################################################
fp <- round(fitted(update(prenatal.glm,.~clinic)),4)


###################################################
### chunk number 89: 
###################################################
chd <- read.csv('Data/chd.csv')


###################################################
### chunk number 90:  eval=FALSE
###################################################
## data(chd,package="SMIR")


###################################################
### chunk number 91: chd.plots
###################################################
chd <- transform(chd,
                 p= chd$r/chd$t,
                 bp=factor(chd$bp,labels=c("<127","127-146","147-166",">166")),
                 chol=factor(chd$chol,labels=c("<200","200-219","220-259",">259")))
#chd$p <- chd$r/chd$t
#chd$bp <- factor(chd$bp,labels=c("<127","127-146","147-166",">166"))
#chd$chol <- factor(chd$chol,labels=c("<200","200-219","220-259",">259"))
print(xyplot(p~bp|chol,data=chd,type="b",xlab="blood pressure",
scales=list(x=list(rot=45))))


###################################################
### chunk number 92: 
###################################################
chd <- transform(chd,
                 p= chd$r/chd$t,
                 bp=factor(chd$bp,labels=c("<127","127-146","147-166",">166")),
                 chol=factor(chd$chol,labels=c("<200","200-219","220-259",">259")))
#chd$p <- chd$r/chd$t
#chd$bp <- factor(chd$bp,labels=c("<127","127-146","147-166",">166"))
#chd$chol <- factor(chd$chol,labels=c("<200","200-219","220-259",">259"))
print(xyplot(p~bp|chol,data=chd,type="b",xlab="blood pressure",
scales=list(x=list(rot=45))))


###################################################
### chunk number 93: chd.logit
###################################################
chd.glm1 <- glm(cbind(r,(t-r))~1,data=chd,family=binomial)
chd.glm2 <- update(chd.glm1,.~.+bp)
chd.glm3 <- update(chd.glm2,.~.+chol)
anova(chd.glm1,chd.glm2,chd.glm3)


###################################################
### chunk number 94: 
###################################################
round(resid(chd.glm3,type="pearson"),4)


###################################################
### chunk number 95: 
###################################################
round(summary(chd.glm3)$coeff,4)


###################################################
### chunk number 96: chd.glm4
###################################################
#chd <- transform(chd,ch=factor(chol,levels=levels(chol)[c(1,1,2,3)]),
#                       b=factor(bp,levels=levels(bp)[c(1,1,2,3)]))
chd$ch <- chd$chol
levels(chd$ch) <- c(1,1,2,3)
chd$b <- chd$bp
levels(chd$b) <- c(1,1,2,3)
chd.glm4 <- update(chd.glm3,.~ch+b)
anova(chd.glm3,chd.glm4)
round(summary(chd.glm4)$coeff,4)
round(resid(chd.glm4,type="pearson"),4)


###################################################
### chunk number 97: chd.glm5
###################################################
chd.glm5 <- update(chd.glm4,.~unclass(ch)+unclass(b))
anova(chd.glm5,chd.glm4)
round(summary(chd.glm5)$coeff,4)
round(resid(chd.glm5,type="pearson"),4)


###################################################
### chunk number 98: chd.glm6
###################################################
chd <- transform(chd,score=unclass(b)+unclass(ch))
chd.glm6 <- update(chd.glm5,.~score)
anova(chd.glm6,chd.glm5)
round(summary(chd.glm6)$coeff,4)
round(resid(chd.glm6,type="pearson"),4)


###################################################
### chunk number 99: chd.glm7
###################################################
chd$ofs <- log(2)*chd$score
chd.glm7 <- update(chd.glm6,.~offset(ofs))
anova(chd.glm7,chd.glm6)
round(summary(chd.glm7)$coeff,4)


###################################################
### chunk number 100: chd.glm8
###################################################
chd$ofs2 <- -5.075+chd$ofs
chd.glm8 <- update(chd.glm7,.~offset(ofs2)-1)
anova(chd.glm8,chd.glm7)
#round(summary(chd.glm8)$coeff,4)


###################################################
### chunk number 101: 
###################################################
newdf <- expand.grid(bp=seq(1,4,len=100),
                     chol=levels(chd$chol))
newdf <- transform(newdf, 
                   b = ifelse(bp>=2, bp - 1, 1),
                   ch = ifelse(unclass(chol)>=2, 
                     unclass(chol)-1, 2))
newdf <- transform(newdf,
                   score = b + ch)
newdf <- transform(newdf,
                   ofs = log(2)*score)
newdf <- transform(newdf,
                   ofs2 = -5.075 + ofs)
newdf <- transform(newdf,
                   fitted = exp(ofs2)/(1+exp(ofs2)))


###################################################
### chunk number 102: chd.prob eval=FALSE
###################################################
## #chd$fval <- fitted(chd.glm8)
##                 
## #
## newdf <- newdf[order(newdf$score),]
## print(xyplot(fitted~(score-1),data=newdf,
## type="l",col="black",xlab="score",ylab="p",
## panel=function(x,y){
## panel.xyplot(x,y,type="l")
## lpoints((chd$score-1),chd$p,pch=chd$chol)}))
## #
## #print(xyplot(p~bp|chol,data=chd,
## #xlab="bp",ylab="p",groups=newdf$chol,
## #scales=list(x=list(at=c(1,2,3,4),lab=1:4)),
## #panel=function(x,y,subscripts,...){
## #panel.xyplot(x[subscripts],y[subscripts],...)
## #llines(newdf$bp[subscripts],newdf$p[subscripts])}))
## #print(xyplot(p+fitted(chd.glm8)~bp|chol,
## #             data=chd,xlab="blood pressure",
## #             scales=list(x=list(rot=45)),
## #             ylab="p",
## #             panel=panel.superpose.2,
## #             type=c("p","l")))


###################################################
### chunk number 103: 
###################################################
#chd$fval <- fitted(chd.glm8)
                
#
newdf <- newdf[order(newdf$score),]
print(xyplot(fitted~(score-1),data=newdf,
type="l",col="black",xlab="score",ylab="p",
panel=function(x,y){
panel.xyplot(x,y,type="l")
lpoints((chd$score-1),chd$p,pch=chd$chol)}))
#
#print(xyplot(p~bp|chol,data=chd,
#xlab="bp",ylab="p",groups=newdf$chol,
#scales=list(x=list(at=c(1,2,3,4),lab=1:4)),
#panel=function(x,y,subscripts,...){
#panel.xyplot(x[subscripts],y[subscripts],...)
#llines(newdf$bp[subscripts],newdf$p[subscripts])}))
#print(xyplot(p+fitted(chd.glm8)~bp|chol,
#             data=chd,xlab="blood pressure",
#             scales=list(x=list(rot=45)),
#             ylab="p",
#             panel=panel.superpose.2,
#             type=c("p","l")))


###################################################
### chunk number 104: byssinosis
###################################################
byssinosis <- read.csv('byssinosis.csv')


###################################################
### chunk number 105: 
###################################################
byssinosis <- transform(byssinosis,
dust=factor(dust,labels=c("most","less","least")),
race=factor(race,labels=c("white","non-white")),
emp =factor(emp,labels=c("<10","10-20",">20")),
sex=factor(sex,labels=c("male","female")),
smok=factor(smok,labels=c("smoker","non")),
n=yes+no,
p=yes/(yes+no))


###################################################
### chunk number 106: byssinosis eval=FALSE
###################################################
## data(byssinosis,package="SMIR")
## byssinosis <- transform(byssinosis,
##                         n = yes + no,
##                         p = yes / (yes + no))


###################################################
### chunk number 107: byssinosis.table
###################################################
yestot <- with(byssinosis,tapply(yes,
list(smok=smok,sex=sex,emp=emp,race=race,dust=dust),sum))
ntot <- with(byssinosis,tapply(n,
list(smok=smok,sex=sex,emp=emp,race=race,dust=dust),sum))
round(with(byssinosis,ftable(yestot/ntot,row.vars=c(1,2,3),
col.vars=c(4,5))),4)


###################################################
### chunk number 108: 
###################################################
for (i in 1:5){
cat("\t",names(byssinosis)[i],"\n")
print(round(tapply(byssinosis$yes,byssinosis[,i],sum,na.rm=TRUE)/
tapply(byssinosis$n,byssinosis[,i],sum,na.rm=TRUE),4))}
with(byssinosis,round(sum(yes)/sum(n),4))


###################################################
### chunk number 109: 
###################################################
with(byssinosis,round(tapply(yes,list(race=race,dust=dust),sum)/
tapply(n,list(race=race,dust=dust),sum),4))


###################################################
### chunk number 110: 
###################################################
with(byssinosis,tapply(n,list(race=race,dust=dust),sum))


###################################################
### chunk number 111: byssinosis.anova
###################################################
bys.glm1 <- glm(cbind(yes,no)~(dust+emp+smok+sex+race)^3,data=byssinosis,
family=binomial)


###################################################
### chunk number 112: 
###################################################
library(xtable)
xtable(anova(bys.glm1),caption="Analysis of deviance table, byssinosis",
       label="tab:byssinosis.anova")


###################################################
### chunk number 113: byssinosis.glm2
###################################################
bys.glm2a <- glm(cbind(yes,no)~(dust+emp+smok+sex+race)^2,
data=byssinosis,family=binomial)
anova(bys.glm2a,bys.glm2 <- glm(cbind(yes,no)~(dust+emp+smok+sex+race)^2+dust:emp:sex,
data=byssinosis,family=binomial))


###################################################
### chunk number 114: byssinosis.glm3
###################################################
bys.glm3 <- update(bys.glm2,.~.-dust:emp:sex)
#round(summary(bys.glm3)$coeff,4)[9:27,]
drop1(bys.glm3,test="Chisq")


###################################################
### chunk number 115: byssinosis.glm4
###################################################
bys.glm4 <- update(bys.glm3,.~.-smok:(emp+sex+race)-sex:race)
#drop1(bys.glm4,test="Chisq")
round(summary(bys.glm4)$coeff,4)[9:22,]


###################################################
### chunk number 116: byssinosis.glm5
###################################################
bys.glm5 <- update(bys.glm4,.~.-emp:sex)
#drop1(bys.glm5,test="Chisq")
round(summary(bys.glm5)$coeff,4)[9:20,]


###################################################
### chunk number 117: byssinosis.glm6
###################################################
bys.glm6 <- update(bys.glm5,.~.-dust:race)
#drop1(bys.glm6,test="Chisq")
round(summary(bys.glm6)$coeff,4)[9:18,]


###################################################
### chunk number 118: byssinosis.glm7
###################################################
bys.glm7 <- update(bys.glm6,.~.-dust:smok)
#drop1(bys.glm7,test="Chisq")
round(summary(bys.glm7)$coeff,4)[9:16,]


###################################################
### chunk number 119: byssinosis.glm8
###################################################
byssinosis$d3e2 <- ifelse((byssinosis$dust=="least")&(byssinosis$emp=="10-20"),
1,0)
bys.glm8 <- update(bys.glm7,.~.-dust:emp+d3e2)
anova(bys.glm8,bys.glm7)
bys.glm9 <- update(bys.glm8,.~.-d3e2)
anova(bys.glm9,bys.glm8)
round(summary(bys.glm9)$coeff,4)[9:12,]


###################################################
### chunk number 120: byssinosis.glm10
###################################################
bys.glm10 <- update(bys.glm9,.~.-emp:race)
anova(bys.glm10,bys.glm9)
summary(bys.glm10)


###################################################
### chunk number 121: byssinosis.glm11
###################################################
bys.glm11 <- update(bys.glm10,.~.-race)
anova(bys.glm11,bys.glm10)
summary(bys.glm11)


###################################################
### chunk number 122: byssinosis.glm12
###################################################
byssinosis$d2 <- byssinosis$dust
levels(byssinosis$d2) <- c(1,2,2)
bys.glm12 <- update(bys.glm11,.~.-dust:sex+d2+d2:sex)
anova(bys.glm12,bys.glm11)
summary(bys.glm12)


###################################################
### chunk number 123: byssinosis.glm13
###################################################
bys.glm13 <- update(bys.glm12,.~.-dust+d2)
anova(bys.glm13,bys.glm12)
summary(bys.glm13)


###################################################
### chunk number 124: byssinosis.glm14
###################################################
byssinosis$e2 <- byssinosis$emp
levels(byssinosis$e2) <- c(1,2,2)
bys.glm14 <- update(bys.glm13,.~.-emp+e2)
anova(bys.glm14,bys.glm13)
summary(bys.glm14)


###################################################
### chunk number 125: byssinosis.glm15
###################################################
byssinosis$line <- as.numeric(unclass(byssinosis$emp))
bys.glm15 <- update(bys.glm14,.~.-e2+line)
anova(bys.glm14,bys.glm15)
summary(bys.glm15)


###################################################
### chunk number 126: byssinosis.fitted.proportion eval=FALSE
###################################################
## fp <- fitted(bys.glm15)
## round(with(byssinosis,ftable(tapply(fp,
## list(smok=smok,emp=emp,dust=d2,sex=sex),
## mean),row.vars=c(1,2),
## col.vars=c(3,4))),3)
## round(with(byssinosis,ftable(
##                              tapply(yes,list(smok=smok,emp=emp,dust=d2,sex=sex),
##                                     sum)/
##                              tapply((yes+no),list(smok=smok,emp=emp,dust=d2,sex=sex),
##                                     sum),
##                              row.vars=c(1,2),
##                              col.vars=c(3,4))),3)


###################################################
### chunk number 127: byssinosis.fitted.counts eval=FALSE
###################################################
## with(byssinosis,ftable(tapply(n,
## list(smok=smok,emp=emp,dust=d2,sex=sex),sum),row.vars=c(1,2)))