## ----surv-setup, message=FALSE, warning=FALSE---------------------------------
library(gtregression)
library(dplyr)

data("data_lungcancer", package = "gtregression")

lung_data <- data_lungcancer |>
  mutate(
    trt = factor(trt, levels = c(1, 2),
                 labels = c("Standard treatment", "Test treatment")),
    prior = factor(prior, levels = c(0, 10), labels = c("No", "Yes")),
    celltype = factor(
      celltype,
      levels = c("squamous", "smallcell", "adeno", "large"),
      labels = c("Squamous", "Small cell", "Adenocarcinoma", "Large cell")
    )
  )

attr(lung_data$time, "label") <- "Survival time"
attr(lung_data$status, "label") <- "Death status"
attr(lung_data$trt, "label") <- "Treatment group"
attr(lung_data$celltype, "label") <- "Cancer cell type"
attr(lung_data$karno, "label") <- "Karnofsky performance score"
attr(lung_data$age, "label") <- "Age"
attr(lung_data$prior, "label") <- "Prior therapy"

surv_exposures <- c("trt", "celltype", "karno", "age", "prior")

## ----surv-describe, message=FALSE, warning=FALSE------------------------------
lung_summary <- descriptive_table(
  data = lung_data,
  exposures = c("time", "status", "celltype", "karno", "age", "prior"),
  by = trt,
  statistic = c(time = "median", karno = "mean", age = "mean"),
  percent = column,
  show_overall = last
)

lung_summary$table

## ----surv-km-plot, message=FALSE, warning=FALSE-------------------------------
km_curve <- km_plot(
  data = lung_data,
  time = time,
  event = status,
  by = trt,
  break_time_by = 200,
  ylim = c(50, 100),
  title = "Kaplan-Meier Survival by Treatment"
)

km_curve

## ----surv-km-panel, fig.width=7, fig.height=5, message=FALSE, warning=FALSE----
km_trt_panel <- km_plot(
  data = lung_data,
  time = time,
  event = status,
  by = trt,
  risk_table = FALSE,
  break_time_by = 200,
  ylim = c(50, 100),
  title = "A. Treatment group",
  title_size = 10,
  title_face = plain,
  legend_position = bottom,
  base_size = 10
)

km_prior_panel <- km_plot(
  data = lung_data,
  time = time,
  event = status,
  by = prior,
  risk_table = FALSE,
  break_time_by = 200,
  ylim = c(50, 100),
  title = "B. Prior therapy",
  title_size = 10,
  title_face = plain,
  legend_position = bottom,
  base_size = 10
)

patchwork::wrap_plots(km_trt_panel, km_prior_panel, ncol = 2) +
  patchwork::plot_layout(guides = "collect") &
  ggplot2::theme(legend.position = "bottom")

## ----surv-observed-tables, message=FALSE, warning=FALSE-----------------------
survival_summary(
  data = lung_data,
  time = time,
  event = status,
  by = trt
)$table

survival_prob(
  data = lung_data,
  time = time,
  event = status,
  by = trt,
  times = c(90, 180, 365)
)$table

## ----surv-rmst, message=FALSE, warning=FALSE----------------------------------
rmst_table(
  data = lung_data,
  time = time,
  event = status,
  by = trt,
  tau = 365
)$table

## ----surv-logrank, message=FALSE, warning=FALSE-------------------------------
logrank_test(
  data = lung_data,
  time = time,
  event = status,
  by = trt
)$table

## ----surv-cox, message=FALSE, warning=FALSE-----------------------------------
cox_crude <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = surv_exposures
)

cox_adjusted <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = c(trt, celltype, prior),
  adjust_for = c(age, karno)
)

cox_adjusted$table

## ----surv-cox-interaction, message=FALSE, warning=FALSE-----------------------
cox_interaction <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = trt,
  adjust_for = c(age, karno),
  interaction = trt*prior
)

cox_interaction$table

## ----surv-check-ph, message=FALSE, warning=FALSE------------------------------
check_ph(cox_adjusted)$table

## ----surv-parametric-checks, message=FALSE, warning=FALSE---------------------
surv_model_compare(
  data = lung_data,
  time = time,
  event = status,
  exposures = c(trt, celltype, prior),
  adjust_for = c(age, karno),
  distributions = c(weibull, exponential, "log-normal", "log-logistic")
)$table

plot_surv_fit(
  data = lung_data,
  time = time,
  event = status,
  by = trt,
  adjust_for = c(age, karno),
  distributions = c(weibull, "log-logistic"),
  break_time_by = 200
)

## ----surv-parametric, message=FALSE, warning=FALSE----------------------------
surv_crude <- surv_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = surv_exposures,
  distribution = loglogistic
)

surv_adjusted <- surv_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = c(trt, celltype, prior),
  adjust_for = c(age, karno),
  distribution = loglogistic,
  model_stats = TRUE
)

surv_adjusted$table
surv_adjusted$model_stats

## ----surv-parametric-interaction, message=FALSE, warning=FALSE----------------
surv_interaction <- surv_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = trt,
  adjust_for = c(age, karno),
  interaction = trt*prior,
  distribution = loglogistic
)

surv_interaction$table

## ----surv-predict, message=FALSE, warning=FALSE-------------------------------
surv_predict(
  model = surv_adjusted$models$trt,
  newdata = data.frame(
    trt = factor("Test treatment", levels = levels(lung_data$trt)),
    age = 60,
    karno = 70
  ),
  times = c(90, 180, 365)
)$table

## ----surv-visualise-export, message=FALSE, warning=FALSE----------------------
plot_reg_combine(
  cox_crude,
  cox_adjusted,
  show_ref = FALSE,
  title_uni = "Crude HR",
  title_multi = "Adjusted HR"
)

surv_forest_data <- forest_df(cox_crude, cox_adjusted, desc = lung_summary)

forest_reg(
  surv_forest_data,
  xlim = list(c(0.25, 8), c(0.25, 8)),
  ticks_at = list(c(0.5, 1, 2, 4, 8), c(0.5, 1, 2, 4, 8)),
  quiet = TRUE
)