1 生存曲线
1.1 Kaplan-Meier curve
Code
### 数据 三列 “group”分组列， “time”生存时间列，“status”事件状态列
d<-read_excel("data/01source.xlsx",sheet=1,range = "A19:E75") |>
rename(
LR_DPVB=`LR-DPVB`
) |>
pivot_longer(
cols = CON:LR_DPVB,
names_to = "treatment",
values_to = "status"
) |> drop_na() |>
mutate(
treatment=factor(treatment,levels=c('CON','LDRT','DPVB','LR_DPVB'))
)
# 加载包
d
#> # A tibble: 56 × 3
#> Days treatment status
#> <dbl> <fct> <dbl>
#> 1 11 CON 1
#> 2 13 CON 1
#> 3 14 CON 1
#> 4 14 CON 1
#> 5 15 CON 1
#> 6 17 CON 1
#> 7 17 CON 1
#> 8 17 CON 1
#> 9 20 CON 1
#> 10 20 CON 1
#> # ℹ 46 more rows
# 使用 Kaplan-Meier 方法拟合生存曲线
km_fit<-survfit2(Surv(Days,status)~treatment,data=d)
summary(km_fit)
#> Call: survfit(formula = Surv(Days, status) ~ treatment, data = d)
#>
#> treatment=CON
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 11 14 1 0.9286 0.0688 0.8030 1.000
#> 13 13 1 0.8571 0.0935 0.6921 1.000
#> 14 12 2 0.7143 0.1207 0.5129 0.995
#> 15 10 1 0.6429 0.1281 0.4351 0.950
#> 17 9 3 0.4286 0.1323 0.2341 0.785
#> 20 6 2 0.2857 0.1207 0.1248 0.654
#> 21 4 2 0.1429 0.0935 0.0396 0.515
#> 25 2 1 0.0714 0.0688 0.0108 0.472
#> 27 1 1 0.0000 NaN NA NA
#>
#> treatment=LDRT
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 13 14 2 0.8571 0.0935 0.6921 1.000
#> 15 12 1 0.7857 0.1097 0.5977 1.000
#> 16 11 1 0.7143 0.1207 0.5129 0.995
#> 18 10 1 0.6429 0.1281 0.4351 0.950
#> 19 9 2 0.5000 0.1336 0.2961 0.844
#> 20 7 3 0.2857 0.1207 0.1248 0.654
#> 24 4 1 0.2143 0.1097 0.0786 0.584
#> 25 3 1 0.1429 0.0935 0.0396 0.515
#> 27 2 1 0.0714 0.0688 0.0108 0.472
#> 30 1 1 0.0000 NaN NA NA
#>
#> treatment=DPVB
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 20 14 1 0.929 0.0688 0.803 1.000
#> 23 13 1 0.857 0.0935 0.692 1.000
#> 27 12 1 0.786 0.1097 0.598 1.000
#> 28 11 1 0.714 0.1207 0.513 0.995
#> 30 10 1 0.643 0.1281 0.435 0.950
#> 32 9 1 0.571 0.1323 0.363 0.899
#> 38 8 1 0.500 0.1336 0.296 0.844
#> 39 7 1 0.429 0.1323 0.234 0.785
#> 45 6 1 0.357 0.1281 0.177 0.721
#>
#> treatment=LR_DPVB
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 30 14 1 0.929 0.0688 0.803 1
#> 40 13 1 0.857 0.0935 0.692 1
1.2 plot()
1.3 survminer::ggsurvplot()
Code
# 使用 ggplot2 和 survminer 绘制生存曲线图
library(survminer)
p <- ggsurvplot(
fit = km_fit,
data = d,
fun = "pct", #event',累计死亡比列，'cumhaz':cumulative hazard累计风险
palette = "jco",
linetype = c(1,5,4,6),
surv.median.line = 'none',
conf.int = F,
pval = TRUE,# "log rank p=0.0053"
log.rank.weights='1',
pval.method = TRUE,
pval.coord=c(40,75),
pval.method.coord=c(32,75),
risk.table = "abs_pct", # "nrisk_cumcensor" and "nrisk_cumevents"
risk.table.col = "treatment",
risk.table.height=0.3,
tables.theme = theme_cleantable(),
ncensor.plot=T,
axes.offset=F, #(0,0)pos
title = "Hepa1-6 tumor",
xlab="Days",
ylab="Percent survival (%)",
legend.title = "",
legend = c(0.2,0.35),
legend.labs=c('CON','LDRT','DPVB','LR_DPVB'),
ggtheme = theme_survminer(),
)
p
Code
fit<- survfit(Surv(time, status) ~ sex, data = lung)
# Basic survival curves
ggsurvplot(fit, data = lung)
Code
# Customized survival curves
p <- ggsurvplot(fit, data = lung,
surv.median.line = "hv", # Add medians survival
# Change legends: title & labels
legend.title = "Sex",
legend.labs = c("Male", "Female"),
# Add p-value and tervals
pval = TRUE,
conf.int = TRUE,
# Add risk table
risk.table = TRUE,
tables.height = 0.2,
tables.theme = theme_cleantable(),
# Color palettes. Use custom color: c("#E7B800", "#2E9FDF"),
# or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco")
palette = c("#E7B800", "#2E9FDF"),
ggtheme = theme_bw() # Change ggplot2 theme
)
#> Warning in geom_segment(aes(x = 0, y = max(y2), xend = max(x1), yend = max(y2)), : All aesthetics have length 1, but the data has 2 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#> a single row.Code
p <- ggsurvplot(fit, data = lung,
surv.median.line = "hv", # Add medians survival
# Change legends: title & labels
legend.title = "Sex",
legend.labs = c("Male", "Female"),
# Add p-value and tervals
pval = F, ##不默认显示P
conf.int = TRUE,
# Add risk table
risk.table = TRUE,
tables.height = 0.2,
tables.theme = theme_cleantable(),
# Color palettes. Use custom color: c("#E7B800", "#2E9FDF"),
# or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco")
palette = c("#E7B800", "#2E9FDF"),
ggtheme = theme_bw() # Change ggplot2 theme
)
#> Warning in geom_segment(aes(x = 0, y = max(y2), xend = max(x1), yend = max(y2)), : All aesthetics have length 1, but the data has 2 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#> a single row.
##计算公式
data.survdiff <- survdiff(Surv(time, status) ~ sex,data = lung)
p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
ci <- paste0(sprintf("%.3f",HR)," [",sprintf("%.3f",low95),", ",sprintf("%.3f",up95),"]")
p$plot <- p$plot +
annotate("text", x = 0, y = 0.05,
label = paste0("P value = ",sprintf("%.3f",p.val),"\n HR (95% CI) = ",ci), ###添加P和HR 95%CI
size = 5, color = "black", hjust = 0)+
theme(text = element_text(size = 15))
p$plot
#> Warning in geom_segment(aes(x = 0, y = max(y2), xend = max(x1), yend = max(y2)), : All aesthetics have length 1, but the data has 2 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#> a single row.
1.4 ggsurvfit
Code
library(ggsurvfit)
km_fit
#> Call: survfit(formula = Surv(Days, status) ~ treatment, data = d)
#>
#> n events median 0.95LCL 0.95UCL
#> treatment=CON 14 14 17.0 15 25
#> treatment=LDRT 14 14 19.5 18 27
#> treatment=DPVB 14 9 38.5 30 NA
#> treatment=LR_DPVB 14 2 NA NA NA
km_fit |> tidy_survfit() |>
ggplot(aes(time,estimate,min=conf.low,ymax=conf.low,
color=strata,fill=strata))+
scale_color_jco()+
geom_step()
Code
km_fit |> ggsurvfit()+
add_confidence_interval()+
# 添加删失点
add_censor_mark() +
add_quantile()+
add_pvalue("annotation",size=3,x=5,y=0.25)+
#add_pvalue(caption = "Log-rank {p.value}")+
add_risktable(size = 3,risktable_height = 0.3)+
add_legend_title("treatment")