### R code from vignette source 'age-coh-sol.rnw'

###################################################
### code chunk number 1: age-coh-sol.rnw:10-15
###################################################
library(Epi)
lung <- read.table( "../data/lung5-M.txt", header=T )
lung$C <- lung$P - lung$A
attach( lung )
table( C )


###################################################
### code chunk number 2: age-coh-sol.rnw:27-30
###################################################
ac.1 <- glm( D ~ factor(A) + factor(C) + offset(log(Y)),
             family=poisson, data=lung )
summary( ac.1 )


###################################################
### code chunk number 3: age-coh-sol.rnw:47-49
###################################################
ac.2 <- glm( D ~ factor(A) - 1 + relevel(factor(C),"1908") + offset(log(Y)),
                 family=poisson, data=lung )


###################################################
### code chunk number 4: agerates
###################################################
age.cf <- ci.lin( ac.2, subset="A", Exp=TRUE)[,5:7]
matplot( as.numeric(levels(factor(A)))+2.5, age.cf,
         log="y", type="l", lty=1, lwd=c(3,1,1), col="black" )


###################################################
### code chunk number 5: RR
###################################################
RR.cf <- ci.lin( ac.2, subset="C", Exp=TRUE )[,5:7]
wh <- grep( "1908", levels(factor(C)) ) - 1
RR.cf <- rbind( RR.cf[1:wh,], c(1,1,1), RR.cf[-(1:wh),] )
RR.cf
matplot( as.numeric(levels(factor(C))), RR.cf,
         type="l", log="y", lwd=c(3,1,1), lty=1, col="black" )


###################################################
### code chunk number 6: ageRR
###################################################
alim <- range( A ) + c(0,5)
clim <- range( C ) + c(-2.5,2.5)
# Compute limits explicitly
rlim <- range(age.cf*10^5)*c(1/1.05,1.05)
RRlim <- 10^(log10(rlim)-ceiling(mean(log10(rlim)))) / 2
# Determine relative width of plots
layout( rbind( c(1,2) ), widths=c(diff(alim),diff(clim)) )
# No space on the sides of the plots, only outer space
par( mar=c(4,0,1,0), oma=c(0,4,0,4), mgp=c(3,1,0)/1.5, las=1 )
matplot( as.numeric(levels(factor(A)))+2.5, age.cf*10^5,
         type="l", lwd=c(3,1,1), lty=1, col="black",
         log="y", xaxs="i", xlim=alim, xlab="Age", ylim=rlim )
mtext( "Male lung cancer per 100,000", las=0, side=2, outer=T, line=2.5 )
matplot( as.numeric(levels(factor(C))), RR.cf,
         type="l", lwd=c(3,1,1), lty=1, col="black",
         log="y", xlab="Date of birth", xlim=clim, yaxt="n", ylim=RRlim, ylab="" )
abline( h=1 )
points( 1908, 1, pch=1, lwd=3 )
axis( side=4 )
mtext( "Rate ratio", side=4, outer=T, las=0, line=2.5 )


###################################################
### code chunk number 7: ageinc-comp
###################################################
load( file = "../data/age-per-est.Rdata" )
matplot( as.numeric(levels(factor(A)))+2.5, age.cf,
         log="y", type="l", lty=1, lwd=c(3,1,1), col="black" )
matlines( age.pt, arates,
          type="l", lty=1, lwd=c(3,1,1), col="blue" )