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

###################################################
### code chunk number 1: surv-lung.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: surv-lung.rnw:27-33
###################################################
library(survival)
library(Epi)
library(popEpi)
# popEpi::splitMulti returns a data.frame rather than a data.table
options("popEpi.datatable" = FALSE)
clear()


###################################################
### code chunk number 3: surv-lung.rnw:35-37
###################################################
noquote(version$version.string)
noquote(installed.packages()[c("Epi","popEpi"),c("Version","Built")])


###################################################
### code chunk number 4: surv-lung.rnw:46-52
###################################################
data(lung)
lung$sex <- factor(lung$sex,
                   levels = 1:2,
                   labels = c("M", "W"))
lung$time <- round(lung$time / (365.25/12), 3)
head(lung)


###################################################
### code chunk number 5: surv-lung.rnw:57-59 (eval = FALSE)
###################################################
## ?Surv
## ?survfit


###################################################
### code chunk number 6: surv-lung.rnw:61-64
###################################################
km  <- survfit(Surv(time, status == 2) ~ 1, data = lung)
km
# summary(km) # very long output


###################################################
### code chunk number 7: surv-lung.rnw:72-73
###################################################
plot(km)


###################################################
### code chunk number 8: surv-lung.rnw:78-80
###################################################
kms <- survfit(Surv(time, status == 2) ~ sex, data = lung)
kms


###################################################
### code chunk number 9: kms
###################################################
 plot(kms, col = c("blue", "red"), lwd = 1, conf.int = TRUE)
lines(kms, col = c("blue", "red"), lwd = 3)


###################################################
### code chunk number 10: surv-lung.rnw:91-92
###################################################
with(lung, tapply(age, sex, mean))


###################################################
### code chunk number 11: surv-lung.rnw:95-97 (eval = FALSE)
###################################################
## ?survdiff
## survdiff(Surv(time, status==2) ~ sex, data = lung)


###################################################
### code chunk number 12: surv-lung.rnw:112-117
###################################################
c0 <- coxph(Surv(time, status == 2) ~ sex            , data = lung)
c1 <- coxph(Surv(time, status == 2) ~ sex + I(age/10), data = lung)
summary(c1)
ci.exp(c0)
ci.exp(c1)


###################################################
### code chunk number 13: surv-lung.rnw:132-133 (eval = FALSE)
###################################################
## ?cox.zph


###################################################
### code chunk number 14: surv-lung.rnw:135-138
###################################################
cox.zph(c0)
(z1 <- cox.zph(c1))
par(mfrow = c(1, 2)) ; plot(z1)


###################################################
### code chunk number 15: surv-lung.rnw:160-162
###################################################
prs <- survfit(c0, newdata = data.frame(sex = c("M","W")))
plot(prs, col = c("blue", "red"))


###################################################
### code chunk number 16: surv-lung.rnw:169-172
###################################################
 plot(prs, col = c("blue", "red"), lwd = 1, lty = 1, conf.int = FALSE)
lines(prs, col = c("blue", "red"), lty = 1, lwd=3)
lines(kms, col = c("blue", "red"), lty = "11", lwd = 2, lend = "butt")


###################################################
### code chunk number 17: surv-lung.rnw:198-199 (eval = FALSE)
###################################################
## ?poisreg


###################################################
### code chunk number 18: surv-lung.rnw:201-207
###################################################
with(lung, cbind(status == 2, time)[1:10,])
p1 <- glm(cbind(status == 2, time) ~ sex + I(age/10),
          family = poisreg,
            data = lung)
ci.exp(p1) # estimates form Poisson
ci.exp(c1) # estimates from Cox


###################################################
### code chunk number 19: surv-lung.rnw:222-231 (eval = FALSE)
###################################################
## px <- glm(status == 2  ~ sex + age + offset(log(time)),
##           family = poisson,
##             data = lung)
## ## or:
## px <- glm(status == 2  ~ sex + age,
##           offset = log(time),
##           family = poisson,
##             data = lung)
## ci.exp(px)


###################################################
### code chunk number 20: surv-lung.rnw:262-263 (eval = FALSE)
###################################################
## ?Lexis


###################################################
### code chunk number 21: surv-lung.rnw:266-272
###################################################
Ll <- Lexis(exit = list(tfl = time),
            exit.status = factor(status,
                                 levels = 1:2,
                                 labels = c("Alive","Dead")),
            data = lung)
head(Ll)


###################################################
### code chunk number 22: surv-lung.rnw:296-297
###################################################
summary(Ll)


###################################################
### code chunk number 23: surv-lung.rnw:302-303 (eval = FALSE)
###################################################
## ?boxes


###################################################
### code chunk number 24: surv-lung.rnw:306-307 (eval = FALSE)
###################################################
## boxes(Ll, boxpos = TRUE)


###################################################
### code chunk number 25: surv-lung.rnw:309-310 (eval = FALSE)
###################################################
## boxes(Ll, boxpos = TRUE, scale.Y = 12, digits.R = 2)


###################################################
### code chunk number 26: surv-lung.rnw:317-318 (eval = FALSE)
###################################################
## ?Surv


###################################################
### code chunk number 27: surv-lung.rnw:320-324
###################################################
cl <- coxph(Surv(tfl,
                 tfl + lex.dur,
                 lex.Xst == "Dead") ~ sex + age,
            data = Ll)


###################################################
### code chunk number 28: surv-lung.rnw:327-328 (eval = FALSE)
###################################################
## ?coxph.Lexis


###################################################
### code chunk number 29: surv-lung.rnw:331-334
###################################################
cL <- coxph.Lexis(Ll, tfl ~ sex + age)
ci.exp(cL)
ci.exp(cl)


###################################################
### code chunk number 30: surv-lung.rnw:338-341
###################################################
pc <- glm(cbind(lex.Xst == "Dead", lex.dur) ~ sex + age,
          family = poisreg,
            data = Ll)


###################################################
### code chunk number 31: surv-lung.rnw:344-347
###################################################
pL <- glm.Lexis(Ll, ~ sex + age)
ci.exp(pL)
ci.exp(pc)


###################################################
### code chunk number 32: surv-lung.rnw:371-372
###################################################
Sl <- splitMulti(Ll, tfl = 0:36)


###################################################
### code chunk number 33: surv-lung.rnw:377-379
###################################################
summary(Ll)
summary(Sl)


###################################################
### code chunk number 34: surv-lung.rnw:383-386
###################################################
wh <- names(Ll)[1:10] # names of variables in some order
subset(Ll, lex.id == 10)[,wh]
subset(Sl, lex.id == 10)[,wh]


###################################################
### code chunk number 35: surv-lung.rnw:401-406
###################################################
ps <- glm(cbind(lex.Xst == "Dead", lex.dur)
          ~ Ns(tfl, knots = seq(0, 36, 12)) + sex + age,
          family = poisreg,
            data = Sl)
ci.exp(ps)


###################################################
### code chunk number 36: surv-lung.rnw:409-412 (eval = FALSE)
###################################################
## ?glm.Lexis
## ps <- glm.Lexis(Sl, ~ Ns(tfl, knots = seq(0, 36, 12)) + sex + age)
## ci.exp(ps)


###################################################
### code chunk number 37: surv-lung.rnw:417-422
###################################################
ests <- cbind(ci.exp(cl),
              ci.exp(ps, subset = c("sex","age")),
              ci.exp(pc, subset = c("sex","age")))
colnames(ests)[1:3*3-2] <- c("Cox","Pois-spline","Pois-const")
round(ests, 3)


###################################################
### code chunk number 38: surv-lung.rnw:449-452
###################################################
prf <- data.frame(tfl = seq(0, 30, 0.2),
                  sex = "W",
                  age = 60)


###################################################
### code chunk number 39: surv-lung.rnw:460-462
###################################################
matshade(prf$tfl, ci.pred(ps, prf), lwd = 3, plot = TRUE, log = "y")
matshade(prf$tfl, ci.pred(pc, prf), lwd = 3, lty = "21", lend = "butt")


###################################################
### code chunk number 40: surv-lung.rnw:478-480
###################################################
matshade(prf$tfl, ci.surv(ps, prf, intl = 0.2),
         plot = TRUE, ylim = 0:1, lwd = 3)


###################################################
### code chunk number 41: surv-lung.rnw:487-492
###################################################
matshade(prf$tfl, ci.surv(ps, prf, intl = 0.2),
         plot = TRUE, ylim = 0:1, lwd = 3)
lines(prf$tfl, ci.surv(pc, prf, intl = 0.2)[,1], lwd = 2, col = gray(0.5))
lines(survfit(c1, newdata = data.frame(sex = "W", age = 60)),
      lwd = 2, lty = 1)


###################################################
### code chunk number 42: ratesurv
###################################################
par(mfrow = c(1,2))
# hazard scale
matshade(prf$tfl, ci.pred(ps, prf),
         plot = TRUE, log = "y", lwd = 3, xlim = c(0,30))
matshade(prf$tfl, ci.pred(pc, prf), lty = 3, lwd = 3)
#
# survival
matshade(prf$tfl, ci.surv(ps, prf, intl = 0.2),
         plot = TRUE, ylim = 0:1, lwd = 3, xlim = c(0,30))
matshade(prf$tfl, ci.surv(pc, prf, intl = 0.2),
         lty = 3, alpha = 0, lwd = 3)
lines(survfit(c1, newdata = data.frame(sex = "W", age = 60)),
      col = "forestgreen", lwd = 3)


###################################################
### code chunk number 43: parkm
###################################################
par(mfrow=c(1,2))
pk <- glm(cbind(lex.Xst == "Dead",
                lex.dur) ~ Ns(tfl, knots = seq(0, 36, 12)),
          family = poisreg,
            data = Sl)
#or
pk <- glm.Lexis(Sl, ~ Ns(tfl, knots = seq(0, 36, 12)))

# hazard
matshade(prf$tfl, ci.pred(pk, prf),
         plot = TRUE, log = "y", lwd = 3, ylim = c(0.01,1))

# survival from smooth model
matshade(prf$tfl, ci.surv(pk, prf, intl = 0.2) ,
         plot = TRUE, lwd = 3, ylim = 0:1)

# K-M estimator
lines(km, lwd = 2)


###################################################
### code chunk number 44: surv-lung.rnw:552-566
###################################################
zz <-
function(dk)
{
kn <<- seq(0, 36, dk) # must be in global environment...
print(kn)
pk <<- glm.Lexis(Sl, ~ Ns(tfl, knots = kn))
matshade(prf$tfl, ci.pred(pk, prf),
         plot = TRUE, log = "y", lwd = 2, ylim = c(0.01,1), xlim = c(0,30))
rug(kn, lwd=3)

plot(km, lwd = 2, col = "limegreen", xlim = c(0,30))
matshade(prf$tfl, ci.surv(pk, prf, intl = 0.2) ,
         lwd = 2, ylim = 0:1, xlim = c(0,30))
}


###################################################
### code chunk number 45: surv-lung.rnw:568-574 (eval = FALSE)
###################################################
## par(mfrow=c(1,2))
## zz(18)
## zz(12)
## zz(6)
## zz(4)
## zz(3)


###################################################
### code chunk number 46: knots2
###################################################
par(mfrow=c(5,2), mar = rep(1,4))
for (nk in c(18, 12, 6, 4, 3)) zz(nk)