### R code from vignette source 'lung-ex.rnw'

###################################################
### code chunk number 1: lung-ex.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) ) )


###################################################
### code chunk number 2: lung-ex.rnw:12-17
###################################################
options( width=120 )
library( survival )
library( Epi )
data( lung )
head( lung )


###################################################
### code chunk number 3: lung-ex.rnw:24-29
###################################################
Lx <- Lexis( exit=list( tfd=time+runif(nrow(lung),-0.5,0.5)),
             exit.status=(status==2),
             data=lung )
summary( Lx, scale=365.25 )
head( Lx )


###################################################
### code chunk number 4: lung-ex.rnw:36-39
###################################################
Sx <- splitLexis( Lx, "tfd", breaks=c(0,unique(exit(Lx))) )
summary( Sx, scale=365.25 )
subset( Sx, lex.id==96 )


###################################################
### code chunk number 5: lung-ex.rnw:46-57
###################################################
c1 <- coxph( Surv(time,status==2) ~ sex + pat.karno, data=lung )
c2 <- coxph( Surv(tfd,tfd+lex.dur,lex.Xst==TRUE) ~ sex + pat.karno, data=Lx )
p1 <- glm( lex.Xst ~ factor(tfd) + sex + pat.karno,
           offset = log(lex.dur), family=poisson,
           data=Sx )
p2 <- glm( lex.Xst ~ ns(tfd,df=6) + sex + pat.karno,
           offset = log(lex.dur), family=poisson,
           data=Sx )
p3 <- glm( lex.Xst ~ ns(tfd,df=2) + sex + pat.karno,
           offset = log(lex.dur), family=poisson,
           data=Sx )


###################################################
### code chunk number 6: lung-ex.rnw:65-75
###################################################
k7 <- c( 0, quantile( rep(Sx$tfd,Sx$lex.Xst), (1:7-0.5)/7 ) )
k3 <- c( 0, quantile( rep(Sx$tfd,Sx$lex.Xst), (1:3-0.5)/3 ) )
xtabs( lex.Xst ~ cut(tfd,breaks=c(k7,Inf)), data=Sx )
xtabs( lex.Xst ~ cut(tfd,breaks=c(k3,Inf)), data=Sx )
p2 <- glm( lex.Xst ~ Ns(tfd,knots=k7) + sex + pat.karno,
           offset = log(lex.dur), family=poisson,
           data=Sx )
p3 <- glm( lex.Xst ~ Ns(tfd,knots=k3) + sex + pat.karno,
           offset = log(lex.dur), family=poisson,
           data=Sx )


###################################################
### code chunk number 7: lung-ex.rnw:82-88
###################################################
ee <- rbind( ci.exp( c1 ), ci.exp( c2 ),
             ci.exp( p1, subset=c("sex","pat") ),
             ci.exp( p2, subset=c("sex","pat") ),
             ci.exp( p3, subset=c("sex","pat") ) )
wh <- 1:5*2
round( cbind( ee[wh-1,], ee[wh,] ), 4 )


###################################################
### code chunk number 8: rates
###################################################
range( Sx$tfd )
nd <- data.frame( tfd=0:1000, lex.dur=36525,
                  pat.karno=80, sex=1 )
pr2 <- predict( p2, newdata=nd, se.fit=TRUE, type="link" )
pr3 <- predict( p3, newdata=nd, se.fit=TRUE, type="link" )
pr2 <- exp( cbind(pr2$fit,pr2$se.fit) %*% ci.mat() )
pr3 <- exp( cbind(pr3$fit,pr3$se.fit) %*% ci.mat() )
matplot( nd$tfd, cbind( pr2, pr3 ),
         type="l", lty=1, lwd=c(4,1,1), col=rep(c("blue","red"),each=3),
         log="y", xlab="Days since diagnosis",
         ylab="Rate per 100 PY (Man, Karnofsky 80)" )
rug( k7, lwd=2, col="blue", ticksize=0.04 )
rug( k3, lwd=4, col="red" , ticksize=0.02 )


###################################################
### code chunk number 9: lung-ex.rnw:110-119
###################################################
pr.exp <-
function ( object, newdata, alpha=0.05, ... )
{
pr <- predict.glm( object,
                  newdata = newdata,
                     type = "link",
                   se.fit = TRUE, ... )
exp( cbind(pr$fit,pr$se.fit) %*% ci.mat(alpha=alpha) )
}