yll.R
In Epi: Statistical Analysis in Epidemiology

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

###################################################
### code chunk number 1: yll.rnw:21-27
###################################################
options(width=90,
         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: yll.rnw:31-33
###################################################
anfang <- Sys.time()
cat("Start time:", format(anfang, "%F, %T"), "\n")


###################################################
### code chunk number 3: states
###################################################
getOption("SweaveHooks")[["fig"]]()
library(Epi)
TM <- matrix(NA, 4, 4)
rownames(TM) <-
colnames(TM) <- c("Well", "DM", "Dead", "Dead(DM)")
TM[1, 2:3] <- TM[2, 4] <- 1
TM
zz <- boxes(TM, boxpos = list(x = c(20, 80, 20, 80),
                              y = c(80, 80, 20, 20)),
                wm = 1.5,
                hm = 4)


###################################################
### code chunk number 4: states
###################################################
getOption("SweaveHooks")[["fig"]]()
zz$Arrowtext <- c(expression(lambda(a)),
                  expression(mu[W](a)),
                  expression(mu[D][M](a,d)))
boxes.MS(zz)


###################################################
### code chunk number 5: yll.rnw:280-281
###################################################
data(DMepi)


###################################################
### code chunk number 6: yll.rnw:287-289
###################################################
str(DMepi)
head(DMepi)


###################################################
### code chunk number 7: yll.rnw:309-315
###################################################
DMepi <- transform(subset(DMepi, A > 30),
                   A = A + 0.5,
                   P = P + 0.5,
                 D.T = D.nD + D.DM,
                 Y.T = Y.nD + Y.DM)
head(DMepi)


###################################################
### code chunk number 8: yll.rnw:321-346
###################################################
# Knots used in all models
(a.kn <- seq(40, 95, , 6))
(p.kn <- seq(1997, 2015, , 4))
(c.kn <- seq(1910, 1976, , 6))
# Check the number of events between knots
ae <- xtabs(cbind(D.nD, D.DM, X) ~ cut(A, c(30, a.kn, Inf)) + sex, data=DMepi)
ftable(addmargins(ae, 1), col.vars=3:2)
pe <- xtabs(cbind(D.nD, D.DM, X) ~ cut(P, c(1990, p.kn, Inf)) + sex, data=DMepi)
ftable(addmargins(pe, 1), col.vars=3:2)
ce <- xtabs(cbind(D.nD, D.DM, X) ~ cut(P-A, c(-Inf, c.kn, Inf)) + sex, data=DMepi)
ftable(addmargins(ce, 1), col.vars=3:2)
# Fit an APC-model for all transitions, separately for men and women
mW.m <- glm(cbind(D.nD, Y.nD) ~ -1 + Ns(    A, knots=a.kn, int=TRUE) +
                                     Ns(P    , knots=p.kn, ref=2005) +
                                     Ns(P - A, knots=c.kn, ref=1950),
            family = poisreg,
              data = subset(DMepi, sex=="M"))
mD.m <- update(mW.m, cbind(D.DM, Y.DM) ~ .)
mT.m <- update(mW.m, cbind(D.T , Y.T ) ~ .)
lW.m <- update(mW.m, cbind(X   , Y.nD) ~ .)
# Model for women
mW.f <- update(mW.m, data = subset(DMepi, sex == "F"))
mD.f <- update(mD.m, data = subset(DMepi, sex == "F"))
mT.f <- update(mT.m, data = subset(DMepi, sex == "F"))
lW.f <- update(lW.m, data = subset(DMepi, sex == "F"))


###################################################
### code chunk number 9: yll.rnw:353-390
###################################################
a.ref <- 30:90
p.ref <- 1996:2016
aYLL <- NArray(list(type = c("Imm", "Tot", "Sus"),
                       sex = levels(DMepi$sex),
                       age = a.ref,
                      date = p.ref))
str(aYLL)
system.time(
for(ip in p.ref)
   {
   nd <- data.frame(A = seq(30, 90, 0.2)+0.1,
                     P = ip,
                  Y.nD = 1,
                  Y.DM = 1,
                  Y.T  = 1)
   muW.m <- ci.pred(mW.m, nd)[, 1]
   muD.m <- ci.pred(mD.m, nd)[, 1]
   muT.m <- ci.pred(mT.m, nd)[, 1]
   lam.m <- ci.pred(lW.m, nd)[, 1]
   muW.f <- ci.pred(mW.f, nd)[, 1]
   muD.f <- ci.pred(mD.f, nd)[, 1]
   muT.f <- ci.pred(mT.f, nd)[, 1]
   lam.f <- ci.pred(lW.f, nd)[, 1]
   aYLL["Imm", "M", , paste(ip)] <- yll(int=0.2, muW.m, muD.m, lam=NULL,
                                      A=a.ref, age.in=30, note=FALSE)[-1]
   aYLL["Imm", "F", , paste(ip)] <- yll(int=0.2, muW.f, muD.f, lam=NULL,
                                      A=a.ref, age.in=30, note=FALSE)[-1]
   aYLL["Tot", "M", , paste(ip)] <- yll(int=0.2, muT.m, muD.m, lam=NULL,
                                      A=a.ref, age.in=30, note=FALSE)[-1]
   aYLL["Tot", "F", , paste(ip)] <- yll(int=0.2, muT.f, muD.f, lam=NULL,
                                      A=a.ref, age.in=30, note=FALSE)[-1]
   aYLL["Sus", "M", , paste(ip)] <- yll(int=0.2, muW.m, muD.m, lam=lam.m,
                                      A=a.ref, age.in=30, note=FALSE)[-1]
   aYLL["Sus", "F", , paste(ip)] <- yll(int=0.2, muW.f, muD.f, lam=lam.f,
                                      A=a.ref, age.in=30, note=FALSE)[-1]
   })
round(ftable(aYLL[, , seq(1, 61, 10), ], col.vars=c(3, 2)), 1)


###################################################
### code chunk number 10: imm
###################################################
getOption("SweaveHooks")[["fig"]]()
plyll <- function(wh, xtxt){
par(mfrow = c(1, 2),
      mar = c(3, 3, 1, 1),
      mgp = c(3, 1, 0) / 1.6,
      bty = "n",
      las = 1)
matplot(a.ref, aYLL[wh, "M", , ],
         type="l", lty=1, col="blue", lwd=1:2,
         ylim=c(0, 12), xlab="Age",
         ylab=paste0("Years lost to DM", xtxt),
         yaxs="i")
abline(v=50, h=1:11, col=gray(0.7))
text(90, 11.5, "Men", col="blue", adj=1)
text(40, aYLL[wh, "M", "40", "1996"], "1996", adj=c(0, 0), col="blue")
text(43, aYLL[wh, "M", "44", "2016"], "2016", adj=c(1, 1), col="blue")

matplot(a.ref, aYLL[wh, "F", , ],
         type="l", lty=1, col="red", lwd=1:2,
         ylim=c(0, 12), xlab="Age",
         ylab=paste0("Years lost to DM", xtxt),
         yaxs="i")
abline(v=50, h=1:11, col=gray(0.7))
text(90, 11.5, "Women", col="red", adj=1)
text(40, aYLL[wh, "F", "40", "1996"], "1996", adj=c(0, 0), col="red")
text(43, aYLL[wh, "F", "44", "2016"], "2016", adj=c(1, 1), col="red")
}
plyll("Imm", " - immunity assumption")


###################################################
### code chunk number 11: tot
###################################################
getOption("SweaveHooks")[["fig"]]()
plyll("Tot", " - total mortality refernce")


###################################################
### code chunk number 12: sus
###################################################
getOption("SweaveHooks")[["fig"]]()
plyll("Sus", " - susceptibility assumed")


###################################################
### code chunk number 13: yll.rnw:473-477
###################################################
ende <- Sys.time()
cat("  Start time:", format(anfang, "%F, %T"), "\n")
cat("    End time:", format(  ende, "%F, %T"), "\n")
cat("Elapsed time:", round(difftime(ende, anfang, units = "mins"), 2), "minutes\n")