R语言高效数据处理-3个自定义函数笔记

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 =TRUEtables.y.text.col=TRUE,#风险表标记图标的颜色是否跟随strata# tables.col='strata',#风险数量颜色调整surv.scale='percent',#生存概率y轴的刻度显示default，percentrisk.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"))
}