yll.R
In Epi: Statistical Analysis in Epidemiology

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

###################################################
### code chunk number 1: yll.rnw:21-24
###################################################
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: 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
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 3: states
###################################################
getOption("SweaveHooks")[["fig"]]()
zz$Arrowtext <- c(expression(lambda),   
                  expression(mu[W]),    
                  expression(mu[D][M]))
boxes(zz)


###################################################
### code chunk number 4: yll.rnw:269-270
###################################################
data(DMepi)


###################################################
### code chunk number 5: yll.rnw:275-277
###################################################
str(DMepi)
head(DMepi)


###################################################
### code chunk number 6: yll.rnw:297-301
###################################################
DMepi <- transform(subset(DMepi, A>30), 
                    D.T = D.nD + D.DM, 
                    Y.T = Y.nD + Y.DM)
head(DMepi)


###################################################
### code chunk number 7: yll.rnw:307-333
###################################################
# 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, seperately for men and women
mW.m <- glm(D.nD ~ -1 + Ns(A  , knots=a.kn, int=TRUE) +
                         Ns( P, knots=p.kn, ref=2005) +
                         Ns(P-A, knots=c.kn, ref=1950), 
           offset = log(Y.nD), 
           family = poisson, 
             data = subset(DMepi, sex=="M"))
mD.m <- update(mW.m,  D.DM ~ . , offset=log(Y.DM))
mT.m <- update(mW.m,  D.T  ~ . , offset=log(Y.T))
lW.m <- update(mW.m,  X ~ .)
# 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 8: yll.rnw:340-377
###################################################
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 9: imm
###################################################
getOption("SweaveHooks")[["fig"]]()
plyll <- function(wh){
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="Years lost to DM", yaxs="i")
abline(v=50, h=1:10, col=gray(0.7))
text(90, 11, "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="Years lost to DM", yaxs="i")
abline(v=50, h=1:10, col=gray(0.7))
text(90, 11, "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")


###################################################
### code chunk number 10: tot
###################################################
getOption("SweaveHooks")[["fig"]]()
plyll("Tot")


###################################################
### code chunk number 11: sus
###################################################
getOption("SweaveHooks")[["fig"]]()
plyll("Sus")