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

###################################################
### code chunk number 1: APC-s.rnw:3-8
###################################################
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: APC-s.rnw:24-32
###################################################
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: APC-s.rnw:47-55
###################################################
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: APC-s.rnw:61-62
###################################################
anova( m.A, m.Ad, m.AP, m.APC, m.AC, m.Ad, test="Chisq" )


###################################################
### code chunk number 5: APC-s.rnw:83-85
###################################################
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: APC-s.rnw:89-91
###################################################
with( lung, table(P,Pr) )
with( lung, table(P-A,Cr) )


###################################################
### code chunk number 7: APC-s.rnw:96-101
###################################################
m.APC1 <- glm( D ~ -1 + factor(A) + factor(Pr) + factor(Cr),
               offset = log(Y),
               family = poisson, 
                 data = lung )
round( m.APC1$coef, 3 )


###################################################
### code chunk number 8: APC-s.rnw:111-119
###################################################
A.eff <- ci.exp( m.APC1, subset="A" )
P.eff <- rbind( c(1,1,1),
                ci.exp( m.APC1, subset="P" ),
                c(1,1,1) )
C.ref <- match( "1908", levels( with(lung,factor(P-A)) ) )
C.eff <- rbind( ci.exp( m.APC1, subset="C" )[1:(C.ref-1),],
                c(1,1,1),
                ci.exp( m.APC1, subset="C" )[C.ref:(nlevels(lung$Cr)-1),] )


###################################################
### code chunk number 9: APC-s.rnw:122-125
###################################################
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 )
matshade( A.pt, A.eff, plot=TRUE,
          xlab="Age", ylab="Rates", log="y" )
matshade( P.pt, P.eff, plot=TRUE,
          xlab="Period", ylab="RR", log="y" )
abline( h=1 )
matshade( C.pt, C.eff, plot=TRUE,
          xlab="Cohort", ylab="RR", log="y" )
abline( h=1 )


###################################################
### code chunk number 11: parm2
###################################################
par( mfrow=c(1,3), las=2 )
matshade( A.pt, A.eff*1000, plot=TRUE,
          xlab="Age", ylab="Rates", ylim=c(0.1,4), log="y"  )
matshade( P.pt, P.eff, plot=TRUE,
          xlab="Period", ylab="RR", ylim=c(0.1,4)/2, log="y"  )
abline( h=1 )
matshade( C.pt, C.eff, plot=TRUE,
          xlab="Cohort", ylab="RR", ylim=c(0.1,4)/2, log="y"  )
abline( h=1 )


###################################################
### code chunk number 12: APC-s.rnw:191-193
###################################################
with( lung, table(P-A) )
with( lung, tapply(D,list(P-A),sum) )


###################################################
### code chunk number 13: APC-s.rnw:197-199
###################################################
( 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: APC-s.rnw:201-203
###################################################
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: APC-s.rnw:207-209
###################################################
m.APC2 <- glm( D ~ -1 + factor(A) + factor(Px) + factor(Cx) + offset( log(Y) ),
               family=poisson, data=lung )


###################################################
### code chunk number 16: APC-s.rnw:213-216
###################################################
c(summary( m.APC  )$deviance, 
  summary( m.APC1 )$deviance, 
  summary( m.APC2 )$deviance )


###################################################
### code chunk number 17: APC-s.rnw:224-235
###################################################
A.Eff <- ci.exp( m.APC2, subset="A" )
P.Eff <- ci.exp( m.APC2, subset="P" )
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.exp( m.APC2, subset="C" )
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, mar=c(4,3,0.5,0.5), mgp=c(3,1,0)/1.6 )
matshade( A.pt, cbind(A.eff,A.Eff)*1000, plot=TRUE,
          xlab="Age", ylab="Rates", ylim=c(0.1,4), 
          log="y", col=c("black","blue") )
matshade( P.pt, cbind(P.eff,P.Eff), plot=TRUE,
          xlab="Period", ylab="RR", ylim=c(0.1,4)/2,
          log="y", col=c("black","blue") )
abline( h=1 )
points( c(1943,1993,1973)+2.5, rep(1,3), pch=16, col=c("black","blue")[c(1,1,2)])
matshade( C.pt, cbind(C.eff,C.Eff), plot=TRUE,
          xlab="Cohort", ylab="RR", ylim=c(0.1,4)/2,
          log="y", col=c("black","blue") )
points( c(1878,1933,1908), rep(1,3), pch=16, col=c("black","blue")[c(2,2,1)])
abline( h=1 )


###################################################
### code chunk number 19: APC-s.rnw:293-296
###################################################
f.cp <- apc.fit( lung, model = "factor", parm = "ACP", ref.c=1908, scale=1000 )
f.pc <- apc.fit( lung, model = "factor", parm = "APC", ref.p=1968, scale=1000 )
names( f.pc )


###################################################
### code chunk number 20: APC-s.rnw:300-302
###################################################
f.cp$Drift
f.pc$Drift


###################################################
### code chunk number 21: APC-s.rnw:308-313
###################################################
( drifts <- rbind(
apc.fit( lung, model="factor", dr="d", pr=FALSE )$Drift,
apc.fit( lung, model="factor", dr="r", pr=FALSE )$Drift,
apc.fit( lung, model="factor", dr="y", pr=FALSE )$Drift,
apc.fit( lung, model="factor", dr="n", pr=FALSE )$Drift)[c(2,1,3,5,7),] )


###################################################
### code chunk number 22: pc-cp
###################################################
par( mar=c(3,4,0,4), las=1 )
plot( f.cp, lwd=1, r.txt="Male lungcancer incidence in Denmark, per 1000 PY" )
   matshade( f.cp$Age[,1], f.cp$Age[,-1] )
pc.matshade( f.cp$Per[,1], f.cp$Per[,-1] )
pc.matshade( f.cp$Coh[,1], f.cp$Coh[,-1] )
   matshade( f.pc$Age[,1], f.pc$Age[,-1], col="blue" )
pc.matshade( f.pc$Per[,1], f.pc$Per[,-1], col="blue" )
pc.matshade( f.pc$Coh[,1], f.pc$Coh[,-1], col="blue" )


###################################################
### code chunk number 23: pc-cp-sp
###################################################
s.cp <- apc.fit( lung, parm = "ACP", ref.c=1908, scale=1000 )
par( mar=c(3,4,0,4), las=1 )
plot( f.cp, lwd=1, r.txt="Male lungcancer incidence in Denmark, per 1000 PY" )
   matshade( f.cp$Age[,1], f.cp$Age[,-1] )
pc.matshade( f.cp$Per[,1], f.cp$Per[,-1] )
pc.matshade( f.cp$Coh[,1], f.cp$Coh[,-1] )
   matshade( f.pc$Age[,1], f.pc$Age[,-1], col="blue" )
pc.matshade( f.pc$Per[,1], f.pc$Per[,-1], col="blue" )
pc.matshade( f.pc$Coh[,1], f.pc$Coh[,-1], col="blue" )
   matshade( s.cp$Age[,1], s.cp$Age[,-1], col="forestgreen" )
pc.matshade( s.cp$Per[,1], s.cp$Per[,-1], col="forestgreen" )
pc.matshade( s.cp$Coh[,1], s.cp$Coh[,-1], col="forestgreen" )


###################################################
### code chunk number 24: APC-s.rnw:379-386
###################################################
Dr <- cbind( drifts, rbind(
apc.fit( lung, dr="d", parm="APC", pr=FALSE )$Drift,
apc.fit( lung, dr="r", parm="APC", pr=FALSE )$Drift,
apc.fit( lung, dr="y", parm="APC", pr=FALSE )$Drift,
apc.fit( lung, dr="n", parm="APC", pr=FALSE )$Drift)[c(2,1,3,5,7),] )
colnames( Dr )[c(1,4)] <- c("Factor","Spline") 
round( (Dr-1)*100, 2 )