### R code from vignette source 'pop-PY.rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: startup ################################################### options( width=90 ) library( Epi ) library( splines ) print( sessionInfo(), loc=F ) clear() ## system.time( ## load( file=url("http://bendixcarstensen.com/DMreg/DMTB/TubDM.Rda") ) ## ) load( file="../data/TubDM.Rda") Ecol <- ecol bwcol <- gray(0:3/4) names(bwcol) <- names(ecol) ################################################### ### code chunk number 2: list-Agg ################################################### # load( file="../data/Agg.Rda" ) str(Agg) round( ftable( xtabs( Y/1000 ~ region + state, data = Agg ), row.vars=c(1) ) , 1 ) ftable( xtabs( cbind( D.tb, D.dm, D.dd ) ~ region + state, data = Agg ), row.vars=c(3,1) ) ################################################### ### code chunk number 3: non-DKWell-agg ################################################### system.time( Cgg <- with( subset( Agg, A<100 & P>1994 & P<2010 & !(region=="" & state=="Well") ), aggregate( cbind( X = Y ), list( A = A, P = P, upper = U, sex = sex ), FUN = sum ) ) ) str( Cgg ) summary( Cgg ) ################################################### ### code chunk number 4: pop-Y ################################################### data( Y.dk ) Y.dk$sex <- factor( Y.dk$sex, labels=c("M","F") ) Y.dk <- subset( Y.dk, A<100 & P>1994 & P<2010, select=c("sex","A","P","upper","Y") ) ################################################### ### code chunk number 5: make-DKWell-PY ################################################### Y.rev <- merge( Cgg, Y.dk, all.y=TRUE ) summary( Y.rev ) Y.rev <- transform( Y.rev, Y.pop = Y-pmax(X,0,na.rm=TRUE), state = "Well", region = "", U = upper )[,c("A","P","U","sex","state","region","Y.pop")] str( Y.rev ) ################################################### ### code chunk number 6: merge-DKWell-Y ################################################### Afu <- merge( subset( Agg, A<100 & P>1994 & P<2010 ), Y.rev, all=TRUE ) Afu <- transform( Afu, Y = pmax( Y,Y.pop,na.rm=TRUE), D.tb = pmax(D.tb, 0,na.rm=TRUE), D.dm = pmax(D.dm, 0,na.rm=TRUE), D.dd = pmax(D.dd, 0,na.rm=TRUE) )[, c("sex","A","P","U","state","region","Y","D.tb","D.dm","D.dd")] str( Afu ) ################################################### ### code chunk number 7: showFU ################################################### round( addmargins( xtabs( Y ~ region + state, data=Afu )/1000 ), 1 ) ftable( addmargins( xtabs( cbind(D.tb,D.dm) ~ region + state, data=Afu ), margin = 1:2 ), row.vars=c(3,1) ) save( Afu, file="../data/Afu.Rda" ) ################################################### ### code chunk number 8: get-Dur-dat ################################################### # load( file="../data/Dgg.Rda" ) with( Dgg, table( state, D.tb ) ) with( Dgg, table( state, D.dm ) ) with( Dgg, table( state, D.dd ) ) ################################################### ### code chunk number 9: append-nonDM ################################################### Dfu <- rbind( subset( Dgg, A<100 ), cbind( subset( Afu, state=="Well" ), dur=NA ) ) str( Dfu ) ################################################### ### code chunk number 10: show-Dur-dat ################################################### round( addmargins( xtabs( Y ~ region + state, data=Dfu ) )/1000, 1 ) ftable( addmargins( xtabs( cbind(D.tb,D.dm,D.dd) ~ region + state, data=Dfu ), margin = 1:2 ), row.vars=c(3,1) ) save( Dfu, file="../data/Dfu.Rda" ) ################################################### ### code chunk number 11: get-xLx ################################################### # load( file="../data/Afu.Rda" ) addmargins( xtabs( D.tb ~ region + state, data=Dfu ) ) formatC( at <- xtabs( cbind( D.tb, D.dm, D.dd, Y ) ~ state, data=Afu ), format="f", digits=0, big.mark=",", preserve.width=NULL ) # load( file="../data/xLx.Rda" ) xLx <- subset( xLx, age-lex.dur <= 100 ) ( tt <- tmat( xLx , Y=TRUE ) ) ( ti <- tmat( subset(xLx,region!=""), Y=TRUE ) ) tt["Well","Well"] <- at["Well","Y"] ( td <- abs( tt - ti ) ) td["Well","Dead"] <- 0 tt["Well","Dead"] <- 0 ################################################### ### code chunk number 12: boxes-final ################################################### aclr <- rep("black",9) aclr[5] <- "red" aclr[3] <- "forestgreen" ## tmpl <- function(){ ## par( mfrow=c(2,1) ) ## boxes.Lexis( td, ## boxpos=list( x=c(10,90,10,50,50,90), ## y=c(65,35,35,90,10,65) ), ## hmult=1.5, col.arr=aclr, ## show=TRUE, scale.Y=1000, digits.R=2 ) ## text( 3, 95, "Danish born", adj=c(0,1), font=2, cex=1.5 ) ## boxes.Lexis( ti, ## boxpos=list( x=c(10,90,10,50,50,90), ## y=c(65,35,35,90,10,65) ), ## hmult=1.5, col.arr=aclr, ## show=TRUE, scale.Y=1000 ) ## text( 3, 95, "Foreign born", adj=c(0,1), font=2, cex=1.5 ) ## } ## tmpl() tmpl <- function(){ par( mfrow=c(2,1) ) bx.td <- boxes.Lexis( td, boxpos=list( x=c(10,90,10,50,50,90), y=c(75,35,35,90,10,75) ), hmult=1.5, col.arr=aclr, #DR.sep=c(" (",")"), show=TRUE, scale.Y=1000, digits.R=2 ) text( 3, 95, "Danish born", adj=c(0,1), font=2, cex=1.5 ) bx.ti <- boxes.Lexis( ti, boxpos=list( x=c(10,90,10,50,50,90), y=c(75,35,35,90,10,75) ), hmult=1.5, col.arr=aclr, #DR.sep=c(" (",")"), show=TRUE, scale.Y=1000 ) text( 3, 95, "Foreign born", adj=c(0,1), font=2, cex=1.5 ) list( bx.td, bx.ti ) } # bfig1 <- tmpl() for( i in 1:2 ) { bfig1[[i]]$State.names <- gsub("\\("," \\(",bfig1[[i]]$State.names ) bfig1[[i]]$State.names <- gsub( ",", " ",bfig1[[i]]$State.names ) bfig1[[i]]$Boxes$wd <- 17 bfig1[[i]]$Boxes$ht <- 12 bfig1[[i]]$Boxes$font <- 1 bfig1[[i]]$Arrows$font.arr <- 1 bfig1[[i]]$Arrows$pos.arr[c(2,9)] <- 0.60 } bfig1[[1]]$State.names[3] <- gsub(" ","" ,bfig1[[1]]$State.names )[3] bfig1[[2]]$State.names[1] <- gsub(" ","" ,bfig1[[2]]$State.names )[1] bfig1[[1]]$Arrowtext[ c(2,4)] <- gsub(","," ",bfig1[[1]]$Arrowtext )[ c(2,4)] bfig1[[1]]$Arrowtext[-c(2,4)] <- gsub(",", "",bfig1[[1]]$Arrowtext )[-c(2,4)] bfig1[[2]]$Arrowtext[ 2] <- gsub(","," ",bfig1[[2]]$Arrowtext )[ 2] bfig1[[2]]$Arrowtext[-2] <- gsub(",", "",bfig1[[2]]$Arrowtext )[-2] ################################################### ### code chunk number 13: Fig1 ################################################### pdf( "IJTLD-Fig1.pdf", height=13, width=9 ) par( mfrow=c(2,1) ) boxes(bfig1[[1]]) boxes(bfig1[[2]]) dev.off() for( i in 1:2 ) { bfig1[[i]]$Boxes$col.txt <- bfig1[[i]]$Boxes$col.border <- c("black",gray(0.6))[c(1,2,1,1,1,2)] bfig1[[i]]$Arrows$col.txt.arr <- bfig1[[i]]$Arrows$col.arr <- c("black",gray(0.6))[c(2,2,1,2,1,2,2,2,2)] } pdf( "IJTLD-Fig1-bw.pdf", height=13, width=9 ) par( mfrow=c(2,1) ) boxes(bfig1[[1]]) boxes(bfig1[[2]]) dev.off() # tiff( "IJTLD-Fig1.tif", height=14, width=10, units="in", res=500 ) # tmpl() # dev.off() # tiff( "IJTLD-Fig1-bw.tif", height=14, width=10, units="in", res=500 ) # tmpl() # dev.off() ################################################### ### code chunk number 2: read-Afu ################################################### # load( file="../data/Afu.Rda" ) str( Afu ) ################################################### ### code chunk number 3: libs ################################################### cnr <- function( xf, yf ) { # A function that gives the coordinates of the # point (xf,yf) from ll corner in the current plot. # if xf or yf are > 1 they are considered percentages # cn <- par()$usr xf <- ifelse( xf>1, xf/100, xf ) yf <- ifelse( yf>1, yf/100, yf ) xx <- ( 1 - xf ) * cn[1] + xf * cn[2] yy <- ( 1 - yf ) * cn[3] + yf * cn[4] if ( par()$xlog ) xx <- 10^xx if ( par()$ylog ) yy <- 10^yy list( x=xx, y=yy ) } rect <- function( x1, y1, x2, y2, ... ) { # A wrapper for rect, that allows the posiiton of the rectangle to be # passed as either two vectors of length 2 or a list of two such # vectors if( is.list(x1) ) { y1 <- x1$y x1 <- x1$x } if( length(x1) > 1 & length(y1) > 1 ) graphics::rect( x1[1], y1[1], x1[2], y1[2], ... ) else graphics::rect( x1 , y1 , x2 , y2 , ... ) } ################################################### ### code chunk number 4: tr2Atb ################################################### Atb <- transform( subset( Afu, state %in% c("Well","DM") ), state = factor(state), Region = Relevel( region, list(Asia=c(4,7),Other=c(3,5,6,8)), first=FALSE ), ax = A+(1+U)/3, px = P+(2-U)/3 ) str( Atb ) levels( Atb$region )[1] <- levels( Atb$Region )[1] <- "DK" ( atab <- addmargins( xtabs( D.tb ~ region + state, data=Atb ) ) ) ( aTab <- addmargins( xtabs( D.tb ~ Region + state, data=Atb ) ) ) atab <- rbind(atab,aTab[4,])[c(1,2,4,10,6,5,3,8,9),1:2] rownames( atab )[4] <- "Remain" atab save( atab, file="atab.Rda" ) ################################################### ### code chunk number 5: TBtabs ################################################### addmargins( xtabs( D.tb ~ region + state, data=Atb ) ) ################################################### ### code chunk number 6: TBtab-grp ################################################### Atb$Region <- Relevel( Atb$region, list(Asia=c(4,7),Other=c(3,5,6,8)), first=FALSE ) formatC( ftable( addmargins( xtabs( cbind( D.tb, Y=Y/1000 ) ~ sex + Region + state, data=Atb ), margin = 1:2 ), row.vars=2:1, col.vars=4:3 ), format="f", digits=1, big.mark=",", pr="c" ) ################################################### ### code chunk number 7: ap-knots ################################################### nk <- 4 ( a.kn <- with( Atb, quantile( rep(ax,D.tb), (1:nk-0.5)/nk ) ) ) nk <- 3 ( p.kn <- with( Atb, quantile( rep(px,D.tb), (1:nk-0.5)/nk ) ) ) ################################################### ### code chunk number 8: ecol ################################################### scol <- c("blue","red") names(scol) <- levels( Atb$sex ) ecol <- c("black","orange","magenta","forestgreen") names(ecol) <- levels( Atb$Region ) c( ecol, scol ) save( ecol, scol, file="../data/clrs.Rda" ) par( mar=c(0,0,0,0) ) plot(1:4,1:4,axes=F, xlab="", ylab="",xlim=c(0,5),ylim=c(0,5),type="n") text(rep(1,4),4:1, names(ecol), col=ecol, font=2, cex=2,adj=0 ) text(rep(4,2),3:2, names(scol), col=scol, font=2, cex=2,adj=0 ) ################################################### ### code chunk number 9: main models ################################################### m1 <- glm( D.tb ~ Ns( ax, kn=a.kn ) + Ns( px, kn=p.kn ) + sex + state, offset = log(Y/10^5), family = poisson, data = Atb ) m2 <- update( m1, . ~ . + Region ) m3 <- update( m2, . ~ . - state + Region:state ) anova( m3, m2, m1, test="Chisq" ) round( ci.exp( m1 ), 3 ) round( ci.exp( m2 ), 3 ) round( ci.exp( m3 ), 3 ) ################################################### ### code chunk number 10: cmats ################################################### n.pt <- 200 a.pt <- seq(0,90,,n.pt) Ca <- Ns( a.pt, kn=a.kn ) p.pt <- seq(1995,2010,,n.pt) Cp <- Ns( p.pt, kn=p.kn ) p.ref <- 2005 Cpr <- Ns( rep(p.ref,n.pt), kn=p.kn ) ################################################### ### code chunk number 11: age-eff ################################################### m.eff <- ci.exp( m3, ctr.mat=cbind(1,Ca,Cpr), subset=c("Int","ax","px") ) f.eff <- ci.exp( m3, ctr.mat=cbind(1,Ca,Cpr,1), subset=c("Int","ax","px","sex") ) tmpl <- function(){ par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( a.pt, m.eff, type="n", log="y",las=1, ylim=c(0.7,7), xlab="Age at follow-up", ylab="TB incidence rate per 100,000 PY in 2005") abline( v=seq(0,90,5), h=c(5:15/10,2:10), col=gray(0.8) ) matlines( a.pt, cbind( m.eff, f.eff), lty=1, lwd=c(3,1,1), col=rep(scol,each=3) ) rect( cnr(c(0,10),c(85,100)), col="white", border=gray(0.8) ) text( cnr(rep(5,2),c(95,90)), levels(Atb$sex), col=scol, cex=1.1 ) box() } tmpl() ################################################### ### code chunk number 12: age-eff-emf ################################################### ################################################### ### code chunk number 13: per-eff ################################################### p.rr <- ci.exp( m3, ctr.mat=Cp-Cpr, subset="px" ) tmpl <- function() { par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( p.pt, p.rr, type="n", log="y",las=1, ylim=c(0.5,2), xlab="Date of follow-up", ylab="RR of TB") abline( v=1995:2010, h=c(1:15/10,2:10), col=gray(0.8) ) abline( h=1 ) matlines( p.pt, p.rr, lty=1, lwd=c(3,1,1), col="black" ) points( c(p.ref,p.ref), c(1,1), cex=1.3, pch=c(16,1), lwd=3, col=c("white","black") ) } tmpl() ################################################### ### code chunk number 14: per-eff-emf ################################################### ################################################### ### code chunk number 15: TB-ana.rnw:196-199 ################################################### m4 <- update( m3, . ~ . - Ns(px,kn=p.kn) + px ) anova( m3, m4, test="Chisq" ) round((1-ci.exp( m4, subset="px" )[,c(1,3,2)])*100,2) ################################################### ### code chunk number 16: TB-ana.rnw:203-204 ################################################### round((ci.exp( m4, subset="px", ctr.mat=matrix(15,1,1) ))*100,2) ################################################### ### code chunk number 17: cim3 ################################################### round( ci.exp( m3, subset="Region" ), 2 ) ################################################### ### code chunk number 18: get-fixeff ################################################### round( e1 <- ci.exp( m1, subset="state" ), 3 ) round( e2 <- ci.exp( m2, subset="state" ), 3 ) round( e3 <- ci.exp( m3, subset="state" ), 3 ) rownames( e3 ) <- gsub( "Region","", gsub( ":stateDM", "", rownames( e3 ) ) ) ee <- rbind( e1, e2, e3 ) rownames( ee )[1:2] <- c("Raw","Region-adj") round( ee, 2 ) ################################################### ### code chunk number 19: DMRR-est ################################################### par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) irr <- function(){ plotEst( ee[-(1:2),], lwd=2, vref=1, cex=1.1, grid=c(3:15/10,2,2.5,3,4), xtic=c(0.3,0.5,1,2,3,4), xlab="TB RR: DM vs. non-DM", xlog=TRUE, col=ecol, y=4:1 )} irr() ################################################### ### code chunk number 20: all-int ################################################### # Set up the relevant contrast matrix nr <- nlevels( Atb$Region ) CRR <- diag( 2*nr ) CRR[nr+1:nr,1:nr] <- diag(nr) CRR[1,] <- 0 CRR <- CRR[,-1] rownames(CRR) <- t( outer( c("","DM "), levels(Atb$Region), paste, sep="" ) ) CRR ci.exp( m3, subset="Region" ) round( e3 <- ci.exp( m3, subset=c("Region"), ctr.mat=CRR ), 2 ) rownames( e3 )[nr+1:nr] <- NA par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) arr <- function(){ plotEst( e3, y=c(nr:1+0.1,nr:1-0.1), txtpos=rep(4:1,2), lwd=2, vref=1, cex=1.1, grid=outer(1:9,10^(0:2)), xtic=c(1,2,5,10,20,50,80,100),xlim=c(0.9,100), xlab="TB RR vs. non-DM, DK born", xlog=TRUE, col=c(rep(gray(0.5),4),ecol) )} arr() ################################################### ### code chunk number 21: all-int-emf ################################################### ################################################### ### code chunk number 22: all-int ################################################### # PLot the two sets of RRs next to each other # First x-axis is from 0.3->4, the second from 1->100 # The second is actually plotted from 10->1000 rr2 <- function(){ par( mar=c(3,1,2,1), mgp=c(3,1,0)/1.6 ) plotEst( ee[-(1:2),], ylim=c(0,4.1), xlim=c(0.3,1000), xtic=c(0.3,0.5,1,2,3,4), lwd=2, vref=1, cex=1.1, #grid=outer(1:9,10^(0:2)), xlog=TRUE, xlab="", grid=c(3:9/10,1.5,2:4), col=ecol ) abline( v=c(1:9*10,1:9*100,1000,15,150), col=gray(0.9) ) abline( v=10 ) axis( side=1, at=c(1,2,5,c(1,2,5)*10,80,100)*10, labels=c(1,2,5,c(1,2,5)*10,80,100) ) linesEst( e3*10, y=c(nr:1+0.1,nr:1-0.1), txtpos=rep(4:1,2), lwd=2, col=c(rep(gray(0.5),4),ecol) ) mtext( "TB RR, DM vs. non-DM", side=1, at=1.0, line=2 ) mtext( "TB RR vs. DK born non-DM", side=1, at=100, line=2 ) mtext( c("A)","B)"), at=c(0.3,10)*1.1, side=3, line=1, font=2 ) } rr2() ################################################### ### code chunk number 23: TB-ana.rnw:301-310 ################################################### pdf("IJTLD-Fig2.pdf",height=2.2,width=7.5) rr2() dev.off() ecol <- bwcol pdf("IJTLD-Fig2-bw.pdf",height=2.2,width=7.5) rr2() dev.off() ecol <- Ecol ################################################### ### code chunk number 24: inter-age ################################################### mi <- update( m3 , . ~ . - Ns(px,knots=p.kn) + I(px-2005) ) mi.s <- update( mi , . ~ . + sex:Ns(ax,knots=a.kn) ) mi.sr <- update( mi.s , . ~ . + Region:Ns(ax,knots=a.kn) ) mi.srs <- update( mi.sr , . ~ . + state:Ns(ax,knots=a.kn) ) mi.Srs <- update( mi.srs , . ~ . + sex:I(px-2005) ) mi.SRs <- update( mi.Srs , . ~ . + Region:I(px-2005) ) mi.SRS <- update( mi.SRs , . ~ . + state:I(px-2005) ) it <- as.matrix( anova( mi, mi.s , mi.sr , mi.srs, mi.Srs, mi.SRs, mi.SRS, test="Chisq" ) )[-1,3:5] rownames( it ) <- c( "sex:age", "Region:age", "DM:age", "sex:per", "Region:per", "DM:per" ) round( it, 3 ) ################################################### ### code chunk number 25: TB-ana.rnw:361-362 ################################################### round( ci.exp( mi.SRS ), 3 ) ################################################### ### code chunk number 26: TB-ana.rnw:366-369 ################################################### m.int <- update( mi.SRS , . ~ . - state:Ns(ax,knots=a.kn) - sex:I(px-2005) ) round( ci.exp(m.int), 3 ) ################################################### ### code chunk number 27: age-rates ################################################### pp <- Atb[1:length(a.pt),c("sex","ax","px","Region","state")] str(pp) pp <- transform( pp, sex = "M", ax = a.pt, px = 2005, Region = "DK", state = "Well", Y = 10^5 ) head(pp) ptb <- function( pp ) { exp( do.call( cbind, predict( m.int, newdata=pp, type="link", se.fit=TRUE )[1:2] ) %*% ci.mat() ) } r.dk <- ptb( transform(pp,state="Well",Region="DK" ) ) d.dk <- ptb( transform(pp,state="DM" ,Region="DK" ) ) r.af <- ptb( transform(pp,state="Well",Region="Africa") ) d.af <- ptb( transform(pp,state="DM" ,Region="Africa") ) r.as <- ptb( transform(pp,state="Well",Region="Asia" ) ) d.as <- ptb( transform(pp,state="DM" ,Region="Asia" ) ) r.ot <- ptb( transform(pp,state="Well",Region="Other" ) ) d.ot <- ptb( transform(pp,state="DM" ,Region="Other" ) ) tmpf <- function(rrcol="red"){ matplot( a.pt, cbind(r.dk,r.af,r.as,r.ot), log="y", xlab="Age (years)", ylab="", ylim=c(0.5,1800), xlim=c(10,85), las=1, yaxt="n", col="transparent" ) axis( side=2, at=outer(c(1,2,5),10^(-1:4)), labels=formatC( outer(c(1,2,5),10^(-1:4)), format="f", digits=1, drop=T ), las=1 ) mtext( "TB incidence per 100 000 PY in 2005", side=2, line=2.5 ) abline( v=1:9*10, h=outer(c(1.5,1:9),10^c(-1:4)), col=gray(0.9) ) matlines( a.pt, cbind(r.dk,r.af,r.as,r.ot), type="l", lty=1, lwd=c(4,1,1), col=rep(ecol,each=3) ) matlines( a.pt, cbind(d.dk,d.af,d.as,d.ot), type="l", lty=rep(c("22","99"),c(1,2)), lwd=c(4,1,1), lend=2, col=rep(ecol,each=3) ) text( cnr(2,101-c(4,1:3)*3), levels(Atb$Region), col=ecol, font=1, adj=c(0,1) ) abline( h=1 ) axis( side=4, at=c(5,7,10,15)/10, las=1 ) matlines( a.pt, ci.exp( m.int, subset="sex", ctr.mat=cbind(1,Ca) ), type="l", lty=1, lwd=c(4,1,1), col=rrcol ) mtext( "Female/Male RR", side=4, at=1, line=2.0, col=rrcol ) box() } par( mar=c(3,4,1,4), mgp=c(3,1,0)/1.6 ) tmpf() ################################################### ### code chunk number 28: TB-ana.rnw:435-447 ################################################### pdf( "IJTLD-Fig3.pdf", width=7, height=7, pointsize=16 ) par( mar=c(3,4,1,4), mgp=c(3,1,0)/1.6 ) tmpf() dev.off() ecol <- bwcol pdf( "IJTLD-Fig3-bw.pdf", width=7, height=7, pointsize=16 ) par( mar=c(3,4,1,4), mgp=c(3,1,0)/1.6 ) tmpf(rrcol="black") dev.off() ecol <- Ecol ################################################### ### code chunk number 29: TB-ana.rnw:460-461 ################################################### round( cfp <- ci.exp( m.int, subset="px" ), 3 ) ################################################### ### code chunk number 30: TB-ana.rnw:465-475 ################################################### rn <- outer( levels(Atb$state), levels(Atb$Region), paste ) CM <- cbind( rep(1,8), rep(c(0,1,0,0),each=2), rep(c(0,0,1,0),each=2), rep(c(0,0,0,1),each=2), rep(0:1,4) ) rownames( CM ) <- rn colnames( CM ) <- rownames( cfp ) CM round( 100*(1-ci.exp( m.int, subset="px", ctr.mat=CM ))[,c(1,3,2)], 1 ) ### R code from vignette source 'TB-dur.rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: TB-dur.rnw:2-8 ################################################### ################################################### ### code chunk number 2: load dfu ################################################### # load( file="../data/Dfu.Rda" ) str( Dfu ) with( Dfu, ftable( state, durOK=!is.na(dur), D.tb, row.vars=1:2 ) ) ################################################### ### code chunk number 3: trans2dtb ################################################### Dtb <- transform( subset( Dfu, state %in% c("Well","DM") ), state = factor(state), Region = Relevel( region, list(Asia=c(4,7),Other=c(3,5,6,8)), first=FALSE ), ax = A+(1+U)/3, px = P+(2-U)/3, dur = pmax( Dfu$dur, 0, na.rm=TRUE ) ) ################################################### ### code chunk number 4: TB-dur.rnw:34-43 ################################################### ( dtab <- addmargins( xtabs( D.tb ~ region + state, data=Dtb ) ) ) ( dTab <- addmargins( xtabs( D.tb ~ Region + state, data=Dtb ) ) ) dtab <- rbind(dtab,dTab[4,])[c(1,2,4,10,6,5,3,8,9),1:2] rownames( dtab )[c(1,4)] <- c("DK","Remain") colnames( dtab )[2] <- "DM(dur)" # load( file="atab.Rda" ) cbind( atab, dtab ) ( tt <- cbind( atab, dtab )[,-3] ) cbind( tt, round( sweep( tt, 2, tt[9,]/100, "/" ), 1 ) )[,c(1,4,2,5,3,6)] ################################################### ### code chunk number 5: knots-def ################################################### nk <- 4 p.dur <- seq(0,16,,100) ( d.kn <- with( subset(Dtb,state=="DM"), c( 0, quantile( rep(dur,D.tb), 1:nk/(nk+0.5) ) ) ) ) ( dmd <- with( subset(Dtb,state=="DM"), median( rep(dur,D.tb) ) ) ) ################################################### ### code chunk number 6: dur-splines ################################################### matplot( p.dur, Ns(p.dur,knots=d.kn), type="l", lwd=3, lty=1 ) rug( d.kn, lwd=3 ) ################################################### ### code chunk number 7: knots-def ################################################### nk <- 4 ( a.kn <- with( Dtb, quantile( rep(ax,D.tb), (1:nk-0.5)/nk ) ) ) nk <- 3 ( p.kn <- with( Dtb, quantile( rep(px,D.tb), (1:nk-0.5)/nk ) ) ) ################################################### ### code chunk number 8: md ################################################### md <- glm( D.tb ~ Ns( ax , kn=a.kn ) + sex + Region*state + Ns( dur, kn=d.kn ) + Region:I(px-2005), offset = log(Y/10^5), family = poisson, data = Dtb ) summary( md ) round( ci.exp( md ), 3 ) ################################################### ### code chunk number 9: dur1 ################################################### # load( file="../data/clrs.Rda" ) d.pt <- seq(0,15,,200) Cd <- Ns( d.pt, kn=d.kn ) round( ci.exp( md, subset=c("dur","DM") ), 3 ) d.eff <- ci.exp( md, subset=c("dur","DM"), ctr.mat=cbind(Cd,1,0,0,0) ) d.afr <- ci.exp( md, subset=c("dur","DM"), ctr.mat=cbind(Cd,1,1,0,0) ) d.asi <- ci.exp( md, subset=c("dur","DM"), ctr.mat=cbind(Cd,1,0,1,0) ) d.oth <- ci.exp( md, subset=c("dur","DM"), ctr.mat=cbind(Cd,1,0,0,1) ) tmpl <- function() { par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( d.pt, d.eff, type="n", log="y",las=1, ylim=c(0.2,5), xlab="DM duration (years)", ylab="RR of TB versus non-DM persons") abline( v=seq(0,16,2), h=c(5:15/10,2:10), col=gray(0.8) ) abline( h= 1 ) matlines( d.pt, cbind(d.eff,d.afr,d.asi,d.oth), type="l", lty=1, lwd=c(5,2,2), col=rep(ecol,each=3) ) rect( cnr(c(80,100),c(80,100)), col="white", border=gray(0.8) ) cc <- cnr(82,97) text( cc$x, cc$y*0.85^(0:3), names(ecol)[c(4,1,3,2)], col=ecol[c(4,1,3,2)], adj=0 ) box() } tmpl() ################################################### ### code chunk number 10: TB-dur.rnw:128-131 ################################################### ################################################### ### code chunk number 11: xtabs ################################################### xtabs( D.tb ~ Region + state, data=Dtb ) ################################################### ### code chunk number 12: dur2 ################################################### nk <- 3 d.kn <- with( subset(Dfu,state=="DM"), c( 0, quantile( rep(dur,D.tb), 1:nk/(nk+0.5) ) ) ) md <- update( md, . ~ . ) Cd <- Ns( d.pt, kn=d.kn ) d.eff <- ci.exp( md, subset=c("dur","DM"), ctr.mat=cbind(Cd,1,0,0,0) ) d.afr <- ci.exp( md, subset=c("dur","DM"), ctr.mat=cbind(Cd,1,1,0,0) ) d.asi <- ci.exp( md, subset=c("dur","DM"), ctr.mat=cbind(Cd,1,0,1,0) ) d.oth <- ci.exp( md, subset=c("dur","DM"), ctr.mat=cbind(Cd,1,0,0,1) ) tmpl <- function(){ par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( d.pt, d.eff, type="n", log="y",las=1, ylim=c(0.2,5), xlim=c(0,13), xlab="DM duration (years)", ylab="RR of TB versus non-DM persons") abline( v=seq(0,16,1), h=c(5:10/10,1.5,2:10), col=gray(0.8) ) abline( h= 1 ) matlines( d.pt, cbind(d.eff,d.afr,d.asi,d.oth), type="l", lty=1, lwd=c(5,2,2), col=rep(ecol,each=3) ) # rect( cnr(c(86,100),c(80,100)), col="white", border=gray(0.8) ) cc <- cnr(3,3) text( cc$x, cc$y*1.15^(0:3), names(ecol)[c(2,3,1,4)], col=ecol[c(2,3,1,4)], adj=c(0,0) ) box() } tmpl() ################################################### ### code chunk number 13: dur2-emf ################################################### pdf( "IJTLD-Fig4.pdf", width=7, height=7, pointsize=16 ) tmpl() dev.off() ecol <- bwcol pdf( "IJTLD-Fig4-bw.pdf", width=7, height=7, pointsize=16 ) tmpl() dev.off() ecol <- Ecol ################################################### ### code chunk number 14: test-int ################################################### gc() mdi <- update( md, . ~ . + Region:dur ) round( ci.exp( mdi ), 3 ) 1-pchisq( md$deviance - mdi$deviance, md$df.res - mdi$df.res ) ################################################### ### code chunk number 15: dur3 ################################################### d.eff <- ci.exp( mdi, subset=c("dur","DM"), ctr.mat=cbind(Cd,d.pt, 0 ,0,0,1,0,0,0) ) d.afr <- ci.exp( mdi, subset=c("dur","DM"), ctr.mat=cbind(Cd,d.pt,d.pt,0,0,1,1,0,0) ) d.asi <- ci.exp( mdi, subset=c("dur","DM"), ctr.mat=cbind(Cd,d.pt,0,d.pt,0,1,0,1,0) ) d.oth <- ci.exp( mdi, subset=c("dur","DM"), ctr.mat=cbind(Cd,d.pt,0,0,d.pt,1,0,0,1) ) par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( d.pt, d.eff, type="n", log="y",las=1, ylim=c(0.2,5), xlim=c(0,13), xlab="DM duration (years)", ylab="RR of TB versus non-DM persons") abline( v=seq(0,16,2), h=c(5:10/10,1.5,2:10), col=gray(0.8) ) abline( h= 1 ) matlines( d.pt, cbind(d.eff,d.afr,d.asi,d.oth), type="l", lty=1, lwd=c(3,1,1), col=rep(ecol,each=3) ) box() ################################################### ### code chunk number 16: TB-dur.rnw:234-239 (eval = FALSE) ################################################### ## md <- glm( D.tb ~ Ns( ax , kn=a.kn ) + sex + Region*state + ## Ns( dur, kn=d.kn ), ## offset = log(Y/10^5), ## family = poisson, ## data = Dtb ) ################################################### ### code chunk number 17: dur-per ################################################### md.dev <- md$deviance md.df <- md$df.res rm( md ) gc() mdp <- glm( D.tb ~ Ns( ax , kn=a.kn ) + sex + Region*state + state:Ns( px, kn=p.kn ) + Ns( dur, kn=d.kn ):Ns( px, kn=p.kn ), offset = log(Y/10^5), family = poisson, data = Dtb ) round( ci.exp( mdp ), 3 ) ci.exp( mdp, subset=c("DM","dur","Well:Ns") ) 1-pchisq(abs(md.dev-mdp$deviance), abs(md.df -mdp$df.res) ) ################################################### ### code chunk number 18: int ################################################### tmpl <- function(ci=FALSE){ par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( 0:16, 0:16, type="n", log="y",las=1, ylim=c(0.5,5), xlab="DM duration (years)", ylab="RR of TB versus non-DM persons") abline( v=seq(0,16,2), h=c(1:15/10,2:10), col=gray(0.8) ) abline( h= 1 ) box() round( ci.exp( mdp ), 3 ) round( ci.exp( mdp, subset=c("DM","dur") ), 3 ) for( yod in 1995:2009 ) { d.pt <- seq( 0,2010-yod,0.1) p.pt <- seq(yod,2010 ,0.1) Cp <- Ns( p.pt, kn=p.kn) CM <- model.matrix( ~ Ns( p.pt, kn=p.kn ):Ns( d.pt, kn=d.kn) )[,-1] d.eff <- ci.exp( mdp, subset=c("DM","dur"), ctr.mat=cbind(1,0,0,0,Cp,CM) ) matlines( d.pt, d.eff, type="l", lty=if(!ci) c(1,0,0) else 1, lwd=c(3,1,1), col=gray(1-(2015-yod)/22) ) np <- length( d.pt ) if( (yod %% 2) == 1 ) text( d.pt[np], d.eff[np,1], paste(yod), adj=c(0,1), font=1, col=gray(1-(2015-yod)/22) ) } } tmpl() ################################################### ### code chunk number 19: int-ci ################################################### tmpl(ci=TRUE) ################################################### ### code chunk number 20: Fig5 ################################################### pdf( "IJTLD-Fig5.pdf", width=7, height=7, pointsize=16 ) tmpl() dev.off() ecol <- bwcol pdf( "IJTLD-Fig5-bw.pdf", width=7, height=7, pointsize=16 ) tmpl() dev.off() ecol <- Ecol