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

###################################################
### code chunk number 1: age-per-coh-sol.rnw:2-3
###################################################
library(Epi)


###################################################
### code chunk number 2: age-per-coh-sol.rnw:14-22
###################################################
lung <- read.table( "../data/lung5-M.txt", header=T )
str( lung )
m.AP <- glm( D ~ factor(A) + factor(P) + offset( log(Y) ),
             family=poisson, data=lung )
m.AC <- glm( D ~ factor(A) + factor(P-A) + offset( log(Y) ),
             family=poisson, data=lung )
m.Ad <- glm( D ~ factor(A) + P + offset( log(Y) ),
             family=poisson, data=lung )


###################################################
### code chunk number 3: age-per-coh-sol.rnw:37-41
###################################################
m.APC <- glm( D ~ factor(A) + factor(P) + factor(P-A) + offset( log(Y) ),
              family=poisson, data=lung )
m.A   <- glm( D ~ factor(A) + offset( log(Y) ),
              family=poisson, data=lung )


###################################################
### code chunk number 4: age-per-coh-sol.rnw:50-51
###################################################
anova( m.A, m.Ad, m.AP, m.APC, m.AC, m.Ad, test="Chisq" )


###################################################
### code chunk number 5: age-per-coh-sol.rnw:87-89
###################################################
lung$Pr <- Relevel( factor(lung$P), list("first & last"=c("1943","1993") ) )
lung$Cr <- Relevel( factor(lung$P-lung$A), "1908" )


###################################################
### code chunk number 6: age-per-coh-sol.rnw:93-95
###################################################
with( lung, table(P,Pr) )
with( lung, table(P-A,Cr) )


###################################################
### code chunk number 7: age-per-coh-sol.rnw:102-105
###################################################
m.APC1 <- glm( D ~ -1 + factor(A) + factor(Pr) + factor(Cr) + offset( log(Y) ),
               family=poisson, data=lung )
m.APC1$coef


###################################################
### code chunk number 8: age-per-coh-sol.rnw:119-127
###################################################
A.eff <- ci.lin( m.APC1, subset="A", Exp=TRUE )[,5:7]
P.eff <- rbind( c(1,1,1),
         ci.lin( m.APC1, subset="P", Exp=TRUE )[,5:7],
                c(1,1,1) )
C.ref <- match( "1908", levels( with(lung,factor(P-A)) ) )
C.eff <- rbind( c(1,1,1),
         ci.lin( m.APC1, subset="C",
                 Exp=TRUE )[,5:7] )[c(2:C.ref,1,C.ref:(nlevels(lung$Cr)-1)),]


###################################################
### code chunk number 9: age-per-coh-sol.rnw:130-133
###################################################
A.pt <- sort( unique( lung$A ) ) + 2.5
P.pt <- sort( unique( lung$P ) ) + 2.5
C.pt <- sort( unique( lung$P-lung$A ) )


###################################################
### code chunk number 10: parm1
###################################################
par( mfrow=c(1,3), las=2 )
matplot( A.pt, A.eff,
         xlab="Age", ylab="Rates",
         log="y", type="l", lty=1, lwd=c(3,1,1), col="black" )
matplot( P.pt, P.eff,
         xlab="Period", ylab="RR",
         log="y", type="l", lty=1, lwd=c(3,1,1), col="black" )
abline( h=1 )
matplot( C.pt, C.eff,
         xlab="Cohort", ylab="RR",
         log="y", type="l", lty=1, lwd=c(3,1,1), col="black" )
abline( h=1 )


###################################################
### code chunk number 11: parm2
###################################################
par( mfrow=c(1,3), las=2 )
matplot( A.pt, A.eff*1000,
         xlab="Age", ylab="Rates", ylim=c(0.05,5),
         log="y", type="l", lty=1, lwd=c(3,1,1), col="black" )
matplot( P.pt, P.eff,
         xlab="Period", ylab="RR", ylim=c(0.05,5),
         log="y", type="l", lty=1, lwd=c(3,1,1), col="black" )
abline( h=1 )
matplot( C.pt, C.eff,
         xlab="Cohort", ylab="RR", ylim=c(0.05,5),
         log="y", type="l", lty=1, lwd=c(3,1,1), col="black" )
abline( h=1 )


###################################################
### code chunk number 12: age-per-coh-sol.rnw:205-207
###################################################
with( lung, table(P-A) )
with( lung, tapply(D,list(P-A),sum) )


###################################################
### code chunk number 13: age-per-coh-sol.rnw:211-213
###################################################
C.ref.pos <- with( lung, match( c("1878","1933"), levels( factor(P-A) ) ) )
P.ref.pos <- with( lung, match( "1973", levels( factor(P) ) ) )


###################################################
### code chunk number 14: age-per-coh-sol.rnw:215-217
###################################################
lung$Cx <- Relevel( factor(lung$P-lung$A), list("first-last"=c("1878","1933") ) )
lung$Px <- Relevel( factor(lung$P), "1973" )


###################################################
### code chunk number 15: age-per-coh-sol.rnw:221-224
###################################################
m.APC2 <- glm( D ~ -1 + factor(A) + factor(Px) + factor(Cx) + offset( log(Y) ),
               family=poisson, data=lung )
m.APC2$coef


###################################################
### code chunk number 16: age-per-coh-sol.rnw:228-231
###################################################
summary( m.APC  )$deviance
summary( m.APC1 )$deviance
summary( m.APC2 )$deviance


###################################################
### code chunk number 17: age-per-coh-sol.rnw:240-251
###################################################
A.Eff <- ci.lin( m.APC2, subset="A", Exp=TRUE )[,5:7]
P.Eff <- ci.lin( m.APC2, subset="P", Exp=TRUE )[,5:7]
nP <- nrow(P.Eff)
P.Eff <- rbind( P.Eff[1:(P.ref.pos-1),],c(1,1,1),P.Eff[P.ref.pos:nP,])
C.Eff <- ci.lin( m.APC2, subset="C",Exp=TRUE )[,5:7]
nC <- nrow(C.Eff)
C.Eff <- rbind(C.Eff[1:(C.ref.pos[1]-1),],
               c(1,1,1),
               C.Eff[(C.ref.pos[1]):(C.ref.pos[2]-2),],
               c(1,1,1),
               C.Eff[(C.ref.pos[2]-1):nC,] )


###################################################
### code chunk number 18: parm3
###################################################
par( mfrow=c(1,3), las=2 )
matplot( A.pt, cbind(A.eff,A.Eff)*1000,
         xlab="Age", ylab="Rates", ylim=c(0.05,5),
         log="y", type="l", lty=1, lwd=c(3,1,1), col=rep(c("black","blue"),each=3) )
matplot( P.pt, cbind(P.eff,P.Eff),
         xlab="Period", ylab="RR", ylim=c(0.05,5),
         log="y", type="l", lty=1, lwd=c(3,1,1), col=rep(c("black","blue"),each=3) )
abline( h=1 )
matplot( C.pt, cbind(C.eff,C.Eff),
         xlab="Cohort", ylab="RR", ylim=c(0.05,5),
         log="y", type="l", lty=1, lwd=c(3,1,1), col=rep(c("black","blue"),each=3) )
abline( h=1 )