### R code from vignette source 'curve-sol.rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: curve-sol.rnw:11-14 ################################################### options( width=90, 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: curve-sol.rnw:60-65 ################################################### library( Epi ) sessionInfo() data( testisDK ) str( testisDK ) head( testisDK ) ################################################### ### code chunk number 3: curve-sol.rnw:70-74 ################################################### 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 ) ################################################### ### code chunk number 4: curve-sol.rnw:89-91 ################################################### ml <- glm( D ~ A, offset=log(Y), family=poisson, data=testisDK ) ci.exp( ml ) ################################################### ### code chunk number 5: curve-sol.rnw:110-111 ################################################### ( cf <- coef( ml ) ) ################################################### ### code chunk number 6: curve-sol.rnw:116-117 ################################################### round( cbind( 25:45, exp( cf[1] + cf[2]*(25:45) )*10^5 ), 3 ) ################################################### ### code chunk number 7: curve-sol.rnw:141-143 ################################################### ( CM <- cbind( 1, 25:45 ) ) round( ci.exp( ml, ctr.mat=CM )*10^5, 3 ) ################################################### ### code chunk number 8: lin ################################################### C1 <- cbind( 1, 15:65 ) matplot( 15:65, ci.exp( ml, ctr.mat=C1 )*10^5, log="y", xlab="Age", ylab="Testis cancer incidence rate per 100,000 PY", type="l", lty=1, lwd=c(3,1,1), col="black" ) ################################################### ### code chunk number 9: curve-sol.rnw:188-190 ################################################### mq <- glm( D ~ A + I(A^2), offset=log(Y), family=poisson, data=testisDK ) ci.exp( mq, Exp=F ) ################################################### ### code chunk number 10: qdr ################################################### aa <- 15:65 C2 <- cbind( 1, aa, aa^2 ) matplot( aa, ci.exp( mq, ctr.mat=C2 )*10^5, log="y", xlab="Age", ylab="Testis cancer incidence rate per 100,000 PY", type="l", lty=1, lwd=c(3,1,1), col="black" ) matlines( aa, ci.exp( ml, ctr.mat=C1 )*10^5, type="l", lty=1, lwd=c(3,1,1), col="blue" ) ################################################### ### code chunk number 11: curve-sol.rnw:211-212 ################################################### 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 ) C3 <- cbind( 1, aa, aa^2, aa^3 ) matplot( aa, ci.exp( mc, ctr.mat=C3 )*10^5, log="y", xlab="Age", ylab="Testis cancer incidence rate per 100,000 PY", type="l", lty=1, lwd=c(3,1,1), col="black" ) ################################################### ### code chunk number 13: spl ################################################### library( splines ) ms <- glm( D ~ Ns(A,knots=seq(15,65,10)), offset=log(Y), family=poisson, data=testisDK ) As <- Ns( aa, knots=seq(15,65,10) ) matplot( aa, ci.exp( ms, ctr.mat=cbind(1,As) )*10^5, log="y", xlab="Age", ylab="Testis cancer incidence rate per 100,000 PY", type="l", lty=1, lwd=c(3,1,1), col="black" ) ################################################### ### code chunk number 14: spl-P ################################################### msp <- glm( D ~ Ns(A,knots=seq(15,65,10)) + P, offset=log(Y), family=poisson, data=testisDK ) CM <- cbind( 1, Ns( aa, knots=seq(15,65,10) ), 1970 ) matplot( aa, ci.exp( msp, ctr.mat=CM )*10^5, log="y", xlab="Age", ylab="Testis cancer incidence rate per 100,000 PY", type="l", lty=1, lwd=c(3,1,1), col="black" ) ################################################### ### code chunk number 15: curve-sol.rnw:322-323 ################################################### ci.exp( msp, subset="P" ) ################################################### ### code chunk number 16: spl-P1 ################################################### yy <- 1943:1996 Cp1 <- cbind( yy - 1970 ) matplot( yy, ci.exp( msp, ctr.mat=Cp1, subset="P" ), log="y", xlab="Date", ylab="RR of Testis cancer", type="l", lty=1, lwd=c(3,1,1), col="black" ) abline( h=1 ) ################################################### ### code chunk number 17: spl-P2 ################################################### msp <- glm( D ~ Ns(A,knots=seq(15,65,10)) + P + I(P^2), offset=log(Y), family=poisson, data=testisDK ) Cq <- cbind( yy, yy^2 ) - cbind( rep(1970,length(yy)), 1970^2 ) matplot( yy, ci.exp( msp, ctr.mat=Cq, subset="P" ), log="y", xlab="Age", ylab="Testis cancer incidence RR", type="l", lty=1, lwd=c(3,1,1), col="black" ) abline( h=1, v=1970 ) ################################################### ### code chunk number 18: curve-sol.rnw:380-384 ################################################### 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( mssp, msp, test="Chisq" ) ################################################### ### code chunk number 19: splA-spl ################################################### Ps <- Ns( yy , knots=seq(1950,1990,10) ) Pr <- Ns( rep(1970,length(yy)), knots=seq(1950,1990,10) ) matplot( yy, ci.exp( mssp, ctr.mat=Ps-Pr, subset="P" ), log="y", xlab="Age", ylab="Testis cancer incidence RR", type="l", lty=1, lwd=c(3,1,1), col="black" ) ################################################### ### code chunk number 20: spl-splP ################################################### Ar <- Ns( rep(1970,length(aa)), knots=seq(1950,1990,10) ) matplot( aa, ci.exp( mssp, ctr.mat=cbind(1,As,Ar) )*10^5, log="y", xlab="Age", ylab="Testis cancer incidence RR", type="l", lty=1, lwd=c(3,1,1), col="black" ) ################################################### ### code chunk number 21: finP ################################################### a.kn <- seq(15,65,10) p.kn <- seq(1950,1990,10) a.pt <- 10:65 p.pt <- 1945:1993 p.ref <- 1970 na <- length(a.pt) np <- length(p.pt) As <- Ns( a.pt, knots=a.kn ) Ps <- Ns( p.pt, knots=p.kn ) Prp <- Ns( rep(p.ref,np), knots=p.kn ) Pra <- Ns( rep(p.ref,na), knots=p.kn ) mAP <- glm( D ~ Ns(A,knots=a.kn) + Ns(P,knots=p.kn), offset=log(Y), family=poisson, data=testisDK ) par( mfrow=c(1,2) ) matplot( a.pt, ci.exp( mAP, ctr.mat=cbind(1,As,Pra) )*10^5, log="y", xlab="Age", ylab=paste( "Testis cancer incidence per 100,000 PY, in", p.ref ), type="l", lty=1, lwd=c(3,1,1), col="black", ylim=c(1,20) ) matplot( p.pt, ci.exp( mAP, ctr.mat=Ps-Prp, subset="P" ), 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=p.ref ) ################################################### ### code chunk number 22: finB ################################################### testisDK <- transform( testisDK, B = P-A ) # with( testisDK, hist( rep(B,D), breaks=100, col="black" ) ) a.kn <- seq(15,65,5) b.kn <- seq(1900,1970,5) a.pt <- 10:65 b.pt <- 1890:1970 b.ref <- 1940 na <- length(a.pt) nb <- length(b.pt) As <- Ns( a.pt, knots=a.kn ) Bs <- Ns( b.pt, knots=b.kn ) Brb <- Ns( rep(b.ref,nb), knots=b.kn ) Bra <- Ns( rep(b.ref,na), knots=b.kn ) mAB <- glm( D ~ Ns(A,knots=a.kn) + Ns(B,knots=b.kn), offset=log(Y), family=poisson, data=testisDK ) par( mfrow=c(1,2) ) matplot( a.pt, ci.exp( mAB, ctr.mat=cbind(1,As,Bra) )*10^5, log="y", xlab="Age", ylab=paste( "Testis cancer incidence per 100,000 PY, in", b.ref, "birth cohort"), type="l", lty=1, lwd=c(3,1,1), col="black", ylim=c(1,20) ) matplot( b.pt, ci.exp( mAB, ctr.mat=Bs-Brb, subset="B" ), 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=b.ref )