### R code from vignette source 'apc-tri.rnw'

###################################################
### code chunk number 1: apc-tri.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: apc-tri.rnw:28-33
###################################################
library( Epi)
ltri <- read.table( "../data/lung5-Mc.txt", header = T)
# ltri <- read.table( "http://bendixcarstensen.com/APC/KEA-2023/data/lung5-Mc.txt", header = T)
head(ltri)
with(ltri, table(P5 - A5 - C5))


###################################################
### code chunk number 3: apc-tri.rnw:37-38
###################################################
ltri$S5 <- ltri$P5 - ltri$A5


###################################################
### code chunk number 4: Lexis-ill
###################################################
par( mar = c(3,3,1,1), mgp = c(3,1,0)/1.6)
Lexis.diagram(age = 30 + c(0,65),
             date = 1938 + c(0,65),
         coh.grid = TRUE)
with(ltri, text(Px, Ax, paste(D), cex = 0.55))
box()


###################################################
### code chunk number 5: apc-tri.rnw:65-68
###################################################
ms <- glm(cbind(D, Y) ~ -1 + factor(A5) + factor(P5) + factor(S5),
           family = poisreg, data = ltri)
summary(ms)


###################################################
### code chunk number 6: apc-tri.rnw:73-77
###################################################
mc <- glm(cbind(D, Y) ~ -1 + factor(A5) + factor(P5) + factor(C5),
           family = poisreg, data = ltri)
summary(mc)$df
summary(ms)$df


###################################################
### code chunk number 7: parmest1
###################################################
par(mfrow = c(1,3))
(a.pt <- as.numeric(levels(factor(ltri$A5))))
(p.pt <- as.numeric(levels(factor(ltri$P5))))
(s.pt <- as.numeric(levels(factor(ltri$S5))))
(c.pt <- as.numeric(levels(factor(ltri$C5))))

matshade(a.pt, cbind(ci.exp(ms, subset = "A5"),
                     ci.exp(mc, subset = "A5")), plot = TRUE,
         lty = 1, lwd = 2, col = c("black", "blue"),
         xlab = "Age", ylab = "Rates", log = "y")

matshade(p.pt, rbind(1,
               cbind(ci.exp(ms, subset = "P5"),
                     ci.exp(mc, subset = "P5"))), plot = TRUE,
         lty = 1, lwd = 2, col = c("black", "blue"),
         xlab = "Period", ylab = "Rates", log = "y")

matshade(s.pt, rbind(1, ci.exp(ms, subset = "S5")), plot = TRUE,
         lty = 1, lwd = 2, col = "black",
         xlab = "cohort", ylab = "Rates", log = "y")

matshade(c.pt, rbind(1, ci.exp(mc, subset = "C5")),
         lty = 1, lwd = 2, col = "blue",
         xlab = "cohort", ylab = "Rates", log = "y")


###################################################
### code chunk number 8: apc-tri.rnw:120-137 (eval = FALSE)
###################################################
## matplot(a.pt, ci.lin(ms, subset = "A5", Exp = TRUE)[,5:7]/10^5,
##          type = "l", lty = 1, lwd = c(3,1,1), col = "black",
##          xlab = "Age", ylab = "Rates", log = "y")
## matlines(a.pt, ci.lin(mc, subset = "A5", Exp = TRUE)[,5:7]/10^5,
##           type = "l", lty = 1, lwd = c(3,1,1), col = "blue")
## 
## matplot(p.pt, rbind(c(1,1,1), ci.lin(ms, subset = "P5",Exp = TRUE)[,5:7]),
##          type = "l", lty = 1, lwd = c(3,1,1), col = "black",
##          xlab = "Period", ylab = "RR", log = "y")
## matlines(p.pt, rbind(c(1,1,1), ci.lin(mc, subset = "P5",Exp = TRUE)[,5:7]),
##           type = "l", lty = 1, lwd = c(3,1,1), col = "blue")
## 
## matplot(s.pt, rbind(c(1,1,1),ci.lin(ms, subset = "S5", Exp = TRUE)[,5:7]),
##          type = "l", lty = 1, lwd = c(3,1,1), col = "black",
##          xlab = "Cohort", ylab = "RR", log = "y")
## matlines(c.pt, rbind(c(1,1,1),ci.lin(mc, subset = "C5", Exp = TRUE)[,5:7]),
##           type = "l", lty = 1, lwd = c(3,1,1), col = "blue")


###################################################
### code chunk number 9: apc-tri.rnw:151-158
###################################################
mt <- glm(cbind(D, Y) ~ -1 + factor(Ax) + factor(Px) + factor(Cx),
           family = poisreg, data = ltri)
summary(mt)$df
sum(!is.na(coef(mt)))
length(coef(mt))
nrow(ltri)
nrow(ltri) - sum(!is.na(coef(mt)))


###################################################
### code chunk number 10: parmest2
###################################################
par(mfrow = c(1,3))
(a.pt <- as.numeric(levels(factor(ltri$Ax))))
(p.pt <- as.numeric(levels(factor(ltri$Px))))
(c.pt <- as.numeric(levels(factor(ltri$Cx))))

matshade(a.pt, ci.exp(mt, subset = "Ax") * 10^5,
         lty = 1, lwd = 2, plot = TRUE,
         xlab = "Age", ylab = "Rates per 10^5", log = "y")

matshade(p.pt, rbind(c(1,1,1), ci.exp(mt, subset = "Px")),
         lty = 1, lwd = 2, plot = TRUE,
        xlab = "Period", ylab = "RR", log = "y")

matshade(c.pt, rbind(c(1,1,1), ci.exp(mt, subset = "Cx")),
        lty = 1, lwd = 2, plot = TRUE,
        xlab = "Cohort", ylab = "RR", log = "y")


###################################################
### code chunk number 11: apc-tri.rnw:192-193
###################################################
summary(mt)$deviance


###################################################
### code chunk number 12: apc-tri.rnw:200-202
###################################################
with(ltri, table(up, P5 - A5 - C5))
with(ltri, table(up, P5 - A5 - S5))


###################################################
### code chunk number 13: apc-tri.rnw:210-218
###################################################
m.up <- glm(cbind(D, Y) ~ -1 + factor(A5) + factor(P5) + factor(S5),
             family = poisreg, data = subset(ltri,up == 1))
summary(m.up)$deviance
m.lo <- glm(cbind(D, Y) ~ -1 + factor(A5) + factor(P5) + factor(S5),
             family = poisreg, data = subset(ltri,up == 0))
summary(m.lo)$deviance
summary(m.lo)$deviance + summary(m.up)$deviance
summary(mt)$deviance


###################################################
### code chunk number 14: parmest3
###################################################
par(mfrow = c(1,3))
a.pt <- as.numeric(levels(factor(ltri$Ax)))
p.pt <- as.numeric(levels(factor(ltri$Px)))
c.pt <- as.numeric(levels(factor(ltri$Cx)))

a5.pt <- as.numeric(levels(factor(ltri$A5)))
p5.pt <- as.numeric(levels(factor(ltri$P5)))
s5.pt <- as.numeric(levels(factor(ltri$S5)))

matplot(a.pt, ci.lin(mt, subset = "Ax", Exp = TRUE)[,5:7]/10^5,
         type = "l", lty = 1, lwd = c(2,1,1), col = gray(0.7),
         xlab = "Age", ylab = "Rates", log = "y")
matpoints(a5.pt, ci.lin(m.up, subset = "A5", Exp = TRUE)[,5:7]/10^5,
           pch = c(16,3,3), col = "blue")
matpoints(a5.pt, ci.lin(m.lo, subset = "A5", Exp = TRUE)[,5:7]/10^5,
           pch = c(16,3,3), col = "red")

matplot(p.pt, rbind(c(1,1,1), ci.lin(mt, subset = "Px",Exp = TRUE)[,5:7]),
         type = "l", lty = 1, lwd = c(2,1,1), col = gray(0.7),
         xlab = "Period", ylab = "RR", log = "y")
matpoints(p5.pt[-1], ci.lin(m.up, subset = "P5", Exp = TRUE)[,5:7],
           pch = c(16,3,3), col = "blue")
matpoints(p5.pt[-1], ci.lin(m.lo, subset = "P5", Exp = TRUE)[,5:7],
           pch = c(16,3,3), col = "red")

matplot(c.pt, rbind(c(1,1,1),ci.lin(mt, subset = "Cx", Exp = TRUE)[,5:7]),
         type = "l", lty = 1, lwd = c(2,1,1), col = gray(0.7),
         xlab = "Cohort", ylab = "RR", log = "y")
matpoints(s5.pt[-1], ci.lin(m.up, subset = "S5", Exp = TRUE)[,5:7],
           pch = c(16,3,3), col = "blue")
matpoints(s5.pt[-1], ci.lin(m.lo, subset = "S5", Exp = TRUE)[,5:7],
           pch = c(16,3,3), col = "red")


###################################################
### code chunk number 15: apc-tri.rnw:286-296
###################################################
library(splines)
mspl <- glm(cbind(D, Y) ~ -1 + ns(Ax,df = 7,intercept = T)
                             + ns(Px,df = 6,intercept = F)
                             + ns(Cx,df = 6,intercept = F),
             family = poisreg, data = ltri)
summary(mspl)
summary(mt)$deviance
summary(mspl)$deviance
summary(mt)$deviance - summary(mspl)$deviance
summary(mt)$df       - summary(mspl)$df


###################################################
### code chunk number 16: parmest4
###################################################
pspl <- predict(mspl, type = "terms", se.fit = TRUE)
str(pspl)
a.ord <- order(ltri$Ax)
p.ord <- order(ltri$Px)
c.ord <- order(ltri$Cx)

par(mfrow = c(1,3))
matplot(ltri$Ax[a.ord], exp(cbind(pspl$fit[,1],
                                  pspl$se.fit[,1])[a.ord,]
                            %*% ci.mat())*10^5,
        type = "l", lty = 1, lwd = c(2,1,1), col = gray(0.2),
        xlab = "Age", ylab = "Rates", log = "y")

matplot(ltri$Px[p.ord], exp(cbind(pspl$fit[,2],
                                  pspl$se.fit[,2])[p.ord,]
                            %*% ci.mat()),
        type = "l", lty = 1, lwd = c(2,1,1), col = gray(0.2),
        xlab = "Period", ylab = "RR", log = "y")

matplot(ltri$Cx[c.ord], exp(cbind(pspl$fit[,3],
                                  pspl$se.fit[,3])[c.ord,]
                            %*% ci.mat()),
        type = "l", lty = 1, lwd = c(2,1,1), col = gray(0.2),
        xlab = "Cohort", ylab = "RR", log = "y" )


###################################################
### code chunk number 17: lungACP
###################################################
lACP <- apc.fit(A = ltri$Ax,
                P = ltri$Px,
                D = ltri$D,
                Y = ltri$Y,
            ref.c = 1900)
class(lACP) ; names(lACP)
plot(lACP)


###################################################
### code chunk number 18: apc-tri.rnw:354-355
###################################################
plot(lACP)


###################################################
### code chunk number 19: lungAPC
###################################################
lAPC <- apc.fit(A = ltri$Ax,
                P = ltri$Px,
                D = ltri$D,
                Y = ltri$Y / 1000,
             parm = "APC",
            ref.p = 1950)
plot(lAPC)