# Cox proportional hazard models
#
# License: GPL version 3
# Jens Mathis Sauer (c) 2020

source("utils.R")
sma_init()

covariates <- c("dx_age", "asct_age", "time_dx_to_asct", "sex",
		"asct_n", "relapse_status", "nicotin_history","pre_rt", "post_rt",
		"pre_alkylating_agents", "pre_alkylating_agents_n",
		"pre_topoisomerase", "pre_topoisomerase_n",
		"pre_anthracyclines", "pre_anthracyclines_n",
		"peri_alkylating_agents", "peri_alkylating_agents_n",
		"peri_topoisomerase", "peri_topoisomerase_n",
		"peri_anthracyclines", "peri_anthracyclines_n",
		"post_alkylating_agents", "post_alkylating_agents_n",
		"post_topoisomerase", "post_topoisomerase_n",
		"post_anthracyclines", "post_anthracyclines_n")

#
# Extract data from Cox PH results
#
sma_cox_res <- function(m) {
	x <- summary(m)
	p.value<-signif(x$wald["pvalue"], digits=2)
	wald.test<-signif(x$wald["test"], digits=2)
	beta<-signif(x$coef[1], digits=2);#coeficient beta
	HR <-signif(x$coef[2], digits=2);#exp(beta)
	HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
	HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
	HR.confint <- paste0(HR.confint.lower, "-", HR.confint.upper)
	logtest <- signif(x$logtest["test"], digits = 2)
	logtest.pvalue <- x$logtest["pvalue"]

	res<-c(beta, HR, HR.confint, wald.test, p.value, logtest, logtest.pvalue)
	names(res)<-c("beta", "HR", "95% CI for HR", "wald test",
		"p value", "logtest", "logtest p value")
	return(res)
}

#
# Univariate Cox ph analysis with multiple covariates.
#
sma_coxph_univ <- function() {
	univ_formulas <- sapply(covariates, function(x) as.formula(
				paste('Surv(time_at_risk, ifelse(time_at_risk_status == 1, 1, 0)) ~ ',
				x)))
	univ_models <- lapply(univ_formulas,
			      function(x) {coxph(x, data = secmal)})
	univ_results <- lapply(univ_models, sma_cox_res)
	res <- t(as.data.frame(univ_results, check.names = FALSE))

	cat("\nUnivariate coxph with multiple covariates:\n")
	cat("##########################################\n\n")
	print(as.data.frame(res))

	cat("\nUnivariat coxph with multiple covariates (HR over 1):\n")
	cat("#####################################################\n\n")
	# Convert from factor to double then filter
	resgt <- as.data.frame(res)
	resgt$HR <- as.double(as.character(resgt$HR))
	resgt <- filter(resgt, HR > 1.0)
	resgt$HR <- sprintf("%0.2f", resgt$HR)
	print(resgt)

	cat("\nUnivariat coxph with multiple covariates (HR lower or 1):\n")
	cat("#########################################################\n\n")
	# Convert from factor to double then filter
	resle <- as.data.frame(res)
	resle$HR <- as.double(as.character(resle$HR))
	resle <- filter(resle, HR <= 1.0)
	resle$HR <- sprintf("%0.2f", resle$HR)
	print(resle)

	cat("\nCoxph model assumptions (check for non-significance)\n")
	cat("###################################################\n\n")
	coxzph <- lapply(univ_models, cox.zph)
	print(coxzph)
}

#
# Univariate Cox ph analysis with multiple covariates.
#
sma_coxph_univ_hr <- function() {
	univ_formulas <- sapply(covariates, function(x) as.formula(
				paste('Surv(time_at_risk, time_at_risk_status) ~ ',
				x)))
	univ_models <- lapply(univ_formulas,
			      function(x) {coxph(x, data = secmal)})
	univ_results <- lapply(univ_models, sma_cox_res)
	res <- t(as.data.frame(univ_results, check.names = FALSE))

	print("Univariate coxph with multiple covariates:")
	print(as.data.frame(res))
}

#
# Stratified univariate Cox ph analysis with multiple covariates.
# The strata is "diagnosis"
#
sma_coxph_univ_strata <- function() {
	univ_formulas <- sapply(covariates, function(x) as.formula(
				paste('Surv(time_at_risk, time_at_risk_status) ~ ',
				x, ' + strata(diagnosis)')))
	univ_models <- lapply(univ_formulas,
			      function(x) {coxph(x, data = secmal)})
	univ_results <- lapply(univ_models, sma_cox_res)
	res <- t(as.data.frame(univ_results, check.names = FALSE))

	print("Univariate coxph with multiple covariates, stratified by diagnosis:")
	print(as.data.frame(res))
}