R语言生存分析函数笔记
A：以下自定义函数为借鉴团队的多个笔记，经过反复调试可以使用的，但部分标签位置可能不恰当，需要自定调整
B：这些函数可能用不到，里面的一些数据处理函数可以借鉴
C：如有疑问随时留言，看到了会解答，如果有用就很荣幸
1、单因素逻辑回归数据结果提取：
传入参数：数据集，变量名称列表，输出：变量回归结果的OR值、95%置信区间、P值
参数解释：
var_list：变量名称列表
var_result：因变量
glm_univar_res=function(data_base,var_list){
glm_model_res <- map_dfr(var_list,~ {
fml <- as.formula(paste0("var_result~", .x))
glm(fml,
family=binomial,data=data_base) %>%
# 提取OR值、置信区间
tidy( exponentiate = TRUE, conf.int = TRUE) %>%
filter(term != "(Intercept)") %>% # 去掉截距
mutate(variable = .x) %>% # 记录变量名
#选择并重命名变量
select(variable, term, OR = estimate, `95% CI lower` = conf.low,
`95% CI upper` = conf.high, p.value)
})
write.xlsx(glm_model_res,'glm_model_res.xlsx')
}
2、生存曲线绘制通用函数-不分组：
参数解释：
data_base：数据集
time_base：生存时间变量
status_base：生存状态变量
time_tag：指定时间生存概率的纵向标记线，接受的值可以是单个数值或者向量
#KM曲线不分组-----------------------------------------------------
km_plot_notgroup=function(data_base,time_base,status_base,time_tag){
data_base$surv_time <- data_base %>% pull({{time_base}})
data_base$surv_status <- data_base %>% pull({{status_base}})
#在survival包中先使用Surv()函数创建生存对象，生存对象是将事件、时间和删失信息合并在一起的数据结构。
data_km_fit_res <- survfit( Surv(surv_time,surv_status) ~ 1, data =data_base)
# ggsurvplot(data_km_fit_res,
# 获取生存函数的摘要信息
data_summ_km_fit_res <- summary(data_km_fit_res)
#在临床上有意义的生存时间点
#计算并提取指定时间下生存数据
data_summary_specify_time_res <- summary(data_km_fit_res,time=time_tag)
data_summary_specify_time_base_res <- data_summary_specify_time_res[c("time", "surv", "lower", "upper")] %>%
data.frame() %>%
mutate(time_tag1=str_c('(95% CI ',str_c(round(lower*100,digits=2),round(upper*100,digits=2),sep='-'),')'),
time_tag2=str_c(percent(surv,digits=2),time_tag1,sep='\n'))
#KM图开始绘制
data_sur_plot_base_res <- ggsurvplot(data_km_fit_res,
data =data_base,
# pval = sprintf("Log-rank test p = %.3f", log_rank_test$pvalue), # 显示log-rank检验p值
pval = FALSE,#组间比较多P值显示
pval.size=3,#组间比较多P值大小
pval.coord=c(0.5,0.2),#组间比较多P值坐标
conf.int = FALSE,# 生存曲线图上显示置信区间
risk.table = 'nrisk_cumcensor',# 显示风险表，可使用的参数TRUE, FALSE, 'absolute', 'percentAge', 'nrisk_cumcensor', 'nrisk_cumevents'
risk.table.pos='out',#风险表在图中的位置
risk.table.title='Number at risk(number censored)',#风险表的标题
tables.y.text=FALSE,#风险表分组标签是否显示在左侧
tables.y.text.col=TRUE,#风险表标记图标的颜色是否跟随strata
# tables.col='strata',#风险数量颜色调整
# tables.theme = theme_void(),
tables.theme = theme_void() +#去除风险表的背景、网格线、坐标
theme(
plot.margin = margin(5, 5, 5, l=15, "mm"),
axis.title.y = element_blank()
),#调整风险表与图形两侧的距离
# ggtheme = theme(axis.line = element_line(size=1),
# axis.title.y = element_text(size=15),
# panel.background = element_blank(), # 去除绘图区背景
# panel.grid.major = element_blank(), # 去除主网格线
# panel.grid.minor = element_blank()),# 去除次网格线
palette = '#4877B9',# JCO杂志配色
# surv.median.line = 'v',
legend.title = "",#将starta图例设置为空
# legend.labs =c("PN","non-PN"),#更改主图例
# legend.labs ="",#更改主图例
legend='none',
# title = 'Kaplan-Meier Survival Curve by Treatment',
# censor.size=6,
ylab='Event-Free Survival(%)',
xlab = 'Time (months)',
break.x.by=3,#x轴的间隔长度
break.y.by=0.25,#y轴的间隔长度，这是默认值，可以不设置，后面还要调整轴的标签
# ylim=c(0,1),#设置y轴显示范围，可以不设置，后面还要调整轴的标签
size=1.1,#生存曲线的粗细
fontsize=4,#风险表数据的大小
# cumcensor =TRUE,#显示累及删失数量表
surv.scale='percent'#生存概率y轴的刻度显示default，percent，但后面会调整刻度值的标签，所以这里设置无效
)
#拼图scale_y_continuous更改y轴坐标显示，geom_segment添加指示线，geom_label添加坐标标签，annotate指定点添加标签不需要数据集的考虑
data_sur_plot_con_res <- data_sur_plot_base_res$plot+
scale_y_continuous(
labels =seq(0, 100, 25),
position = 'left',#y轴的位置
limits = c(0, NA), # 最小值=0，最大值自动计算
expand = expansion(mult = c(0, 0))# 去除上下扩展空白
)+
#指定时间点线条
geom_segment(
data = data_summary_specify_time_base_res,
aes(x = time, xend =time,y = 0, yend = surv),linetype = "dashed")+
theme(plot.margin = margin(1, 0.5, 1, 0.5, "cm"),
axis.line = element_line(size=1),
axis.ticks = element_line(size = 1,color = "black"),# 可选：显示刻度线
axis.ticks.length = unit(0.2, "cm"),
# axis.title.y = element_text(size=15),
panel.background = element_blank(), # 去除绘图区背景
panel.grid.major = element_blank(), # 去除主网格线
panel.grid.minor = element_blank())+# 去除次网格线
# 指定时间生存概率，标签避让函数，但效果一般
geom_text_repel(
data=data_summary_specify_time_base_res,
aes(x = time,y = surv,label = time_tag2),
direction = "both", # 优先垂直方向调整
# nudge_y = 0.1, # 初始偏移量
# nudge_x = 0.1, # 初始偏移量
max.overlaps = 20,
min.segment.length = Inf # 最短引线长度
)+
#中位PFS或res标签添加
annotate(geom='text',x = 3,
y = 0.3,size=4,#注意x、y值的调整，影响中位res值的展示
# y = if_else(is.na(data_summ_km_fit$table['0.95LCL']/100)==TRUE,0.15,data_summ_km_fit$table['0.95LCL']/100),size=3,
label = str_c('m',str_sub(deparse(substitute(time_base)),1,3),' (95%CI):',coalesce(as.character(data_summ_km_fit_res$table['median']),"NR"),'month','(',
coalesce(as.character(data_summ_km_fit_res$table['0.95LCL']),"NR"),'-',
coalesce(as.character(data_summ_km_fit_res$table['0.95UCL']),"NR"),')'),na.rm =FALSE)
#拼接生存曲线调整图和调整后的风险表
grid.arrange(
data_sur_plot_con_res,
data_sur_plot_base_res$table,
ncol = 1,
heights = c(4, 1),
padding = unit(2, "line")
)
}
3、生存曲线绘制通用函数-分组：
参数解释：
data_base：数据集
time_base：生存时间变量
status_base：生存状态变量
group_var：分组变量
time_tag：指定时间生存概率的纵向标记线，接受的值可以是单个数值或者向量
groupA、groupB：分组变量唯一值自定义标签
km_plot_group=function(data_base,time_base,status_base,group_var,time_tag,groupA,groupB){
# data_obj_res <- Surv(time = pull(data_base,{{time_base}}),event = pull(data_base,{{status_base}}))
data_base$surv_time <- data_base %>% pull({{time_base}})
data_base$surv_status <- data_base %>% pull({{status_base}})
data_base$surv_group <- data_base %>% pull({{group_var}})
data_survobj_kmfit_res <- survfit(Surv(surv_time,surv_status) ~ surv_group, data =data_base)
# 获取生存函数的摘要信息
data_survobj_kmfit_summ_res <- summary(data_survobj_kmfit_res)
data_survobj_kmfit_summ_res$table
#在临床上有意义的生存时间点
#计算并提取指定时间下生存数据
data_kmfit_specify_time_res <- summary(data_survobj_kmfit_res,time=time_tag)
data_kmfit_specify_time_res_extract <- data_kmfit_specify_time_res[c("time", "surv", "lower", "upper")] %>%
data.frame() %>%
mutate(time_tag1=str_c('(95% CI ',str_c(round(lower*100,digits=2),round(upper*100,digits=2),sep='-'),')'),
time_tag2=str_c(percent(surv,digits=2),time_tag1,sep='\n'))
#风险比HR计算
data_surv_grouphr <-coxph(Surv(surv_time,surv_status) ~surv_group,
data = data_base)
data_surv_grouphr_ci <- summary(data_surv_grouphr)
data_surv_grouphr_ci_pvalue=if_else(data_surv_grouphr_ci$coefficients[,'Pr(>|z|)']<0.001,'),p <0.001',
str_c('),p =',sprintf("%.3f",data_surv_grouphr_ci$coefficients[,'Pr(>|z|)'])))
data_grouphr_ci <- str_c('Hazard Ratio:',
sprintf("%.2f",data_surv_grouphr_ci$coefficients[,'exp(coef)']),
'(95% CI,',
str_c(sprintf("%.2f",data_surv_grouphr_ci$conf.int[,'lower .95']),',',
sprintf("%.2f",data_surv_grouphr_ci$conf.int[,'upper .95'])),data_surv_grouphr_ci_pvalue)
#log-rank对数秩检验
data_log_rank_test <- survdiff(Surv(surv_time,surv_status) ~surv_group, data=data_base)
data_survobj_kmfit_plot_main_res <- ggsurvplot(data_survobj_kmfit_res,
data =data_base,
# pval = TRUE, # 显示log-rank检验p值
# pval = sprintf("Log-rank test p = %.3f", data_log_rank_test$pvalue), # 显示log-rank检验p值
# pval = "Log-rank test p <0.001",
pval.size=4,
pval.coord=c(1.5,0.25),
conf.int = FALSE,# 显示置信区间
risk.table = 'nrisk_cumcensor',# 显示风险表，可使用的参数TRUE, FALSE, 'absolute', 'percentage', 'nrisk_cumcensor', 'nrisk_cumevents'
tables.y.text=FALSE,
tables.theme = theme_void(),
risk.table.pos='out',#风险表在图中的位置
# palette = "hue",# JCO杂志配色
palette = c('#009FC3', '#B30437'),
# surv.median.line = 'v',
legend.title = "",#将starta图例设置为空
# legend.labs =c("0","1"),#更改主图例
legend.labs =c(groupA,groupB),#更改主图例
legend= c(0.9,0.9),
# cumcensor=TRUE,
# title = 'Kaplan-Meier Survival Curve by Treatment',
ylab='Progression-free survival(%)',
xlab = 'Time (months)',
break.x.by=3,#x轴的间隔长度
break.y.by=0.25,
ylim=c(0,1),
fontsize=4,
# cumcensor =TRUE
tables.y.text.col=TRUE,#风险表标记图标的颜色是否跟随strata
# tables.col='strata',#风险数量颜色调整
surv.scale='percent',#生存概率y轴的刻度显示default，percent
risk.table.title='Number at risk(number censored)')
#拼图scale_y_continuous更改y轴坐标显示，geom_segment添加指示线，geom_label添加坐标标签，annotate指定点添加标签不需要数据集的考虑
data_sur_plot_con_res <- data_survobj_kmfit_plot_main_res$plot+
scale_y_continuous(
labels =seq(0, 100, 25),
position = 'left',
limits = c(0, NA), # 最小值=0，最大值自动计算
expand = expansion(mult = c(0, 0))# 去除上下扩展空白
)+
geom_segment(
data = data_kmfit_specify_time_res_extract,
aes(x = time, xend =time,y = 0, yend = surv),linetype = "dashed")+
theme(plot.margin = margin(1, 0.5, 1, 0.5, "cm"))+ #调整图形与图片的边距
# 这个标签避让函数，默认参数比手动调整设置的效果好一些
geom_text_repel(
data=data_kmfit_specify_time_res_extract,
aes(x = time,y = surv,label = time_tag2),
direction = "both", # 优先垂直方向调整
max.overlaps = 20,
min.segment.length = Inf # 最短引线长度
)+
annotate(geom='text',x = 3,
y = 0.2,size=4,label=data_grouphr_ci)+
annotate(geom='text',x = 3,
y = 0.15,size=4,
label = str_c(groupA,' :m',str_sub(deparse(substitute(time_base)),1,3),' (95%CI):',coalesce(as.character(data_survobj_kmfit_summ_res$table[,'median'][1]),"NR"),'month','(',
coalesce(as.character(data_survobj_kmfit_summ_res$table[,'0.95LCL' ][1]),"NR"),'-',
coalesce(as.character(data_survobj_kmfit_summ_res$table[,'0.95UCL' ][1]),"NR"),')'),na.rm =FALSE)+
annotate(geom='text',x = 3,
y = 0.1,size=4,
label = str_c(groupB,' :m',str_sub(deparse(substitute(time_base)),1,3),' (95%CI):',coalesce(as.character(data_survobj_kmfit_summ_res$table['median'][2]),"NR"),'month','(',
coalesce(as.character(data_survobj_kmfit_summ_res$table[,'0.95LCL' ][2]),"NR"),'-',
coalesce(as.character(data_survobj_kmfit_summ_res$table[,'0.95UCL' ][2]),"NR"),')'),na.rm =FALSE)
#拼接生存曲线调整图和调整后的风险表
grid.arrange(
data_sur_plot_con_res,
data_survobj_kmfit_plot_main_res$table,
ncol = 1,
heights = c(4, 1),
padding = unit(2, "line")
)
}
4、疗效评价统计通用函数：
data：数据集
group_var：分组变量，如果需要可以传入，不需要可以不传
A：针对数据集中包含NE、missing均有判断，可以直接放入数据集中，若不需要将该数据纳入统计，可以在传入数据集前过滤
response_statistic_group=function(data,group_var){
#
data_response <- data %>%
mutate(
#此处调整增加变量的因子水平
best_response= {
factor(best_response, levels =c("CR", "PR", "SD", "PD") %>%
c(if("NE" %in% best_response) "NE",
if("missing" %in% best_response) "missing"))
}) %>%
#原始数据没有相关评价水平时，这里设置水平后可以展示相应情况，对于全部都有时没有影响
tbl_summary(statistic =list(all_categorical() ~ "{n}({p}%)"),
type=list(best_response~'categorical'),
digits =~c(p=2),
by={{group_var}},
missing_text='missing') %>%
modify_header(all_stat_cols() ~ "Total(N = {n})") %>%
remove_footnote_header(columns = all_stat_cols()) %>%
add_stat_label() %>%
add_overall()
if(missing(group_var)){
data_ci_list <- data_response$table_body %>%
mutate(stat_0_split=stri_extract_first_regex(stat_0,'\\d+') %>% as.integer())
data_ci_list_orr <-data_ci_list %>%filter(label %in% c("CR","PR"))%>% summarise(across(ends_with('split'),.fns=sum))
data_ci_list_dcr <-data_ci_list %>% filter(label %in% c("CR","PR","SD"))%>% summarise(across(ends_with('split'),.fns=sum))
data_ci_list_all <-data_ci_list %>% tail(-1)%>% summarise(across(ends_with('split'),.fns=sum))
data_ci_orr_res <- BinomCI(data_ci_list_orr[[1]],data_ci_list_all[[1]], method = 'clopper-pearson') %>% t()%>%
data.frame() %>%
tail(-1) %>%
mutate(across(everything(),~percent(.x,2))) %>%
sapply(., function(x) paste0(x,collapse= "-")) %>%
upData(rename = c(.='stat_0'))
data_ci_dcr_res <- BinomCI(data_ci_list_dcr[[1]],data_ci_list_all[[1]], method = 'clopper-pearson') %>% t() %>%
data.frame() %>%
tail(-1) %>%
mutate(across(everything(),~percent(.x,2))) %>%
sapply(., function(x) paste0(x,collapse= "-")) %>%
upData(rename = c(.='stat_0'))
#将orr、dcr数据合并进影像学评价总表table_body里面
data_response$table_body <- bind_rows( data_response$table_body,
mapply(`str_c`, data_ci_list_orr,'(',percent(data_ci_list_orr/data_ci_list_all,2),')') %>%
upData(rename = c(stat_0_split='stat_0')),
data_ci_orr_res,
mapply(`str_c`, data_ci_list_dcr,'(',percent(data_ci_list_dcr/data_ci_list_all,2),')') %>%
upData(rename = c(stat_0_split='stat_0')),
data_ci_dcr_res) %>%
mutate(across(everything(), ~str_remove_all(.,"%"))) %>%
rownames_to_column("row_num") %>%
mutate(
label = replace(label,(n_distinct(label, na.rm = TRUE) + 1):n(),values = c("ORR", "95% CI", "DCR", "95% CI"))) %>%
select(-row_num)
}else{
#计算ORR、DCR的95%置信区间
data_ci_list <- data_response$table_body %>%
mutate(stat_0_split=stri_extract_first_regex(stat_0,'\\d+') %>% as.integer(),
stat_1_split=stri_extract_first_regex(stat_1,'\\d+') %>% as.integer(),
stat_2_split=stri_extract_first_regex(stat_2,'\\d+') %>% as.integer()
)
data_ci_list_orr <-data_ci_list %>%filter(label %in% c("CR","PR"))%>% summarise(across(ends_with('split'),.fns=sum))
data_ci_list_dcr <-data_ci_list %>% filter(label %in% c("CR","PR","SD"))%>% summarise(across(ends_with('split'),.fns=sum))
data_ci_list_group <-data_ci_list %>% tail(-1)%>% summarise(across(ends_with('split'),.fns=sum))
data_ci_orr_res <- do.call(cbind,lapply(1:3, function(x) BinomCI(data_ci_list_orr[[x]],data_ci_list_group[[x]], method = 'clopper-pearson') %>% t())) %>%
data.frame() %>%
tail(-1) %>%
mutate(across(everything(),~percent(.x,2))) %>%
sapply(., function(x) paste0(x,collapse= "-")) %>%
upData(rename = c(X1='stat_0', X2='stat_1',X3='stat_2'))
data_ci_dcr_res <- do.call(cbind,lapply(1:3, function(x) BinomCI(data_ci_list_dcr[[x]],data_ci_list_group[[x]], method = 'clopper-pearson') %>% t())) %>%
data.frame() %>%
tail(-1) %>%
mutate(across(everything(),~percent(.x,2))) %>%
sapply(., function(x) paste0(x,collapse= "-")) %>%
upData(rename = c(X1='stat_0', X2='stat_1',X3='stat_2'))
#将orr、dcr数据合并进影像学评价总表table_body里面
data_response$table_body <- bind_rows( data_response$table_body,
mapply(`str_c`, data_ci_list_orr,'(',percent(data_ci_list_orr/data_ci_list_group,2),')') %>%
upData(rename = c(stat_0_split='stat_0',
stat_1_split='stat_1',
stat_2_split='stat_2')),
data_ci_orr_res,
mapply(`str_c`, data_ci_list_dcr,'(',percent(data_ci_list_dcr/data_ci_list_group,2),')') %>%
upData(rename = c(stat_0_split='stat_0',
stat_1_split='stat_1',
stat_2_split='stat_2')),
data_ci_dcr_res) %>%
mutate(across(everything(), ~str_remove_all(.,"%"))) %>%
rownames_to_column("row_num") %>%
mutate(
label = replace(label,(n_distinct(label, na.rm = TRUE) + 1):n(),values = c("ORR", "95% CI", "DCR", "95% CI"))) %>%
select(-row_num)
}
#将标签的首字母都调整为大写
data_response$table_body <- data_response$table_body %>%
mutate(var_label=upFirst(data_response$table_body$var_label),
label=upFirst(data_response$table_body$label))
#将数据转换为三线表格式，可用在学术内容中的格式
data_response_flex <- data_response %>%
as_flex_table()%>%
set_table_properties(layout = "autofit") %>%
style(j=1:ncol(data_response),pr_t = fp_text(font.family = "Times New Roman", # 字体
bold = FALSE # 不加粗
)) %>%
border(part = "header",border.top= fp_border(color = "black", width = 1.5)) %>%
hline_bottom(border= fp_border(color = "black", width = 1.5))
save_as_docx(data_response_flex, path = "data_response_flex.docx")
}
扫一扫
453
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
