### R code from vignette source 'param'
### Encoding: UTF-8

###################################################
### code chunk number 1: param.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,las=1,bty="n") ) )


###################################################
### code chunk number 2: categ.rnw:16-25
###################################################
library( Epi )
data( diet )
attach( diet )                  
y <- cal.yr(dox) - cal.yr(doe)  # Follow-up in years for each person
Y <- sum( y )                   # Total no. person-years
D <- sum( chd )                 # Total no. of CHD events
rate <- D / Y
erf <- exp( 1.96 / sqrt( D ) )
c( rate, rate/erf, rate * erf )


###################################################
### code chunk number 3: categ.rnw:37-40
###################################################
m1 <- glm( chd ~ 1, offset=log(y), family=poisson )
round( ci.lin( m1, Exp=TRUE ), 3 )
round( ci.exp( m1 ) )


###################################################
### code chunk number 4: categ.rnw:56-58
###################################################
m1 <- glm( chd ~ 1, offset=log(y/1000), family=poisson )
round( ci.exp( m1 ), 2 )


###################################################
### code chunk number 5: categ.rnw:64-66
###################################################
m2 <- glm( cbind(chd,y/1000) ~ 1, family=poisreg )
round( ci.exp( m2 ), 2 )


###################################################
### code chunk number 6: categ.rnw:81-94
###################################################
st <- stat.table( index = job,
               contents = list(D=sum(chd),
                               Y=sum(y),
                            rate=ratio(chd,y,1000)),
                 margins = TRUE,
                    data = diet )
print( st, digits=c(sum=0,ratio=2) )
str(st)
D <- st[1,]
Y <- st[2,]
( rate <- D/Y*1000 )
erf <- exp( 1.96 / sqrt( D ) )
round( cbind( D, Y, rate, rate/erf, rate * erf ), 2 )


###################################################
### code chunk number 7: categ.rnw:120-124
###################################################
RR <- rate[2] / rate[3]
erf.RR <- exp( 1.96 * sqrt( sum(1/D[2:3]) ) )
round( setNames( c( RR, RR/erf.RR, RR*erf.RR), 
                 c("RR",     "lo",      "hi") ), 2 )


###################################################
### code chunk number 8: categ.rnw:147-151
###################################################
m2a <- glm( cbind(chd,y/1000) ~ job    , family=poisreg )
round( ci.exp( m2a ), 2 )
m2b <- glm( cbind(chd,y/1000) ~ job - 1, family=poisreg )
round( ci.exp( m2b ), 2 )


###################################################
### code chunk number 9: categ.rnw:178-189
###################################################
CM <- rbind( c(1, 0, 0),
             c(0, 1, 0),
             c(0, 0, 1),
             c(1,-1, 0),
             c(1, 0,-1),
             c(0, 1,-1) )
rownames( CM ) <- c( levels(job), paste( levels(job)[c(1,1,2)],'vs.',
                                         levels(job)[c(2,3,3)] ) )
colnames( CM ) <- gsub( 'job', '', names(coef(m2b)) )
CM
round( ci.exp( m2b, ctr.mat=CM ), 2 )


###################################################
### code chunk number 10: categ.rnw:226-227
###################################################
round( ci.exp( m2b, subset='job', diff=TRUE ), 2 )


###################################################
### code chunk number 11: categ.rnw:231-233
###################################################
round( ci.exp( m2b, subset=1:3, diff=TRUE ), 2 )
round( ci.exp( m2b, subset=3:1, diff=TRUE ), 2 )


###################################################
### code chunk number 12: quant.rnw:65-70
###################################################
library( Epi )
data( testisDK )
tdk <- subset( testisDK, A>14 & A<65 )
tdk <- transform( tdk, A = A+1/2, P=P+1/2 )
head( tdk, 4 )


###################################################
### code chunk number 13: quant.rnw:72-74
###################################################
m1 <- glm( cbind(D,Y/10^5) ~ A, family=poisreg, data=tdk )
round( ci.exp( m1, pval=TRUE ), 3 )


###################################################
### code chunk number 14: quant.rnw:112-114
###################################################
m1c <- glm( cbind(D,Y/10^5) ~ I(A-30), family=poisreg, data=tdk )
round( ci.exp( m1c ), 3 )


###################################################
### code chunk number 15: lin
###################################################
with( tdk, plot( A, D/Y*10^5, pch=16, cex=0.3, col="gray", log="y" ) )
nd <- data.frame( A=seq(15,65,0.5) )
pr <- ci.pred( m1c, newdata=nd )
matshade( nd$A, pr, plot=FALSE, lwd=2 )


###################################################
### code chunk number 16: quant.rnw:216-218
###################################################
m2 <- glm( cbind(D,Y/10^5) ~ A + I(A^2), family=poisreg, data=tdk )
round( ci.lin( m2 ), 4 )


###################################################
### code chunk number 17: pr2
###################################################
pr2 <- ci.pred( m2, newdata=nd )
matshade( nd$A, cbind(pr2,pr), plot=TRUE )


###################################################
### code chunk number 18: inc123
###################################################
m3 <- glm( cbind(D,Y/10^5) ~ A + I(A^2) + I(A^3), family=poisreg, data=tdk )
pr3 <- ci.pred( m3, newdata=nd )
matshade( nd$A, cbind(pr3,pr2,pr), plot=TRUE, 
          lty=c("solid","32","32"), lend="butt", lwd=2, col=1, 
          log="y", ylim=c(1,20), xlab="Age at FU (years)", 
          ylab="Testis canxer incidence per 100,000 PY" )


###################################################
### code chunk number 19: s32
###################################################
kk <- quantile( rep(tdk$A,tdk$D), (1:5-0.5)/5 )
xtabs( D ~ cut(A,c(14,kk,65)), data=tdk )
ms <- glm( cbind(D,Y/10^5) ~ Ns(A,knots=kk ), family=poisreg, data=tdk )
round( ci.lin( ms ), 3 )
prs <- ci.pred( ms, newdata=nd )
matshade( nd$A, cbind(prs,pr3,pr2,pr), plot=TRUE, 
          col=gray(c(4,0)/10), lty=c("solid","32","32","32"), lend="butt", 
          log="y", ylim=c(1,20), lwd=c(4,2,2,2) )
rug( kk, lwd=2 )


###################################################
### code chunk number 20: gam-glm
###################################################
library( mgcv )
mg <- gam( cbind(D,Y/10^5) ~ s(A,k=15), family=poisreg, data=tdk )
summary( mg )
prg <- ci.pred( mg, newdata=nd )
matshade( nd$A, cbind(prg,prs), plot=TRUE,  
          col=1, lwd=2, lty=c("solid","31"), lend="butt", 
          log="y", ylim=c(1,20), xlab="Age", 
          ylab="Testis cancer incidence per 100,000 PY" )


###################################################
### code chunk number 21: gam-lin
###################################################
mgl <- gam( cbind(D,Y/10^5) ~ s(A,k=15) + P, family=poisreg, data=tdk )
summary( mgl )
round( ci.exp( mgl, subset="P" ), 3 )
nd <- cbind( nd, P=1970 )
prgl <- ci.pred( mgl, newdata=nd )


###################################################
### code chunk number 22: gam-ss
###################################################
mgp <- gam( cbind(D,Y/10^5) ~ s(A,k=15) + s(P,k=5), family=poisreg, data=tdk )
summary( mgp )
anova( mgl, mgp, test="Chisq" )
prgs <- ci.pred( mgp, newdata=nd )
matshade( nd$A, cbind(prgs,prgl), 
          col=gray(c(0,4)/10), lty=c("solid","31"), lend="butt",
          xlab="Age", ylab="Testis cancer incidence per 100,000 PY",
          plot=TRUE, log="y", lwd=c(3,2) )


###################################################
### code chunk number 23: gam-def
###################################################
par( mfrow=c(1,2), mar=c(4,4,0.1,0.1) )
plot( mgp )


###################################################
### code chunk number 24: incRR
###################################################
par( mfrow=c(1,2) )
matshade( nd$A, ci.pred( mgp, nd ), plot=TRUE,
          log="y",  ylim=c(1,20), xlim=13+c(0,53),
          xlab="Age", ylab="Testis cancer incidence 1970 (per 100,00 PY)" )
np <- data.frame( A=30, P=1943:1996 )
nr <- data.frame( A=30, P=1970 )
RR <- ci.exp( mgp, ctr.mat=list(np,nr) )
matshade( np$P, RR, plot=TRUE,
          log="y",  ylim=c(1,20)/4, xlim=1943+c(0,53),
          xlab="Date of FU", ylab="Incidence RR relative to 1970" )
abline( h=1, v=1970, lty=3 )


###################################################
### code chunk number 25: incRRx
###################################################
par( mfrow=c(1,2), mar=c(3,3,0.1,0.1), mgp=c(3,1,0)/1.6 )
matshade( nd$A, ci.pred( mgp, nd ), plot=TRUE,
          log="y",  ylim=c(1,20), xlim=13+c(0,53),
          xlab="Age at FU", ylab="Testis cancer incidence 1970 (per 100,00 PY)" )
np <- data.frame( A=30, P=1943:1996 )
nr <- data.frame( A=30, P=1970 )
RR <- ci.exp( mgp, ctr.mat=list(np,nr) )
matshade( np$P, RR, plot=TRUE,
          log="y",  ylim=c(1,20)/5, xlim=1943+c(0,53),
          xlab="Date of FU", ylab="Incidence RR relative to 1970" )
abline( h=1, v=1970, lty=3 )


###################################################
### code chunk number 26: quant.rnw:632-634
###################################################
tdk$B <- tdk$P - tdk$A
range( tdk$B )


###################################################
### code chunk number 27: quant.rnw:640-644
###################################################
mgc <- gam( cbind(D,Y/10^5) ~ s(A,k=15) + s(B,k=40), family=poisreg, data=tdk )
summary( mgc )
anova( mgp, mgl, mgc, test="Chisq" )
c( mgp$aic, mgl$aic, mgc$aic )


###################################################
### code chunk number 28: incRRc
###################################################
par( mar=c(3,3,0.1,0.1), mgp=c(3,1,0)/1.6 )
layout( rbind(1:2), widths=c(1,1.7) )
na  <- data.frame( A=15:65, B=1930 )
matshade( na$A, ci.pred( mgc, na ), plot=TRUE,
          log="y",  ylim=c(1,20), xlim=13+c(0,53), xlab="Age at FU", 
          ylab="Testis cancer incidence: 1930 cohort (per 100,000 PY)" )
nc  <- data.frame( A=30, B=1880:1980 )
ncr <- data.frame( A=30, B=1930 )
RRc <- ci.exp( mgc, ctr.mat=list(nc,ncr) )
matshade( nc$B, RRc, plot=TRUE,
          log="y",  ylim=c(1,20)/5, xlim=1878+c(0,106),
          xlab="Date of birth", ylab="Incidence RR relative to 1930 cohort" )
abline( h=1, v=c(1930,1918+10/12,1945+5/12), lty=3 )


###################################################
### code chunk number 29: m-cont
###################################################
nd <- expand.grid( A=15:65, P=1943:1998 )
nd$B <- nd$P-nd$A
nd$Prate <- ci.pred( mgp, nd )[,1]
nd$Brate <- ci.pred( mgc, nd )[,1]
head( nd, 4 )
library( ggplot2 )
library( gridExtra )
library( ggthemes )
cp <- ggplot( nd, aes(P, A, z=log(Prate)) ) +
              geom_contour( binwidth=log(1.1), col=1 ) +
              theme_tufte( base_size  =24, 
                           base_family="Helvetica")
cb <- ggplot( nd, aes(P, A, z=log(Brate)) ) +
              geom_contour( binwidth=log(1.1), col=1 ) +
              theme_tufte( base_size  =24, 
                           base_family="Helvetica") +
              geom_abline( intercept=-c(1945+5/12,1918+10/12),
                               slope=1,
                               color=gray(0.4)) + 
              geom_label( aes(label=c("1918","1945")),
                          data=data.frame(P=c(1983,1998),
                                          A=c(  65,  53),
                                      Brate=NA) )
grid.arrange(cp, cb, nrow=1)


###################################################
### code chunk number 30: quant.rnw:798-802
###################################################
miap <- gam( cbind(D,Y/10^5) ~ s(A,P,k=50), family=poisreg, data=tdk )
summary(miap)
anova( mgp, miap, mgc, test="Chisq" )
c( ap=mgp$aic, api=miap$aic, ac=mgc$aic )


###################################################
### code chunk number 31: intp1
###################################################
yr <- seq(1950,1990,5) 
cl <- gray(seq(.7,.1,length=9))
np <- data.frame( A=seq(15,65,0.5), P=yr[1] )
plot( NA, log="y", ylim=c(0.5,30), xlim=c(15,65),
          xlab="Age at FU", ylab="Incidence rate per 100,000 PY" )
for( i in 1:length(yr) )
   {
np$P <- yr[i]
matshade( np$A, ci.pred(miap,np), col=cl[i] )
   }
text(65,30*0.85^(8:0),yr,adj=1,col=cl)
wh.a <- c(30:34,30:33)
text(wh.a,ci.pred(miap,data.frame(A=wh.a,P=yr))[,1],yr,cex=0.6)


###################################################
### code chunk number 32: intp2
###################################################
ag <- seq(15,65,10) 
cl <- gray(seq(.7,.1,length=6))
np <- data.frame( P=seq(1945,1995,0.5), A=ag[1] )
plot( NA, log="y", ylim=c(0.5,30), xlim=c(1945,1995),
          xlab="Date of FU", ylab="Incidence rate per 100,000 PY" )
for( i in 1:6 )
   {
np$A <- ag[i]
matshade( np$P, ci.pred(miap,np), col=cl[i] )
   }
text(1945,30*0.85^(5:0),ag,adj=0,col=cl)
wh.p <- 1984
text(wh.p,ci.pred(miap,data.frame(A=ag,P=wh.p))[,1],ag,cex=0.5)


###################################################
### code chunk number 33: quant.rnw:962-966
###################################################
miac <- gam( cbind(D,Y/10^5) ~ s(A,B,k=80), family=poisreg, data=tdk )
summary(miac)
anova( mgp, miac, mgc, test="Chisq" )
c( ap=mgp$aic, iap=miap$aic, iac=miac$aic, ac=mgc$aic )


###################################################
### code chunk number 34: intc1
###################################################
bc <- seq(1900,1970,10) 
cl <- gray(seq(.7,.1,length=8))
plot( NA, log="y",  ylim=c(0.5,30), xlim=c(15,65),
          xlab="Age at FU", ylab="Incidence rate per 100,000 PY" )
for( i in 1:8 )
   {
nc <- subset( data.frame( A=seq(15,65,0.5), B=bc[i] ),
              A+B>1943 & A+B<1997 )
nc$P <- nc$A + nc$B
matshade( nc$A, ci.pred(miac,nc), col=cl[i] )
matshade( nc$A, ci.pred(miap,nc), col=cl[i], lty="44",lend="butt" )
   }
text(65,30,"Date of birth",adj=1)
text(65,30*0.85^(8:1),bc,adj=1,col=cl)


###################################################
### code chunk number 35: intc2
###################################################
ag <- seq(15,65,10) 
cl <- gray(seq(.7,.1,length=6))
plot( NA, log="y", ylim=c(0.5,30), xlim=c(1880,1980),
          xlab="Date of birth", ylab="Incidence rate per 100,000 PY" )
for( i in 1:6 )
   {
nb <- subset( data.frame( B=1880:1980, A=ag[i] ),
              A+B>1943 & A+B<1998 )
matshade( nb$B, ci.pred(miac,nb), col=cl[i] )
   }
text(1885,30,"Age",adj=1)
text(1885,30*0.85^(6:1),ag,adj=1,col=cl)
wh.p <- 1988
text(wh.p-ag,ci.pred(miac,data.frame(A=ag,B=wh.p-ag))[,1],ag,cex=0.5)


###################################################
### code chunk number 36: bi-cont
###################################################
nd$Irate <- ci.pred( miac, nd )[,1]
head( nd, 4 )
ci <- ggplot( nd, aes(P, A, z=log(Irate)) ) +
              geom_contour( binwidth=log(1.1), col=1 ) +
              theme_tufte( base_size=24, 
                         base_family="Helvetica") +
              geom_abline( intercept=-c(1945+5/12,1918+10/12),
                               slope=1,
                               color=gray(0.4)) + 
              geom_label( aes(label=c("1918","1945")),
                          data=data.frame(P=c(1983,1998),
                                          A=c(  65,  53),
                                      Irate=NA) ) 
grid.arrange(cb, ci, nrow = 1)


###################################################
### code chunk number 37: quant.rnw:1210-1215
###################################################
mgpv  <- gam( cbind(D,Y/10^5) ~ s(A,k=15) + 
                                s(P,k=15) + 
                                s(A,k=15,by=P), family=poisreg, data=tdk )
summary( mgpv )
anova( mgp, mgpv, miap, test="Chisq" )


###################################################
### code chunk number 38: int1v
###################################################
yr <- seq(1950,1990,5) 
cl <- gray(seq(.8,.1,length=9))
nd <- data.frame( A=seq(15,65,0.5), P=yr[1] )
plot( NA, log="y", ylim=c(0.5,30), xlim=c(15,65),
          xlab="Age at FU", ylab="Incidence rate per 100,000 PY" )
for( i in 1:9 )
   {
nd$P <- yr[i]
matshade( nd$A, ci.pred(mgpv,nd), col=cl[i] )
   }
text(65,30*0.85^(8:0),yr,adj=1,col=cl)
wh.a <- c(30:34,30:33)
text(wh.a,ci.pred(mgpv,data.frame(A=wh.a,P=yr))[,1],yr,cex=0.5)


###################################################
### code chunk number 39: int2v
###################################################
ag <- seq(15,65,10) 
cl <- gray(seq(.7,.1,length=6))
nd <- data.frame( P=seq(1945,1995,0.5), A=ag[1] )
plot( NA, log="y", ylim=c(0.5,30), xlim=c(1945,1995),
          xlab="Date of FU", ylab="Incidence rate per 100,000 PY" )
for( i in 1:6 )
   {
nd$A <- ag[i]
matshade( nd$P, ci.pred(mgpv,nd), col=cl[i], lwd=2 )
   }
text(1945,30*0.85^(5:0),ag,adj=0,col=cl)
wh.p <- 1984
text(wh.p,ci.pred(mgpv,data.frame(A=ag,P=wh.p))[,1],ag,cex=0.5)