## Anemia data
require(survival)
anemia <- read.delim("http://myweb.uiowa.edu/pbreheny/data/anemia.txt")
head(anemia)

## Surv objects
S.gvhd <- with(anemia, Surv(Time_gvhd, Status_gvhd))
S.mort <- with(anemia, Surv(Time, Status))
head(S.mort)

## survfit
Trt <- anemia$Trt
fit <- survfit(S.mort ~ Trt)
fit

## survfit plots
plot(fit, xlab="Days on study", ylab="Survival")
plot(fit, xlab="Days on study", ylab="Survival", mark.time=FALSE)
plot(fit, xlab="Days on study", ylab="Survival", mark.time=FALSE, conf.int=TRUE,
     col=c("red", "blue"))
legend("bottomleft", legend=unique(Trt), col=c("red", "blue"), lty=1)

## survdiff
survdiff(S.mort ~ Trt)
survdiff(S.gvhd ~ Trt)

## npsurv
require(rms)
fit.gvhd <- npsurv(S.gvhd ~ Trt)
fit.mort <- npsurv(S.mort ~ Trt)

## survplot
col <- c("gray","slateblue")
survplot(fit.gvhd, mark.time=FALSE, col=col, xlab="Days on study", ylab="Pr(GVHD-free)",
         xlim=c(0,80), lty=1, conf="none", lwd=3, n.risk=TRUE, y.n.risk=1.05, cex.n.risk=0.7,
         sep.n.risk=0.03, levels.only=TRUE, time.inc=20)
col.fill <- c(rgb(.5,.5,.5,alpha=0.25), rgb(0,0,1,alpha=0.25))
survplot(fit.gvhd, mark.time=FALSE, col=col, col.fill=col.fill, xlab="Days on study",
         ylab="Pr(GVHD-free)", xlim=c(0,80), lty=1, lwd=3, levels.only=TRUE, time.inc=20)
survplot(fit.mort, mark.time=FALSE, col=col, xlab="Days on study", ylab="Survival", lty=1,
         conf="none", lwd=3, n.risk=TRUE, y.n.risk=1.05, cex.n.risk=0.7, sep.n.risk=0.03,
         levels.only=TRUE)
survplot(fit.mort, mark.time=FALSE, col=col, col.fill=col.fill, xlab="Days on study",
         ylab="Survival", lty=1, lwd=3, levels.only=TRUE)

## survdiffplot
survdiffplot(fit.gvhd, xlim=c(0,80), time.inc=20, xlab="Days on Study",
             ylab="Difference in GVHD probability", order=2:1)
survdiffplot(fit.mort, xlab="Days on Study", order=2:1)

