### R code from vignette source 'surv-lung.rnw'
### Encoding: UTF-8

###################################################
### 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:42-48
###################################################
data(lung)
lung$sex <- factor(lung$sex,
                   levels = 1:2,
                   labels = c("M", "W"))
lung$time <- lung$time / (365.25/12)
head(lung)


###################################################
### code chunk number 4: surv-lung.rnw:53-55 (eval = FALSE)
###################################################
## ?Surv
## ?survfit


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


###################################################
### code chunk number 6: surv-lung.rnw:68-69
###################################################
plot(km)


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


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


###################################################
### code chunk number 9: surv-lung.rnw:87-88
###################################################
with(lung, tapply(age, sex, mean))


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


###################################################
### code chunk number 11: surv-lung.rnw:107-112
###################################################
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 12: surv-lung.rnw:129-130 (eval = FALSE)
###################################################
## ?cox.zph


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


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


###################################################
### code chunk number 15: surv-lung.rnw:162-165
###################################################
 plot(prs, col = c("blue", "red"), lwd = 1, lty = 1, conf.int = TRUE)
lines(prs, col = c("blue", "red"), lty = 1, lwd=3)
lines(kms, col = c("blue", "red"), lty = 2, lwd=2)


###################################################
### code chunk number 16: surv-lung.rnw:183-184 (eval = FALSE)
###################################################
## ?poisreg


###################################################
### code chunk number 17: surv-lung.rnw:186-191
###################################################
p1 <- glm(cbind(status == 2, time) ~ sex + age, 
          family = poisreg,
            data = lung)
ci.exp(p1) # estimates form Poisson
ci.exp(c1) # estimates from Cox


###################################################
### code chunk number 18: surv-lung.rnw:205-214 (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 19: surv-lung.rnw:243-244 (eval = FALSE)
###################################################
## ?Lexis


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


###################################################
### code chunk number 21: surv-lung.rnw:274-275
###################################################
summary(Ll)


###################################################
### code chunk number 22: surv-lung.rnw:280-282 (eval = FALSE)
###################################################
## ?boxes
## boxes(Ll, boxpos = TRUE)


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


###################################################
### code chunk number 24: surv-lung.rnw:292-293 (eval = FALSE)
###################################################
## ?Surv


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


###################################################
### code chunk number 26: surv-lung.rnw:302-303 (eval = FALSE)
###################################################
## ?coxph.Lexis


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


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


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


###################################################
### code chunk number 30: surv-lung.rnw:345-346
###################################################
Sl <- splitMulti(Ll, tfl = 0:36)


###################################################
### code chunk number 31: surv-lung.rnw:351-353
###################################################
summary(Ll)
summary(Sl)


###################################################
### code chunk number 32: surv-lung.rnw:357-360
###################################################
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 33: surv-lung.rnw:371-376
###################################################
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 34: surv-lung.rnw:379-382 (eval = FALSE)
###################################################
## ?glm.Lexis
## ps <- glm.Lexis(Sl, ~ Ns(tfl, knots = seq(0, 36, 12)) + sex + age)
## ci.exp(ps)


###################################################
### code chunk number 35: surv-lung.rnw:387-390
###################################################
round(cbind(ci.exp(cl), 
            ci.exp(ps, subset = c("sex","age")), 
            ci.exp(pc, subset = c("sex","age"))), 3)


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


###################################################
### code chunk number 37: surv-lung.rnw:422-425
###################################################
matshade(prf$tfl, ci.pred(ps, prf), 
         plot = TRUE, log = "y", lwd = 3)
matshade(prf$tfl, ci.pred(pc, prf), lty = 2, lwd = 3)


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


###################################################
### code chunk number 39: surv-lung.rnw:447-452
###################################################
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])
lines(survfit(c1, newdata = data.frame(sex = "W", age = 60)), 
      lwd = 2, lty = 1)


###################################################
### code chunk number 40: ratesurv
###################################################
par(mfrow = c(1,2))

# hazard scale
matshade(prf$tfl, ci.pred(ps, prf), 
         plot = TRUE, log = "y", lwd = 3)
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)
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 41: 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)

# 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 42: surv-lung.rnw:511-527
###################################################
zz <- 
function(dk)
{
kn <- seq(0, 36, dk)
pk <- glm(cbind(lex.Xst == "Dead",
                lex.dur) ~ Ns(tfl, knots = kn),
          family = poisreg,
            data = Sl)
matshade(prf$tfl, ci.pred(pk, prf), 
         plot = TRUE, log = "y", lwd = 3, ylim = c(0.01,1))
rug(kn, lwd=3)

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


###################################################
### code chunk number 43: surv-lung.rnw:529-531 (eval = FALSE)
###################################################
## par(mfrow=c(1,2))
## zz(12)


###################################################
### code chunk number 44: knots2
###################################################
par(mfrow=c(4,2))
for (nk in c(6, 4, 3, 2)) zz(nk)