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

###################################################
### code chunk number 1: APC.rnw:2-10
###################################################
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: APC.rnw:24-40
###################################################
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 <- transform(lung, A = A + 2.5,
                        P = P + 2.5)
str(lung)
head(lung)
m.AP <- glm(cbind(D, Y / 1000) ~ factor(A) + factor(P),
            family = poisreg,
              data = lung)
m.AC <- glm(cbind(D, Y / 1000) ~ factor(A) + factor(P-A),
            family = poisreg,
              data = lung )
m.Ad <- glm(cbind(D, Y / 1000) ~ factor(A) + P,
            family = poisreg,
              data = lung )


###################################################
### code chunk number 3: APC.rnw:55-61
###################################################
m.APC <- glm(cbind(D, Y / 1000) ~ factor(A) + factor(P) + factor(P-A),
             family = poisreg,
               data = lung )
m.A   <- glm(cbind(D, Y / 1000) ~ factor(A),
             family = poisreg,
               data = lung )


###################################################
### code chunk number 4: APC.rnw:66-67
###################################################
anova(m.A, m.Ad, m.AP, m.APC, m.AC, m.Ad, test = "Chisq")


###################################################
### code chunk number 5: APC.rnw:84-86
###################################################
summary(m.APC)
ci.exp(m.APC)


###################################################
### code chunk number 6: APC.rnw:99-101
###################################################
lung$Pr <- Relevel(factor(lung$P), list("first & last"=c("1945.5", "1995.5")))
lung$Cr <- Relevel(factor(lung$P - lung$A), "1908")


###################################################
### code chunk number 7: APC.rnw:105-107
###################################################
with(lung, table(P  , Pr))
with(lung, table(P-A, Cr))


###################################################
### code chunk number 8: APC.rnw:111-115
###################################################
m.APC1 <- glm(cbind(D, Y / 1000) ~ -1 + factor(A) + factor(Pr) + factor(Cr),
              family = poisreg,
                data = lung)
coef( m.APC1 )


###################################################
### code chunk number 9: APC.rnw:123-132
###################################################
A.eff <- ci.exp(m.APC1, subset = "A")
P.eff <- rbind(c(1,1,1),
               ci.exp( m.APC1, subset="P" ),
               c(1,1,1))
(C.ref <- match("1908", levels(with(lung, factor(P - A)))))
(C.nlv <- nlevels(with(lung,factor(P-A))))
C.eff <- rbind(ci.exp(m.APC1, subset="C")[1:10,],
               c(1,1,1),
               ci.exp(m.APC1, subset="C")[11:19,] )


###################################################
### code chunk number 10: APC.rnw:135-138
###################################################
A.pt <- sort(unique(lung$A) )
P.pt <- sort(unique(lung$P) )
C.pt <- sort(unique(lung$P - lung$A ))


###################################################
### code chunk number 11: parm1
###################################################
par(mfrow = c(1, 3), las = 1)
matshade(A.pt, A.eff, plot = TRUE,
         xlab = "Age", ylab = "Lung cancer rate per 1000 PY", log="y" )
matshade(P.pt, P.eff, plot = TRUE,
         xlab = "Period", ylab = "RR", log="y")
abline(h = 1, lty = 3)
matshade(C.pt, C.eff, plot=TRUE,
         xlab = "Cohort", ylab = "RR", log = "y")
abline(h = 1, v = 1908, lty = 3)


###################################################
### code chunk number 12: parm2
###################################################
par(mfrow = c(1, 3), las = 1)
matshade(A.pt, A.eff, plot=TRUE,
         xlab="Age", ylab="Rates", ylim=c(0.1,4), log="y"  )
matshade(P.pt, P.eff, plot=TRUE,
         xlab="Period", ylab="RR", ylim=c(0.1,4)/2, log="y"  )
abline(h = 1)
matshade(C.pt, C.eff, plot=TRUE,
         xlab="Cohort", ylab="RR", ylim=c(0.1,4)/2, log="y"  )
abline( h=1 )


###################################################
### code chunk number 13: APC.rnw:204-206
###################################################
with(lung, table(P - A))
with(lung, tapply(D, list(P - A), sum))


###################################################
### code chunk number 14: APC.rnw:210-212
###################################################
(C.ref.pos <- with(lung, match(c("1878", "1933"), levels(factor(P - A)))))
(P.ref.pos <- with(lung, match("1975.5", levels(factor(P)))))


###################################################
### code chunk number 15: APC.rnw:214-216
###################################################
lung$Cx <- Relevel(factor(lung$P - lung$A), list("refCoh"=c("1878", "1933")))
lung$Px <- Relevel(factor(lung$P), "1975.5" )


###################################################
### code chunk number 16: APC.rnw:220-223
###################################################
m.APC2 <- glm(cbind(D, Y / 1000) ~ -1 + factor(A) + factor(Px) + factor(Cx),
              family = poisreg,
                data = lung)


###################################################
### code chunk number 17: APC.rnw:227-230
###################################################
c(deviance(m.APC ),
  deviance(m.APC1),
  deviance(m.APC2))


###################################################
### code chunk number 18: APC.rnw:238-251
###################################################
A.Eff <- ci.exp(m.APC2, subset = "A")
P.Eff <- ci.exp(m.APC2, subset = "P")
nP <- nrow(P.Eff)
P.Eff <- rbind(P.Eff[1:(P.ref.pos-1),],
               c(1,1,1),
               P.Eff[P.ref.pos:nP,])
C.Eff <- ci.exp(m.APC2, subset = "C")
nC <- nrow(C.Eff)
C.Eff <- rbind(C.Eff[1:(C.ref.pos[1]-1),],
               c(1,1,1),
               C.Eff[(C.ref.pos[1]):(C.ref.pos[2]-2),],
               c(1,1,1),
               C.Eff[(C.ref.pos[2]-1):nC,] )


###################################################
### code chunk number 19: parm3
###################################################
par(mfrow = c(1, 3), las = 2, mar = c(4, 3, 0.5, 0.5), mgp = c(3, 1, 0) / 1.6)
matshade(A.pt, cbind(A.eff, A.Eff), plot = TRUE,
         xlab = "Age", ylab = "Rates", ylim = c(0.1, 4),
         log = "y", col = c("black", "blue"))
matshade(P.pt, cbind(P.eff, P.Eff), plot = TRUE,
         xlab = "Period", ylab = "RR", ylim = c(0.1, 4)/2,
         log = "y", col = c("black", "blue"))
abline(h = 1)
points(c(1943, 1993, 1973)+2.5, rep(1, 3), pch = 16, col = c("black", "blue")[c(1, 1, 2)])
matshade(C.pt, cbind(C.eff, C.Eff), plot = TRUE,
          xlab = "Cohort", ylab = "RR", ylim = c(0.1, 4)/2,
          log = "y", col = c("black", "blue"))
points(c(1878, 1933, 1908), rep(1, 3), pch = 16, col = c("black", "blue")[c(2, 2, 1)])
abline(h = 1)


###################################################
### code chunk number 20: APC.rnw:306-309
###################################################
f.cp <- apc.fit(lung, model = "factor", parm = "ACP", ref.c = 1908, scale = 1000)
f.pc <- apc.fit(lung, model = "factor", parm = "APC", ref.p = 1968, scale = 1000)
names(f.pc)


###################################################
### code chunk number 21: APC.rnw:313-315
###################################################
f.cp$Drift
f.pc$Drift


###################################################
### code chunk number 22: APC.rnw:321-326
###################################################
(drifts <- rbind(
apc.fit(lung, model = "factor", dr = "d", pr = FALSE)$Drift,
apc.fit(lung, model = "factor", dr = "r", pr = FALSE)$Drift,
apc.fit(lung, model = "factor", dr = "y", pr = FALSE)$Drift,
apc.fit(lung, model = "factor", dr = "n", pr = FALSE)$Drift)[c(2, 1, 3, 5, 7), ])


###################################################
### code chunk number 23: pc-cp
###################################################
par(mar = c(3, 4, 0, 4), las = 1)
plot(f.cp, lwd = 1, r.txt = "Male lungcancer incidence in Denmark, per 1000 PY",
     col = "transparent")
pc.points(1908, 1, lwd = 2)
   matshade(f.cp$Age[, 1], f.cp$Age[, -1], lwd = 2)
pc.matshade(f.cp$Per[, 1], f.cp$Per[, -1], lwd = 2)
pc.matshade(f.cp$Coh[, 1], f.cp$Coh[, -1], lty = "21", lwd = 2)
pc.points(1965.5, 1, col = "blue", lwd = 2)
   matshade(f.pc$Age[, 1], f.pc$Age[, -1], col = "blue")
pc.matshade(f.pc$Per[, 1], f.pc$Per[, -1], col = "blue")
pc.matshade(f.pc$Coh[, 1], f.pc$Coh[, -1], col = "blue", lty = "21", lwd = 2)


###################################################
### code chunk number 24: pc-cp-sp
###################################################
s.cp <- apc.fit(lung, parm = "ACP", ref.c = 1908, scale = 1000)
par(mar = c(3, 4, 0, 4), las = 1)
plot(f.cp, lwd = 1,
     r.txt = "Male lungcancer incidence in Denmark, per 1000 PY")
   matshade(f.cp$Age[, 1], f.cp$Age[, -1], col = "red")
pc.matshade(f.cp$Per[, 1], f.cp$Per[, -1], col = "red")
pc.matshade(f.cp$Coh[, 1], f.cp$Coh[, -1], col = "red", lty = "21")
   matshade(f.pc$Age[, 1], f.pc$Age[, -1], col = "blue")
pc.matshade(f.pc$Per[, 1], f.pc$Per[, -1], col = "blue")
pc.matshade(f.pc$Coh[, 1], f.pc$Coh[, -1], col = "blue", lty = "21")
   matshade(s.cp$Age[, 1], s.cp$Age[, -1], col = "forestgreen", lwd =2)
pc.matshade(s.cp$Per[, 1], s.cp$Per[, -1], col = "forestgreen", lwd =2)
pc.matshade(s.cp$Coh[, 1], s.cp$Coh[, -1], col = "forestgreen", lwd =2, lty = "21")


###################################################
### code chunk number 25: APC.rnw:394-401
###################################################
Dr <- cbind(drifts, rbind(
apc.fit(lung, dr = "d", parm = "APC", pr = FALSE)$Drift,
apc.fit(lung, dr = "r", parm = "APC", pr = FALSE)$Drift,
apc.fit(lung, dr = "y", parm = "APC", pr = FALSE)$Drift,
apc.fit(lung, dr = "n", parm = "APC", pr = FALSE)$Drift)[c(2, 1, 3, 5, 7), ])
colnames(Dr)[c(1, 4)] <- c("Factor", "Spline")
round((Dr-1)*100, 2)