library(survival)
source("http://myweb.uiowa.edu/pbreheny/7210/f19/notes/fun.R")
Podz <- read.delim("https://s3.amazonaws.com/pbreheny-data-sets/Podzielinski2012.txt")

# Slide 14
table(Podz$sDeath)
fit <- survfit(Surv(tDeath/12, sDeath) ~ 1, Podz)
head(fit$pstate)
P <- fit$pstate
dimnames(P) <- list(round(fit$time, 3), c(levels(Podz$sDeath)[-1], 'Alive'))
head(P)
head(fit$lower)

# Slide 16
Plot(fit, xlab='Time (Years)', ylab='Cumulative incidence', ylim=0:1, legend="none")
lines(fit$time, 1-fit$pstate[,1], type='s', lwd=3, col='gray70')
toplegend(legend=c('Cancer', 'Other', 'Total'), col=c(pal(2), 'gray70'), lwd=3)

# Slide 17
Plot(fit, xlab='Time (Years)', ylab='Cumulative incidence', ylim=0:1, legend="none")
ciband(fit)
toplegend(legend=c('Cancer', 'Other'), col=pal(2), lwd=3)

# Estimating conditional prevalence: Initial work
OS <- with(Podz, Surv(tDeath/12, sDeath!='Censored'))
PFS <- OS
ind <- with(Podz, Podz$tRecur <= Podz$tDeath & Podz$sRecur==1)
PFS[ind,1] <- Podz$tRecur[ind]/12
PFS[ind,2] <- 1
fit.os <- survfit(OS ~ 1)
fit.pfs <- survfit(PFS ~ 1)

# Slide 20
t <- unique(sort(c(fit.os$time,fit.pfs$time)))
Time <- t[t <= 7]
s.os <- summary(fit.os, Time)$surv
s.pfs <- summary(fit.pfs, Time)$surv
plot(Time, 1-s.pfs/s.os, type="s", ylim=c(0,1), las=1, bty='n', xlim=c(0,5), lwd=2,
     ylab="Prevalence of recurrence among surviving patients", xlab="Years")

# Bootstrap (Slide 21)
B <- 1000
Q <- matrix(NA, nrow=B, ncol=length(Time))
n <- nrow(OS)
for (i in 1:B) {
  ind <- sample(n, replace=TRUE)
  S1 <- OS[ind,]
  S2 <- PFS[ind,]
  s1 <- summary(survfit(S1~1), Time, extend=TRUE)$surv
  s2 <- summary(survfit(S2~1), Time, extend=TRUE)$surv
  Q[i,] <- 1-s2/s1
}

# Slide 22
lower <- apply(Q, 2, quantile, p=.025, na.rm=TRUE)
upper <- apply(Q, 2, quantile, p=.975, na.rm=TRUE)
plot(Time, 1-s.pfs/s.os, type="s", ylim=c(0,1), las=1, bty='n', xlim=c(0,5),
     ylab="Prevalence of recurrence among surviving patients", xlab="Years")
polygon.step(Time, lower[-n], upper[-n], col="gray85")
lines(Time, 1-s.pfs/s.os, type="s", lwd=2)