### R code from vignette source 'DMDK-s' ### Encoding: UTF-8 ################################################### ### code chunk number 1: DMDK-s.rnw:2-7 ################################################### options( width=90, # prompt=" ", continue=" ", # Absence of prompts makes it easier for # students to copy from the final pdf document 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: DMDK-s.rnw:43-49 ################################################### options( width=90 ) library( Epi ) data( DMlate ) str( DMlate ) head( DMlate ) summary( DMlate ) ################################################### ### code chunk number 3: DMDK-s.rnw:73-75 ################################################### with( DMlate, table( dead=!is.na(dodth), same=(dodth==dox), exclude=NULL ) ) ################################################### ### code chunk number 4: DMDK-s.rnw:79-86 ################################################### 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 5: DMDK-s.rnw:93-95 ################################################### summary( LL ) head( LL ) ################################################### ### code chunk number 6: DMDK-s.rnw:104-109 ################################################### stat.table( sex, list( D=sum( lex.Xst=="Dead" ), Y=sum( lex.dur ), rate=ratio( lex.Xst=="Dead", lex.dur, 1000 ) ), data=LL ) ################################################### ### code chunk number 7: DMDK-s.rnw:137-140 ################################################### SL <- splitLexis( LL, breaks=seq(0,125,1/2), time.scale="A" ) summary( SL ) summary( LL ) ################################################### ### code chunk number 8: DMDK-s.rnw:162-168 ################################################### 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 9: DMDK-s.rnw:201-206 ################################################### 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 ) str( p.m ) ################################################### ### code chunk number 10: a-rates ################################################### matplot( nd$A, cbind(p.m,p.f), type="l", col=rep(c("blue","red"),each=3), lwd=c(3,1,1), lty=1, log="y", xlab="Age", ylab="Mortality of DM ptt per 1000 PY") ################################################### ### code chunk number 11: DMDK-s.rnw:317-323 ################################################### ( kn.A <- with( subset( SL, lex.Xst=="Dead" ), quantile( A+lex.dur, probs=seq(5,95,20)/100 ) ) ) ( kn.P <- with( subset( SL, lex.Xst=="Dead" ), quantile( P+lex.dur, probs=seq(5,95,20)/100 ) ) ) ( kn.dur <- c(0,with( subset( SL, lex.Xst=="Dead" ), quantile( dur+lex.dur, probs=seq(10,90,20)/100 ) )) ) ################################################### ### code chunk number 12: DMDK-s.rnw:343-352 ################################################### Mm <- glm( (lex.Xst=="Dead") ~ Ns( A, kn=kn.A ) + Ns( P, kn=kn.P ) + Ns( dur, kn=kn.dur ), offset = log( lex.dur ), family = poisson, data = subset( SL, sex=="M" ) ) summary( Mm ) Mf <- update( Mm, data = subset( SL, sex=="F" ) ) round( cbind( ci.exp(Mm), ci.exp(Mf) ), 3 ) ################################################### ### code chunk number 13: DMDK-s.rnw:366-368 ################################################### anova( Mm, r.m, test="Chisq" ) anova( Mf, r.f, test="Chisq" ) ################################################### ### code chunk number 14: DMDK-s.rnw:414-421 ################################################### mm <- glm( (lex.Xst=="Dead") ~ Ns( A, kn=kn.A , intercept=TRUE ) - 1 + Ns( P, kn=kn.P , ref=2000 ) + Ns( dur, kn=kn.dur, ref=5 ), offset = log( lex.dur/100 ), family = poisson, data = subset( SL, sex=="M" ) ) mf <- update( mm, data = subset( SL, sex=="F" ) ) ################################################### ### code chunk number 15: DMDK-s.rnw:425-426 ################################################### c( deviance(Mm), deviance(mm) ) ################################################### ### code chunk number 16: Tpl-m ################################################### Termplot( mm ) ################################################### ### code chunk number 17: Tpl-f ################################################### Termplot( mf ) ################################################### ### code chunk number 18: DMDK-s.rnw:499-505 ################################################### pts <- seq(0,20,2) nd <- data.frame( A= 50+pts, P=1995+pts, dur= pts, lex.dur=1000 ) cbind( nd$A, ci.pred( mm, newdata=nd ) ) ################################################### ### code chunk number 19: DMDK-s.rnw:513-526 ################################################### mpr <- fpr <- NULL pts <- seq(0,20,0.1) for( ip in c(1995,2005) ) for( ia in c(50,60,70) ) { nd <- data.frame( A=ia+pts, P=ip+pts, dur= pts, lex.dur=1000 ) mpr <- cbind( mpr, ci.pred( mm, nd) ) fpr <- cbind( fpr, ci.pred( mf, nd) ) } str( fpr ) ################################################### ### code chunk number 20: rates ################################################### par( mfrow=c(1,2) ) matplot( cbind(50+pts,60+pts,70+pts)[,rep(1:3,2,each=3)], cbind( mpr[,1:9], fpr[,1:9] ), ylim=c(5,500), log="y", xlab="Age", ylab="Mortality, diagnosed 1995", type="l", lwd=c(4,1,1), lty=1, col=rep(c("blue","red"),each=9) ) matplot( cbind(50+pts,60+pts,70+pts)[,rep(1:3,2,each=3)], cbind( mpr[,1:9+9], fpr[,1:9+9] ), ylim=c(5,500), log="y", xlab="Age", ylab="Mortality, diagnosed 2005", type="l", lwd=c(4,1,1), lty=1, col=rep(c("blue","red"),each=9) ) ################################################### ### code chunk number 21: DMDK-s.rnw:631-640 ################################################### str( SL ) SL$Am <- floor( SL$A+0.25 ) SL$Pm <- floor( SL$P+0.25 ) data( M.dk ) str( M.dk ) M.dk <- transform( M.dk, Am = A, Pm = P, sex = factor( sex, labels=c("M","F") ) ) str( M.dk ) ################################################### ### code chunk number 22: DMDK-s.rnw:645-648 ################################################### SLr <- merge( SL, M.dk[,c("sex","Am","Pm","rate")] ) dim( SL ) dim( SLr ) ################################################### ### code chunk number 23: DMDK-s.rnw:664-677 ################################################### SLr$E <- SLr$lex.dur * SLr$rate / 1000 stat.table( sex, list( D = sum(lex.Xst=="Dead"), Y = sum(lex.dur), E = sum(E), SMR = ratio(lex.Xst=="Dead",E) ), data = SLr ) stat.table( list( Age = floor(pmax(A,39)/10)*10 ), list( D = sum(lex.Xst=="Dead"), Y = sum(lex.dur), E = sum(E), SMR = ratio(lex.Xst=="Dead",E) ), data = SLr ) ################################################### ### code chunk number 24: DMDK-s.rnw:700-707 ################################################### sm <- glm( (lex.Xst=="Dead") ~ Ns( A, kn=kn.A , intercept=TRUE ) - 1 + Ns( P, kn=kn.P , ref=2000 ) + Ns( dur, kn=kn.dur, ref=5 ), offset = log( E ), family = poisson, data = subset( SLr, E>0 & sex=="M" ) ) sf <- update( mm, data = subset( SLr, E>0 & sex=="F" ) ) ################################################### ### code chunk number 25: Tpl-sm ################################################### Termplot( sm ) ################################################### ### code chunk number 26: Tpl-sf ################################################### Termplot( sf ) ################################################### ### code chunk number 27: DMDK-s.rnw:728-739 ################################################### s0 <- glm( (lex.Xst=="Dead") ~ Ns( A, kn=kn.A , intercept=TRUE ) - 1 + Ns( P, kn=kn.P , ref=2000 ) + Ns( dur, kn=kn.dur, ref=5 ), offset = log( E ), family = poisson, data = subset( SLr, E>0 ) ) s1 <- update( s0, . ~ . + sex ) sA <- update( s1, . ~ . + sex:A ) sAP <- update( sA, . ~ . + sex:P ) sAPd <- update( sAP, . ~ . + sex:dur ) anova( s0, s1, sA, sAP, sAPd, test="Chisq" ) ################################################### ### code chunk number 28: DMDK-s.rnw:745-749 ################################################### sA <- update( s1, . ~ . + sex:A + sex:I(A^2) ) sAP <- update( sA, . ~ . + sex:P + sex:I(P^2) ) sAPd <- update( sAP, . ~ . + sex:dur + sex:I(dur^2) ) anova( s0, s1, sA, sAP, sAPd, test="Chisq" ) ################################################### ### code chunk number 29: Tpl-s0 ################################################### Termplot( s0 ) ################################################### ### code chunk number 30: DMDK-s.rnw:767-779 ################################################### psmr <- NULL pts <- seq(0,20,0.1) for( ip in c(1995,2005) ) for( ia in c(50,60,70) ) { nd <- data.frame( A=ia+pts, P=ip+pts, dur= pts, E=1 ) psmr <- cbind( psmr, ci.pred( s0, nd) ) } str( psmr ) ################################################### ### code chunk number 31: smr ################################################### par( mfrow=c(1,2) ) matplot( cbind(50+pts,60+pts,70+pts)[,rep(1:3,each=3)], psmr[,1:9], ylim=c(0.7,7), log="y", xlab="Age", ylab="SMR, diagnosed 1995", type="l", lwd=c(4,1,1), lty=1, col="black" ) abline( h=1 ) matplot( cbind(50+pts,60+pts,70+pts)[,rep(1:3,each=3)], psmr[,1:9+9], ylim=c(0.7,7), log="y", xlab="Age", ylab="SMR, diagnosed 2005", type="l", lwd=c(4,1,1), lty=1, col="black" ) abline( h=1 ) ################################################### ### code chunk number 32: DMDK-s.rnw:821-828 ################################################### sx <- glm( (lex.Xst=="Dead") ~ I(A-dur) + I(P-2000) + Ns( dur, kn=c(0,1,3,6,10), ref=5 ), offset = log( E ), family = poisson, data = subset( SLr, E>0 ) ) anova( s0, sx, test="Chisq" ) ################################################### ### code chunk number 33: DMDK-s.rnw:834-836 ################################################### round( ci.exp(sx), 4 ) round( ci.exp( sx, subset=c("A","P"), ctr.mat=10*diag(2) ), 4 ) ################################################### ### code chunk number 34: xsmr ################################################### xsmr <- NULL for( ip in c(1995,2005) ) for( ia in c(50,60,70) ) { nd <- data.frame( A=ia+pts, P=ip+pts, dur= pts, E=1 ) xsmr <- cbind( xsmr, ci.pred( sx, nd) ) } par( mfrow=c(1,2) ) matplot( cbind(50+pts,60+pts,70+pts)[,rep(1:3,each=3)], xsmr[,1:9], ylim=c(0.7,7), log="y", xlab="Age", ylab="SMR, diagnosed 1995", type="l", lwd=c(4,1,1), lty=1, col="black" ) abline( h=1 ) matplot( cbind(50+pts,60+pts,70+pts)[,rep(1:3,each=3)], xsmr[,1:9+9], ylim=c(0.7,7), log="y", xlab="Age", ylab="SMR, diagnosed 2005", type="l", lwd=c(4,1,1), lty=1, col="black" ) abline( h=1 )