Jens Sauer
63c2883c7e
Since the data specifies death as a competing risk, the cox analysis must be adjusted. Death will now, as it was before, treated as a censor.
114 lines
3.9 KiB
R
114 lines
3.9 KiB
R
# 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))
|
|
}
|