### R code from vignette source 'age-drift-sol.rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: age-drift-sol.rnw:2-7
###################################################
options( width=90,
         prompt=" ", continue=" ",
         SweaveHooks=list( fig=function()
         par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6), bty="n", las=1 ) )
library( Epi )


###################################################
### code chunk number 2: age-drift-sol.rnw:22-25
###################################################
lung <- read.table( "../data/lung5-M.txt", header=T )
lung$C <- lung$P - lung$A
table( lung$C )


###################################################
### code chunk number 3: age-drift-sol.rnw:42-47
###################################################
mp <- glm( D ~ -1 + factor(A) + I(P-1968),
           offset = log(Y),
           family = poisson, 
             data = lung )
round( ci.lin( mp ), 3 )


###################################################
### code chunk number 4: age-drift-sol.rnw:55-59
###################################################
lung <- transform( lung, A=A+2.5, P=P+2.5 )
mp <- glm( D ~ -1 + factor(A) + I(P-1970.5) + offset( log(Y) ),
           family=poisson, data=lung )
ci.lin( mp )[,1:2]


###################################################
### code chunk number 5: age-drift-sol.rnw:69-72
###################################################
mc <- glm( D ~ -1 + factor(A) + I(C-1908) + offset( log(Y) ),
           family=poisson, data=lung )
ci.lin( mc )[,1:2]


###################################################
### code chunk number 6: age-drift-sol.rnw:83-85
###################################################
c( summary( mp )$deviance,
   summary( mc )$deviance )


###################################################
### code chunk number 7: age-drift-sol.rnw:103-108
###################################################
ap <- ci.lin( mp )[1:10,1]
ac <- ci.lin( mc )[1:10,1]
c.sl <- ci.lin( mc )[11,1]
a.pt <- seq(40,85,5)
cbind( ap, ac + c.sl*(62.5-a.pt) )


###################################################
### code chunk number 8: rates
###################################################
matshade( a.pt + 2.5, cbind( ci.exp( mp, subset="A" ),
                             ci.exp( mc, subset="A" ) ) * 10^5, plot=TRUE,
          log="y", xlab="Age", ylab="Lung cancer incidence rates / 100,000",
          lty=1, lwd=1, col=c("black","blue") )


###################################################
### code chunk number 9: RRs
###################################################
p.pt <- seq(min(lung$P),max(lung$P),,10)+2.5
c.pt <- seq(min(lung$C),max(lung$C),,10)
ctr.p <- cbind( p.pt - 1970.5 )
ctr.c <- cbind( c.pt - 1908   )
matshade( c.pt, ci.exp( mc, subset="C", ctr.mat=ctr.c ), plot=TRUE,
          log="y", xlab="Calendar time", ylab="Rate ratio", xlim=c(1850,2000),
          type="l", lty=1, lwd=1, col="blue" )
matshade( p.pt, ci.exp( mp, subset="P", ctr.mat=ctr.p ),
          type="l", lty=1, lwd=1, col="black" )
abline( h=1, lty=3 )
points( c(1908,1970.5), c(1,1), pch=16 )