### R code from vignette source 'brcapr.rnw'

###################################################
### code chunk number 1: brcapr.rnw:16-21
###################################################
options(width=90,
        prompt=" ", continue=" ", # Absence of prompts makes it easier for
                                  # students to copy from the final pdf document
        SweaveHooks=list( fig=function()
        par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n", lend = "butt") ))


###################################################
### code chunk number 2: brcapr.rnw:38-42
###################################################
library( Epi )
breast <- read.table("../data/breast.txt", header=T )
str( breast )
summary( breast )


###################################################
### code chunk number 3: brcapr.rnw:55-61
###################################################
breast <- transform(breast, up = P - A - C )
breast <- transform(breast, A = A + (1 + up) / 3,
                            P = P + (2 - up) / 3,
                            C = C + (1 + up) / 3 )
with(breast, summary(P - A - C))
head(breast)


###################################################
### code chunk number 4: ratetab
###################################################
ti <- 4
rt <- with(subset(breast, A>30 ),
           tapply(D, list(floor( A      /ti)*ti+ti/2,
                          floor((P-1943)/ti)*ti+ti/2+1943), sum ) /
           tapply(Y, list(floor( A      /ti)*ti+ti/2,
                          floor((P-1943)/ti)*ti+ti/2+1943), sum ) * 10^5 )
par( mfrow=c(2,2), mar=c(3,3,0,0), oma=c(0,0,1,1), mgp=c(3,1,0)/1.6 )
rateplot( rt, which= c( "ap", "ac", "pa", "ca"),
          col=heat.colors(22), ann=TRUE )


###################################################
### code chunk number 5: apcfit-1
###################################################
par( mfrow=c(1,1), mar=c(3,4,1,3) )
m1 <- apc.fit(subset(breast, A > 30),
              npar = c(8, 6, 10),
             scale = 10^5,
             ref.c = 1920)
names(m1)
m1$Knots
plot(m1)


###################################################
### code chunk number 6: apcfit-2
###################################################
par( las=1, mar=c(3,4,1,4), mgp=c(3,1,0)/1.5 )
fp <-  apc.frame( a.lab = seq(30,90,10),
                 cp.lab = seq(1860,2005,20),
                  r.lab = c(c(1,2,5)*10,c(1,2,5)*100),
#                 rr.lab = r.lab / rr.ref,
                 rr.ref = 100,
                  a.tic = seq(30,90,5),
                 cp.tic = seq(1855,2005,5),
                  r.tic = c(9,1:9*10,1:5*100),
#                 rr.tic = r.tic / rr.ref,
                tic.fac = 1.3,
                  a.txt = "Age",
                 cp.txt = "Calendar time",
                  r.txt = "Rate per 100,000 person-years",
                 rr.txt = "Rate ratio",
                    gap = 8,
               col.grid = gray(0.85),
                  sides = c(1,2,4) )
# lines( m1, ci=T, col="red" )
   matshade(m1$Age[,1], m1$Age[,-1], col="red", lwd = 2 )
pc.matshade(m1$Per[,1], m1$Per[,-1], col="red", lwd = 2 )
pc.matshade(m1$Coh[,1], m1$Coh[,-1], col="red", lwd = 2, lty = "21")
pc.points( 1920, 1, pch=16, col="red" )


###################################################
### code chunk number 7: brcapr.rnw:170-178
###################################################
# Last knot and last point on period scale
( P.rf <- c( max(m1$Knots$Per), max(m1$Per[,1]) ) )
# Last point plus one 20 years later
( P.pt <- P.rf[2] + 0:1*20 )
# Linear interpolation of log-rates at the two reference points
( Pp <- approx( m1$Per[,1], log(m1$Per[,2]), P.rf )$y )
# Liner extrapoltion through these two points to the future points
( P.eff <- Pp[2] + (Pp[2]-Pp[1])/diff(P.rf)*(P.pt-P.rf[2]) )


###################################################
### code chunk number 8: brcapr.rnw:181-185
###################################################
( C.rf <- c( max( m1$Knots$Coh ), max( m1$Coh[,1] ) ) )
( C.pt <- C.rf[2] + 0:1*20 )
( Cp <- approx( m1$Coh[,1], log(m1$Coh[,2]), C.rf )$y )
( C.eff <- Cp[2] + (Cp[2]-Cp[1])/diff(C.rf)*(C.pt-C.rf[2]) )


###################################################
### code chunk number 9: apcfit-3
###################################################
par( las=1, mar=c(3,4,1,4), mgp=c(3,1,0)/1.5 )
fp <-  apc.frame( a.lab = seq(30,90,10),
                 cp.lab = seq(1860,2020,20),
                  r.lab = c(c(1,2,5)*10,c(1,2,5)*100),
#                 rr.lab = r.lab / rr.ref,
                 rr.ref = 100,
                  a.tic = seq(30,90,5),
                 cp.tic = seq(1855,2025,5),
                  r.tic = c(9,1:9*10,1:5*100),
#                 rr.tic = r.tic / rr.ref,
                tic.fac = 1.3,
                  a.txt = "Age",
                 cp.txt = "Calendar time",
                  r.txt = "Rate per 100,000 person-years",
                 rr.txt = "Rate ratio",
                    gap = 8,
               col.grid = gray(0.85),
                  sides = c(1,2,4) )
   matshade(m1$Age[,1], m1$Age[,-1], col="red", lwd = 2 )
pc.matshade(m1$Per[,1], m1$Per[,-1], col="red", lwd = 2 )
pc.matshade(m1$Coh[,1], m1$Coh[,-1], col="red", lwd = 2, lty = "21")
pc.points( 1920, 1, pch=16, col="red" )

lines( P.pt-fp[1], exp(P.eff)*fp[2], col=gray(0.0), lty="11", lwd=2 )
lines( C.pt-fp[1], exp(C.eff)*fp[2], col=gray(0.0), lty="11", lwd=2 )


###################################################
### code chunk number 10: brcapr.rnw:249-255
###################################################
M1 <- glm(D ~ Ns(  A, knots = m1$Knots$Age) +
              Ns(P  , knots = m1$Knots$Per) +
              Ns(P-A, knots = m1$Knots$Coh)[,-1],
          offset = log(Y),
          family = poisson,
            data = subset(breast, A > 30))


###################################################
### code chunk number 11: brcapr.rnw:263-265
###################################################
c(M1$deviance, m1$Model$deviance)
summary(predict(M1) - predict(m1$Model))


###################################################
### code chunk number 12: brcapr.rnw:290-297
###################################################
a.pt <- seq(30, 90, 1/10)
Pfr <- rbind(data.frame(A = a.pt, P = 2020, Y = 1000), NA,
             data.frame(A = a.pt, P = 2030, Y = 1000) )
Cfr <- rbind(data.frame(A = a.pt, P = a.pt + 1960, Y = 1000), NA,
             data.frame(A = a.pt, P = a.pt + 1970, Y = 1000) )
prP <- ci.pred(M1, Pfr)
prC <- ci.pred(M1, Cfr)


###################################################
### code chunk number 13: pred1
###################################################
(ct <- c(0, which(is.na(prP[,1])), nrow(prP) + 1))
for( i in 1:2 )
   {
wh <- (ct[i]+1):(ct[i+1]-1)
matshade(Pfr$A[wh], cbind(prP, prC)[wh,], plot = (i == 1),
         log = "y", las = 1, xlim = c(27, 90), xlab = "",
         ylab = "Predicted breast cancer incidence per 1000 PY",
         type = "l", lwd = 2, lty = 1, col = c("red", "forestgreen"))
   }
text(rep(29.5,2), prP[c(1,603),1], paste(c(2020,2030)),
     col = "red", adj = 1, cex = 0.8 )
text(rep(29.5,2), prC[c(1,603),1], paste(c(1960,1970)),
     col = "forestgreen", adj = 1, cex = 0.8 )
mtext(side = 1, at = 50, "Age at 2020/2030",    col = "red", line = 2)
mtext(side = 1, at = 70, "Current age", col = "forestgreen", line = 2)


###################################################
### code chunk number 14: brcapr.rnw:336-346
###################################################
mp <- apc.fit(subset(breast, A > 30),
              npar = list(A = m1$Knots$Age,
                          P = m1$Knots$Per[-length(m1$Knots$Per)],
                          C = m1$Knots$Coh),
              ref.c = 1920, scale = 10^5)
mpc <- apc.fit(subset(breast, A>30),
               npar = list(A = m1$Knots$Age,
                           P = m1$Knots$Per[-length(m1$Knots$Per)],
                           C = m1$Knots$Coh[-length(m1$Knots$Coh)]),
               ref.c = 1920, scale = 10^5)


###################################################
### code chunk number 15: brcapr.rnw:351-351
###################################################



###################################################
### code chunk number 16: apcfit-4
###################################################
par( las=1, mar=c(3,4,1,4), mgp=c(3,1,0)/1.5 )
fp <-  apc.frame( a.lab = seq(30,90,10),
                 cp.lab = seq(1860,2020,20),
                  r.lab = c(c(1,2,5)*10,c(1,2,5)*100),
#                 rr.lab = r.lab / rr.ref,
                 rr.ref = 100,
                  a.tic = seq(30,90,5),
                 cp.tic = seq(1855,2025,5),
                  r.tic = c(9,1:9*10,1:5*100),
#                 rr.tic = r.tic / rr.ref,
                tic.fac = 1.3,
                  a.txt = "Age",
                 cp.txt = "Calendar time",
                  r.txt = "Rate per 100,000 person-years",
                 rr.txt = "Rate ratio",
                    gap = 8,
               col.grid = gray(0.85),
                  sides = c(1,2,4) )
lines( m1 , frame.par=fp, ci=T, col="black", lwd=c(3,1,1), knots=TRUE )
lines( mp , frame.par=fp, ci=T, col="red"      , lty=1, lwd=c(3,1,1) )
lines( mpc, frame.par=fp, ci=T, col="limegreen", lty="21", lwd=c(3,1,1) )


###################################################
### code chunk number 17: brcapr.rnw:384-389
###################################################
pr.slopes <- matrix( NA, 3, 3 )
rownames( pr.slopes ) <- c("Org","-lastP","-lastPC")
colnames( pr.slopes ) <- c("P-slope","C-slope","P-C-slope")
pr.slopes["Org","P-slope"] <- diff(Pp)/diff(P.rf)
pr.slopes["Org","C-slope"] <- diff(Cp)/diff(C.rf)


###################################################
### code chunk number 18: brcapr.rnw:393-418
###################################################
( P.rf <- c( max( mp$Knots$Per ), max( mp$Per[,1] ) ) )
P.pt <- P.rf[2] + 0:20
Pp <- approx( mp$Per[,1], log(mp$Per[,2]), P.rf )$y
P.eff <- Pp[2] + (Pp[2]-Pp[1])/diff(P.rf)*(P.pt-P.rf[2])
( C.rf <- c( max( mp$Knots$Coh ), max( mp$Coh[,1] ) ) )
C.pt <- C.rf[2] + 0:20
Cp <- approx( mp$Coh[,1], log(mp$Coh[,2]), C.rf )$y
C.eff <- Cp[2] + (Cp[2]-Cp[1])/diff(C.rf)*(C.pt-C.rf[2])
pr.slopes["-lastP","P-slope"] <- diff(Pp)/diff(P.rf)
pr.slopes["-lastP","C-slope"] <- diff(Cp)/diff(C.rf)

( P.rf <- c( max( mpc$Knots$Per ), max( mpc$Per[,1] ) ) )
P.pt <- P.rf[2] + 0:20
Pp <- approx( mpc$Per[,1], log(mpc$Per[,2]), P.rf )$y
P.eff <- Pp[2] + (Pp[2]-Pp[1])/diff(P.rf)*(P.pt-P.rf[2])
( C.rf <- c( max( mpc$Knots$Coh ), max( mpc$Coh[,1] ) ) )
C.pt <- C.rf[2] + 0:20
Cp <- approx( mpc$Coh[,1], log(mpc$Coh[,2]), C.rf )$y
C.eff <- Cp[2] + (Cp[2]-Cp[1])/diff(C.rf)*(C.pt-C.rf[2])
pr.slopes["-lastPC","P-slope"] <- diff(Pp)/diff(P.rf)
pr.slopes["-lastPC","C-slope"] <- diff(Cp)/diff(C.rf)

pr.slopes[,3] <- pr.slopes[,1] + pr.slopes[,2]
round( pr.slopes, 4 )
round( 100*(exp(pr.slopes)-1), 4 )


###################################################
### code chunk number 19: brcapr.rnw:431-443
###################################################
Mp <- glm( D ~ Ns(   A, knots=mp$Knots$Age ) +
               Ns( P  , knots=mp$Knots$Per ) +
               Ns( P-A, knots=mp$Knots$Coh )[,-1],
               family = poisson,
               offset = log(Y),
                 data =  subset( breast, A>30 ) )
Mpc <- glm( D ~ Ns(   A, knots=mpc$Knots$Age ) +
                Ns( P  , knots=mpc$Knots$Per ) +
                Ns( P-A, knots=mpc$Knots$Coh )[,-1],
                family = poisson,
                offset = log(Y),
                  data =  subset( breast, A>30 ) )


###################################################
### code chunk number 20: brcapr.rnw:449-453
###################################################
prPp <- ci.pred( Mp, Pfr )
prCp <- ci.pred( Mp, Cfr )
prPpc <- ci.pred( Mpc, Pfr )
prCpc <- ci.pred( Mpc, Cfr )


###################################################
### code chunk number 21: predx
###################################################
par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
matplot( Pfr$A, cbind( prP[,1], prPp[,1], prPpc[,1] ),
         log="y", las=1, xlim=c(25,90), xlab="Age", ylim=c(0.1,6),
         ylab="Predicted breast cancer incidence per 100,000 PY",
         type="l", lwd=3, lty=1, col=c("gray","limegreen","red") )
matplot( Pfr$A, cbind( prC[,1], prCp[,1], prCpc[,1] ),
         log="y", las=1, xlim=c(25,90), xlab="Age", ylim=c(0.1,6),
         ylab="Predicted breast cancer incidence per 100,000 PY",
         type="l", lwd=3, lty=1, col=c("gray","limegreen","red") )