### R code from vignette source 'LCa-lungF-s' ### Encoding: UTF-8 ################################################### ### code chunk number 1: LCa-lungF-s.rnw:4-9 ################################################### 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: LCa-lungF-s.rnw:27-31 ################################################### library( Epi ) lC <- read.table( "../data/lung-mf.txt", header=TRUE ) lF <- subset( lC, sex==2 ) head( lF ) ################################################### ### code chunk number 3: LCa-lungF-s.rnw:37-40 ################################################### t5 <- xtabs( cbind(D,Y=Y/1000) ~ A5 + P5, data=lF ) str( t5 ) r5 <- t5[,,"D"]/t5[,,"Y"] ################################################### ### code chunk number 4: rates ################################################### par( mfrow=c(2,2),mar=c(3,3,0,0),oma=c(0,0,1,0),mgp=c(3,1,0)/1.6,bty="n",las=1 ) rateplot( r5*100, ylab="", col=heat.colors(20)[1:20], lwd=3 ) mtext( "Lung cancer rates per 100,000 PY in Danish women", outer=TRUE ) ################################################### ### code chunk number 5: LCa-lungF-s.rnw:63-67 ################################################### nk <- 6 ( a.kn <- with( lF, quantile( rep( A,D), probs=(1:nk-0.5)/nk ) ) ) ( p.kn <- with( lF, quantile( rep(P ,D), probs=(1:nk-0.5)/nk ) ) ) ( c.kn <- with( lF, quantile( rep(P-A,D), probs=(1:nk-0.5)/nk ) ) ) ################################################### ### code chunk number 6: LCa-lungF-s.rnw:73-75 ################################################### APC <- apc.fit( lF, npar=list(A=a.kn,P=p.kn,C=c.kn), ref.p=1980, ref.c=1930, scale=10^3 ) ################################################### ### code chunk number 7: plotAPC ################################################### par( mar=c(3,4,0,4), mgp=c(3,1,0)/1.6, las=1, bty="n") plot( APC, "Female lung cancer in Danmark per 100,000 PY", col="red" ) ################################################### ### code chunk number 8: LCa-lungF-s.rnw:94-100 ################################################### LCaP <- LCa.fit( lF, npar=list(a=a.kn,p=p.kn,pi=a.kn,c=c.kn,ci=a.kn), a.ref=60, p.ref=1980, model="APa", VC=TRUE, quiet=FALSE ) LCaC <- LCa.fit( lF, npar=list(a=a.kn,p=p.kn,pi=a.kn,c=c.kn,ci=a.kn), a.ref=60, p.ref=1930, model="ACa", VC=T, quiet=FALSE ) ################################################### ### code chunk number 9: LCa-lungF-s.rnw:104-107 ################################################### round( rbind( c( LCaP$df, LCaP$dev ), c( LCaC$df, LCaC$dev ) ), 1 ) APC$Anova ################################################### ### code chunk number 10: LCaplot ################################################### par( mfrow=c(2,3) ) plot( LCaP ) plot( LCaC ) ################################################### ### code chunk number 11: LCa-lungF-s.rnw:140-145 ################################################### p.pt <- 1950:1997 ; np <- length(p.pt) a.pt <- 5:8*10 ; na <- length(a.pt) nd <- data.frame( A = rep(a.pt,each=np+1), P = rep(c(NA,p.pt), na), Y = 1000 )[-1,] ################################################### ### code chunk number 12: LCa-lungF-s.rnw:151-157 ################################################### AP <- glm( D ~ Ns(A,knots=a.kn)+Ns(P,knots=p.kn), offset=log(Y), family=poisson, data=lF ) AC <- glm( D ~ Ns(A,knots=a.kn)+ Ns(P-A,knots=c.kn), offset=log(Y), family=poisson, data=lF ) APC <- glm( D ~ Ns(A,knots=a.kn)+Ns(P,knots=p.kn)+Ns(P-A,knots=c.kn), offset=log(Y), family=poisson, data=lF ) ################################################### ### code chunk number 13: LCa-lungF-s.rnw:161-166 ################################################### fAP <- ci.pred( AP , nd ) fAC <- ci.pred( AC , nd ) fAPC <- ci.pred( APC, nd ) fLCaP <- predict( LCaP, nd, sim=10000 )*1000 fLCaC <- predict( LCaC, nd, sim=10000 )*1000 ################################################### ### code chunk number 14: fL-cmpt ################################################### ppm <- function( prd, mod ) { matplot( nd$P-nd$A, prd, type="l", lwd=c(2,1,1), lty=c(1,3,3), col="black", log="y", ylim=c(0.05,3), xlim=1860+c(0,90), ylab="", xlab="Date of birth" ) text( 1860, 2, mod, adj=c(0,1) ) matplot( nd$P , prd, type="l", lwd=c(2,1,1), lty=c(1,3,3), col="black", log="y", ylim=c(0.05,3), xlim=1920+c(0,90), ylab="", xlab="Date of event" ) text( 1920, 2, mod, adj=c(0,1) ) } par( mfcol=c(2,5), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n", las=1 ) ppm( fAP , "Age-Period") ppm( fLCaP, "Lee-Carter-Period") ppm( fAPC , "Age-Period-Cohort") ppm( fLCaC, "Lee-Carter Cohort") ppm( fAC , "Age-Cohort") ################################################### ### code chunk number 15: LCa-lungF-s.rnw:195-233 ################################################### # Age-specific rates by period p.pt <- 1950+0:4*10 ; np <- length(p.pt) a.pt <- 40:90 ; na <- length(a.pt) nd <- data.frame( A = rep(a.pt, np), P = rep(p.pt,each=na), Y = 1000 ) nd <- rbind( nd[ 1:na,], NA, nd[1*na+1:na,], NA, nd[2*na+1:na,], NA, nd[3*na+1:na,], NA, nd[4*na+1:na,] ) pAP <- ci.pred( AP , nd ) pAC <- ci.pred( AC , nd ) pAPC <- ci.pred( APC, nd ) pLCP <- predict( LCaP, nd, sim=10000 )*1000 pLCC <- predict( LCaC, nd, sim=10000 )*1000 # Age-specific rates by cohort c.pt <- 1870+0:8*10 ; nc <- length(c.pt) a.pt <- 40:90 ; na <- length(a.pt) nc <- data.frame( A = rep(a.pt, nc), C = rep(c.pt,each=na), Y = 1000 ) nc <- rbind( nc[ 1:na,], NA, nc[1*na+1:na,], NA, nc[2*na+1:na,], NA, nc[3*na+1:na,], NA, nc[4*na+1:na,], NA, nc[5*na+1:na,], NA, nc[6*na+1:na,], NA, nc[7*na+1:na,], NA, nc[8*na+1:na,] ) nc$P <- nc$C + nc$A nc <- subset( nc, (P>1943 & P<2000) | is.na(A) ) cAP <- ci.pred( AP , nc ) cAC <- ci.pred( AC , nc ) cAPC <- ci.pred( APC , nc ) cLCP <- predict( LCaP, nc, sim=10000 )*1000 cLCC <- predict( LCaC, nc, sim=10000 )*1000 ################################################### ### code chunk number 16: fL-cmpa ################################################### ppm <- function( prp, prc, mod ) { matplot( nd$A, prp, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(0.02,2), xlim=c(30,90), ylab="", xlab="Age" ) text( 30, 2, mod, adj=c(0,1) ) matplot( nc$A, prc, type="l", lty=1, lwd=c(3,1,1), col="black", log="y", ylim=c(0.02,2), xlim=c(30,90), ylab="", xlab="Age" ) text( 30, 2, mod, adj=c(0,1) ) } par( mfcol=c(2,5), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n", las=1 ) ppm( pAP , cAP , "Age-Period") ppm( pLCP, cLCP, "Lee-Carter-Period") ppm( pAPC, cAPC, "Age-Period-Cohort") ppm( pLCC, cLCC, "Lee-Carter Cohort") ppm( pAC , cAC , "Age-Cohort")