################ Question 2 ##################

m<-matrix(scan("c:/keith/teaching/smalldatasets/plague.dat"),10,2,T)
all<-m[1:9,1]
plague<-m[1:9,2]
p<-plague/all
week<-1:9
postscript("c:/keith/teaching/523/2008/fig1.eps",horizontal=FALSE)
par(mfrow=c(2,2))
plot(week,p,main="Fig 2.1 Plague mortality")
glm1<-glm(p~week,family=binomial,weight=all)
summary(glm1)
week2<-week^2
glm2<-glm(p~week+week2,family=binomial,weight=all)
summary(glm2)
week3<-week^3
glm3<-glm(p~week+week2+week3,family=binomial,weight=all)
summary(glm3)
lines(week,fitted(glm2))

################ Question 3 ##################

m<-matrix(scan("c:/keith/teaching/smalldatasets/hillrace.dat"),35,3,T)
dist<-m[,1]
height<-m[,2]
time<-m[,3]
plot(dist,time,main="Fig 3.1: Time vs dist")
plot(height,time,main="Fig 3.2: Time vs height")
glm1<-glm(time~dist+height)
summary(glm1)
plot(fitted(glm1),resid(glm1),main="Fig 3.3: Residuals")
abline(0,0)
dev.off()
cbind(time,fitted(glm1),resid(glm1))
time0<-time[-18]
dist0<-dist[-18]
height0<-height[-18]
glm2<-glm(time0~dist0+height0)
summary(glm2)
cbind(25:35, -qt(0.05/2/33,25:35), -qt(0.05/2/34,25:35), -qt(0.05/2/35,25:35))

postscript("c:/keith/teaching/523/2008/fig2.eps",horizontal=FALSE)
par(mfrow=c(2,2))
plot(fitted(glm2),resid(glm2),main="Fig 3.4: Omitting obs 18")
abline(0,0)
glm3<-glm(time0~dist0+height0+height0:dist0)
summary(glm3)
plot(fitted(glm3),resid(glm3),main="Fig 3.5: With interaction")
abline(0,0)

vl<-predict.lm(glm3,se.fit=T)$se.fit^2
sc<-summary(glm3)$dispersion
z<-resid(glm3)/sqrt(sc-vl)
z
hist(z,main="Figure 3.6 Normalized residuals, z")
u<-pnorm(z)
hist(u,main="Figure 3.7 Uniforms, u")
dev.off()
uhat<-((1:34)-0.5)/34
u<-sort(u)
postscript("c:/keith/teaching/523/2008/fig3.eps",horizontal=FALSE)
par(mfrow=c(2,2))
plot(uhat,u,xaxs="i",xlim=c(0,1),yaxs="i",
ylim=c(0,1),main="Figure 3.8 P-P plot")
abline(0,1)
zhat<-qnorm(uhat)
z<-sort(z)
plot(zhat,z,main="Figure 3.9 Q-Q plot")
abline(0,1)
u-uhat

################ Question 4 ##################

m<-matrix(scan("c:/keith/teaching/smalldatasets/remiss.dat"),42,3,T)
pair<-factor(m[,1])
time<-m[,2]
censoring<-m[,3]
treatment<-factor(c(rep(c(1,2,2,1),10),c(1,2)))
rate<-censoring/time
summary(glm(rate~treatment,family=poisson,weight=time))
summary(glm(rate~treatment+pair,family=poisson,weight=time))

################ Question 5 ##################

m<-matrix(scan("c:/keith/teaching/smalldatasets/ear.dat"),287,5,T)
frq<-factor(m[,1])
loc<-factor(m[,2])
Age<-factor(m[,3])
sex<-factor(m[,4])
inf<-m[,5]
summary(glm(inf~frq+loc+Age+sex,family=poisson))
hist(fitted(glm(inf~frq+loc+Age+sex,family=poisson)),main="Figure 5.1 Fitted values")
dev.off()
age<-m[,3]
summary(glm(inf~frq+loc+age+sex,family=poisson))
summary(glm(inf~frq+loc+frq:loc,family=poisson))
summary(glm(inf~frq+loc,family=poisson))

inf1<-inf+1
y<-1/inf1
summary(glm(y~frq+loc+Age+sex,family=binomial,weight=inf1))
summary(glm(y~frq+loc+age+sex,family=binomial,weight=inf1))
summary(glm(y~frq+loc+frq:loc,family=binomial,weight=inf1))
summary(glm(y~frq+loc,family=binomial,weight=inf1))