### R code from vignette source 'DMDK-s.rnw' ################################################### ### code chunk number 1: DMDK-s.rnw:34-40 ################################################### options( width=120 ) library( Epi ) data( DMlate ) str( DMlate ) head( DMlate ) summary( DMlate ) ################################################### ### code chunk number 2: DMDK-s.rnw:64-66 ################################################### with( DMlate, table( dead=!is.na(dodth), same=(dodth==dox), exclude=NULL ) ) ################################################### ### code chunk number 3: DMDK-s.rnw:70-77 ################################################### 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: DMDK-s.rnw:81-82 ################################################### summary( LL ) ################################################### ### code chunk number 5: DMDK-s.rnw:91-96 ################################################### 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 6: DMDK-s.rnw:123-125 ################################################### SL <- splitLexis( LL, breaks=seq(0,125,1), time.scale="A" ) summary( SL ) ################################################### ### code chunk number 7: DMDK-s.rnw:164-170 ################################################### 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 8: DMDK-s.rnw:174-183 ################################################### nd <- data.frame( A = seq(10,90,0.5), lex.dur = 1000) p.m <- predict.glm( r.m, type = "link", newdata = nd, se.fit = TRUE ) p.f <- predict.glm( r.f, type = "link", newdata = nd, se.fit = TRUE ) str( p.m ) ################################################### ### code chunk number 9: DMDK-s.rnw:188-191 ################################################### ci.mat() lr.m <- cbind(p.m$fit,p.m$se.fit) %*% ci.mat() lr.f <- cbind(p.f$fit,p.f$se.fit) %*% ci.mat() ################################################### ### code chunk number 10: a-rates ################################################### matplot( seq(10,90,0.5), exp( cbind(lr.m,lr.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-s.rnw:223-225 ################################################### data( M.dk ) head( M.dk ) ################################################### ### code chunk number 12: DMDK-s.rnw:228-230 (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), exp( cbind(lr.m,lr.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-s.rnw:264-273 ################################################### 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 <- predict.glm( R.m, type = "link", newdata = nd ) P.f <- predict.glm( R.f, type = "link", newdata = nd ) ################################################### ### code chunk number 15: a-rates-s ################################################### matplot( seq(10,90,0.5), exp( cbind(lr.m,lr.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), exp(cbind( P.m,P.f )), lty="12", col=c("blue","red"), lwd=3 ) ################################################### ### code chunk number 16: DMDK-s.rnw:336-342 ################################################### 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-s.rnw:379-389 ################################################### mm <- glm( (lex.Xst=="Dead") ~ ns( A, kn=kn.A[-c(1,length(kn.A))], Bo=kn.A[ c(1,length(kn.A))] ) + ns( P, kn=kn.P[-c(1,length(kn.P))], Bo=kn.P[ c(1,length(kn.P))] ) + ns( dur, kn=kn.dur[-c(1,length(kn.dur))], Bo=kn.dur[ c(1,length(kn.dur))] ), offset = log( lex.dur ), family = poisson, data = subset( SL, sex=="M" ) ) summary( mm ) ################################################### ### code chunk number 18: DMDK-s.rnw:395-397 ################################################### source( "http://BendixCarstensen.com/SPE/R/Ns.r" ) Ns ################################################### ### code chunk number 19: DMDK-s.rnw:400-408 ################################################### 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 20: DMDK-s.rnw:422-424 ################################################### anova( mm, r.m, test="Chisq" ) anova( mf, r.f, test="Chisq" ) ################################################### ### code chunk number 21: DMDK-s.rnw:528-534 ################################################### 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 22: DMDK-s.rnw:538-543 ################################################### 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 23: DMDK-s.rnw:560-566 ################################################### m.A <- ci.exp( mm, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 m.P <- ci.exp( mm, subset="P" , ctr.mat=PC-PR ) m.d <- ci.exp( mm, subset="dur", ctr.mat=dC-dR ) f.A <- ci.exp( mf, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 f.P <- ci.exp( mf, subset="P" , ctr.mat=PC-PR ) f.d <- ci.exp( mf, subset="dur", ctr.mat=dC-dR ) ################################################### ### code chunk number 24: 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 25: 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 m.P <- ci.exp( mm, subset="P" , ctr.mat=PC-PR ) m.d <- ci.exp( mm, subset="dur", ctr.mat=dC-dR ) f.A <- ci.exp( mf, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 f.P <- ci.exp( mf, subset="P" , ctr.mat=PC-PR ) 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 26: 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 27: DMDK-s.rnw:692-700 ################################################### 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 28: DMDK-s.rnw:719-722 ################################################### j.dev <- mm$dev + mf$dev j.df <- mm$df.r + mf$df.r 1 - pchisq( m2$dev - j.dev, m2$df.r - j.df ) ################################################### ### code chunk number 29: DMDK-s.rnw:750-752 ################################################### ci.exp( m2, subset=c("Int","sexM","P","dur") ) ci.exp( m2, subset=c("Int","sexF","P","dur") ) ################################################### ### code chunk number 30: DMDK-s.rnw:756-760 ################################################### 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 31: 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 32: DMDK-s.rnw:824-830 ################################################### pts <- seq(0,20,1) nd <- data.frame( A= 50+pts, P=1995+pts, dur= pts, lex.dur=1000 ) predict( mm, newdata=nd, se.fit=TRUE ) ################################################### ### code chunk number 33: DMDK-s.rnw:833-834 ################################################### sapply(predict( mm, newdata=nd, se.fit=TRUE )[1:2],cbind) %*% ci.mat() ################################################### ### code chunk number 34: DMDK-s.rnw:839-852 ################################################### DMm <- function( A, P, pts=seq(0,15,0.1) ) { nd <- data.frame( A=A+pts, P=P+pts, dur= pts, lex.dur=1000 ) cbind( nd$A, exp(sapply(predict( mm, newdata=nd, se.fit=TRUE )[1:2],cbind) %*% ci.mat()), exp(sapply(predict( mf, newdata=nd, se.fit=TRUE )[1:2],cbind) %*% ci.mat()) ) } DMm( 50, 1996, pts=0:10 ) ################################################### ### code chunk number 35: 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 36: DMDK-s.rnw:941-961 ################################################### 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 ) ) 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 37: DMDK-s.rnw:980-998 ################################################### 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 M.P <- ci.exp( Mm, subset="P" , ctr.mat=PC-PR ) M.d <- ci.exp( Mm, subset="kn.dur", ctr.mat=dC-dR ) F.A <- ci.exp( Mf, ctr.mat=cbind(1,AC,PR,dR) ) * 1000 F.P <- ci.exp( Mf, subset="P" , ctr.mat=PC-PR ) F.d <- ci.exp( Mf, subset="kn.dur", ctr.mat=dC-dR ) ################################################### ### code chunk number 38: 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 39: 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 40: DMDK-s.rnw:1101-1110 ################################################### 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 41: DMDK-s.rnw:1113-1116 ################################################### SLr <- merge( SL, M.dk[,c("Am","Pm","sex","rate")] ) dim( SL ) dim( SLr ) ################################################### ### code chunk number 42: DMDK-s.rnw:1128-1135 ################################################### 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 43: DMDK-s.rnw:1153-1162 ################################################### 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 44: DMDK-s.rnw:1166-1172 ################################################### sM.A <- ci.exp( Sm, ctr.mat=cbind(1,AC,PR,dR) ) sM.P <- ci.exp( Sm, subset="P" , ctr.mat=PC-PR ) sM.d <- ci.exp( Sm, subset="kn.dur", ctr.mat=dC-dR ) sF.A <- ci.exp( Sf, ctr.mat=cbind(1,AC,PR,dR) ) sF.P <- ci.exp( Sf, subset="P" , ctr.mat=PC-PR ) sF.d <- ci.exp( Sf, subset="kn.dur", ctr.mat=dC-dR ) ################################################### ### code chunk number 45: 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 46: DMDK-s.rnw:1204-1210 ################################################### 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 47: 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 48: DMDK-s.rnw:1256-1269 ################################################### 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 49: 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 50: DMDK-s.rnw:1303-1304 ################################################### 100*( 1 - ci.exp( Sx, subset="P" ) ) ################################################### ### code chunk number 51: DMDK-s.rnw:1311-1314 ################################################### 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 52: DMDK-s.rnw:1330-1333 ################################################### Six <- update( Sx, . ~. + I(A-dur):Ns(dur,knots=kn.dur) ) anova( Six, Sx, test="Chisq" ) ci.exp( Six ) ################################################### ### code chunk number 53: SMR-x ################################################### Six.A <- ci.exp( Six, ctr.mat=cbind(1,AC,rf.P,dR,dR*pr.A) ) Six.P <- ci.exp( Six, subset="P" , ctr.mat=cbind(pr.P-rf.P) ) Six.d <- ci.exp( Six, subset="kn.dur", ctr.mat=cbind(dC-dR,(dC-dR)*50) ) for( a in seq(55,90,5) ) Six.d <- cbind( Six.d, ci.exp( Six, subset="kn.dur", ctr.mat=cbind(dC-dR,(dC-dR)*a) ) ) dim( Six.d ) par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( pr.A, Six.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, Six.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, Six.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 54: DMDK-s.rnw:1371-1379 ################################################### 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.si <- exp(sapply(predict( Six, newdata=nd, se.fit=TRUE )[1:2],cbind) %*% ci.mat()) A.sm <- exp(sapply(predict( Sx , newdata=nd, se.fit=TRUE )[1:2],cbind) %*% ci.mat()) ################################################### ### code chunk number 55: 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.si,A.sm), type="l", lty=rep(c(1,3),each=3), lwd=c(3,1,1), col="forestgreen" ) abline( h=1 ) ################################################### ### code chunk number 56: DMDK-s.rnw:1404-1406 ################################################### SiX <- update( Sx, . ~. + Ns(A-dur,knots=kn.Ad):Ns(dur,knots=kn.dur) ) anova( SiX, Six, Sx, test="Chisq" ) ################################################### ### code chunk number 57: SMR-xAx ################################################### A.sX <- exp(sapply(predict( SiX, newdata=nd, se.fit=TRUE )[1:2],cbind) %*% ci.mat()) 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.si,A.sm), type="l", lty=rep(c(1,3),c(6,3)), lwd=c(3,1,1), col=rep(c("magenta","forestgreen"),c(3,6)) ) abline( h=1 )