---
title: "Survival Data"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Survival Data}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r chunkname, echo=-1}
data.table::setDTthreads(2)
```


```{r, echo = FALSE, message = FALSE}
library(simstudy)
# library(ggplot2)
library(scales)
library(grid)
library(gridExtra)
library(survival)
library(gee)
library(data.table)
```

```{r, echo = FALSE}
plotcolors <- c("#B84226", "#1B8445", "#1C5974")

cbbPalette <- c("#B84226","#B88F26", "#A5B435", "#1B8446",
                "#B87326","#B8A526", "#6CA723", "#1C5974") 

ggtheme <- function(panelback = "white") {
  
  ggplot2::theme(
    panel.background = ggplot2::element_rect(fill = panelback),
    panel.grid = ggplot2::element_blank(),
    axis.ticks =  ggplot2::element_line(colour = "black"),
    panel.spacing = ggplot2::unit(0.25, "lines"),  # requires package grid
    panel.border = ggplot2::element_rect(fill = NA, colour="gray90"), 
    plot.title = ggplot2::element_text(size = 8,vjust=.5,hjust=0),
    axis.text = ggplot2::element_text(size=8),
    axis.title = ggplot2::element_text(size = 8)
  )  
  
}
```

Time-to-event data, including both survival and censoring times, are created using functions `defSurv` and `genSurv`. The survival data definitions require a variable name as well as a specification of a scale value, which determines the mean survival time at a baseline level of covariates (i.e. all covariates set to 0). The Weibull distribution is used to generate these survival times. In addition, covariates (which have been defined previously) that influence survival time can be included in the `formula` field. Positive coefficients are associated with longer survival times (and lower hazard rates). Finally, the *shape* of the distribution can be specified. A `shape` value of 1 reflects the *exponential* distribution. As of `simstudy` version 0.5.0, it is also possible to generate survival data that violate a proportional hazards assumption. In addition, data with two or more competing risks can be generated.

### Weibull distribution

The density, mean, and variance of the Weibull distribution that is used in the data generation process are defined by the parameters $\lambda$ (scale) and $\nu$ (shape) as shown below.

\begin{aligned}
f(t) &= \frac{t^{\frac{1}{\nu}-1}}{\lambda \nu} exp\left(-\frac{t^\frac{1}{\nu}}{\lambda}\right) \\
E(T) &= \lambda ^ \nu \Gamma(\nu + 1) \\
Var(T) &= (\lambda^2)^\nu \left( \Gamma(2 \nu + 1) - \Gamma^2(\nu + 1) \right) \\
\end{aligned}

<br>

The survival time $T$ data are generated based on this formula:

$$
T = \left( -\frac{log(U) \lambda}{exp(\beta ^ \prime x)}  \right)^\nu,
$$

where $U$ is a uniform random variable between 0 and 1, $\beta$ is a vector of parameters in a Cox proportional hazard model, and $x$ is a vector of covariates that impact survival time. $\lambda$ and $\nu$ can also vary by covariates.

### Generating standard survival data with censoring

Here is an example showing how to generate data with covariates. In this case the scale and shape parameters will vary by group membership.

```{r, tidy = TRUE}

# Baseline data definitions

def <- defData(varname = "x1", formula = .5, dist = "binary")
def <- defData(def,varname = "grp", formula = .5, dist = "binary")

# Survival data definitions

set.seed(282716)

sdef <- defSurv(varname = "survTime", formula = "1.5*x1", scale = "grp*50 + (1-grp)*25", shape = "grp*1 + (1-grp)*1.5")
sdef <- defSurv(sdef, varname = "censorTime", scale = 80, shape = 1)

sdef

```

The data are generated with calls to `genData` and `genSurv`:

```{r, tidy = TRUE}

# Baseline data definitions

dtSurv <- genData(300, def)
dtSurv <- genSurv(dtSurv, sdef)

head(dtSurv)

# A comparison of survival by group and x1

dtSurv[,round(mean(survTime),1), keyby = .(grp,x1)]

```

Observed survival times and censoring indicators can be generated using the competing risk functionality and specifying a censoring variable:

```{r, tidy = TRUE}
dtSurv <- genData(300, def)
dtSurv <- genSurv(dtSurv, sdef, timeName = "obsTime", 
            censorName = "censorTime", eventName = "status", 
            keepEvents = TRUE)

head(dtSurv)

# estimate proportion of censoring by x1 and group

dtSurv[,round(1-mean(status),2), keyby = .(grp,x1)]
```


Here is a Kaplan-Meier plot of the data by the four groups:


```{r, tidy = TRUE, echo = FALSE, fig.width = 6.5, fig.height = 3.5, warning=FALSE}
fit <- survfit(Surv(obsTime, status) ~ x1+grp, data=dtSurv)

survminer::ggsurvplot(fit, data = dtSurv,
  palette = cbbPalette,
  size = .5,
  ggtheme = ggtheme("grey94")
    # ggplot2::theme(axis.title = ggplot2::element_text(size = 9),
    #                        panel.grid = ggplot2::element_blank())
)
```

Here is a survival analysis (using a Cox proportional hazard model) of a slightly simplified data set with two baseline covariates only:

```{r, tidy = TRUE}

# Baseline data definitions

def <- defData(varname = "x1", formula = .5, dist = "binary")
def <- defData(def,varname = "x2", formula = .5, dist = "binary")

# Survival data definitions

sdef <- defSurv(varname = "survTime", formula = "1.5*x1 - .8*x2", scale = 50, shape = 1/2)
sdef <- defSurv(sdef, varname = "censorTime", scale = 80, shape = 1)

dtSurv <- genData(300, def)
dtSurv <- genSurv(dtSurv, sdef, timeName = "obsTime", 
            censorName = "censorTime", eventName = "status")

coxfit <- survival::coxph(Surv(obsTime, status) ~ x1 + x2, data = dtSurv)
```

The 95\% confidence intervals of the parameter estimates include the values used to generate the data:

```{r, echo=FALSE}
gtsummary::tbl_regression(coxfit)
```

### Competing risks

In the previous example, we actually used the competing risk mechanism in `genSurv` to generate an observed time variable (which was the earliest of the censoring and event time). This is done by specifying a *timeName* argument that will represent the observed time value. The event status is indicated in the field set by the *eventName* argument (which defaults to "event"). If a variable name is indicated in the *censorName* argument, the censored events automatically have a value of 0. As we saw above, competing risk information can be generated as part of `genSurv`. However, there is an additional function `addCompRisk` that will generate the competing risk information using an existing data set. The example here will take that approach.

```{r}
d1 <- defData(varname = "x1", formula = .5, dist = "binary")
d1 <- defData(d1, "x2", .5, dist = "binary")

dS <- defSurv(varname = "event_1", formula = "-10 - 0.6*x1 + 0.4*x2", shape = 0.3)
dS <- defSurv(dS, "event_2", "-6.5 + 0.3*x1 - 0.5*x2", shape = 0.5)
dS <- defSurv(dS, "censor", "-7", shape = 0.55)

dtSurv <- genData(1001, d1)
dtSurv <- genSurv(dtSurv, dS)

dtSurv
```

```{r}
dtSurv <- addCompRisk(dtSurv, events = c("event_1", "event_2", "censor"), 
            timeName = "time", censorName = "censor")
dtSurv
```

The competing risk data can be plotted using the cumulative incidence functions (rather than the survival curves):

```{r, tidy = TRUE, echo = FALSE, fig.width = 6.5, fig.height = 3.5, warning=FALSE}
fit <- survfit(Surv(time, event, type="mstate") ~ 1, data=dtSurv)
survminer::ggcompetingrisks(fit, ggtheme = ggtheme("grey94"))  + 
  ggplot2::scale_fill_manual(values = cbbPalette)
```

The data generation can all be done in two (instead of three) steps:

```{r}
dtSurv <- genData(101, d1)
dtSurv <- genSurv(dtSurv, dS, timeName = "time", censorName = "censor")
dtSurv
```

### Introducing non-proportional hazards

In the standard `simstudy` data generation process for survival/time-to-event outcomes that includes covariates that effect the hazard rate at various time points, the ratio of hazards comparing different levels of a covariate are constant across all time points. For example, if we have a single binary covariate $x$, the hazard $\lambda(t)$ at time $t$ is

$$\lambda(t|x) = \lambda_0(t) e ^ {\beta x}$$ 
where $\lambda_0(t)$ is a baseline hazard when $x=0$. The ratio of the hazards for $x=1$ compared to $x=0$ is

$$\frac{\lambda_0(t) e ^ {\beta}}{\lambda_0(t)} = e ^ \beta,$$

so the log of the hazard ratio is a constant $\beta$, and the hazard ratio is always $e^\beta$.

However, we may not always want to make the assumption that the hazard ratio is constant over all time periods. To facilitate this, it is possible to specify two different data definitions for the same outcome, using the *transition* field to specify when the second definition replaces the first. (While it would theoretically be possible to generate data for more than two periods, the process is more involved, and has not been implemented at this time.)

<br>

**Constant/proportional hazard ratio**

To start, here is an example assuming a constant log hazard ratio of -0.7:

```{r, tidy = TRUE}
def <- defData(varname = "x", formula = .4, dist="binary")

defS <- defSurv(varname = "death", formula = "-14.6 - 0.7*x", shape = .35)
defS <- defSurv(defS, varname = "censor", scale = exp(13), shape = .5)

dd <- genData(500, def)
dd <- genSurv(dd, defS, digits = 2, timeName = "time", censorName = "censor")

fit <- survfit( Surv(time, event) ~ x, data = dd )
```

```{r, tidy = TRUE, echo = FALSE, fig.width = 6.5, fig.height = 3.5, warning=FALSE}
survminer::ggsurvplot(fit, data = dd, 
                      ggtheme = ggtheme("grey94"),
                      palette = cbbPalette
)
```

The Cox proportional hazards model recovers the correct log hazards rate:

```{r}
coxfit <- coxph(formula = Surv(time, event) ~ x, data = dd)
```

```{r, echo=FALSE}
gtsummary::tbl_regression(coxfit)
```

We can test the assumption of proportional hazards using weighted residuals. If the $\text{p-value} < 0.05$, then we would conclude that the assumption of proportional hazards is not warranted. In this case $p = 0.22$, so the model is apparently reasonable:

```{r}
cox.zph(coxfit)
```

<br>

**Non-constant/non-proportional hazard ratio**

In this next case, the risk of death when $x=1$ is lower at all time points compared to when $x=0$, but the relative risk (or hazard ratio) changes at 150 days:

```{r, tidy = TRUE}
def <- defData(varname = "x", formula = .4, dist="binary")

defS <- defSurv(varname = "death", formula = "-14.6 - 1.3*x", shape = .35, transition = 0)
defS <- defSurv(defS, varname = "death", formula = "-14.6 - 0.4*x", shape = .35, transition = 150)
defS <- defSurv(defS, varname = "censor", scale = exp(13), shape = .5)

dd <- genData(500, def)
dd <- genSurv(dd, defS, digits = 2, timeName = "time", censorName = "censor")

fit <- survfit( Surv(time, event) ~ x, data = dd )
```

The survival curve for the sample with $x=1$ has a slightly different shape under this data generation process compared to the previous example under a constant hazard ratio assumption; there is more separation early on (prior to day 150), and then the gap is closed at a quicker rate.

```{r, tidy = TRUE, echo = FALSE, fig.width = 6.5, fig.height = 3.5, warning=FALSE}
survminer::ggsurvplot(fit, data = dd, 
                      ggtheme = ggtheme("grey94"),
                      palette = cbbPalette
)
```

If we ignore the possibility that there might be a different relationship over time, the Cox proportional hazards model gives an estimate of the log hazard ratio quite close to -0.70:

```{r}
coxfit <- survival::coxph(formula = Surv(time, event) ~ x, data = dd)
```


```{r, echo=FALSE}
gtsummary::tbl_regression(coxfit)
```

However, further inspection of the proportionality assumption should make us question the appropriateness of the model. Since $p<0.05$, we would be wise to see if we can improve on the model.

```{r}
cox.zph(coxfit)
```

We might be able to see from the plot where proportionality diverges, in which case we can split the data set into two parts at the identified time point. (In many cases, the transition point or points won't be so obvious, in which case the investigation might be more involved.) By splitting the data at day 150, we get the desired estimates:

```{r}
dd2 <- survSplit(Surv(time, event) ~ ., data= dd, cut=c(150),
                 episode= "tgroup", id="newid")

coxfit2 <- survival::coxph(Surv(tstart, time, event) ~ x:strata(tgroup), data=dd2)
```

```{r, echo=FALSE}
gtsummary::tbl_regression(coxfit2)
```

And the diagnostic test of proportionality confirms the appropriateness of the model:

```{r}
cox.zph(coxfit2)
```

### Generating parameters for survival distribution

Throughout this vignette, I have been using various assumptions for the parameters - *formula*, *scale*, and *shape* - that define the Weibull-based survival distribution. Where do these assumptions come from and how can we determine what is appropriate to use in our simulations? That will depend, of course, on each specific application and use of the simulation, but there are two helper functions in `simstudy`, `survGetParams` and `survParamPlot`, that are intended to guide the process.

`survGetParams` will provide the *formula* and *shape* parameters (the *scale* parameter will always be set to 1) that define a curve close to points provided as inputs. For example, if we would like to find the parameters for a distribution where 80% survive until day 100, and 10% survive until day 200 (any number of points may be provided):

```{r}
points <- list(c(100, 0.80), c(200, 0.10))
r <- survGetParams(points)
r
```

We can visualize the curve that is defined by these parameters:

```{r, tidy = TRUE, fig.width = 6.5, fig.height = 3.5, warning=FALSE}
survParamPlot(f = r[1], shape = r[2], points)
```

And we can generate data based on these parameters:

```{r}

defS <- defSurv(varname = "death", formula = -17, scale = 1, shape = 0.3)
defS <- defSurv(defS, varname = "censor", formula = 0, scale = exp(18.5), shape = 0.3)

dd <- genData(500)
dd <- genSurv(dd, defS, timeName = "time", censorName = "censor")
```

```{r, tidy = TRUE, echo = FALSE, fig.width = 6.5, fig.height = 3.5, warning=FALSE}
fit <- survfit( Surv(time, event) ~ 1, data = dd )

survminer::ggsurvplot(fit, data = dd, 
                      ggtheme = ggtheme("grey94"),
                      palette = cbbPalette,
                      legend = "none"
)
```