── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.2 ✔ tibble 3.3.0
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.4
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(survival)
library(survminer)
Loading required package: ggpubr
Attaching package: 'survminer'
The following object is masked from 'package:survival':
myeloma
km_all <- survfit(Surv(time, status) ~ 1, data = OVC_Data)
s <- summary(km_all)
life_table <- data.frame(
time = s$time,
n_risk = s$n.risk,
n_event = s$n.event,
n_censor = s$n.censor,
survival = s$surv,
std_err = s$std.err,
lower = s$lower,
upper = s$upper
) %>%
mutate(cum_hazard = -log(survival))
print(head(life_table, 20))
time n_risk n_event n_censor survival std_err lower upper
1 59 26 1 0 0.9615385 0.03771464 0.8903890 1.0000000
2 115 25 1 0 0.9230769 0.05225894 0.8261294 1.0000000
3 156 24 1 0 0.8846154 0.06265627 0.7699542 1.0000000
4 268 23 1 0 0.8461538 0.07075894 0.7182378 0.9968513
5 329 22 1 0 0.8076923 0.07729201 0.6695613 0.9743199
6 353 21 1 0 0.7692308 0.08262864 0.6231935 0.9494899
7 365 20 1 0 0.7307692 0.08698929 0.5787019 0.9227957
8 431 17 1 2 0.6877828 0.09188148 0.5293447 0.8936430
9 464 15 1 1 0.6419306 0.09652130 0.4780800 0.8619371
10 475 14 1 0 0.5960784 0.09992615 0.4291495 0.8279387
11 563 12 1 1 0.5464052 0.10320939 0.3773402 0.7912188
12 638 11 1 0 0.4967320 0.10510266 0.3281088 0.7520148
cum_hazard
1 0.03922071
2 0.08004271
3 0.12260232
4 0.16705408
5 0.21357410
6 0.26236426
7 0.31365756
8 0.37428218
9 0.44327505
10 0.51738302
11 0.60439440
12 0.69970458
ggsurvplot(km_all,
data = OVC_Data,
risk.table = TRUE,
surv.median.line = "hv",
xlab = "Time",
ylab = "Survival probability",
title = "Overall Kaplan-Meier Survival Curve")
km_by_trt <- survfit(Surv(time, status) ~ treatment, data = OVC_Data)
logrank <- survdiff(Surv(time, status) ~ treatment, data = OVC_Data)
logrank
Call:
survdiff(formula = Surv(time, status) ~ treatment, data = OVC_Data)
N Observed Expected (O-E)^2/E (O-E)^2/V
treatment=1 13 7 5.23 0.596 1.06
treatment=2 13 5 6.77 0.461 1.06
Chisq= 1.1 on 1 degrees of freedom, p= 0.3
p_val <- 1 - pchisq(logrank$chisq, df = length(logrank$n) - 1)
cat("Log-rank test p-value:", signif(p_val, 4), "\n")
Log-rank test p-value: 0.3026
ggsurvplot(km_by_trt,
data = OVC_Data,
pval = TRUE,
pval.method = TRUE, # shows log-rank method _label_
risk.table = TRUE,
conf.int = FALSE,
legend.title = "Group",
legend.labs = levels(OVC_Data$treatment),
xlab = "Time",
ylab = "Survival probability",
title = "Kaplan-Meier by Treatment")
cox_mod <- coxph(Surv(time, status) ~ age + treatment, data = OVC_Data)
summary(cox_mod)
Call:
coxph(formula = Surv(time, status) ~ age + treatment, data = OVC_Data)
n= 26, number of events= 12
coef exp(coef) se(coef) z Pr(>|z|)
age 0.14657 1.15786 0.04585 3.196 0.00139 **
treatment -0.79593 0.45116 0.63294 -1.258 0.20857
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.1579 0.8637 1.0583 1.267
treatment 0.4512 2.2165 0.1305 1.560
Concordance= 0.794 (se = 0.079 )
Likelihood ratio test= 15.82 on 2 df, p=4e-04
Wald test = 13.55 on 2 df, p=0.001
Score (logrank) test = 18.61 on 2 df, p=9e-05
tidy_cox <- broom::tidy(cox_mod, exponentiate = TRUE, conf.int = TRUE)
print(tidy_cox)
# A tibble: 2 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 age 1.16 0.0459 3.20 0.00139 1.06 1.27
2 treatment 0.451 0.633 -1.26 0.209 0.130 1.56
zph <- cox.zph(cox_mod)
print(zph)
chisq df p
age 0.276 1 0.60
treatment 0.776 1 0.38
GLOBAL 0.959 2 0.62
ggadjustedcurves(cox_mod, data = OVC_Data, variable = "treatment",
legend.title = "Group",
title = "Cox-adjusted survival curves by Treatment",
xlab = "Time", ylab = "Adjusted survival")
