### R code from vignette source 'brcapr-s' ### Encoding: UTF-8 ################################################### ### code chunk number 1: brcapr-s.rnw:5-10 ################################################### 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: brcapr-s.rnw:20-24 ################################################### library( Epi ) breast <- read.table("../data/breast.txt", header=T ) str( breast ) summary( breast ) ################################################### ### code chunk number 3: brcapr-s.rnw:33-39 ################################################### breast <- transform( breast, up = P-A-C ) breast <- transform( breast, A = A+(1+up)/3, P = P+(2-up)/3, C = C+(1+up)/3 ) with( breast, summary( P-A-C ) ) head( breast ) ################################################### ### code chunk number 4: ratetab ################################################### ti <- 4 rt <- with( subset( breast, A>30 ), tapply( D, list(floor( A /ti)*ti+ti/2, floor((P-1943)/ti)*ti+ti/2+1943), sum ) / tapply( Y, list(floor( A /ti)*ti+ti/2, floor((P-1943)/ti)*ti+ti/2+1943), sum ) * 10^5 ) par( mfrow=c(2,2), mar=c(3,3,0,0), oma=c(0,0,1,1), mgp=c(3,1,0)/1.6 ) rateplot( rt, which= c( "ap", "ac", "pa", "ca"), col=heat.colors(22), ann=TRUE ) ################################################### ### code chunk number 5: apcfit-1 ################################################### par( mfrow=c(1,1), mar=c(3,3,1,3) ) m1 <- apc.fit( subset( breast, A>30 ), npar = c(8,6,10), ref.c = 1920, scale = 10^5 ) plot( m1 ) ################################################### ### code chunk number 6: apcfit-2 ################################################### par( las=1, mar=c(3,4,1,4), mgp=c(3,1,0)/1.5 ) fp <- apc.frame( a.lab = seq(30,90,10), cp.lab = seq(1860,2005,20), r.lab = c(c(1,2,5)*10,c(1,2,5)*100), # rr.lab = r.lab / rr.ref, rr.ref = 100, a.tic = seq(30,90,5), cp.tic = seq(1855,2005,5), r.tic = c(9,1:9*10,1:5*100), # rr.tic = r.tic / rr.ref, tic.fac = 1.3, a.txt = "Age", cp.txt = "Calendar time", r.txt = "Rate per 100,000 person-years", rr.txt = "Rate ratio", gap = 8, col.grid = gray(0.85), sides = c(1,2,4) ) lines( m1, ci=T, col="red" ) matshade( m1$Age[,1], m1$Age[,-1], col="red" ) pc.matshade( m1$Per[,1], m1$Per[,-1], col="red" ) pc.matshade( m1$Coh[,1], m1$Coh[,-1], col="red" ) pc.points( 1920, 1, pch=16, col="red" ) ################################################### ### code chunk number 7: brcapr-s.rnw:132-140 ################################################### # Last knot and last point on period scale ( P.rf <- c( max(m1$Knots$Per), max(m1$Per[,1]) ) ) # Last point plus one 20 years ago ( P.pt <- P.rf[2] + 0:1*20 ) # Linear interpolation of log-rates at the two reference points ( Pp <- approx( m1$Per[,1], log(m1$Per[,2]), P.rf )$y ) # Liner extrapoltion throug these two points to the future points ( P.eff <- Pp[2] + (Pp[2]-Pp[1])/diff(P.rf)*(P.pt-P.rf[2]) ) ################################################### ### code chunk number 8: brcapr-s.rnw:143-147 ################################################### ( C.rf <- c( max( m1$Knots$Coh ), max( m1$Coh[,1] ) ) ) ( C.pt <- C.rf[2] + 0:1*20 ) ( Cp <- approx( m1$Coh[,1], log(m1$Coh[,2]), C.rf )$y ) ( C.eff <- Cp[2] + (Cp[2]-Cp[1])/diff(C.rf)*(C.pt-C.rf[2]) ) ################################################### ### code chunk number 9: apcfit-3 ################################################### par( las=1, mar=c(3,4,1,4), mgp=c(3,1,0)/1.5 ) fp <- apc.frame( a.lab = seq(30,90,10), cp.lab = seq(1860,2020,20), r.lab = c(c(1,2,5)*10,c(1,2,5)*100), # rr.lab = r.lab / rr.ref, rr.ref = 100, a.tic = seq(30,90,5), cp.tic = seq(1855,2025,5), r.tic = c(9,1:9*10,1:5*100), # rr.tic = r.tic / rr.ref, tic.fac = 1.3, a.txt = "Age", cp.txt = "Calendar time", r.txt = "Rate per 100,000 person-years", rr.txt = "Rate ratio", gap = 8, col.grid = gray(0.85), sides = c(1,2,4) ) lines( m1, frame.par=fp, ci=T, col="red", lwd=c(4,1,1), knots=TRUE ) lines( P.pt-fp[1], exp(P.eff)*fp[2], col=gray(0.0), lty="11", lwd=2 ) lines( C.pt-fp[1], exp(C.eff)*fp[2], col=gray(0.0), lty="11", lwd=2 ) ################################################### ### code chunk number 10: brcapr-s.rnw:196-202 ################################################### M1 <- glm( D ~ Ns( A, knots=m1$Knots$Age ) + Ns( P , knots=m1$Knots$Per ) + Ns( P-A, knots=m1$Knots$Coh )[,-1], family = poisson, offset = log(Y), data = subset( breast, A>30 ) ) ################################################### ### code chunk number 11: brcapr-s.rnw:210-212 ################################################### c( M1$deviance, m1$Model$deviance ) summary( fitted(M1) - fitted(m1$Model) ) ################################################### ### code chunk number 12: brcapr-s.rnw:218-225 ################################################### a.pt <- seq(30,90,1/10) Pfr <- rbind( data.frame(A=a.pt,P=2020,Y=1000), NA, data.frame(A=a.pt,P=2030,Y=1000) ) Cfr <- rbind( data.frame(A=a.pt,P=a.pt+1960,Y=1000), NA, data.frame(A=a.pt,P=a.pt+1970,Y=1000) ) prP <- ci.pred( M1, Pfr ) prC <- ci.pred( M1, Cfr ) ################################################### ### code chunk number 13: pred1 ################################################### ( ct <- c(0,which( is.na( prP[,1] ) ),nrow(prP)+1 ) ) for( i in 1:2 ) { wh <- (ct[i]+1):(ct[i+1]-1) matshade( Pfr$A[wh], cbind( prP, prC )[wh,], plot=(i==1), log="y", las=1, xlim=c(25,90), xlab="Age", ylab="Predicted breast cancer incidence per 1000 PY", type="l", lwd=1, lty=1, col=c("red","forestgreen") ) } text( rep(29.5,2), prP[c(1,603),1], paste(c(2020,2030)), col="red", adj=1, cex=0.8 ) text( rep(29.5,2), prC[c(1,603),1], paste(c(1960,1970)), col="forestgreen", adj=1, cex=0.8 ) ################################################### ### code chunk number 14: brcapr-s.rnw:251-261 ################################################### mp <- apc.fit( subset(breast, A>30), npar=list(A=m1$Knots$Age, P=m1$Knots$Per[-length(m1$Knots$Per)], C=m1$Knots$Coh), ref.c=1920, scale=10^5 ) mpc <- apc.fit( subset(breast, A>30), npar=list(A=m1$Knots$Age, P=m1$Knots$Per[-length(m1$Knots$Per)], C=m1$Knots$Coh[-length(m1$Knots$Coh)]), ref.c=1920, scale=10^5 ) ################################################### ### code chunk number 15: brcapr-s.rnw:266-266 ################################################### ################################################### ### code chunk number 16: apcfit-4 ################################################### par( las=1, mar=c(3,4,1,4), mgp=c(3,1,0)/1.5 ) fp <- apc.frame( a.lab = seq(30,90,10), cp.lab = seq(1860,2020,20), r.lab = c(c(1,2,5)*10,c(1,2,5)*100), # rr.lab = r.lab / rr.ref, rr.ref = 100, a.tic = seq(30,90,5), cp.tic = seq(1855,2025,5), r.tic = c(9,1:9*10,1:5*100), # rr.tic = r.tic / rr.ref, tic.fac = 1.3, a.txt = "Age", cp.txt = "Calendar time", r.txt = "Rate per 100,000 person-years", rr.txt = "Rate ratio", gap = 8, col.grid = gray(0.85), sides = c(1,2,4) ) lines( m1 , frame.par=fp, ci=T, col="black", lwd=c(3,1,1), knots=TRUE ) lines( mp , frame.par=fp, ci=T, col="red" , lty=1, lwd=c(3,1,1) ) lines( mpc, frame.par=fp, ci=T, col="limegreen", lty=3, lwd=c(3,1,1) ) ################################################### ### code chunk number 17: brcapr-s.rnw:299-304 ################################################### pr.slopes <- matrix( NA, 3, 3 ) rownames( pr.slopes ) <- c("Org","-lastP","-lastPC") colnames( pr.slopes ) <- c("P-slope","C-slope","P-C-slope") pr.slopes["Org","P-slope"] <- diff(Pp)/diff(P.rf) pr.slopes["Org","C-slope"] <- diff(Cp)/diff(C.rf) ################################################### ### code chunk number 18: brcapr-s.rnw:308-333 ################################################### ( P.rf <- c( max( mp$Knots$Per ), max( mp$Per[,1] ) ) ) P.pt <- P.rf[2] + 0:20 Pp <- approx( mp$Per[,1], log(mp$Per[,2]), P.rf )$y P.eff <- Pp[2] + (Pp[2]-Pp[1])/diff(P.rf)*(P.pt-P.rf[2]) ( C.rf <- c( max( mp$Knots$Coh ), max( mp$Coh[,1] ) ) ) C.pt <- C.rf[2] + 0:20 Cp <- approx( mp$Coh[,1], log(mp$Coh[,2]), C.rf )$y C.eff <- Cp[2] + (Cp[2]-Cp[1])/diff(C.rf)*(C.pt-C.rf[2]) pr.slopes["-lastP","P-slope"] <- diff(Pp)/diff(P.rf) pr.slopes["-lastP","C-slope"] <- diff(Cp)/diff(C.rf) ( P.rf <- c( max( mpc$Knots$Per ), max( mpc$Per[,1] ) ) ) P.pt <- P.rf[2] + 0:20 Pp <- approx( mpc$Per[,1], log(mpc$Per[,2]), P.rf )$y P.eff <- Pp[2] + (Pp[2]-Pp[1])/diff(P.rf)*(P.pt-P.rf[2]) ( C.rf <- c( max( mpc$Knots$Coh ), max( mpc$Coh[,1] ) ) ) C.pt <- C.rf[2] + 0:20 Cp <- approx( mpc$Coh[,1], log(mpc$Coh[,2]), C.rf )$y C.eff <- Cp[2] + (Cp[2]-Cp[1])/diff(C.rf)*(C.pt-C.rf[2]) pr.slopes["-lastPC","P-slope"] <- diff(Pp)/diff(P.rf) pr.slopes["-lastPC","C-slope"] <- diff(Cp)/diff(C.rf) pr.slopes[,3] <- pr.slopes[,1] + pr.slopes[,2] round( pr.slopes, 4 ) round( 100*(exp(pr.slopes)-1), 4 ) ################################################### ### code chunk number 19: brcapr-s.rnw:343-355 ################################################### Mp <- glm( D ~ Ns( A, knots=mp$Knots$Age ) + Ns( P , knots=mp$Knots$Per ) + Ns( P-A, knots=mp$Knots$Coh )[,-1], family = poisson, offset = log(Y), data = subset( breast, A>30 ) ) Mpc <- glm( D ~ Ns( A, knots=mpc$Knots$Age ) + Ns( P , knots=mpc$Knots$Per ) + Ns( P-A, knots=mpc$Knots$Coh )[,-1], family = poisson, offset = log(Y), data = subset( breast, A>30 ) ) ################################################### ### code chunk number 20: brcapr-s.rnw:361-365 ################################################### prPp <- ci.pred( Mp, Pfr ) prCp <- ci.pred( Mp, Cfr ) prPpc <- ci.pred( Mpc, Pfr ) prCpc <- ci.pred( Mpc, Cfr ) ################################################### ### code chunk number 21: predx ################################################### par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) matplot( Pfr$A, cbind( prP[,1], prPp[,1], prPpc[,1] ), log="y", las=1, xlim=c(25,90), xlab="Age", ylim=c(0.1,6), ylab="Predicted breast cancer incidence per 100,000 PY", type="l", lwd=3, lty=1, col=c("gray","limegreen","red") ) matplot( Pfr$A, cbind( prC[,1], prCp[,1], prCpc[,1] ), log="y", las=1, xlim=c(25,90), xlab="Age", ylim=c(0.1,6), ylab="Predicted breast cancer incidence per 100,000 PY", type="l", lwd=3, lty=1, col=c("gray","limegreen","red") )