###################################################
### chunk number 1: 
###################################################
library( Epi )
library( splines )


###################################################
### chunk number 2: 
###################################################
lung <- read.table( "../data/apc-Lung.txt", header=T )
head( lung)


###################################################
### chunk number 3: 
###################################################
lung <- transform( lung, up = P-A-C, At = A, Pt = P, Ct = C )
lung <- transform( lung, A = At + 1/3 + up/3, 
                         P = Pt + 2/3 - up/3 )
lung <- transform( lung, C = P - A )
head( lung )


###################################################
### chunk number 4: 
###################################################
lrate <- with( subset( lung, A>40 & A<90 ),
               tapply( D, list(sex,
                               floor(A/5)*5+2.5,
                               floor((P-1943)/5)*5+1943+2.5), 
                          sum ) /
               tapply( Y, list(sex,
                               floor(A/5)*5+2.5,
                               floor((P-1943)/5)*5+1943+2.5), 
                          sum ) * 10^5 )


###################################################
### chunk number 5: rates
###################################################
par( mfrow=c(2,4), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
rateplot( lrate[1,,], col="blue", ylim=range(lrate,na.rm=T) )
rateplot( lrate[2,,], col="red" , ylim=range(lrate,na.rm=T) )


###################################################
### chunk number 6: 
###################################################
apc.m <- apc.fit( subset(lung,sex==1 & A>40), npar=c(8,8,15), ref.c=1930, scale=10^5 )
apc.f <- apc.fit( subset(lung,sex==2 & A>40), npar=c(8,8,15), ref.c=1930, scale=10^5 )


###################################################
### chunk number 7: apc-m
###################################################
apc.plot( apc.m, col="blue" )


###################################################
### chunk number 8: apc-f
###################################################
apc.plot( apc.f, col="red" )


###################################################
### chunk number 9: apc-2
###################################################
par( las=1, mar=c(4,4,1,4), mgp=c(3,1,0)/1.6 )
           r.lab <- c(6,c(1,2,5)*10,c(1,2,5)*100)
          rr.ref <- 200
           r.tic <- c(5:9,1:9*10,1:6*100)
apc.frame( a.lab = seq(40,90,20),
          cp.lab = seq(1880,2000,20),
           r.lab = r.lab,
          rr.ref = rr.ref,
          rr.lab = r.lab / rr.ref,
           a.tic = seq(35,90,5),
          cp.tic = seq(1855,2005,5),
           r.tic = r.tic,
          rr.tic = r.tic / rr.ref,
         tic.fac = 1.3,
           a.txt = "Age",
          cp.txt = "Calendar time",
           r.txt = "Lung cancer rate per 100,000 person-years",
          rr.txt = "Rate ratio",
        ref.line = TRUE,
             gap = 13,
        col.grid = gray(0.85),
           sides = c(1,2,4) )
apc.lines( apc.m, col="blue", ci=T )
apc.lines( apc.f, col="red" , ci=T )


###################################################
### chunk number 10: 
###################################################
rr <- 
function( one, two )
{
one[,-1] <- log(one[,-1] )
two[,-1] <- log(two[,-1] )
sd.dif <- sqrt( ((one[,4]-one[,3])/3.92)^2 +
                ((two[,4]-two[,3])/3.92)^2 )
rat <- one
rat[,-1] <- exp( cbind( one[,2]-two[,2], sd.dif ) %*% 
                 rbind( c(1,1,1), 1.96*c(0,-1,1) ) )
rat               
} 
rr.Age <- rr( apc.m$Age, apc.f$Age )
rr.Per <- rr( apc.m$Per, apc.f$Per )
rr.Coh <- rr( apc.m$Coh, apc.f$Coh )


###################################################
### chunk number 11: rr0
###################################################
( RRr <- range( rbind(rr.Age[,-1],
                      rr.Per[,-1],
                      rr.Coh[,-1]) ) )


###################################################
### chunk number 12: 
###################################################
apc.mf <- apc.m
apc.mf$Age <- rr.Age
apc.mf$Per <- rr.Per
apc.mf$Coh <- rr.Coh


###################################################
### chunk number 13: apc-rr
###################################################
par( las=1, mar=c(4,3,1,2), mgp=c(3,1,0)/1.6 )
fp <- apc.frame( a.lab = seq(40,90,20),
                cp.lab = seq(1880,2000,20),
                 r.lab = c(0.2,0.5,1,2,5),
                rr.ref = 1,
                 a.tic = seq(35,90,5),
                cp.tic = seq(1855,2005,5),
                 r.tic = c(2:9/10,1:5),
               tic.fac = 1.3,
                 a.txt = "Age",
                cp.txt = "Calendar time",
                 r.txt = "M/F Rate ratio of lung cancer",
                rr.txt = "",
              ref.line = TRUE,
                   gap = 13,
              col.grid = gray(0.85),
                 sides = c(1,2,4) )
abline( h=1 )
apc.lines( apc.mf, col="black", ci=T, frame.par=fp )


###################################################
### chunk number 14: 
###################################################
A.kn <- apc.m$Knots$Age
nk.A <- length(A.kn)
MA <- ns( lung$A, knots=A.kn[-c(1,nk.A)], Bo=A.kn[c(1,nk.A)], intercept=TRUE )
P.kn <- apc.m$Knots$Per
nk.P <- length(P.kn)
MP <- ns( lung$P, knots=P.kn[-c(1,nk.P)], Bo=P.kn[c(1,nk.P)], intercept=TRUE )
MP <- detrend( MP, lung$P )                  
C.kn <- apc.m$Knots$Coh
nk.C <- length(C.kn)
MC <- ns( lung$C, knots=C.kn[-c(1,nk.C)], Bo=C.kn[c(1,nk.C)], intercept=TRUE )
MC <- detrend( MC, lung$C )


###################################################
### chunk number 15: 
###################################################
lung$sex <- factor(lung$sex,labels=c("M","F"))
m.int <- glm( D ~ -1 + MA:sex + MP:sex + MC:sex + I(C-1930):sex +
                  offset( log(Y) ), family=poisson, data=lung )


###################################################
### chunk number 16: 
###################################################
m.per <- update( m.int, . ~ . - MP:sex + MP )
m.coh <- update( m.int, . ~ . - MC:sex + MC )
anova( m.coh, m.int, m.per, test="Chisq" )                  


###################################################
### chunk number 17: 
###################################################
m.age <- update( m.int, . ~ . - MA:sex + MA + sex + sex:A )
m.aln <- update( m.age, . ~ . - sex:A )
anova( m.int, m.age, m.aln, test="Chisq" )                  


###################################################
### chunk number 18: 
###################################################
m.RR <- glm( D ~ -1 + MA     + MP + cbind(MC,C-1930) +
                      MA:sex +      cbind(MC,C-1930):sex,
                 offset = log(Y), family=poisson, data=lung )
pr.RR <- predict( m.RR, type="terms", se.fit=TRUE )
str( pr.RR )
dimnames( pr.RR$fit )[[2]]


###################################################
### chunk number 19: 
###################################################
# Unique ages and cohort
au <- match( sort(unique(lung$A)), lung$A)  
cu <- match( sort(unique(lung$C)), lung$C)


###################################################
### chunk number 20: 
###################################################
ci.lin( m.RR )[,1:2]


###################################################
### chunk number 21: 
###################################################
au <- match( sort(unique(lung$A)), lung$A[lung$sex=="M"])
cu <- match( sort(unique(lung$C)), lung$C[lung$sex=="F"])


###################################################
### chunk number 22: 
###################################################
A.term <- exp( cbind(pr.RR$fit   [lung$sex=="M","MA:sex"][au],
                     pr.RR$se.fit[lung$sex=="M","MA:sex"][au]) %*% ci.mat() )
C.term <- exp(-cbind(pr.RR$fit   [lung$sex=="F","cbind(MC, C - 1930):sex"][cu],
                     pr.RR$se.fit[lung$sex=="F","cbind(MC, C - 1930):sex"][cu]) %*% ci.mat() )


###################################################
### chunk number 23: 
###################################################
# Unique ages and cohort
au <- match( sort(unique(lung$A)), lung$A)  
cu <- match( sort(unique(lung$C)), lung$C)
# Corresponding subsets of the design matrices
A.ctr <- MA[au,]
C.ctr <- cbind( MC[cu,], (lung$C-1930)[cu] )
# Paramter names
parnam <- names( coef(m.RR) )
# Have we found the age-paramters we want?
a.par <- intersect( grep("MA",parnam), grep("sexM",parnam) )
parnam[a.par]
# Have we found the acohort-paramters we want?
c.par <- c( grep("MC",parnam), grep("I",parnam) )
c.par <- intersect( c.par, grep("sex",parnam) )
parnam[c.par]
# Then we can extract effects, the parametrization for the cohort 
# effect is for F/M, hence we use -C.ctr
A.eff <- ci.lin( m.RR, subset=a.par, ctr.mat= A.ctr, Exp=TRUE )[,5:7]
C.eff <- ci.lin( m.RR, subset=c.par, ctr.mat=-C.ctr, Exp=TRUE )[,5:7]


###################################################
### chunk number 24: AC-rr0
###################################################
par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
matplot( lung$A[au], A.eff, 
         log="y", ylim=c(0.5,5),
         type="l", lty=1, col="black", lwd=c(3,1,1) )
matlines( lung$A[au], A.term, lty=2, col="red", lwd=c(3,1,1) )
abline(h=1)
matplot( lung$C[cu], C.eff, 
         log="y", ylim=c(0.5,5),
         type="l", lty=1, col="black", lwd=c(3,1,1) )
matlines( lung$C[cu], C.term, lty=2, col="red", lwd=c(3,1,1) )
abline(h=1)


###################################################
### chunk number 25: AC-rr1
###################################################
par( las=1, mar=c(4,3,1,2), mgp=c(3,1,0)/1.6 )
fp <- apc.frame( a.lab = seq(40,90,20),
                cp.lab = seq(1880,2000,20),
                 r.lab = c(0.5,1,2,5),
                rr.ref = 1,
                 a.tic = seq(35,90,5),
                cp.tic = seq(1855,2005,5),
                 r.tic = c(4:9/10,1:6),
               tic.fac = 1.3,
                 a.txt = "Age",
                cp.txt = "Calendar time",
                 r.txt = "M/F Rate ratio of lung cancer",
                rr.txt = "",
              ref.line = TRUE,
                   gap = 13,
              col.grid = gray(0.85),
                 sides = c(1,2,4) )
abline( h=1 )
apc.lines( apc.mf, col="black", ci=F, frame.par=fp, lwd=2 )
matlines( lung$A[au], A.eff, lwd=c(1,1,1), lty=1, col="blue" )
matlines( lung$C[cu]-fp[1], C.eff, lwd=c(1,1,1), lty=1, col="blue" )


###################################################
### chunk number 26: 
###################################################
maleA <- ns( lung$A, knots=A.kn[-c(1,nk.A)], Bo=A.kn[c(1,nk.A)], intercept=TRUE ) *
         (lung$sex=="M")
maleC <- ( ns(               lung$C, knots=C.kn[-c(1,nk.C)], Bo=C.kn[c(1,nk.C)] ) -
           ns( rep(1930,nrow(lung)), knots=C.kn[-c(1,nk.C)], Bo=C.kn[c(1,nk.C)] ) ) *
         (lung$sex=="M")


###################################################
### chunk number 27: 
###################################################
A.pt <- 40:90
C.pt <- 1860:1960
ctr.A <- ns(                 A.pt  , knots=A.kn[-c(1,nk.A)], Bo=A.kn[c(1,nk.A)], 
                                     intercept=TRUE ) 
ctr.C <- ns(                 C.pt  , knots=C.kn[-c(1,nk.C)], Bo=C.kn[c(1,nk.C)] ) -
         ns( rep(1930,length(C.pt)), knots=C.kn[-c(1,nk.C)], Bo=C.kn[c(1,nk.C)] ) 


###################################################
### chunk number 28: 
###################################################
M.RR <- glm( D ~ -1 + MA     + MP + cbind(MC,C-1930) +
                      maleA + maleC,
                 offset = log(Y), family=poisson, data=lung )
A.eff <- ci.lin( M.RR, subset="maleA", ctr.mat=ctr.A, E=T )[,5:7]
C.eff <- ci.lin( M.RR, subset="maleC", ctr.mat=ctr.C, E=T )[,5:7]


###################################################
### chunk number 29: AC-rr2
###################################################
par( las=1, mar=c(4,3,1,2), mgp=c(3,1,0)/1.6 )
fp <- apc.frame( a.lab = seq(40,90,20),
                cp.lab = seq(1880,2000,20),
                 r.lab = c(0.5,1,2,5),
                rr.ref = 1,
                 a.tic = seq(35,90,5),
                cp.tic = seq(1855,2005,5),
                 r.tic = c(4:9/10,1:6),
               tic.fac = 1.3,
                 a.txt = "Age",
                cp.txt = "Calendar time",
                 r.txt = "M/F Rate ratio of lung cancer",
                rr.txt = "",
              ref.line = TRUE,
                   gap = 13,
              col.grid = gray(0.85),
                 sides = c(1,2,4) )
abline( h=1 )
apc.lines( apc.mf, col="black", ci=TRUE, frame.par=fp, lwd=c(2,1,1) )
matlines( A.pt, A.eff, lwd=c(3,1,1), lty=1, col="blue" )
matlines( C.pt-fp[1], C.eff, lwd=c(3,1,1), lty=1, col="blue" )