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

###################################################
### code chunk number 1: Ad.rnw:2-12
###################################################
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)
library(tidyverse)
clear()


###################################################
### code chunk number 2: Ad.rnw:25-30
###################################################
lung <- read.table("http://bendixcarstensen.com/APC/KEA-2023/data/lung5-M.txt",
                   header = T)
# lung <- read.table( "../data/lung5-M.txt", header=T )
lung$C <- lung$P - lung$A
table(lung$C)


###################################################
### code chunk number 3: Ad.rnw:38-42
###################################################
mp <- glm(cbind(D, Y / 1000) ~ -1 + factor(A) + I(P - 1968),
          family = poisreg,
            data = lung )
round(ci.lin(mp), 4)


###################################################
### code chunk number 4: Ad.rnw:52-59
###################################################
lung <- transform(lung, A = A + 2.5,
                        P = P + 2.5)
mp <- glm(cbind(D, Y / 1000) ~ -1 + factor(A) + I(P-1970.5),
          family = poisreg,
            data = lung )
ci.lin(mp)[,1:2]
ci.exp(mp)


###################################################
### code chunk number 5: Ad.rnw:68-73
###################################################
mc <- glm(cbind(D, Y / 1000) ~ -1 + factor(A) + I(C - 1908),
          family = poisreg,
            data = lung)
round(ci.exp(mc), 6)
round(ci.exp(mp), 6)


###################################################
### code chunk number 6: Ad.rnw:85-88
###################################################
oo <- options(digits = 9)
c(deviance(mp), deviance(mc))
options(oo)


###################################################
### code chunk number 7: Ad.rnw:106-112
###################################################
(ap <- ci.lin(mp)[1:10,1])
(ac <- ci.lin(mc)[1:10,1])
cbind(ap, ac, ap - ac, diff(ap - ac))
c.sl <- ci.lin(mc)[11,1]
a.pt <- seq(40, 85, 5) + 2.5
cbind( ap, ac + c.sl*(62.5-a.pt) )


###################################################
### code chunk number 8: rates
###################################################
matshade(a.pt, cbind(ci.exp(mp, subset = "A"),
                     ci.exp(mc, subset = "A")) * 10^5, plot = TRUE,
         log = "y", xlab = "Age", ylab = "Lung cancer incidence rates / 100,000",
         lty = 1, lwd = 1, col = c("black","blue") )


###################################################
### code chunk number 9: RRs
###################################################
p.pt <- seq(min(lung$P), max(lung$P), , 10) + 2.5
c.pt <- seq(min(lung$C), max(lung$C), , 10)
ctr.p <- cbind(p.pt - 1970.5)
ctr.c <- cbind(c.pt - 1908  )
matshade(c.pt, ci.exp(mc, subset = "C", ctr.mat = ctr.c), plot = TRUE,
         log = "y", xlab = "Calendar time", ylab = "Rate ratio", xlim = c(1850,2000),
         type = "l", lty = 1, lwd = 1, col = "blue" )
matshade(p.pt, ci.exp(mp, subset = "P", ctr.mat = ctr.p),
         type = "l", lty = 1, lwd = 1, col = "black" )
abline(h = 1, lty = 3)
points(c(1908, 1970.5), c(1, 1), pch = 16)