# Setup/ read in GVHD data
source("http://myweb.uiowa.edu/pbreheny/7210/f19/notes/fun.R")
Data <- read.delim("https://s3.amazonaws.com/pbreheny-data-sets/gvhd.txt")

# Surv
require(survival)
S <- Surv(Data$Time, Data$Status)
class(S)
head(S)
head(S[,1])
head(S[,2])

# survfit
fit <- survfit(S~Data$Group)
fit
summary(fit, time=40)

# Kaplan-Meier curve
fit <- survfit(Surv(Time, Status) ~ Group, Data)
plot(fit)  # Default: kind of ugly
Plot(fit)  # My version: feel free to use if you want

# Plotting options
Plot(fit, mark.time=TRUE)
Plot(fit, xmax=80)
Plot(fit, xscale=365.25)

# Adding an at-risk table
op <- par(oma=c(0,0,1,0))
Plot(fit, xmax=80, xlab='Time on study (Days)', ylab='Probability (GVHD-free)')
nrisk(fit)
par(op)

# Another way to do it
op <- par(mar=c(5,5,4,8))
Plot(fit, xmax=80, xlab='Time on study (Days)', ylab='Probability (GVHD-free)', legend='right')
nrisk(fit)
par(op)

# Yet another way
Plot(fit, xmax=80, xlab='Time on study (Days)', ylab='Probability (GVHD-free)', legend='none')
nrisk(fit)
text(80, 0.55, 'MTX', xpd=TRUE, cex=0.8)
text(80, 0.87, 'MTX+CSP', xpd=TRUE, cex=0.8)

# Median
median(fit)

# survminer
library(survminer)
ggsurvplot(fit, Data, risk.table=TRUE, xlim=c(0, 60), break.x.by=20)
surv_median(fit)
surv_summary(fit, Data)[1:5,] %>% knitr::kable('latex', booktabs=TRUE, digits=3)