### R code from vignette source 'AP-AC.rnw'

###################################################
### code chunk number 1: AP-AC.rnw:2-9
###################################################
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")))


###################################################
### code chunk number 2: AP-AC.rnw:22-24
###################################################
library(Epi)
library(tidyverse)


###################################################
### code chunk number 3: AP-AC.rnw:30-41
###################################################
lung <- read.table("http://bendixcarstensen.com/APC/KEA-2023/data/lung5-M.txt",
                   header = T)
lung <- read.table("../data/lung5-M.txt", header = T)
head(lung)
with(lung , table(A))
with(lung , table(P))
round(ftable(addmargins(xtabs(cbind(D = D,
                                    Y = Y/1000) ~  A + P,
                              data = lung),
                        margin = 1),
             row.vars = c(3, 1)), 0)


###################################################
### code chunk number 4: AP-AC.rnw:51-56
###################################################
ap.1 <- glm(D ~ factor(A) + factor(P),
            offset = log(Y / 1000),
            family = poisson,
              data = lung)
summary(ap.1)


###################################################
### code chunk number 5: AP-AC.rnw:61-62
###################################################
ap.1<-glm(D~factor(A)+factor(P),offset=log(Y/1000),family=poisson,data=lung)


###################################################
### code chunk number 6: AP-AC.rnw:70-74
###################################################
ap.1 <- glm(cbind(D, Y/1000) ~ factor(A) + factor(P),
            family = poisreg,
              data = lung)
summary(ap.1)


###################################################
### code chunk number 7: AP-AC.rnw:88-89
###################################################
round(ci.exp(ap.1), 2)


###################################################
### code chunk number 8: AP-AC.rnw:94-98
###################################################
ap.0 <- glm(cbind(D, Y / 1000) ~ -1 + factor(A) + factor(P),
             family = poisreg,
               data = lung)
round(ci.exp(ap.0), 3)


###################################################
### code chunk number 9: AP-AC.rnw:108-111
###################################################
ap.3 <- glm(cbind(D, Y / 1000) ~ factor(A) - 1 + relevel(factor(P), "1968"),
            family = poisreg,
            data = lung)


###################################################
### code chunk number 10: AP-AC.rnw:114-115
###################################################
round(ci.exp(ap.3), 3)


###################################################
### code chunk number 11: AP-AC.rnw:122-123
###################################################
(ap.cf <- ci.exp(ap.3, subset = "A"))


###################################################
### code chunk number 12: AP-AC.rnw:131-134
###################################################
nd <- data.frame(A = seq(40, 85, 5), P = 1968)
(ap.rt <- ci.pred(ap.3, nd))
(ap.rt <- ci.pred(ap.0, nd))


###################################################
### code chunk number 13: agesh
###################################################
matshade(seq(40, 85, 5) + 2.5, ci.exp(ap.3, subset = "A"),
         plot = TRUE,
         type = "l", lty = 1, lwd = 1, col = 1,
          log = "y",
         xlab = "Age",
         ylab = "Male lung cancer incidence rate per 1000 PY")


###################################################
### code chunk number 14: AP-AC.rnw:157-158
###################################################
(RR.cf <- ci.exp(ap.3, subset = "P"))


###################################################
### code chunk number 15: AP-AC.rnw:163-164
###################################################
(RR.cf <- rbind(RR.cf[1:5, ], 1, RR.cf[6:10, ]))


###################################################
### code chunk number 16: APrrLung
###################################################
matshade(seq(1943, 1993, 5) + 2.5, RR.cf,
          lty = 1, lwd = 1, log = "y", col = 1, plot = TRUE,
          xlab = "Calendar time", ylab = "Rate ratio rel. to 1968--72")
abline(h = 1, v = 1970.5, lty = 3)


###################################################
### code chunk number 17: AP-AC.rnw:187-191
###################################################
ap.x <- glm(cbind(D, Y / 1000) ~ -1 + A + P,
            family = poisreg,
              data = transform(lung, A = factor(A), P = factor(P)))
summary(ap.x)


###################################################
### code chunk number 18: AP-AC.rnw:198-201
###################################################
nd <- data.frame(P = seq(1943, 1993, 5))
nr <- data.frame(P = 1968)
(rrx <- ci.exp(ap.x, list(nd, nr)))


###################################################
### code chunk number 19: AP-AC.rnw:217-218
###################################################
with(lung, table(P-A))


###################################################
### code chunk number 20: AP-AC.rnw:226-230
###################################################
ac.0 <- glm(cbind(D, Y / 1000) ~ A + C,
            family = poisreg,
              data = transform(lung, A = factor(A), C = factor(P-A)))
summary(ac.0)


###################################################
### code chunk number 21: AP-AC.rnw:239-242
###################################################
table(near(residuals(ac.0), 0))
subset(lung, near(residuals(ac.0), 0))
table(residuals(ac.0) == 0)


###################################################
### code chunk number 22: AP-AC.rnw:250-256
###################################################
ac.r <- glm(cbind(D, Y / 1000) ~ -1 + A + C,
             family = poisreg,
               data = transform(lung,
                                A = factor(A),
                                C = relevel(factor(P - A), "1908")))
round(ci.exp(ac.r), 3)


###################################################
### code chunk number 23: AP-AC.rnw:282-286
###################################################
ndc <- data.frame(C = seq(1858, 1953, 5))
ndr <- data.frame(C = 1908)
try(RR.C <- ci.exp(ac.r, list(ndc, ndr)))
   (RR.C <- ci.exp(ac.0, list(ndc, ndr)))


###################################################
### code chunk number 24: cohRR
###################################################
matshade(ndc$C, RR.C, log = 'y', plot = TRUE,
          xlab = "Date of birth", ylab = "Lung cancer incidence RR")
abline(h = 1, v = 1908, lty = 3)


###################################################
### code chunk number 25: AP-AC.rnw:297-299
###################################################
ai.coh <- ci.pred(ac.0, data.frame(A = factor(seq(40, 85, 5)),
                                   C = '1908'))


###################################################
### code chunk number 26: Aincmp
###################################################
matshade(seq(40, 85, 5), ai.coh, log = "y", plot = TRUE)