### R code from vignette source 'DMDK-demo.rnw' ################################################### ### code chunk number 1: DMDK-demo.rnw:2-7 ################################################### options( width=90, # prompt=" ", continue=" ", # 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-demo.rnw:27-33 ################################################### options( width=120 ) library( Epi ) data( DMlate ) str( DMlate ) head( DMlate ) summary( DMlate ) ################################################### ### code chunk number 3: DMDK-demo.rnw:43-45 ################################################### with( DMlate, table( dead=!is.na(dodth), same=(dodth==dox), exclude=NULL ) ) ################################################### ### code chunk number 4: DMDK-demo.rnw:49-56 ################################################### 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-demo.rnw:60-61 ################################################### summary( LL ) ################################################### ### code chunk number 6: DMDK-demo.rnw:65-70 ################################################### 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-demo.rnw:86-88 ################################################### SL <- splitLexis( LL, breaks=seq(0,125,1), time.scale="A" ) summary( SL ) ################################################### ### code chunk number 8: DMDK-demo.rnw:95-101 ################################################### library( splines ) r.m <- glm( (lex.Xst=="Dead") ~ ns( A, df=10, intercept=TRUE ) - 1, 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-demo.rnw:106-110 ################################################### 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: DMDK-demo.rnw:128-130 ################################################### data( M.dk ) head( M.dk ) ################################################### ### code chunk number 12: DMDK-demo.rnw:133-135 (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.1,200), 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: DMDK-demo.rnw:154-163 ################################################### R.m <- glm( D ~ ns( A, df=10, intercept=TRUE ) - 1, 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: DMDK-demo.rnw:200-206 ################################################### ( kn.A <- with( subset( SL, lex.Xst=="Dead" ), quantile( A+lex.dur, probs=seq(5,95,10)/100 ) ) ) ( kn.P <- with( subset( SL, lex.Xst=="Dead" ), quantile( P+lex.dur, probs=seq(5,95,30)/100 ) ) ) ( kn.dur <- c(0,with( subset( SL, lex.Xst=="Dead" ), quantile( dur+lex.dur, probs=seq(5,95,10)/100 ) ) ) ) ################################################### ### code chunk number 17: DMDK-demo.rnw:210-218 ################################################### 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" ) ) ################################################### ### code chunk number 18: DMDK-demo.rnw:222-224 ################################################### anova( mm, r.m, test="Chisq" ) anova( mf, r.f, test="Chisq" ) ################################################### ### code chunk number 19: DMDK-demo.rnw:260-266 ################################################### N <- 100 pr.A <- seq(10,90,,N) pr.P <- seq(1995,2010,,N) pr.d <- seq(0,15,,N) rf.P <- 2009 rf.d <- 2 ################################################### ### code chunk number 20: DMDK-demo.rnw:270-275 ################################################### AC <- Ns( pr.A, knots=kn.A ) PC <- Ns( pr.P, knots=kn.P ) dC <- Ns( pr.d, knots=kn.dur ) PR <- Ns( rep(rf.P,N), knots=kn.P ) dR <- Ns( rep(rf.d,N), knots=kn.dur ) ################################################### ### code chunk number 21: DMDK-demo.rnw:292-298 ################################################### m.A <- ci.exp( mm, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 f.A <- ci.exp( mf, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 m.P <- ci.exp( mm, subset="P" , ctr.mat=PC-PR ) f.P <- ci.exp( mf, subset="P" , ctr.mat=PC-PR ) m.d <- ci.exp( mm, subset="dur", ctr.mat=dC-dR ) f.d <- ci.exp( mf, subset="dur", ctr.mat=dC-dR ) ################################################### ### code chunk number 22: effects0 ################################################### par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, cbind(m.A,f.A), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", xlab="Age", ylab="Mortality rate per 1000 PY" ) matplot( pr.P, cbind(m.P,f.P), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", xlab="Date of follow-up", ylab="Mortality rate ratio" ) matplot( pr.d, cbind(m.d,f.d), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", xlab="Diabetes duration", ylab="Mortality rate ratio" ) ################################################### ### code chunk number 23: effects1 ################################################### kn.dur <- c(0,with( subset( SL, lex.Xst=="Dead" ), quantile( dur+lex.dur, probs=seq(5,95,30)/100 ) )) dC <- Ns( pr.d, knots=kn.dur ) dR <- Ns( rep(rf.d,N), knots=kn.dur ) 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" ) ) mf <- update( mm, data = subset( SL, sex=="F" ) ) m.A <- ci.exp( mm, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 f.A <- ci.exp( mf, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 m.P <- ci.exp( mm, subset="P" , ctr.mat=PC-PR ) f.P <- ci.exp( mf, subset="P" , ctr.mat=PC-PR ) m.d <- ci.exp( mm, subset="dur", ctr.mat=dC-dR ) f.d <- ci.exp( mf, subset="dur", ctr.mat=dC-dR ) par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, cbind(m.A,f.A), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(0.2,180), xlab="Age", ylab="Mortality rate per 1000 PY" ) matplot( pr.P, cbind(m.P,f.P), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/30,30), xlab="Date of follow-up", ylab="Mortality rate ratio" ) abline( h=1 ) matplot( pr.d, cbind(m.d,f.d), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/30,30), xlab="Diabetes duration", ylab="Mortality rate ratio" ) abline( h=1 ) ################################################### ### code chunk number 24: effects2 ################################################### par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, cbind(m.A,f.A), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(0.2,180), xlab="Age", ylab="Mortality rate per 1000 PY" ) matplot( pr.P, cbind(m.P,f.P), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/3,3), xlab="Date of follow-up", ylab="Mortality rate ratio" ) abline( h=1 ) matplot( pr.d, cbind(m.d,f.d), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/3,3), xlab="Diabetes duration", ylab="Mortality rate ratio" ) abline( h=1 ) ################################################### ### code chunk number 25: DMDK-demo.rnw:399-407 ################################################### m2 <- glm( (lex.Xst=="Dead") ~ sex + sex:Ns( A, kn=kn.A ) + Ns( P, kn=kn.P ) + Ns( dur, kn=kn.dur ), offset = log( lex.dur ), family = poisson, data = SL ) ci.exp(m2) ################################################### ### code chunk number 26: DMDK-demo.rnw:419-422 ################################################### j.dev <- mm$dev + mf$dev j.df <- mm$df.r + mf$df.r + 1 1 - pchisq( m2$dev - j.dev, m2$df.r - j.df ) ################################################### ### code chunk number 27: DMDK-demo.rnw:440-442 ################################################### ci.exp( m2, subset=c("Int","sexM","P","dur") ) ci.exp( m2, subset=c("Int","sexF","P","dur") ) ################################################### ### code chunk number 28: DMDK-demo.rnw:446-450 ################################################### mi.A <- ci.exp( m2, subset=c("Int","sexM","P","dur"), ctr.mat=cbind(1 ,AC,PR,dR) ) * 1000 fi.A <- ci.exp( m2, subset=c("Int","sexF","P","dur"), ctr.mat=cbind(1,1,AC,PR,dR) ) * 1000 b.P <- ci.exp( m2, subset="P" , ctr.mat=PC-PR ) b.d <- ci.exp( m2, subset="dur", ctr.mat=dC-dR ) ################################################### ### code chunk number 29: effects3 ################################################### par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, cbind(m.A,f.A,mi.A,fi.A), type="l", lty=rep(c(3,1),each=6), lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(0.2,180), xlab="Age", ylab="Mortality rate per 1000 PY" ) matplot( pr.P, cbind(m.P,f.P,b.P), type="l", lty=rep(c(3,1),c(6,3)), lwd=c(3,1,1), col=rep(c("blue","red","black"),each=3), log="y", ylim=c(1/3,3), xlab="Date of follow-up", ylab="Mortality rate ratio" ) abline( h=1 ) matplot( pr.d, cbind(m.d,f.d,b.d), type="l", lty=rep(c(3,1),c(6,3)), lwd=c(3,1,1), col=rep(c("blue","red","black"),each=3), log="y", ylim=c(1/3,3), xlab="Diabetes duration", ylab="Mortality rate ratio" ) abline( h=1 ) ################################################### ### code chunk number 30: DMDK-demo.rnw:492-498 ################################################### pts <- seq(0,20,1) nd <- data.frame( A= 50+pts, P=1995+pts, dur= pts, lex.dur=1000 ) ci.pred( mm, newdata=nd ) ################################################### ### code chunk number 31: DMDK-demo.rnw:501-501 ################################################### ################################################### ### code chunk number 32: DMDK-demo.rnw:506-517 ################################################### DMm <- function( A, P, pts=seq(0,25,0.1) ) { nd <- data.frame( A=A+pts, P=P+pts, dur= pts, lex.dur=1000 ) cbind( nd$A, ci.pred( mm, newdata=nd ), ci.pred( mf, newdata=nd ) ) } DMm( 50, 1996, pts=0:10 ) ################################################### ### code chunk number 33: morts ################################################### DMm.1996 <- rbind( DMm( 30, 1996 ), NA, DMm( 40, 1996 ), NA, DMm( 50, 1996 ), NA, DMm( 60, 1996 ), NA, DMm( 70, 1996 ), NA, DMm( 80, 1996 ), NA, DMm( 90, 1996 ) ) DMm.2005 <- rbind( DMm( 30, 2005 ), NA, DMm( 40, 2005 ), NA, DMm( 50, 2005 ), NA, DMm( 60, 2005 ), NA, DMm( 70, 2005 ), NA, DMm( 80, 2005 ), NA, DMm( 90, 2005 ) ) par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( DMm.1996[,1], DMm.1996[,-1], type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1,1000), xlim=c(30,95), las=1, xlab="Age", ylab="Mortality rate per 1000 PY" ) text( 30, 1000, "DM diagnosed 1996", adj=c(0,1) ) matplot( DMm.2005[,1], DMm.2005[,-1], type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1,1000), xlim=c(30,95), las=1, xlab="Age", ylab="Mortality rate per 1000 PY" ) text( 30, 1000, "DM diagnosed 2005", adj=c(0,1) ) ################################################### ### code chunk number 34: DMDK-demo.rnw:581-585 ################################################### kn.Ad <- with( subset( SL, lex.Xst=="Dead" ), quantile( A-dur, probs=seq(5,95,10)/100 ) ) kn.Pd <- with( subset( SL, lex.Xst=="Dead" ), quantile( P-dur, probs=seq(5,95,20)/100 ) ) ################################################### ### code chunk number 35: DMDK-demo.rnw:589-605 ################################################### anova( mm, update( mm, . ~ . + Ns(A-dur,knots=kn.Ad) ), update( mm, . ~ . + Ns(A-dur,knots=kn.Ad) - Ns(A,knots=kn.A) ), test = "Chisq" ) anova( mm, update( mm, . ~ . + Ns(P-dur,knots=kn.Pd) ), update( mm, . ~ . + Ns(P-dur,knots=kn.Pd) - Ns(P,knots=kn.P) ), test = "Chisq" ) anova( mf, update( mf, . ~ . + Ns(A-dur,knots=kn.Ad) ), update( mf, . ~ . + Ns(A-dur,knots=kn.Ad) - Ns(A,knots=kn.A) ), test = "Chisq" ) anova( mf, update( mf, . ~ . + Ns(P-dur,knots=kn.Pd) ), update( mf, . ~ . + Ns(P-dur,knots=kn.Pd) - Ns(P,knots=kn.P) ), test = "Chisq" ) ################################################### ### code chunk number 36: DMDK-demo.rnw:624-642 ################################################### AC <- Ns( pr.A, knots=kn.Ad ) PC <- Ns( pr.P, knots=kn.Pd ) PR <- Ns( rep(rf.P,N), knots=kn.Pd ) Mm <- glm( (lex.Xst=="Dead") ~ Ns( A-dur, kn=kn.Ad ) + Ns( P-dur, kn=kn.Pd ) + Ns( dur, kn=kn.dur ), offset = log( lex.dur ), family = poisson, data = subset( SL, sex=="M" ) ) Mf <- update( Mm, data = subset( SL, sex=="F" ) ) M.A <- ci.exp( Mm, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 F.A <- ci.exp( Mf, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 M.P <- ci.exp( Mm, subset="P" , ctr.mat=PC-PR ) F.P <- ci.exp( Mf, subset="P" , ctr.mat=PC-PR ) M.d <- ci.exp( Mm, subset="kn.dur", ctr.mat=dC-dR ) F.d <- ci.exp( Mf, subset="kn.dur", ctr.mat=dC-dR ) ################################################### ### code chunk number 37: APd2 ################################################### par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, cbind(M.A,F.A), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(0.2,180), xlab="Age at diagnosis", ylab="Mortality rate at 2 years duration per 1000 PY" ) matplot( pr.P, cbind(M.P,F.P), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/3,3), xlab="Date of diagnosis", ylab="Mortality rate ratio" ) abline( h=1 ) matplot( pr.d, cbind(M.d,F.d), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/3,3), xlab="Diabetes duration", ylab="Mortality rate ratio" ) abline( h=1 ) ################################################### ### code chunk number 38: APd-cmp ################################################### par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, cbind(M.A,F.A,m.A,f.A), type="l", lty=rep(c(1,3),each=6), lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(0.2,180), xlab="Age at diagnosis/follow-up", ylab="Mortality rate at 2 years duration per 1000 PY" ) matplot( pr.P, cbind(M.P,F.P,m.P,f.P), type="l", lty=rep(c(1,3),each=6), lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/3,3), xlab="Date of diagnosis/follow-up", ylab="Mortality rate ratio" ) abline( h=1 ) matplot( pr.d, cbind(M.d,F.d,m.d,f.d), type="l", lty=rep(c(1,3),each=6), lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/3,3), xlab="Diabetes duration", ylab="Mortality rate ratio" ) abline( h=1 ) ################################################### ### code chunk number 39: DMDK-demo.rnw:745-754 ################################################### str( SL ) SL$Am <- floor( SL$A+0.5 ) SL$Pm <- floor( SL$P+0.5 ) 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 40: DMDK-demo.rnw:757-760 ################################################### SLr <- merge( SL, M.dk[,c("Am","Pm","sex","rate")] ) dim( SL ) dim( SLr ) ################################################### ### code chunk number 41: DMDK-demo.rnw:770-777 ################################################### stat.table( list( Age=floor(A/10)*10, Sex=sex ), list( D=sum(lex.Xst=="Dead"), E=sum(lex.dur*rate/1000), SMR=ratio(lex.Xst=="Dead",lex.dur*rate/1000) ), margins = TRUE, data = SLr ) ################################################### ### code chunk number 42: DMDK-demo.rnw:787-796 ################################################### SLr <- subset( SLr, rate>0) SLr$E <- SLr$lex.dur * SLr$rate / 1000 Sm <- glm( (lex.Xst=="Dead") ~ Ns( A-dur, kn=kn.Ad ) + Ns( P-dur, kn=kn.Pd ) + Ns( dur, kn=kn.dur ), offset = log( E ), family = poisson, data = subset( SLr, sex=="M" ) ) Sf <- update( Sm, data = subset( SLr, sex=="F" ) ) ################################################### ### code chunk number 43: DMDK-demo.rnw:800-806 ################################################### sM.A <- ci.exp( Sm, ctr.mat=cbind(1,AC,PR,dR) ) sF.A <- ci.exp( Sf, ctr.mat=cbind(1,AC,PR,dR) ) sM.P <- ci.exp( Sm, subset="P" , ctr.mat=PC-PR ) sF.P <- ci.exp( Sf, subset="P" , ctr.mat=PC-PR ) sM.d <- ci.exp( Sm, subset="kn.dur", ctr.mat=dC-dR ) sF.d <- ci.exp( Sf, subset="kn.dur", ctr.mat=dC-dR ) ################################################### ### code chunk number 44: SMR ################################################### par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, cbind(sM.A,sF.A), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/3,3), xlab="Age at follow-up", ylab="SMR" ) abline( h=1 ) matplot( pr.P, cbind(sM.P,sF.P), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/3,3), xlab="Date of follow-up", ylab="SMR ratio" ) abline( h=1 ) matplot( pr.d, cbind(sM.d,sF.d), type="l", lty=1, lwd=c(3,1,1), col=rep(c("blue","red"),each=3), log="y", ylim=c(1/3,3), xlab="Diabetes duration", ylab="SMR ratio" ) abline( h=1 ) ################################################### ### code chunk number 45: DMDK-demo.rnw:836-842 ################################################### Sb <- update( Sm, data = SLr ) Sb.s <- update( Sb, . ~. + sex ) Sb.i <- update( Sb, . ~. + sex:( Ns( A-dur, kn=kn.Ad ) + Ns( P-dur, kn=kn.Pd ) + Ns( dur, kn=kn.dur ) ) ) anova( Sb, Sb.s, Sb.i, test="Chisq" ) ################################################### ### code chunk number 46: SMR-b ################################################### Sb.A <- ci.exp( Sb, ctr.mat=cbind(1,AC,PR,dR) ) Sb.P <- ci.exp( Sb, subset="P" , ctr.mat=PC-PR ) Sb.d <- ci.exp( Sb, subset="kn.dur", ctr.mat=dC-dR ) par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, Sb.A, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(1/3,3), xlab="Age at diagnosis", ylab="SMR" ) abline( h=1 ) matplot( pr.P, Sb.P, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(1/3,3), xlab="Date of diagnosis", ylab="SMR ratio" ) abline( h=1 ) matplot( pr.d, Sb.d, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(1/3,3), xlab="Diabetes duration", ylab="SMR ratio" ) abline( h=1 ) ################################################### ### code chunk number 47: DMDK-demo.rnw:879-892 ################################################### kn.Ad <- with( subset( SL, lex.Xst=="Dead" ), quantile( A-dur, probs=seq(5,95,20)/100 ) ) kn.dur <- 0:2 AC <- Ns( pr.A, knots=kn.Ad ) dC <- Ns( pr.d, knots=kn.dur ) dR <- Ns( rep(rf.d,N), knots=kn.dur ) Sx <- glm( (lex.Xst=="Dead") ~ Ns( A-dur, kn=kn.Ad ) + I( P-dur ) + Ns( dur, kn=kn.dur ), offset = log( E ), family = poisson, data = SLr ) ################################################### ### code chunk number 48: SMR-x ################################################### Sx.A <- ci.exp( Sx, ctr.mat=cbind(1,AC,rf.P,dR) ) Sx.P <- ci.exp( Sx, subset="P" , ctr.mat=cbind(pr.P-rf.P) ) Sx.d <- ci.exp( Sx, subset="kn.dur", ctr.mat=dC-dR ) par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, Sx.A, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(1/2,4), xlab="Age at diagnosis", ylab="SMR" ) abline( h=1 ) abline( v=4:8*10, col="gray" ) matplot( pr.P, Sx.P, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(1/2,4), xlab="Date of diagnosis", ylab="SMR ratio" ) abline( h=1 ) matplot( pr.d, Sx.d, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(1/2,4), xlab="Diabetes duration", ylab="SMR ratio" ) abline( h=1,v=2 ) ################################################### ### code chunk number 49: DMDK-demo.rnw:923-924 ################################################### 100*( 1 - ci.exp( Sx, subset="P" ) ) ################################################### ### code chunk number 50: DMDK-demo.rnw:929-932 ################################################### d6 <- Ns( 6, knots=kn.dur ) d5 <- Ns( 5, knots=kn.dur ) 100*( ci.exp( Sx, subset="kn.dur", ctr.mat=d6-d5 ) - 1 ) ################################################### ### code chunk number 51: DMDK-demo.rnw:949-952 ################################################### Slx <- update( Sx, . ~. + I(A-dur):Ns(dur,knots=kn.dur) ) anova( Slx, Sx, test="Chisq" ) ci.exp( Slx ) ################################################### ### code chunk number 52: SMR-x ################################################### Slx.A <- ci.exp( Slx, ctr.mat=cbind(1,AC,rf.P,dR,dR*pr.A) ) Slx.P <- ci.exp( Slx, subset="P" , ctr.mat=cbind(pr.P-rf.P) ) Slx.d <- ci.exp( Slx, subset="kn.dur", ctr.mat=cbind(dC-dR,(dC-dR)*50) ) for( a in seq(55,90,5) ) Slx.d <- cbind( Slx.d, ci.exp( Slx, subset="kn.dur", ctr.mat=cbind(dC-dR,(dC-dR)*a) ) ) dim( Slx.d ) par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, Slx.A, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(1/2,4), xlab="Age at diagnosis", ylab="SMR" ) abline( h=1 ) abline( v=4:8*10, col="gray" ) matplot( pr.P, Slx.P, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(1/2,4), xlab="Date of diagnosis", ylab="SMR ratio" ) abline( h=1 ) matplot( pr.d, Slx.d, type="l", lty=1, lwd=c(3,1,1), col=rep(heat.colors(9),each=3), log="y", ylim=c(1/2,4), xlab="Diabetes duration", ylab="SMR ratio" ) abline( h=1 ) ################################################### ### code chunk number 53: DMDK-demo.rnw:990-998 ################################################### pts <- c(seq(0,15,0.1),NA) np <- length( pts ) nd <- data.frame( A=rep(seq(50,90,5),each=np)+pts, P=rf.P+pts, dur= pts, E=1 ) A.sx <- ci.pred( Sx , newdata=nd ) A.sl <- ci.pred( Slx, newdata=nd ) ################################################### ### code chunk number 54: SMR-xA ################################################### matplot( NA, NA, log="y", ylim=c(1/2,5), xlim=c(50,100), xlab="Age at follow-up", ylab="SMR" ) abline( h=c(5:19/10,seq(2,5,0.5)), v=seq(50,100,5), col=gray(0.8) ) matlines( nd$A, cbind(A.sx,A.sl), type="l", lty=rep(c(1,3),each=3), lwd=c(3,1,1), col="forestgreen" ) abline( h=1 ) ################################################### ### code chunk number 55: DMDK-demo.rnw:1024-1027 ################################################### Six <- update( Sx, . ~. + Ns(A-dur,knots=kn.Ad):Ns(dur,knots=kn.dur) ) anova( Six, Slx, Sx, test="Chisq" ) A.si <- ci.pred( Six, newdata=nd ) ################################################### ### code chunk number 56: SMR-xAx ################################################### matplot( NA, NA, log="y", ylim=c(1/2,5), xlim=c(50,100), xlab="Age at follow-up", ylab="SMR" ) abline( h=c(5:19/10,seq(2,5,0.5)), v=seq(50,100,5), col=gray(0.8) ) matlines( nd$A, cbind(A.si,A.sl,A.sx), type="l", lty=rep(c(1,3),c(6,3)), lwd=c(3,1,1), col=rep(c("magenta","limegreen"),c(3,6)) ) abline( h=1 )