# Secondary malignancies calculations
#
# License: GPL version 3
# Jens Mathis Sauer (c) 2020

source("utils.R")
sma_init()

# Setup survival object
surv_sec <- Surv(time = secmal$time_at_risk, event = ifelse(secmal$time_at_risk_status == 1, 1, 0))

# plot cummulative events
sma_plot_secmal_event <- function() {
	ggsurvplot(survfit(surv_sec ~ 1), data = secmal, xscale = "d_y",
		   title = "Secondary malignacies",
		   fun = "event",
		   break.time.by = sma_break.time.by,
		   surv.median.line = "hv",
		   risk.table = "nrisk_cumevents",
		   ggtheme = theme_bw())
}

# plot cummulative hazard
sma_plot_secmal_haz <- function() {
	ggsurvplot(survfit(surv_sec ~ 1), data = secmal, xscale = "d_y",
		   title = "Secondary malignacies",
		   fun = "cumhaz",
		   break.time.by = sma_break.time.by,
		   risk.table = "abs_pct",
		   ggtheme = theme_bw())
}

# plot cummulative events per sex
sma_plot_secmal_event_sex <- function() {
	ggsurvplot(survfit(surv_sec ~ sex, data = secmal), data = secmal, xscale = "d_y",
		   title = "Secondary malignacies",
		   fun = "event",
		   break.time.by = sma_break.time.by,
		   surv.median.line = "hv",
		   risk.table = "nrisk_cumevents",
		   pval = TRUE,
		   ggtheme = theme_bw())
}

# plot cummulative hazard per sex
sma_plot_secmal_haz_sex <- function() {
	ggsurvplot(survfit(surv_sec ~ sex, data = secmal), data = secmal, xscale = "d_y",
		   title = "Secondary malignacies",
		   fun = "cumhaz",
		   break.time.by = sma_break.time.by,
		   risk.table = "abs_pct",
		   pval = TRUE,
		   ggtheme = theme_bw())
}

# plot cummulative events per diagnosis
sma_plot_secmal_event_dx <- function() {
	ggsurvplot(survfit(surv_sec ~ diagnosis, data = secmal), data = secmal, xscale = "d_y",
		   title = "Secondary malignacies",
		   fun = "event",
		   break.time.by = sma_break.time.by,
		   risk.table = "nrisk_cumevents",
		   pval = TRUE,
		   ggtheme = theme_bw())
}

# plot cummulative hazard per diagnosis
sma_plot_secmal_haz_dx <- function() {
	ggsurvplot(survfit(surv_sec ~ diagnosis, data = secmal), data = secmal, xscale = "d_y",
		   title = "Secondary malignacies",
		   fun = "cumhaz",
		   break.time.by = sma_break.time.by,
		   risk.table = "abs_pct",
		   pval = TRUE,
		   ggtheme = theme_bw())
}

# plot competing risk of SM
sma_plot_secmal_cmprsk <- function() {
	cpr <- cuminc(secmal$time_at_risk / 365.25, secmal$time_at_risk_status)
	ggcompetingrisks(cpr[1], conf.int = TRUE, multiple_panels = FALSE,
			legend = "none", title = "", ylim = c(0,0.5),
			xlab = "Years", ggtheme = theme_bw(), palette = "lancet")
}

# Calculate cumuative risk with ci from cuminc.
# CI calculation shamelessly stolen from ggcompetingrisks()
sma_secmal_cmprsk_time <- function(c, t) {
	coef <- 1.96
	tp <- timepoints(c[1], t)
	df <- as.data.frame(matrix(nrow = 1, ncol = 6))
	colnames(df) <- c("est", "var", "time", "std", "ci.lower", "ci.upper")
	df$var <- as.numeric(tp$var)
	df$est <- as.numeric(tp$est)
	df$time <- t
	df$std <- sqrt(df$var)
	df$ci.lower <- df$est - coef * df$std
	df$ci.upper <- df$est + coef * df$std

	return(df)
}

# Write secondary malignancy plots to files
sma_plot_file_secmal <- function() {
	sma_plot_file("secmal_event.png", png, sma_plot_secmal_event)
	sma_plot_file("secmal_haz.png", png, sma_plot_secmal_haz)
	sma_plot_file("secmal_event_sex.png", png, sma_plot_secmal_event_sex)
	sma_plot_file("secmal_haz_sex.png", png, sma_plot_secmal_haz_sex)
	sma_plot_file("secmal_event_dx.png", png, sma_plot_secmal_event_dx)
	sma_plot_file("secmal_haz_dx.png", png, sma_plot_secmal_haz_dx)
}