library("survminer")
## Loading required package: ggplot2
## Loading required package: ggpubr
## Loading required package: magrittr
require("survival")
## Loading required package: survival
library(car)
## Loading required package: carData
km<- read.csv("match_all.csv")
fit <- survfit(Surv(sm, death) ~ an_type, data = km)
# KM plot
ggsurvplot(
fit, # survfit object with calculated statistics.
data = km, # data used to fit survival curves.
risk.table = TRUE, # show risk table.
pval = TRUE, pval.coord = c(5, 0.75), # show p-value of log-rank test.
conf.int = TRUE, # show confidence intervals for
# point estimates of survival curves.
xlim = c(0,60),
ylim = c(0.7,1),# present narrower X axis, but not affect
# survival estimates.
xlab = "Time (month)", # customize X axis label.
ylab = "Survival probability",
legend.labs =
c("VIA", "TIVA"),
break.time.by = 12, # break X axis in time intervals by 60.
ggtheme = theme_light(), # customize plot and risk table with a theme.
risk.table.y.text.col = T, # colour risk table text annotations.
risk.table.y.text = FALSE, # show bars instead of names in text annotations
surv.median.line = "hv")
## Warning in .add_surv_median(p, fit, type = surv.median.line, fun = fun, :
## Median survival not reached.
covariates <- c("age","sex","ht","wt","bmi", "asa", "htn" , "dm", "an_time" ,"an_type", "op_time" , "meta" , "tf")
univ_formulas <- sapply(covariates,
function(x) as.formula(paste('Surv(sm,death)~', x)))
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = km)})
# Extract data
univ_results <- lapply(univ_models,
function(x){
x <- summary(x)
p.value<-signif(x$wald["pvalue"], digits=4)
wald.test<-signif(x$wald["test"], digits=4)
beta<-signif(x$coef[1], digits=4);#coeficient beta
HR <-signif(x$coef[2], digits=4);#exp(beta)
HR.confint.lower <- signif(x$conf.int[,"lower .95"],4)
HR.confint.upper <- signif(x$conf.int[,"upper .95"],4)
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)
#return(exp(cbind(coef(x),confint(x))))
})
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)
## beta HR (95% CI for HR) wald.test p.value
## age 0.003897 1.004 (0.9887-1.019) 0.25 0.6167
## sex 0.4713 1.602 (1.114-2.305) 6.45 0.01107
## ht 0.01165 1.012 (0.9904-1.034) 1.15 0.2837
## wt -0.01483 0.9853 (0.9672-1.004) 2.45 0.1172
## bmi -0.07016 0.9322 (0.8833-0.9839) 6.5 0.01081
## asa 0.2215 1.248 (0.8702-1.79) 1.45 0.2285
## htn 0.1827 1.2 (0.7884-1.828) 0.73 0.3944
## dm -0.1015 0.9035 (0.5484-1.489) 0.16 0.6903
## an_time 0.002794 1.003 (1.001-1.005) 9.37 0.002203
## an_type 0.2272 1.255 (0.8822-1.785) 1.59 0.2066
## op_time 0.003236 1.003 (1.001-1.005) 12.03 0.0005237
## meta -2.088 0.1239 (0.08536-0.1798) 120.7 4.538e-28
## tf -0.3789 0.6846 (0.1693-2.769) 0.28 0.5951
# ca만 단변량분석
res.cox <- coxph(Surv(sm, death) ~ ca, data = km)
res.cox
## Call:
## coxph(formula = Surv(sm, death) ~ ca, data = km)
##
## coef exp(coef) se(coef) z p
## caC 1.7242 5.6078 0.4820 3.577 0.000347
## caL 2.3311 10.2891 0.5479 4.254 2.10e-05
## caP 2.0907 8.0909 0.5579 3.747 0.000179
## caS 2.0772 7.9818 0.4630 4.487 7.24e-06
##
## Likelihood ratio test=40.47 on 4 df, p=3.466e-08
## n= 1458, number of events= 125
summary(res.cox)
## Call:
## coxph(formula = Surv(sm, death) ~ ca, data = km)
##
## n= 1458, number of events= 125
##
## coef exp(coef) se(coef) z Pr(>|z|)
## caC 1.7242 5.6078 0.4820 3.577 0.000347 ***
## caL 2.3311 10.2891 0.5479 4.254 2.10e-05 ***
## caP 2.0907 8.0909 0.5579 3.747 0.000179 ***
## caS 2.0772 7.9818 0.4630 4.487 7.24e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## caC 5.608 0.17832 2.180 14.42
## caL 10.289 0.09719 3.515 30.11
## caP 8.091 0.12360 2.711 24.15
## caS 7.982 0.12528 3.221 19.78
##
## Concordance= 0.631 (se = 0.02 )
## Rsquare= 0.027 (max possible= 0.706 )
## Likelihood ratio test= 40.47 on 4 df, p=3e-08
## Wald test = 23.28 on 4 df, p=1e-04
## Score (logrank) test = 31.19 on 4 df, p=3e-06
#다변량분석
res.cox1 <- coxph(Surv(sm, death) ~ age + sex + ht + wt + bmi + asa + htn + dm +an_type+ an_time + meta + tf + ca, data = km)
summary (res.cox1)
## Call:
## coxph(formula = Surv(sm, death) ~ age + sex + ht + wt + bmi +
## asa + htn + dm + an_type + an_time + meta + tf + ca, data = km)
##
## n= 1458, number of events= 125
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age -0.008437 0.991598 0.009563 -0.882 0.377637
## sexM 0.159222 1.172598 0.288591 0.552 0.581138
## ht -0.030469 0.969990 0.033872 -0.900 0.368371
## wt 0.018422 1.018593 0.041261 0.446 0.655247
## bmi -0.086697 0.916955 0.101795 -0.852 0.394392
## asaII 0.303566 1.354681 0.216902 1.400 0.161647
## htnX 0.191567 1.211146 0.252958 0.757 0.448867
## dmX -0.257291 0.773143 0.275894 -0.933 0.351041
## an_typeT 0.277456 1.319768 0.181502 1.529 0.126346
## an_time 0.001160 1.001161 0.001341 0.865 0.387082
## metaX -2.136327 0.118088 0.217330 -9.830 < 2e-16 ***
## tfX 0.403868 1.497607 0.730832 0.553 0.580527
## caC 0.978264 2.659835 0.537981 1.818 0.069003 .
## caL 1.975964 7.213567 0.596595 3.312 0.000926 ***
## caP 1.966594 7.146295 0.588360 3.342 0.000830 ***
## caS 1.936387 6.933653 0.488427 3.965 7.35e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 0.9916 1.0085 0.97319 1.0104
## sexM 1.1726 0.8528 0.66604 2.0644
## ht 0.9700 1.0309 0.90769 1.0366
## wt 1.0186 0.9817 0.93946 1.1044
## bmi 0.9170 1.0906 0.75110 1.1194
## asaII 1.3547 0.7382 0.88554 2.0724
## htnX 1.2111 0.8257 0.73770 1.9885
## dmX 0.7731 1.2934 0.45021 1.3277
## an_typeT 1.3198 0.7577 0.92471 1.8836
## an_time 1.0012 0.9988 0.99853 1.0038
## metaX 0.1181 8.4683 0.07713 0.1808
## tfX 1.4976 0.6677 0.35753 6.2731
## caC 2.6598 0.3760 0.92667 7.6346
## caL 7.2136 0.1386 2.24041 23.2259
## caP 7.1463 0.1399 2.25563 22.6409
## caS 6.9337 0.1442 2.66203 18.0598
##
## Concordance= 0.764 (se = 0.021 )
## Rsquare= 0.092 (max possible= 0.706 )
## Likelihood ratio test= 140.5 on 16 df, p=<2e-16
## Wald test = 151.7 on 16 df, p=<2e-16
## Score (logrank) test = 215.6 on 16 df, p=<2e-16
ggforest(res.cox1, data=km)
## Warning: Removed 7 rows containing missing values (geom_errorbar).
#