<!--START_EXPLANATION-->
Rmarkdown template to generate this page can be found on [Rmarkdown-Template](Rmarkdown/Multiple_Ascending_Dose_PKPD_time_to_event.Rmd).
<!--END_EXPLANATION-->

```{r, error = TRUE, echo = TRUE, warning=FALSE, message=FALSE}
library(ggplot2)
library(dplyr)
library(broom)
library(knitr)
library(tidyr)
library(survival)
library(survminer)
library(xgxr)
set.seed(123456)
```

# Explore Exposure-Response Relationship

## Create a Dataset

```{r, error = TRUE, warning=FALSE}
#Use the lung dataset to create a fake exposure dataset
km_data <- lung %>%
  mutate(fake_Cmax = ph.ecog + age/50 + sex + runif(nrow(lung)),
         fake_AUC  = ph.ecog*3 + (age/50)^sex + runif(nrow(lung))) %>%
  filter(!is.na(fake_Cmax)) %>%
  mutate(Cmax_quartile= cut(fake_Cmax,
                            breaks=quantile(fake_Cmax,c(0,.25,.5,.75,1)),
                            include.lowest=TRUE,
                            labels=paste0("Q",c(4,3,2,1))),
         AUC_quartile = cut(fake_AUC,
                            breaks=quantile(fake_AUC,c(0,.25,.5,.75,1)),
                            include.lowest=TRUE,
                            labels=paste0("Q",c(4,3,2,1))),)
#these columns are required for your dataset
km_data <- km_data %>%
  mutate(exposure = Cmax_quartile, #exposure quantile,
         time   = time,   #time of the event (or censoring)
         event  = status) #status: there are a three options for this column (see ?Surv)
#        a) 0 = censored (alive), 1 = dead (event)
#        b) 1 = censored (alive), 2 = dead (event)
#        c) 0 = right censored, 1 = event at time, 
#           2 = left censored,  3 = interval censored.
km_data_Cmax_AUC <- km_data %>%
  gather(NCAparam, NCAvalue, Cmax_quartile:AUC_quartile) 
```

## Plot Kaplan Meier with Confidence Intervals using only Cmax
<!--START_EXPLANATION-->
Time-to-event plots can be summarized by Kaplan-Meier plots and stratified by exposure quartile to give an overview of the dose-response.  To see if there's an effect on exposure vs response, look to see if there is separation between the quartiles.  However note that it may be that there are factors that impact both exposure and response and so a more careful model-based analysis may be needed to formally assess any causal relationship, as described in the reference below.

Yang, Jun, et al. "The combination of exposure‐response and case‐control analyses in regulatory decision making." The Journal of Clinical Pharmacology 53.2 (2013): 160-166.
<!--END_EXPLANATION-->

```{r, error = TRUE, warning=FALSE, message=FALSE}
km_fit <- survfit(Surv(time, status) ~ exposure, data = km_data, conf.int = 0.95)
gg <- ggsurvplot(km_fit, km_data, conf.int = TRUE, ggtheme = xgx_theme())
gg <- gg + xgx_scale_x_time_units(units_dataset = "day", units_plot = "year")
print(gg)
```

## Plot Kaplan Meier faceted by exposure metric
<!--START_EXPLANATION-->
Using the `ggsurvplot_facet` command, one can look at many exposure metrics at once.  Again, look for separation in the quartiles, but be aware that something like a case-control analysis may be needed to assess whether exposure causes reduced response or whether other factors may affect both exposure and response.
<!--END_EXPLANATION-->

```{r, error = TRUE, warning=FALSE, message=FALSE}
km_fit_Cmax_AUC <- survfit(Surv(time, status) ~ NCAvalue, data = km_data_Cmax_AUC, conf.int = 0.95)
gg <- ggsurvplot_facet(km_fit_Cmax_AUC, 
                       km_data_Cmax_AUC, 
                       facet.by = "NCAparam",
                       conf.int = TRUE, 
                       ggtheme = xgx_theme())
gg <- gg + xgx_scale_x_time_units(units_dataset = "day", units_plot = "year")
print(gg)
```

## Plot Kaplan Meier faceted by exposure quartile for assessing covariates
<!--START_EXPLANATION-->
This plot can help in assessing the impact of a covariate (ECOG) on outcome.

Note, `ph.ecog == 3` patient are removed from dataset because there is only one patient
with ECOG = 3 and this causes an error when generating the confidence interval.  This bug was [reported to the survminer developer and fixed in the github version](https://github.com/kassambara/survminer/issues/443).
<!--END_EXPLANATION-->

```{r, error = TRUE, warning=FALSE, message=FALSE}
km_data_noecog3 <- km_data %>% filter (ph.ecog != 3)
km_fit_ecog <- survfit(Surv(time, status) ~ ph.ecog, 
                       data = km_data_noecog3 , 
                       conf.int = 0.95)
gg <- ggsurvplot_facet(km_fit_ecog, 
                       km_data_noecog3, 
                       facet.by = c("Cmax_quartile"),
                       conf.int = TRUE,
                       ggtheme = xgx_theme(),
                       legend.labs = paste0("ECOG=",sort(unique(km_data_noecog3$ph.ecog))))
gg <- gg + xgx_scale_x_time_units(units_dataset = "day", units_plot = "year")
print(gg)
```

## Performing Cox regression on the above data. 
<!--START_EXPLANATION-->
In this data only 2 variables are available, so variable selection isn't needed. If there are multiple variables, applying a variable selection approach (e.g. stepwiseAIC) would be advisable. Look for small p values or large hazard ratios to identify covariates that have a large effect.
<!--END_EXPLANATION-->

```{r, error = TRUE, warning=FALSE, message=FALSE}
km_cox <- coxph(Surv(time) ~ fake_AUC + fake_Cmax,
             data =km_data)
km_cox_summary <- broom::tidy(km_cox) %>%
  mutate_at(vars(-term),signif,2)
kable(km_cox_summary)  
```

## R Session Info
```{r}
sessionInfo()
```
