From 8f07c45f6367b90683aa84648690a2e4959298d7 Mon Sep 17 00:00:00 2001
From: Jens Sauer <jsauer65@gmail.com>
Date: Wed, 25 Nov 2020 17:26:18 +0100
Subject: [PATCH] coxph: First commit of coxph models

First univariate with covariates coxph models.
---
 coxph.R | 72 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 72 insertions(+)
 create mode 100644 coxph.R

diff --git a/coxph.R b/coxph.R
new file mode 100644
index 0000000..189c691
--- /dev/null
+++ b/coxph.R
@@ -0,0 +1,72 @@
+# 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(x) {
+	x <- summary(x)
+	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 <- paste0(HR, " (",
+	       HR.confint.lower, "-", HR.confint.upper, ")")
+
+	res<-c(beta, HR, wald.test, p.value)
+	names(res)<-c("beta", "HR (95% CI for HR)", "wald.test",
+		"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, 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))
+}