### R code from vignette source 'DMdemo-sol.rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: DMdemo-sol.rnw:9-12 ################################################### options( width=90, SweaveHooks=list( fig=function() par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,bty="n",las=1) ) ) ################################################### ### code chunk number 2: DMdemo-sol.rnw:57-63 ################################################### # options( width=120 ) library( Epi ) data( DMlate ) str( DMlate ) head( DMlate ) summary( DMlate ) ################################################### ### code chunk number 3: DMdemo-sol.rnw:84-91 ################################################### LL <- Lexis( entry = list( A = dodm-dobth, P = dodm, dur = 0 ), exit = list( P = dox ), exit.status = factor( !is.na(dodth), labels=c("Alive","Dead") ), data = DMlate ) ################################################### ### code chunk number 4: DMdemo-sol.rnw:95-96 ################################################### summary( LL ) ################################################### ### code chunk number 5: DMdemo-sol.rnw:122-124 ################################################### SL <- splitLexis( LL, breaks=seq(0,125,1), time.scale="A" ) summary( SL ) ################################################### ### code chunk number 6: DMdemo-sol.rnw:141-147 ################################################### library( splines ) r.m <- glm( (lex.Xst=="Dead") ~ ns( A, df=10 ), offset = log( lex.dur ), family = poisson, data = subset( SL, sex=="M" ) ) r.f <- update( r.m, data = subset( SL, sex=="F" ) ) ################################################### ### code chunk number 7: DMdemo-sol.rnw:178-180 ################################################### ( a.kn <- with( subset(SL,lex.Xst=="Dead"), quantile( A+lex.dur, (1:10-0.5)/10 ) ) ) ################################################### ### code chunk number 8: DMdemo-sol.rnw:187-192 ################################################### r.m <- glm( (lex.Xst=="Dead") ~ Ns( A, knots=a.kn ), offset = log( lex.dur ), family = poisson, data = subset( SL, sex=="M" ) ) r.f <- update( r.m, data = subset( SL, sex=="F" ) ) ################################################### ### code chunk number 9: DMdemo-sol.rnw:213-217 ################################################### nd <- data.frame( A = seq(10,90,0.5), lex.dur = 1000) p.m <- ci.pred( r.m, newdata = nd ) p.f <- ci.pred( r.f, newdata = nd ) ################################################### ### code chunk number 10: a-rates ################################################### matplot( seq(10,90,0.5), cbind(p.m,p.f), type="l", lty=1, lwd=c(3,1,1), las=1, col=rep(c("blue","red"),each=3), log="y", ylim=c(0.1,200), xlab="Age", ylab="Mortality rates per 1000 PY" ) ################################################### ### code chunk number 11: DMdemo-sol.rnw:255-257 ################################################### data( M.dk ) head( M.dk ) ################################################### ### code chunk number 12: DMdemo-sol.rnw:260-262 (eval = FALSE) ################################################### ## with( subset( M.dk, sex==1 & P==2005 ), lines( A, rate, col="blue", lty="12", lwd=3 ) ) ## with( subset( M.dk, sex==2 & P==2005 ), lines( A, rate, col="red" , lty="12", lwd=3 ) ) ################################################### ### code chunk number 13: a-rates-p ################################################### matplot( seq(10,90,0.5), cbind(p.m,p.f), type="l", lty=1, lwd=c(3,1,1), las=1, col=rep(c("blue","red"),each=3), log="y", ylim=c(0.05,400), xlab="Age", ylab="Mortality rates per 1000 PY" ) with( subset( M.dk, sex==1 & P==2005 ), lines( A, rate, col="blue", lty="12", lwd=3 ) ) with( subset( M.dk, sex==2 & P==2005 ), lines( A, rate, col="red" , lty="12", lwd=3 ) ) ################################################### ### code chunk number 14: DMdemo-sol.rnw:301-310 ################################################### R.m <- glm( D ~ Ns( A, knots=seq(10,90,10) ), offset = log( Y ), family = poisson, data = subset( M.dk, sex==1 & P>1994 ) ) R.f <- update( R.m, data = subset( M.dk, sex==2 & P>1994 ) ) nd <- data.frame( A = seq(10,90,0.5), Y = 1000) P.m <- ci.pred( R.m, newdata = nd ) P.f <- ci.pred( R.f, newdata = nd ) ################################################### ### code chunk number 15: a-rates-s ################################################### matplot( seq(10,90,0.5), cbind(p.m,p.f), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(0.1,200), xlab="Age", ylab="Mortality rates per 1000 PY" ) matlines( seq(10,90,0.5), cbind(P.m,P.f), lty="12", col=c("blue","red"), lwd=c(3,1,1) ) ################################################### ### code chunk number 16: DMdemo-sol.rnw:381-384 ################################################### mid.pt <- 0:999/10 + 1/20 mid.pt[1:5] nd <- data.frame( A = mid.pt, Y = 1 ) ################################################### ### code chunk number 17: DMdemo-sol.rnw:394-396 ################################################### S.m <- exp( -cumsum( ci.pred(R.m,newdata=nd)[,1]*1/10 ) ) S.f <- exp( -cumsum( ci.pred(R.f,newdata=nd)[,1]*1/10 ) ) ################################################### ### code chunk number 18: surv-pop ################################################### matplot( mid.pt+1/20, cbind(S.m,S.f), type="l", lty=1, col=c("blue","red"), lwd=3, ylim=c(0,1), xlim=c(0,100), ylab="Survival probability", xlab="Age (years)" ) ################################################### ### code chunk number 19: DMdemo-sol.rnw:426-428 ################################################### ( EL.m <- sum(S.m)/10 ) ( EL.f <- sum(S.f)/10 ) ################################################### ### code chunk number 20: DMdemo-sol.rnw:438-440 ################################################### ( EL.m <- sum(S.m)/10+1/20 ) ( EL.f <- sum(S.f)/10+1/20 ) ################################################### ### code chunk number 21: DMdemo-sol.rnw:453-458 ################################################### nd <- data.frame( A = mid.pt, lex.dur = 1 ) s.m <- exp( -cumsum( ci.pred(r.m,newdata=nd)[,1]*1/10 ) ) s.f <- exp( -cumsum( ci.pred(r.f,newdata=nd)[,1]*1/10 ) ) ( el.m <- sum(s.m)/10 ) ( el.f <- sum(s.f)/10 ) ################################################### ### code chunk number 22: DMdemo-sol.rnw:481-483 ################################################### EL.m - el.m EL.f - el.f ################################################### ### code chunk number 23: DMdemo-sol.rnw:505-509 ################################################### S50.m <- S.m[500:1000]/S.m[500] S50.f <- S.f[500:1000]/S.f[500] s50.m <- s.m[500:1000]/s.m[500] s50.f <- s.f[500:1000]/s.f[500] ################################################### ### code chunk number 24: DMdemo-sol.rnw:513-519 ################################################### ( EL50.m <- sum(S50.m)/10 ) ( EL50.f <- sum(S50.f)/10 ) ( el50.m <- sum(s50.m)/10 ) ( el50.f <- sum(s50.f)/10 ) EL50.m - el50.m EL50.f - el50.f ################################################### ### code chunk number 25: DMdemo-sol.rnw:533-548 ################################################### YLL <- function( A ) { SA.m <- S.m[(A*10):1000]/S.m[(A*10)] SA.f <- S.f[(A*10):1000]/S.f[(A*10)] sA.m <- s.m[(A*10):1000]/s.m[(A*10)] sA.f <- s.f[(A*10):1000]/s.f[(A*10)] ELA.m <- sum(SA.m)/10 ELA.f <- sum(SA.f)/10 elA.m <- sum(sA.m)/10 elA.f <- sum(sA.f)/10 c( Men=ELA.m - elA.m, Women=ELA.f - elA.f ) } YLL( 50 ) YLL( 60 ) ################################################### ### code chunk number 26: DMdemo-sol.rnw:551-557 ################################################### yll <- matrix( NA, 41, 2 ) rownames( yll ) <- 40:80 colnames( yll ) <- c("Men","Women") t(yll) for( a in 40:80 ) yll[a-39,] <- YLL(a) round( t(yll), 1 ) ################################################### ### code chunk number 27: yll ################################################### matplot( 40:80, yll, type="l", lty=1, lwd=3, col=c("blue","red"), ylim=c(0,7), ylab="Years of life lost to DM", xlab="Age at evaluation" ) ################################################### ### code chunk number 28: DMdemo-sol.rnw:591-595 ################################################### rp.m <- update( r.m, . ~ . + P ) rp.f <- update( r.f, . ~ . + P ) Rp.m <- update( R.m, . ~ . + P ) Rp.f <- update( R.f, . ~ . + P ) ################################################### ### code chunk number 29: DMdemo-sol.rnw:604-605 ################################################### ci.exp( rp.m, subset="P" ) ################################################### ### code chunk number 30: DMdemo-sol.rnw:608-616 ################################################### chg <- cbind( rbind( ci.exp(rp.m,subset="P"), ci.exp(rp.f,subset="P") ), rbind( ci.exp(Rp.m,subset="P"), ci.exp(Rp.f,subset="P") ) ) rownames( chg ) <- c("DM","Pop") colnames( chg ) <- paste( c("M:","","","F:","",""), colnames( chg ) ) round( (chg-1)*100, 1 ) ################################################### ### code chunk number 31: DMdemo-sol.rnw:636-642 ################################################### a.pt <- 0:999/10 + 1/20 mort <- NArray( list( A = a.pt, P = seq(1995,2010,5), sex = c("M","W"), type = c("DM","Pop") ) ) str( mort ) ################################################### ### code chunk number 32: DMdemo-sol.rnw:646-660 ################################################### nd <- data.frame( A = a.pt, Y = 1, lex.dur = 1 ) for( ip in seq(1995,2010,5) ) { nd$P <- ip mort[,paste(ip),"M","DM"] <- ci.pred( rp.m, newdata=nd )[,1] mort[,paste(ip),"W","DM"] <- ci.pred( rp.f, newdata=nd )[,1] mort[,paste(ip),"M","Pop"] <- ci.pred( Rp.m, newdata=nd )[,1] mort[,paste(ip),"W","Pop"] <- ci.pred( Rp.f, newdata=nd )[,1] } round( ftable( mort[1:5+600,,,]*1000, row.vars=c(4,1), col.vars=c(3,2) ), 1 ) ################################################### ### code chunk number 33: DMdemo-sol.rnw:662-665 (eval = FALSE) ################################################### ## matplot( a.pt, cbind( mort[,,"M","DM"], ## mort[,,"M","Pop"] ), ## log="y", type="l", lwd=3, lty=1 ) ################################################### ### code chunk number 34: DMdemo-sol.rnw:670-675 ################################################### Mort <- apply( mort, 2:4, cumsum )/10 str( Mort ) round( ftable( Mort[200+1:5,,,], row.vars=c(4,1), col.vars=c(3,2) ), 5 ) ################################################### ### code chunk number 35: M-surv ################################################### Surv <- exp( -Mort ) matplot( a.pt, cbind( Surv[,,"M","DM"], Surv[,,"M","Pop"] ), type="l", lwd=3, lty=1, col=heat.colors(6)[1:4] ) ################################################### ### code chunk number 36: M-surv40 ################################################### S40 <- Surv[400:1000,,,]/Surv[rep(400,601),,,] matplot( a.pt[400:1000], cbind( S40[,,"M","DM"], S40[,,"M","Pop"] ), type="l", lwd=3, lty=1, col=heat.colors(6)[1:4] ) ################################################### ### code chunk number 37: DMdemo-sol.rnw:710-713 ################################################### ERL40 <- apply( S40, 2:4, sum )/10 str( ERL40 ) ftable( ERL40 ) ################################################### ### code chunk number 38: DMdemo-sol.rnw:716-720 ################################################### ERL40 <- ERL40[,,c(1,2,2)] dimnames(ERL40)[[3]][3] <- "YLL" ERL40[,,"YLL"] <- ERL40[,,"Pop"] - ERL40[,,"DM"] ftable( ERL40, row.vars=2:1 ) ################################################### ### code chunk number 39: DMdemo-sol.rnw:735-736 ################################################### musmuc <- function(x) rev(cumsum(rev(x))) ################################################### ### code chunk number 40: DMdemo-sol.rnw:741-746 ################################################### ERL <- ( apply( Surv, 2:4, musmuc )/10 ) / Surv LL <- ERL[,,,c(1,2,2)] dimnames( LL )[[4]][3] <- "YLL" LL[,,,3] <- LL[,,,"Pop"]-LL[,,,"DM"] str( LL ) ################################################### ### code chunk number 41: YLL4 ################################################### matplot( a.pt, cbind( LL[,,"M","YLL"], LL[,,"W","YLL"] ), type="l", lty=1, lwd=1:4, col=rep(c("blue","red"),each=4), xlim=c(40,80), ylim=c(0,10), yaxs="i", xlab="Age alive", ylab="Years of life lost to DM" ) abline( h=0 )