From df588619bbadf9b24241883469ac174c9be132c6 Mon Sep 17 00:00:00 2001 From: seanzheng33 <53749449+seanzheng33@users.noreply.github.com> Date: Wed, 6 Jan 2021 13:30:29 +0000 Subject: [PATCH] Add files via upload --- .../clinical_outcome_summary_code_sz.R | 186 ++++++++++++++++++ 1 file changed, 186 insertions(+) create mode 100644 statistical_analysis/clinical_outcome_summary_code_sz.R diff --git a/statistical_analysis/clinical_outcome_summary_code_sz.R b/statistical_analysis/clinical_outcome_summary_code_sz.R new file mode 100644 index 0000000..e83d3fa --- /dev/null +++ b/statistical_analysis/clinical_outcome_summary_code_sz.R @@ -0,0 +1,186 @@ +#script to run survival analysis - by Sean L Zheng + +#load in required packages +library(data.table) +library(tidyverse) +library("survival") +library("survminer") +library("lessR") +library(RColorBrewer) +library(gtsummary) +library(cmprsk) + +####Load in survival dataset#### +sa <- fread('/Path_to_survival_dataset') + +#filter out those that have withdrawn consent +sa <- sa %>% filter(hasWithdrawn==FALSE) + +#if running incident cases then include line below to remove individuals who have had an event preceding recruitment +#if not then bypass +sa <- sa %>% filter(is.na(pre_death_mace)) + +#####fit survival model#### +fit <- survfit(Surv(time_death_mace, event_death_mace) ~ geno, data=sa) + +#produce CH survival curve +ggsurv<- ggsurvplot(fit, + ylim=c(0,1), + xlab = "Years", + ylab="Probability freedom from events", + legend='bottom', + legend.title = 'Genotype', font.legend=c(14), + legend.labs = c('SARC-NEG', 'SARC-IND', 'SARC-P/LP'), + palette=c("#BDBDBD", "#FFA500",'#A52A2A'), + risk.table = TRUE, risk.table.title = NULL, risk.table.height=0.20, + linetype='strata', + fontsize=4, font.x=14, font.y=14, font.tickslab=12, + break.time.by = 10, + break.y.by=0.1, + fun='cumhaz', + title=NULL) +ggsurv$plot <- ggsurv$plot + theme(legend.text=element_text(size=14)) +ggsurv + +#####Cox PH - compare between pairs of arms #### +res.cox <- coxph(Surv(sa$time_death_mace, event_death_mace) ~ + gene_positive + gene_vus + sexgenetic, + data=sa) +print(res.cox) +coxph(Surv(sa$time_death_mace, event_death_mace) ~ + gene_positive + gene_vus + sexgenetic, + data=sa) %>% gtsummary::tbl_regression(exp=TRUE) + +#if plotting - then separate into 1 vs 1 comparison and then combine later +sa1 <- sa %>% filter(geno!=3) +geno_df <- with(sa1, data.frame(geno=c(1,2), sexgenetic = c(1,1), agerecruit=rep(mean(sa1$agerecruit, na.rm=TRUE), 2))) +res.cox <- coxph(Surv(sa1$time_death_mace, sa1$event_death_mace) ~ geno + agerecruit + sexgenetic, data=sa1) +fit <- survfit(res.cox, data=sa1, newdata=geno_df) +plots <- list() +plots[[1]] <- ggsurvplot(fit, conf.int = TRUE, + legend.labs=c("Sarc-ve", "Sarc Other"), + ylim=c(0,1), + xlim=c(0,80), + xlab = "Years", ylab="Freedom from outcome probability", + legend='bottom', + legend.title = 'Genotype', + palette=c("#BDBDBD", "#FFA500"), + risk.table = FALSE, tables.height=0.2, + linetype = 'strata', + ggtheme=theme_bw(), + fontsize=3, + censor=FALSE, + break.time.by = 20, + break.y.by=0.2, + title=NULL, + size=1) + +sa2 <- sa %>% filter(geno!=2) +geno_df <- with(sa2, data.frame(geno=c(1,3), sexgenetic = c(1,1), agerecruit=rep(mean(sa2$agerecruit, na.rm=TRUE), 2))) +res.cox <- coxph(Surv(sa2$time_death_mace, sa2$event_death_mace) ~ geno + agerecruit + sexgenetic, data=sa2) +fit1 <- survfit(res.cox, data=sa2, newdata=geno_df) +plots[[2]] <- ggsurvplot(fit1, conf.int = TRUE, + legend.labs=c("Sarc-ve", "Sarc+ve"), + ylim=c(0,1), + xlim=c(0,80), + xlab = "Years", ylab="Freedom from outcome probability", + legend='bottom', + legend.title = 'Genotype', + palette=c("#BDBDBD", "#A52A2A"), + risk.table = FALSE, tables.height=0.2, + linetype = 'strata', + ggtheme=theme_bw(), + fontsize=3, + censor=FALSE, + break.time.by = 20, + break.y.by=0.2, + title=NULL, + size=1) + +plots1 <- arrange_ggsurvplots(plots, print=TRUE, nrow=1, ncol=2) + +####Competing risk analysis#### +#labels +sa$cr_death_mace <- factor(sa$cr_death_mace, levels = c(0,1,2), labels = c('Censored', 'MACE', 'Death')) + +sa1 <- sa %>% filter(sexgenetic==0) + +#Event free survival - regression +res <- coxph(formula = Surv(sa1$time_death_mace, cr_death_mace =='MACE') ~ gene_positive + gene_vus, data=sa1, ties = c("efron","breslow","exact")[1]) +print(res) +coxph(Surv(sa1$time_death_mace, cr_death_mace =='MACE') ~ gene_positive + gene_vus, data=sa1, ties = c("efron","breslow","exact")[1]) %>% gtsummary::tbl_regression(exp=TRUE) + +#####Test for assumption of proportional hazards##### +#Using the Schoenfeld residuals +fit <- coxph(Surv(sa$time_death_mace, sa$event_death_mace) ~ geno + sexgenetic, data=sa) +test <- cox.zph(fit) +print(test) +plot(test) + +#####Sex stratified analysis#### +sa$sexgenetic1 <- ifelse(sa$sexgenetic==0, 2, 1) +res.cox <- coxph(Surv(sa$time_death_mace, event_death_mace) ~ + (gene_positive * sexgenetic1) + (gene_vus * sexgenetic1), + data=sa) +print(res.cox) +coxph(Surv(sa$time_death_mace, event_death_mace) ~ + (gene_positive * sexgenetic1) + (gene_vus * sexgenetic1), + data=sa) %>% gtsummary::tbl_regression(exp=TRUE) + +fit2<- survfit(Surv(sa$time_death_mace, event_death_mace) ~ geno + sexgenetic, data=sa) +ggsurv <- ggsurvplot(fit2, + ylim=c(0,0.4), + ggtheme=theme_bw(), + xlab = "Years from recruitment", + ylab="Probability freedom from events", + fun='cumhaz', + legend='bottom', + legend.title = 'Genotype', font.legend=c(14), + legend.labs = c('SARC-NEG (F)', 'SARC-NEG (M)','SARC-IND (F)', 'SARC-IND (M)','SARC-P/LP (F)','SARC-P/LP (M)'), + palette=c("#BDBDBD", "#BDBDBE","#FFA500","#FFA501",'#A52A2A','#A52A2B'), + #palette=c("#BDBDBD", "#FFA500",'#A52A2A'), + risk.table = TRUE, risk.table.title = NULL, + risk.table.height = 0.3, + linetype=c('solid', 'dotted', 'solid', 'dotted','solid', 'dotted'), + censor=FALSE, + fontsize=4, font.x=14, font.y=14, font.tickslab=14) +ggsurv$plot <- ggsurv$plot + theme(legend.text=element_text(size=14)) +ggsurv + +####Wall thickness#### +#note if running FD then use same dataset but change covariate to fd + +wt <- fread('/Path_to_survival_dataset_in_participants_with_cmr') +wt <- wt %>% filter(hasWithdrawn==FALSE) +fit <- survfit(Surv(time_death_mace, event_death_mace) ~ wt, data=wt) + +ggsurv<- ggsurvplot(fit, + ylim=c(0,1), + xlim=c(0,80), + xlab = "Years", + ylab="Probability freedom from events", + legend='bottom', + legend.title = 'LV max wall thickness (mm)', font.legend=c(14), + legend.labs = c('<13', '13 to 15', '>15'), + palette='Blues', + risk.table = TRUE, risk.table.title = NULL, risk.table.height=0.20, + linetype='strata', + fontsize=4, font.x=14, font.y=14, font.tickslab=12, + break.time.by = 10, + break.y.by=0.1, + fun='cumhaz', + title=NULL) +ggsurv$plot <- ggsurv$plot + theme(legend.text=element_text(size=14)) +ggsurv + +#Cox PH - change between using categorical and continuous variable for LV wall thickness +res.cox <- coxph(Surv(time_death_mace, event_death_mace) ~ + medium + thick + # WT_Max..mm. + + autoSysBP + sexgenetic, data=wt) +print(res.cox) +coxph(Surv(time_death_mace, event_death_mace) ~ + medium + thick + # WT_Max..mm. + + autoSysBP + sexgenetic, data=wt) %>% gtsummary::tbl_regression(exp=TRUE) +