### R code from vignette source 'cont-eff-s.rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: cont-eff-s.rnw:2-6 ################################################### options( width=90, prompt=" ", continue=" ", SweaveHooks=list( fig=function() par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6), bty="n", las=1 ) ) ################################################### ### code chunk number 2: cont-eff-s.rnw:22-27 ################################################### library( Epi ) sessionInfo() data( testisDK ) str( testisDK ) head( testisDK ) ################################################### ### code chunk number 3: cont-eff-s.rnw:32-45 ################################################### round( ftable( xtabs( cbind(D,PY=Y/1000) ~ I(floor(A/10)*10) + I(floor(P/10)*10), data=testisDK ), row.vars=c(3,1) ), 1 ) ST <- stat.table( list(A=floor(A/10)*10, P=floor(P/10)*10), list( D=sum(D), Y=sum(Y/1000), rate=ratio(D,Y,10^5) ), margins=TRUE, data=testisDK ) print( ST, digits=c(sum=0,rate=2) ) ################################################### ### code chunk number 4: cont-eff-s.rnw:52-54 ################################################### ml <- glm( D ~ A, offset=log(Y), family=poisson, data=testisDK ) ci.exp( ml ) ################################################### ### code chunk number 5: cont-eff-s.rnw:66-67 ################################################### ( cf <- coef( ml ) ) ################################################### ### code chunk number 6: cont-eff-s.rnw:72-73 ################################################### round( cbind( 25:45, exp( cf[1] + cf[2]*(25:45) )*10^5 ), 3 ) ################################################### ### code chunk number 7: cont-eff-s.rnw:84-86 ################################################### nd <- data.frame( A = 15:65, Y = 10^5 ) head( ci.pred( ml, nd ) ) ################################################### ### code chunk number 8: lin ################################################### matshade( nd$A, ci.pred( ml, nd ), plot=TRUE, lwd=2, col="black", log="y", xlab="Age", ylab="Testis cancer incidence rate per 100,000 PY" ) ################################################### ### code chunk number 9: cont-eff-s.rnw:102-104 ################################################### mq <- glm( D ~ A + I(A^2), offset=log(Y), family=poisson, data=testisDK ) ci.exp( mq, Exp=F ) ################################################### ### code chunk number 10: qdr ################################################### matshade( nd$A, cbind( ci.pred( mq, nd ), ci.pred( ml, nd ) ), plot=TRUE, lwd=2, col=c("black","blue"), log="y", xlab="Age", ylab="Testis cancer incidence rate per 100,000 PY" ) ################################################### ### code chunk number 11: cont-eff-s.rnw:123-124 ################################################### anova( mq, ml, test="Chisq" ) ################################################### ### code chunk number 12: cub ################################################### mc <- glm( D ~ A + I(A^2) + I(A^3), offset = log(Y), family = poisson, data = testisDK ) matshade( nd$A, cbind( ci.pred( mc, nd ), ci.pred( mq, nd ) ), plot=TRUE, lwd=2, col=c("black","blue"), log="y", xlab="Age", ylab="Testis cancer incidence rate per 100,000 PY" ) ################################################### ### code chunk number 13: cont-eff-s.rnw:147-148 ################################################### anova( ml, mq, mc, test="Chisq" ) ################################################### ### code chunk number 14: spl ################################################### library( splines ) ms <- glm( D ~ Ns(A,knots=seq(15,65,10)), offset = log(Y), family = poisson, data = testisDK ) matshade( nd$A, cbind( ci.pred( ms, nd ), ci.pred( mc, nd ) ), plot=TRUE, lwd=2, col=c("black","blue"), log="y", xlab="Age", ylab="Testis cancer incidence rate per 100,000 PY" ) ################################################### ### code chunk number 15: cont-eff-s.rnw:179-184 ################################################### msp <- glm( D ~ Ns(A,knots=seq(15,65,10)) + P, offset = log(Y), family = poisson, data = testisDK ) round( ci.exp( msp ), 3 ) ################################################### ### code chunk number 16: cont-eff-s.rnw:191-193 ################################################### nd <- cbind( nd, P=1970 ) head( nd ) ################################################### ### code chunk number 17: spl-P ################################################### matshade( nd$A, cbind( ci.pred( msp, nd ), ci.pred( ms , nd ) ), plot=TRUE, log="y", xlab="Age", ylab="Testis cancer incidence rate in 1970 per 100,000 PY", type="l", lty=1, lwd=2, col=c("black","blue") ) ################################################### ### code chunk number 18: cont-eff-s.rnw:210-211 ################################################### ci.exp( msp, subset="P" ) ################################################### ### code chunk number 19: spl-P1 ################################################### nl <- list( data.frame(A=50,P=1943:1996), data.frame(A=50,P=1970)) RR <- ci.exp( msp, nl ) matshade( nl[[1]]$P, RR, plot=TRUE, log="y", xlab="Age", ylab="Testis cancer incidence RR", lty=1, lwd=2, col="black" ) abline( h=1, v=1970, lty=3 ) ################################################### ### code chunk number 20: cont-eff-s.rnw:241-245 ################################################### msp2 <- glm( D ~ Ns(A,knots=seq(15,65,10)) + P + I(P^2), offset = log(Y), family = poisson, data = testisDK ) ################################################### ### code chunk number 21: spl-P2 ################################################### RR <- ci.exp( msp2, nl ) matshade( nl[[1]]$P, RR, plot=TRUE, log="y", xlab="Age", ylab="Testis cancer incidence RR", lty=1, lwd=2, col="black" ) abline( h=1, v=1970, lty=3 ) ################################################### ### code chunk number 22: cont-eff-s.rnw:263-267 ################################################### mssp <- glm( D ~ Ns(A,knots=seq(15,65,10)) + Ns(P,knots=seq(1950,1990,10)), offset=log(Y), family=poisson, data=testisDK ) anova( msp, msp2, mssp, test="Chisq" ) ################################################### ### code chunk number 23: splA-Pspl ################################################### matshade( nl[[1]]$P, ci.exp( mssp, nl ), plot=TRUE, log="y", xlab="Date of FU", ylab="Testis cancer incidence RR", lty=1, lwd=2, col="black" ) abline( h=1, v=1970, lty=3 ) ################################################### ### code chunk number 24: spl-splP ################################################### matshade( nd$A, ci.pred( mssp, newdata=nd ), plot=TRUE, log="y", xlab="Age", ylab="Testis cancer incidence in 1970", lty=1, lwd=1, col="black" ) ################################################### ### code chunk number 25: finB ################################################### testisDK <- transform( testisDK, B = P - A ) mAB <- glm( D ~ Ns(A,knots=seq(15,65,10)) + Ns(B,knots=seq(1900,1970,5)), offset=log(Y), family=poisson, data=testisDK ) nd <- data.frame( A=15:65, B=1936, Y=10^5 ) nb <- data.frame( A=40, B=1854:1996, Y=10^5 ) nr <- data.frame( A=40, B=1936, Y=10^5 ) par( mfrow=c(1,2) ) matshade( nd$A, ci.pred( mAB, newdata=nd ), plot=TRUE, log="y", xlab="Age", ylab="Testis cancer incidence per 100,000 PY, in 1908 birth cohort", type="l", lty=1, lwd = 2, col="black", ylim=c(1,20) ) matshade( nb$B, ci.exp( mAB, ctr.mat=list(nb,nr) ), plot=TRUE, log="y", xlab="Age", ylab="Testis cancer incidence RR", type="l", lty=1, lwd=c(3,1,1), col="black", ylim=c(1,20)/4 ) abline( h=1, v=1936, lty=3 ) rect( cal.yr("1914-07-28"), 0.01, cal.yr("1918-11-11"), 10, col="#0000BB44", border="transparent") rect( cal.yr("1939-09-01"), 0.01, cal.yr("1945-05-05"), 10, col="#0000BB44", border="transparent") ################################################### ### code chunk number 26: cont-eff-s.rnw:334-338 ################################################### library( mgcv ) mAB <- gam( D ~ s(A,k=40) + s(B,k=40) + offset(log(Y)), family=poisson, data=testisDK ) gam.check( mAB ) ################################################### ### code chunk number 27: fingam ################################################### par( mfrow=c(1,2) ) matshade( nd$A, ci.pred( mAB, newdata=nd ), plot=TRUE, log="y", xlab="Age", ylab="Testis cancer incidence per 100,000 PY, in 1908 birth cohort", type="l", lty=1, lwd = 2, col="black", ylim=c(1,20) ) matshade( nb$B, ci.exp( mAB, ctr.mat=list(nb,nr) ), plot=TRUE, log="y", xlab="Age", ylab="Testis cancer incidence RR", type="l", lty=1, lwd=c(3,1,1), col="black", ylim=c(1,20)/4 ) abline( h=1, v=1936, lty=3 ) rect( cal.yr("1914-07-28"), 0.01, cal.yr("1918-11-11"), 10, col="#0000BB33", border="transparent") rect( cal.yr("1939-09-01"), 0.01, cal.yr("1945-05-05"), 10, col="#0000BB33", border="transparent")