### R code from vignette source 'lung-sex.rnw'

###################################################
### code chunk number 1: lung-sex.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: lung-sex.rnw:11-13
###################################################
library(Epi)
library(tidyverse)


###################################################
### code chunk number 3: lung-sex.rnw:39-43
###################################################
library( Epi )
library( splines )
lung <- read.table( "../data/apc-Lung.txt", header=T )
head( lung)


###################################################
### code chunk number 4: lung-sex.rnw:56-61
###################################################
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 )


###################################################
### code chunk number 5: lung-sex.rnw:77-86
###################################################
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 )


###################################################
### code chunk number 6: rates
###################################################
# x11(h=18,w=27,p=24)
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) )


###################################################
### code chunk number 7: lung-sex.rnw:112-114
###################################################
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 )


###################################################
### code chunk number 8: apc-m
###################################################
par( mfrow=c(1,1) )
apc.plot( apc.m, col="blue" )


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


###################################################
### code chunk number 10: apc-2
###################################################
    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:7*100)
par( las=1, mar=c(4,3,1,4), mgp=c(3,1,0)/1.6 )
apc.frame( a.lab = seq(40,90,20),
          cp.lab = seq(1880,2000,20),
           r.lab = c(6,c(1,2,5)*10,c(1,2,5)*100),
          rr.lab = r.lab / rr.ref,
          rr.ref = 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 = 7,
        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 )


###################################################
### code chunk number 11: lung-sex.rnw:199-203
###################################################
rr <- function(one, two) cbind(one[,1], ci.ratio(one[,-1], two[,-1]))
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 )


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


###################################################
### code chunk number 13: lung-sex.rnw:220-224
###################################################
apc.mf <- apc.m
apc.mf$Age <- rr.Age
apc.mf$Per <- rr.Per
apc.mf$Coh <- rr.Coh


###################################################
### code chunk number 14: apc-rr
###################################################
par( las=1, mar=c(4,3,1,2), mgp=c(3,1,0)/1.6 )
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 )


###################################################
### code chunk number 15: apc-3
###################################################
    r.lab <- c(6,c(1,2,5)*10,c(1,2,5)*100)
   rr.ref <- 100
    r.tic <- c(5:9,1:9*10,1:6*100)
par( las=1, mar=c(4,3,1,4), mgp=c(3,1,0)/1.6 )
apc.frame( a.lab = seq(40,90,20),
          cp.lab = seq(1880,2000,20),
           r.lab = c(6,c(1,2,5)*10,c(1,2,5)*100),
          rr.lab = r.lab / rr.ref,
          rr.ref = 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 )
   matlines( rr.Age[,1], rr.Age[,-1]*rr.ref, lwd=c(3,1,1), col="black", lty=1 )
pc.matlines( rr.Per[,1], rr.Per[,-1]       , lwd=c(3,1,1), col="black", lty=1 )
pc.matlines( rr.Coh[,1], rr.Coh[,-1]       , lwd=c(3,1,1), col="black", lty=1 )


###################################################
### code chunk number 16: lung-sex.rnw:306-315
###################################################
A.kn <- apc.m$Knots$Age ; nk.A <- length(A.kn)
P.kn <- apc.m$Knots$Per ; nk.P <- length(P.kn)
C.kn <- apc.m$Knots$Coh ; nk.C <- length(C.kn)
MA   <- Ns( lung$A, knots=A.kn, intercept=TRUE )
MP   <- Ns( lung$P, knots=P.kn, intercept=TRUE )
MP   <- detrend( MP, lung$P )
MC   <- Ns( lung$C, knots=C.kn, intercept=TRUE )
MC   <- detrend( MC, lung$C )
lung$sex <- factor(lung$sex,labels=c("M","F"))


###################################################
### code chunk number 17: lung-sex.rnw:320-326
###################################################
mc.int <- glm( D ~ -1 + MA:sex + MP:sex + MC:sex + I(C-1930):sex,
                    offset=log(Y), family=poisson, data=lung )
mp.int <- glm( D ~ -1 + MA:sex + MP:sex + MC:sex + I(P-1980):sex,
                    offset=log(Y), family=poisson, data=lung )
rbind( ci.exp( mc.int, subset="I" ),
       ci.exp( mp.int, subset="I" ) )


###################################################
### code chunk number 18: lung-sex.rnw:331-339
###################################################
mcc.int <- update( mc.int, . ~ . - MC:sex + MC )
mpc.int <- update( mp.int, . ~ . - MC:sex + MC )
rbind( ci.exp( mcc.int, subset="I" ),
       ci.exp( mpc.int, subset="I" ) )
mcp.int <- update( mc.int, . ~ . - MP:sex + MP )
mpp.int <- update( mp.int, . ~ . - MP:sex + MP )
rbind( ci.exp( mcp.int, subset="I" ),
       ci.exp( mpp.int, subset="I" ) )


###################################################
### code chunk number 19: lung-sex.rnw:356-359
###################################################
m.per <- update( mc.int, . ~ . - MP:sex + MP )
m.coh <- update( mc.int, . ~ . - MC:sex + MC )
anova( m.coh, mc.int, m.per, test="Chisq" )


###################################################
### code chunk number 20: lung-sex.rnw:371-374
###################################################
m.age <- update( mc.int, . ~ . - MA:sex + MA + sex + sex:A )
m.aln <- update( m.age, . ~ . - sex:A )
anova( mc.int, m.age, m.aln, test="Chisq" )


###################################################
### code chunk number 21: lung-sex.rnw:393-399
###################################################
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]]


###################################################
### code chunk number 22: lung-sex.rnw:406-409
###################################################
# Unique ages and cohort
au <- match( sort(unique(lung$A)), lung$A)
cu <- match( sort(unique(lung$C)), lung$C)


###################################################
### code chunk number 23: lung-sex.rnw:413-414
###################################################
ci.lin( m.RR )[,1:2]


###################################################
### code chunk number 24: lung-sex.rnw:419-421
###################################################
au <- match( sort(unique(lung$A)), lung$A[lung$sex=="M"])
cu <- match( sort(unique(lung$C)), lung$C[lung$sex=="F"])


###################################################
### code chunk number 25: lung-sex.rnw:426-430
###################################################
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() )


###################################################
### code chunk number 26: lung-sex.rnw:436-455
###################################################
# 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] )
# Parameter names
parnam <- names( coef(m.RR) )
# Have we found the age-parameters we want?
a.par <- intersect( grep("MA",parnam), grep("sexM",parnam) )
parnam[a.par]
# Have we found the cohort-parameters 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]


###################################################
### code chunk number 27: 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)


###################################################
### code chunk number 28: AC-rr1
###################################################
par( las=1, mar=c(4,3,1,2), mgp=c(3,1,0)/1.6 )
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, lwd=2 )
   matlines( lung$A[au], A.eff, lwd=c(1,1,1), lty=1, col="blue" )
pc.matlines( lung$C[cu], C.eff, lwd=c(1,1,1), lty=1, col="blue" )


###################################################
### code chunk number 29: lung-sex.rnw:533-538
###################################################
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")


###################################################
### code chunk number 30: lung-sex.rnw:541-547
###################################################
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)] )


###################################################
### code chunk number 31: lung-sex.rnw:552-557
###################################################
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]


###################################################
### code chunk number 32: AC-rr2
###################################################
par( las=1, mar=c(4,3,1,2), mgp=c(3,1,0)/1.6 )
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, lwd=c(2,1,1) )
   matlines( A.pt, A.eff, lwd=c(3,1,1), lty=1, col="blue" )
pc.matlines( C.pt, C.eff, lwd=c(3,1,1), lty=1, col="blue" )