#Load Packages
library(ggplot2)
library(survival)
library(survminer)
library(broom)
library(knitr)
library(gtsummary)
library(muhaz)
library(flexsurv)
library(data.table)
library(rpart)
library(randomForestSRC)
library(MASS)
#===============================================================================;


  #Example 1.1: Cox Regression, treatment only
head(stagec)
hist(stagec$pgtime) #years
table(stagec$pgstat)
fit1 <- ggsurvplot(
  fit = survfit(Surv(pgtime, pgstat)~1, data=stagec),
  xlab="Years",
  ylab="Overall non-progression probability",
  palette = "orange"
)
fit1

mod1 <-coxph(Surv(pgtime,pgstat)~ploidy, data=stagec)
mod1 %>%
  gtsummary::tbl_regression(exp=TRUE)
fit2 <- ggsurvplot(
  fit = survfit(Surv(pgtime, pgstat)~ploidy, data=stagec),
  xlab="Years",
  ylab="Overall non-progression probability"
)
fit2
#===============================================================================;


  #Example 1.2: Cox Regression, other covariates
mod2.1 <-coxph(Surv(pgtime, pgstat)~ploidy + age + factor(eet) + 
               g2 + factor(grade) + factor(gleason), data=stagec)
mod2.1 %>%
  gtsummary::tbl_regression(exp=TRUE)

mod2.2 <-coxph(Surv(pgtime, pgstat)~ploidy + age + factor(eet) + 
               g2 + grade + gleason, data=stagec)
mod2.2 %>%
  gtsummary::tbl_regression(exp=TRUE)

mod2.3 <-coxph(Surv(pgtime, pgstat)~ploidy + 
                 g2, data=stagec)
mod2.3 %>%
  gtsummary::tbl_regression(exp=TRUE)
#===============================================================================;


  #Example 1.3: Parametric Models
k_haz_est <-muhaz(stagec$pgtime, stagec$pgstat)
k_haz <-data.table(time=k_haz_est$est.grid,
                    est= k_haz_est$haz.est,
                    method="Kernel density")

dists <- c("exp", "weibull", "gompertz", "gamma", 
           "lognormal", "llogis", "gengamma")
dists_long <- c("Exponential", "Weibull (AFT)",
                "Gompertz", "Gamma", "Lognormal", "Log-logistic",
                "Generalized gamma")
parametric_haz <- vector(mode = "list", length = length(dists))
for (i in 1:length(dists)){
  fit <- flexsurvreg(Surv(pgtime, pgstat) ~ 1, data = stagec, dist = dists[i]) 
  parametric_haz[[i]] <- summary(fit, type = "hazard", 
                                 ci = FALSE, tidy = TRUE)
  parametric_haz[[i]]$method <- dists_long[i]
}
parametric_haz <- rbindlist(parametric_haz)

haz <- rbind(k_haz, parametric_haz)
haz[, method := factor(method,
                       levels = c("Kernel density",
                                  dists_long))]
n_dists <- length(dists) 
ggplot(haz, aes(x = time, y = est, col = method, linetype = method)) +
  geom_line(size=2) +
  xlab("Years") + ylab("Hazard") + 
  scale_colour_manual(name = "", 
                      values = c("black", rainbow(n_dists))) +
  scale_linetype_manual(name = "",
                        values = c(1, rep_len(2:6, n_dists)))

mod3 <-flexsurvreg(Surv(pgtime, pgstat)~ploidy, 
                   data=stagec, dist='gengamma')
mod3 %>%
  gtsummary::tbl_regression(exp=TRUE)

mod4 <-flexsurvreg(Surv(pgtime, pgstat)~ploidy + g2 + factor(grade), 
                   data=stagec, dist='gengamma')
mod4 %>%
  gtsummary::tbl_regression(exp=TRUE)
#===============================================================================;


  #Example 2.1: Survival Tree
stagec$progstat <- factor(stagec$pgstat, levels = 0:1, labels = c("No", "Prog"))
tree1 <- rpart(progstat ~ age + factor(eet) + g2 + factor(grade) + factor(gleason) + ploidy,
                data = stagec, method = 'class')
par(mar=rep(0.1, 4))
plot(tree1)
text(tree1, use.n=TRUE, min=2)
par(mar=c(5.1, 4.1, 4.1, 2.1))
#===============================================================================;


  #Example 2.2: Survival Random Forest
v.obj <-rfsrc(Surv(pgtime, pgstat) ~ age + factor(eet) + g2 + factor(grade) 
              + factor(gleason) + ploidy, data=stagec,
              ntree=100, block.size=1)
plot(get.tree.rfsrc(v.obj, 3))
plot(get.tree.rfsrc(v.obj, 33))
plot(get.tree.rfsrc(v.obj, 99))


print(v.obj)
plot(v.obj)

plot.survival(v.obj, subset=1:10)
plot.survival(v.obj)
#===============================================================================;


  #Example 3: Frailty Model
      #https://www.rdocumentation.org/packages/survival/versions/3.2-13/topics/frailty
      
head(lung)

  #no random effect
mod5.1 <-coxph(Surv(time, status) ~ age, 
               data=lung)
mod5.1 %>%
  gtsummary::tbl_regression(exp=TRUE)

  #random institutional effect
mod5.2 <-coxph(Surv(time, status) ~ age 
               + frailty(inst, df=4), 
               data=lung)
mod5.2 %>%
  gtsummary::tbl_regression(exp=TRUE)