---
title: "Introduction to optband"
author: "Tom Chen & Sam Tracy"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{opt.ci}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

````{r, echo=FALSE, message=FALSE}
library(utils)
library(stats)
library(LambertW)
library(survival)
library(km.ci)
library(optband)
````

Classical simultaneous confidence bands for survival functions (Hall and Wellner, 1980; Nair, 1984) are derived from transformations of the convergence to a Weiner process with a strictly-increasing variance function These transformations are often motivated by factors such as:

- Time intervals for which the analyst prefers to be thinner or wider
- Tractability of computing critical values of the asymptotic distribution in attaining the prescribed (nominal) coverage level

This package instead approaches the problem purely from an optimization perspective: given a certain coverage level, obtain bands such that the area between is minimized. While an exact solution cannot be obtained in closed form, `optband` provides an approximate solution based off local time arguments for both the survival and cumulative-hazard functions, generalizing the results specified by Kendall et al. (2007).

## Usage

`opt.ci` takes a `survfit` object from the `survival` package with the desired $1-\alpha$ coverage level, function of interest (either `'surv'` for the survival function or `'cumhaz'` for the cumulative-hazard function), and optional upper or lower bounds for data truncation. Defaults are $\alpha=0.05$, `fun = 'surv'`, `tl = NA`, `tu = NA`.

## Example

First, let's generate some survival data using the `stats` package and estimate the Kaplan-Meier curve:

```{r, echo=TRUE}
set.seed(1990)
N = 200
x1 <- stats::rweibull(N, 1, 1)
x2 <- stats::rweibull(N, 2, 1)
d <- x1 < x2
x <- pmin(x1,x2)
mydata = data.frame(stop = x, event = d)
S = survival::survfit(Surv(x, d) ~ 1, type="kaplan-meier")
```

Now we can estimate the optimized confidence bands and plot the curve:

```{r, fig.show='hold', fig.height=6, fig.width=6, echo=TRUE}
opt_S <- optband::opt.ci(S, conf.level = 0.95, fun = "surv", tl = NA, tu = NA)
plot(opt_S, xlab="time", ylab="KM", mark.time=FALSE)
```

And we can do the same with the estimated cumulative-hazard function:

```{r, fig.show='hold', fig.height=6, fig.width=6, echo=TRUE}
opt_H <- optband::opt.ci(S, conf.level = 0.95, fun = "cumhaz", tl = NA, tu = NA)
plot(opt_H, fun="cumhaz", xlab="time", ylab="CH", mark.time=FALSE)
```

We can further play with the procedure by adjusting the coverage level and also truncating the data:

```{r, fig.show='hold', fig.height=6, fig.width=6, echo=TRUE}
opt_H <- optband::opt.ci(S, conf.level = 0.90, fun = "cumhaz", tl = .1, tu = .9)
plot(opt_H, fun="cumhaz", xlab="time", ylab="CH", mark.time=FALSE)
```

And compare the results of this band to the Equal Precision and Hall-Wellner bands (for this we will use the `km.ci` package):

````{r, fig.height=6, fig.width=6, echo=TRUE}
color <- c("grey", "darkblue", "red", "green")
plot(S, mark.time=FALSE, xlab="time", ylab="KM", col = "white")

lines(opt_S, col=color[2], lty=2, lwd=2, mark.time=F)
e <- km.ci::km.ci(S, conf.level = 0.95, method = "epband")
h <- km.ci::km.ci(S, conf.level = 0.95, method = "hall-wellner")
lines(e, col=color[3], lty=3, lwd=2, mark.time=F)
lines(h, col=color[4], lty=4, lwd=2, mark.time=F)
lines(S,col="grey", lwd=2, mark.time=F)

legend("topright", c("KM", "optband", "epband", "hall-wellner"), 
       lwd=2, lty=1:4, col=color)
````

Lastly, while the 2-sample survival function difference confidence bands are intractable with this method, we can estimate confidence bands for the cumulative hazard function difference. For this, we'll use the `bladder` data set from the `survival` package, first subsetting upon only the times until first recurrence:

````{r, fig.height=6, fig.width=6, echo=TRUE}
dat <- bladder[bladder$enum==1,]
Hdif = survival::survfit(Surv(stop, event) ~ rx, type="kaplan-meier", data=dat)
opt_Hdif <- optband::opt.ci(Hdif, fun="cumhaz", conf.level = 0.95, samples=2)

plot(opt_Hdif$difference~opt_Hdif$time, xlab="t", ylim=c(-1.5,1),
     main=expression(hat(Lambda)(t)[1]-hat(Lambda)(t)[2]), ylab="", col = "white")
lines(opt_Hdif$difference~opt_Hdif$time, col = "grey", lty=1, lwd=2)
lines(-log(opt_Hdif$upper)~opt_Hdif$time, col=4, lty=3, lwd=2)
lines(-log(opt_Hdif$lower)~opt_Hdif$time, col=4, lty=3, lwd=2)
````

## References

Hall, W. and Wellner, J.(1980).Confidence bands for a survival curve from censored data. *Biometrika*, **67**(1):233-143.

Kendall, W., Marin, J.,and Robert, C. (2007).Confidence bands for brownian motion and applications to monte carlo simulation. *Statistics and Computing*, **17**(1):1–10.

Nair, V. (1984). Confidence bands for survival functions with censored data: a comparative study. *Technometrics*, **26**(3):265–275.