### R code from vignette source 'curve-s.rnw'

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


###################################################
### code chunk number 2: curve-s.rnw:15-19
###################################################
library( Epi )
library( splines )
data( DMlate )
head( DMlate )


###################################################
### code chunk number 3: curve-s.rnw:23-29
###################################################
DMlate <- transform( DMlate, D = !is.na(dodth),
                             Y = dox-dodm,
                             A = dodm-dobth,
                             P = dodm )
DMlate <- subset( DMlate, Y>0 & sex=="M" )
str( DMlate )


###################################################
### code chunk number 4: curve-s.rnw:33-34
###################################################
m0 <- glm( D ~ A + P, family=poisson, offset=log(Y), data=DMlate )


###################################################
### code chunk number 5: curve-s.rnw:37-39
###################################################
ci.lin( m0 )
round( ci.exp( m0 )[-1,], 3 )


###################################################
### code chunk number 6: curve-s.rnw:44-46
###################################################
mq <- glm( D ~ A + P + I(P^2), family=poisson, offset=log(Y), data=DMlate )
round( ci.lin( mq ), 3 )


###################################################
### code chunk number 7: rr2005
###################################################
P.pt <- 1995:2010
p1 <- P.pt - 2005
p2 <- P.pt^2 - 2005^2
coef( mq )
lRR <- coef(mq)[3] * p1 + coef(mq)[4] * p2
 RR <- exp( lRR )
plot( P.pt, RR, type="l", lwd=3, log="y" )
abline( h=1,v=2005)


###################################################
### code chunk number 8: rr2000
###################################################
P.pt <- 1995:2010
p1 <- P.pt - 2000
p2 <- P.pt^2 - 2000^2
coef( mq )
lRR <- coef(mq)[3] * p1 + coef(mq)[4] * p2
 RR <- exp( lRR )
plot( P.pt, RR, type="l", lwd=3, log="y" )
abline( h=1,v=2000)


###################################################
### code chunk number 9: curve-s.rnw:103-105
###################################################
ci.lin( mq, subset="P" )
ci.exp( mq, subset="P" )


###################################################
### code chunk number 10: curve-s.rnw:110-112
###################################################
cbind( p1, p2 )
ci.exp( mq, subset="P", ctr.mat=cbind(p1,p2) )


###################################################
### code chunk number 11: RRci
###################################################
matplot( P.pt, ci.exp( mq, subset="P", ctr.mat=cbind(p1,p2) ),
         type="l", col="Black", lty=1, lwd=c(3,1,1), log="y" )
abline( h=1, v=2000 )


###################################################
### code chunk number 12: curve-s.rnw:127-134
###################################################
( MP <- cbind( P.pt, P.pt^2 ) )
( Mr <- cbind( rep(2000,length(P.pt)), 2000^2 ) )
MP-Mr
ci.exp( mq, subset="P", ctr.mat=MP-Mr )
matplot( P.pt, ci.exp( mq, subset="P", ctr.mat=MP-Mr ),
         type="l", col="Black", lty=1, lwd=c(3,1,1), log="y" )
abline( h=1,v=2000)


###################################################
### code chunk number 13: curve-s.rnw:143-145
###################################################
p.kn <- c(1997,2000,2003,2006,2009)
Ns( P.pt, knots=p.kn )


###################################################
### code chunk number 14: curve-s.rnw:149-152
###################################################
ms <- glm( D ~ A + Ns(P,knots=p.kn), family=poisson, offset=log(Y), data=DMlate )
ci.exp( ms )
ci.exp( ms, subset="P" )


###################################################
### code chunk number 15: curve-s.rnw:158-161
###################################################
CP <- Ns( P.pt, knots=p.kn )
Cr <- Ns( rep(2005,length(P.pt)), knots=p.kn )
CP-Cr


###################################################
### code chunk number 16: RRspl
###################################################
RR <- ci.exp( ms, subset="P", ctr.mat=CP-Cr )
matplot( P.pt, RR,
         type="l", col="Black", lty=1, lwd=c(3,1,1), log="y" )
abline( h=1,v=2000)


###################################################
### code chunk number 17: curve-s.rnw:174-176
###################################################
ma <- update( ms, . ~ . + I(A^2) )
round( ci.lin( ma ), 3 )


###################################################
### code chunk number 18: curve-s.rnw:186-187
###################################################
nd <- data.frame( A = 40:85, P=2005, Y=1000 )


###################################################
### code chunk number 19: amort1
###################################################
rate <- ci.pred( ma, newdata=nd )
matplot( nd$A, rate,
         type="l", col="Black", lty=1, lwd=c(3,1,1), log="y" )


###################################################
### code chunk number 20: amort2
###################################################
ci.exp( ma )
Rate <- ci.exp( ma, ctr.mat=cbind(1,40:85,Ns(rep(2005,46),knots=p.kn),(40:85)^2) )*1000
matplot( nd$A, Rate,
         type="l", col="Black", lty=1, lwd=c(3,1,1), log="y" )


###################################################
### code chunk number 21: curve-s.rnw:226-229
###################################################
( a.kn <- 3:9*10 )
mA <- update( ms, . ~ . - A + Ns(A,knots=a.kn) )
round( ci.lin( mA ), 3 )


###################################################
### code chunk number 22: aspl
###################################################
rate <- ci.pred( mA, newdata=nd )
matplot( nd$A, rate,
         type="l", col="Black", lty=1, lwd=c(3,1,1), log="y" )