### R code from vignette source 'AP-AC-s.rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: AP-AC-s.rnw:3-9
###################################################
options( width=90,
         prompt=" ", continue=" ", # Absence of prompts makes it easier for
                                   # students to copy from the final pdf document
         SweaveHooks=list( fig=function()
         par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) )
library( Epi )


###################################################
### code chunk number 2: AP-AC-s.rnw:24-26
###################################################
library( Epi )
print( sessionInfo(), l=F )


###################################################
### code chunk number 3: AP-AC-s.rnw:33-39
###################################################
lung <- read.table( "../data/lung5-M.txt", header=T )
with( lung , table( A ) )
with( lung , table( P ) )
round( ftable( addmargins( xtabs( cbind(D,Y/1000) ~  A + P, data = lung ), 
                           margin = 1 ), 
               row.vars=c(3,1) ), 0 )


###################################################
### code chunk number 4: AP-AC-s.rnw:51-56
###################################################
ap.1 <- glm( D ~ factor(A) + factor(P),
             offset = log(Y/1000),
             family = poisson, 
               data = lung )
summary( ap.1 )


###################################################
### code chunk number 5: AP-AC-s.rnw:69-70
###################################################
round( ci.exp(ap.1), 2 )


###################################################
### code chunk number 6: AP-AC-s.rnw:76-81
###################################################
ap.0 <- glm( D ~ -1 + factor(A) + factor(P),
             offset = log(Y/1000),
             family = poisson, 
               data = lung )
round( ci.exp(ap.0), 3 )


###################################################
### code chunk number 7: AP-AC-s.rnw:92-96
###################################################
ap.3 <- glm( D ~ factor(A) - 1 + relevel(factor(P),"1968"),
             offset = log(Y/1000),
             family = poisson, 
               data = lung )


###################################################
### code chunk number 8: AP-AC-s.rnw:99-100
###################################################
round( ci.exp( ap.3 ), 3 )


###################################################
### code chunk number 9: AP-AC-s.rnw:107-108
###################################################
( ap.cf <- ci.exp( ap.3, subset="A" ) )


###################################################
### code chunk number 10: agesh
###################################################
matshade( seq(40,85,5)+2.5, ci.exp( ap.3, subset="A" ),
          type="l", lty=1, lwd=1, log="y", col=1, plot=TRUE,
          xlab="Age",
          ylab="Male lung cancer incicdence rate per 1000 PY")


###################################################
### code chunk number 11: AP-AC-s.rnw:129-130
###################################################
( RR.cf <- ci.exp( ap.3, subset="P" ) )


###################################################
### code chunk number 12: AP-AC-s.rnw:135-136
###################################################
( RR.cf <- rbind( RR.cf[1:5,], 1, RR.cf[6:10,] ) )


###################################################
### code chunk number 13: APrrLung
###################################################
matshade( seq(1943,1993,5)+2.5, RR.cf,
          lty=1, lwd=1, log="y", col=1, plot=TRUE,
          xlab="Calendar time", ylab="Rate ratio rel. to 1968--72")
abline( h=1, v=1970.5, lty=3 )


###################################################
### code chunk number 14: AP-AC-s.rnw:158-162
###################################################
nd <- data.frame( A = seq(40,85,5),
                  P = 1968,
                  Y = 1000 )
( rt <- ci.pred( ap.3, nd ) )


###################################################
### code chunk number 15: AP-AC-s.rnw:177-182
###################################################
ap.x <- glm( D ~ -1 + A + P,
             offset = log(Y),
             family = poisson, 
               data = transform(lung,A=factor(A),P=factor(P)) )
summary( ap.x )


###################################################
### code chunk number 16: AP-AC-s.rnw:189-192
###################################################
nd <- data.frame( P = seq(1943,1993,5) )
nr <- data.frame( P = 1968 )
( rrx <- ci.exp( ap.x, list(nd,nr) ) )


###################################################
### code chunk number 17: AP-AC-s.rnw:213-214
###################################################
with( lung, table( P-A ) )


###################################################
### code chunk number 18: AP-AC-s.rnw:224-229
###################################################
ac.0 <- glm( D ~ A + C,
             offset = log(Y),
             family = poisson, 
               data = transform(lung,A=factor(A),C=factor(P-A)) )
summary( ac.0 )


###################################################
### code chunk number 19: AP-AC-s.rnw:238-243
###################################################
ac.r <- glm( D ~ -1 + A + relevel(C,'1908'),
             offset = log(Y),
             family = poisson, 
               data = transform(lung,A=factor(A),C=factor(P-A)) )
round( ci.exp( ac.r ), 3 )


###################################################
### code chunk number 20: AP-AC-s.rnw:267-271
###################################################
ndc <- data.frame( C=seq(1858,1953,5) )
ndr <- data.frame( C=1908 )
try( RR.C <- ci.exp( ac.r, list(ndc,ndr) ) )
   ( RR.C <- ci.exp( ac.0, list(ndc,ndr) ) )


###################################################
### code chunk number 21: cohRR
###################################################
matshade( ndc$C, RR.C, log='y', plot=TRUE,
          xlab="Date of birth", ylab="Lung cancer incidence RR" )
abline( h=1, v=1908, lty=3 )


###################################################
### code chunk number 22: AP-AC-s.rnw:284-285
###################################################
ai.coh <- ci.pred( ac.0, data.frame(A=factor(seq(40,85,5)),C='1908',Y=1000) )


###################################################
### code chunk number 23: Aincmp
###################################################
matshade( seq(40,85,5), cbind( rt, ai.coh ), col=c("black","blue"),
          log="y", plot=TRUE )
abline( v=60, lty=3 )