R软件一键绘制Table，解放医生双手

本章思维导图：

前言

上期介绍了SAS macro自动绘制SCI表格，10小时出的表格瞬间逼格提升到10分钟，战力力出众。但是由于SAS越发的小众，显得只有SASer们才自嗨。如今医学研究时代潮流中，R、SPSS、python逐步成为广大医学科研工作者的左右手，尤其是R。全面能打的R堪称医学统计的六边形战士（图1）。

图1 R全能王

如何才能用最低成本，最省时间，输出SCI规范用最美表格？答案只有，R。

本章将陪伴各位小伙伴，开启冒险打怪升级的新一站：利用R一键自动化生成SCI表格，解放各位的双手。

 自动Table包

累赘的、陈旧的、民工手动挡dplyr等方法就不多说了，一切以节约医生时间为金标准。直接上硬菜。

R中自动Table包有几个，包括Tableone包，finalfit包，Hmisc包(大牛Frank Harrel制作)。2014年面试鼎鼎大名的Tableone，BMJ子刊统计方法还为它专门建立shiny网页版（图2）https://esm.uoi.gr/shiny/tableone/，¹ 笔者实操体验卡顿，但偶尔用一次还是不错。上传本地xlsx等文件，网页会绘制好表格你再下载，整体色彩搭配都很美。Tableone做基线表其实还不错，但更绝佳的应用场景是倾向评分分析匹配前后基线均衡性，见后面倾向评分章节。下面直接进入本章核心部分，基线表和主分析表主推的都是finalfit包，2018年末推出的专门为SCI论文而生的超级表格强者。

finalfit包不光可以做所有的单变量，分组，按level构成比排序的基线table，还可以做单因素logistic、cox分析等，甚至多因素logistic、cox、多水平模型均可完成，相乘交互作用、相加交互作用都不是问题，甚至多个图表还能合并，简直是Table变形金刚。虽然相应的建模细节会比专门的package粗糙一些，但是作为课题前期快速出图，最快速度获取项目idea，finalfit包不要说太牛逼。

图2 Tableone shiny网页版

实战R表格

本章依旧使用rhc文件作为讲解，rhc.csv下载链接详见上一章节《SAS论文表格自动宏》。导入R后为rhc.rdata文件，在Rstudio中格式如下图，为了方便大家重现rhc，本章节R中附带了导入程序。R采用最新R4.1，Rstudio 1.4版本撰写。

一键绘表R程序，红色部分为核心语法，灰色部分酌情替代变量即可。

图3 rhc文件Rstudio格式

#需要 dplyr  finalfit knitr 各类package混搭

#rhc数据导入

setwd("D:/写书/R绘制SCI表格")

load("D:/写书/R绘制SCI表格/rhc.rdata")

#可在线下载

#rhc <- read.csv("https://biostat.app.vumc.org/wiki/pub/Main/DataSets/rhc.csv")

#-----------直接撸自动Table

#数据准备，部分连续变量，分类变量

conVars <- c("age","edu","das2d3pc","surv2md1")

catVars <- c("sex","race","income","swang1","death")

Vars<-c(convars,catvars)< span="">

#把结局y拎出来

explanatory <-  c(conVars,"sex","race","income","swang1")

dependent <- "death"

#分类变量转为→factor因子，后面dummy哑变量用

#rhc$death<- as.factor(rhc$death)

rhc[catVars] <- lapply(rhc[catVars], factor)

#--------方法① finalfit包

pacman::p_load(dplyr,knitr,Tableone,finalfit,Hmisc)

#基线;knitr:: kable为快捷应用knitr包对应函数

#亚组分析filter()，特定变量select()

#当存在单元格《5时候,分类fisher检验：p_cat = "fisher"。

#缺失值分析

rhc%>% missing_pattern(dependent, explanatory)

#summary_factorlist 进行绘表；默认正态分布，mean±sd，%;用F 和 X2 分析;digits = c(1,1,3,1)可设置小数位。

rhc %>% summary_factorlist(dependent, explanatory,p=TRUE,p_cat = "fisher") -> baselinenorm

knitr:: kable(baselinenorm, format = "pandoc", row.names=FALSE,

                col.names = c("", "", "Live", "Death","p-value"))

#两组，t检验，p_cont_para = "t.test"

#连续变量均为偏态

library(pastecs)

by(rhc[conVars],rhc$death,function(x)stat.desc(x,norm=TRUE))

#用非参数检验，cont = "median"对应 Kruskal-Wallis/Mann-Whitney U

rhc %>%  summary_factorlist(dependent, explanatory,p=TRUE,cont = "median",p_cat = "fisher") -> baselineskew

  knitr:: kable(baselineskew, format = "pandoc", row.names=FALSE,

                col.names = c("", "", "Live", "Death","p-value"))

#假如连续变量中仅仅edu surv2md1为偏态，cont_nonpara = c(2,4)序号进行非参数描述,。

rhc %>%  summary_factorlist(dependent, explanatory,p=TRUE,p_cat = "fisher",cont_nonpara = c(2,4)) ->  baselinemix

  knitr:: kable(baselinemix, format = "pandoc", row.names=FALSE,

                col.names = c("", "", "Live", "Death","p-value"))

#total_col = TRUE 列出total mean

rhc %>%  summary_factorlist(dependent, explanatory,p=TRUE,p_cat = "fisher",cont_nonpara = c(2,4), column = TRUE, total_col = TRUE) %>%

  knitr:: kable(., format = "pandoc", row.names=FALSE,

                col.names = c("", "", "Live", "Death","Total","p-value"))

#=====基线表最优代码，考虑miss; 分类水平按%大小排序,orderbytotal

#-----直接运行这段就是最完整的基线表

rhc %>% summary_factorlist(dependent, explanatory,p=TRUE,p_cat = "fisher",cont_nonpara = c(2,4), column = TRUE,total_col = TRUE, orderbytotal = TRUE,  na_to_p = TRUE,na_include = TRUE,add_row_total = TRUE, include_row_missing_col = FALSE) -> baseline

 knitr:: kable(baseline, format = "pandoc", row.names=FALSE)

#=====单因素logistic；多因素enter法

rhc %>%  finalfit(dependent, explanatory) -> univariatelogistic

knitr::kable(univariatelogistic, align=c("l", "l", "r", "r", "r"))

#=====多因素logistic reduced model; %>% ff_metrics()或者metrics = TRUE获取模型AIC C HL等参数，通用

#手动输入变量多因素；keep_models = TRUE；暂无stepwise过程

explanatory_reduce <-  c("sex","race","swang1")

rhc %>% finalfit(dependent, explanatory, explanatory_reduce,  metrics = TRUE,digits= c(2,2,3)) -> mutilogistic

knitr::kable(mutilogistic, align=c("l", "l", "r", "r", "r"))

#== OR 多因素enter法森林图形式；or_plot, hr_plot and surv_plot

rhc%>% or_plot(dependent, explanatory)

#========单因素cox表格；enter法多因素cox

explanatory <-  c(conVars,"sex","race","income","swang1")

death01 <- ifelse(rhc$death=="Yes",1,0)

dependent = "Surv(time, death01)"

rhc %>%  finalfit(dependent, explanatory) -> unicox

knitr::kable(unicox, align=c("l", "l", "r", "r", "r"))

#== HR 多因素enter法森林图形式

rhc%>% hr_plot(dependent, explanatory)

#finalfit包还可以做linear regression，Mixed effects random-intercept model，Bayesian logistic；可以finalfit_merge()多表格合并

#还可以分组km曲线带number at risk

explanatory <-c(< span="">"sex")

rhc %>% surv_plot(dependent, explanatory, xlab="Time (days)", pval=TRUE, legend="none")

#---上述结果导出到xlsx文件编辑；同理其他表格；这个从R console无缝链接超级实用。

library('openxlsx')

write.xlsx(baseline, file = "baseline.xlsx", colNames = TRUE, borders = "surrounding")

write.xlsx(univariatelogistic, file = "univariatelogistic.xlsx", colNames = TRUE, borders = "surrounding")

#-----其他细节

#移除p值；%>% ff_remove_p()

#有相乘交互作用，仅需要sex*race，而不需要额外再写sex race；race3分类，这里生产2行交互项

explanatory_interactive <-  c("sex*race","swang1")

rhc %>%  finalfit(dependent, explanatory, explanatory_interactive,  metrics = TRUE) -> logistic_inter

logistic_inter

#测试校正其他变量与sex swang1 交互，假设校正conVars，注意不能包含交互变量

rhc[catVars] <- lapply(rhc[catVars], factor)

dependent <- "death"

inter_cov <-  c(conVar,"race","income")

#叉生分析,相加交互，注意交互项变成sex_swang1

rhc %>%  ff_interaction(sex, swang1) %>%

  finalfit(dependent, c(inter_cov, "sex_swang1"),  metrics = TRUE) -> logistic_inter2

logistic_inter2

#--------方法②tableone 包

#载入包

pacman::p_load(dplyr,knitr,Tableone)

#快速获取变量名向量c("")

dput(names(rhc))

conVars <- c("age","edu","das2d3pc","surv2md1")

catVars <- c("sex","race","income","swang1","death")

Vars<-c(convars,catvars)< span="">

#单变量描述默认正态

tone_univariate <- CreateTableOne(vars = Vars,  data = rhc, test = FALSE)

summary(tone_univariate)

tone_univariate2 <- print(tone_univariate, printToggle = FALSE, noSpaces = TRUE,showAllLevels = TRUE, formatOptions = list(big.mark = ","))

kable(tone_univariate2, format = "pandoc", caption = "tableone Univariate table")

#基线变量按照y分组描述,全正态；默认F X2

tone_bivariate_test <- CreateTableOne(vars = Vars,strata = "death", data = rhc, test = T)

tone_bivariate_test2 <- print(tone_bivariate_test, printToggle = FALSE, noSpaces = TRUE,showAllLevels = TRUE, formatOptions = list(big.mark = ","))

kable(tone_bivariate_test2, format = "pandoc", caption = "tableone Bivariate table with statistical test")

#连续变量正态偏态分布测试

library(pastecs)

by(rhc[conVars],rhc$death,function(x)stat.desc(x,norm=TRUE))

#假设有非正态c("edu","surv2md1")；exact test，nonnormal用KW(MW)，单元格<5 fisher x2。

nonnormal_var <- c("edu","surv2md1")

#quote = TRUE 则需要手动复制console内容到excel

tone_outexcel <- print(tone_bivariate_test, nonnormal = nonnormal_var, exact = "race", quote =F, noSpaces = TRUE, printToggle = FALSE, showAllLevels = TRUE)

tone_outexcel

#表格数据导出csv，不能xlsx会漏列

write.csv(tone_outexcel, file = "Tableone_bivariate.csv")

#--------方法③ Hmisc包

library(Hmisc)

#带检验统计量与p

tab_hmisc<- summaryM(age+edu+das2d3pc+sex+race+swang1~ death, data = rhc,test = T)

print(tab_hmisc, long = T,what = "%",  prmsd = T, round = 2)

由于Hmisc重在线性模型，在一键成表方面逊色很多，不过多介绍。

 小结

一键出SCI基线表，课题快速预览统计表， finalfit包可以说是至强唯一。

解放医生们的双手，让科研更有效率，就请盘它吧。

原创不易，转载请说明来自本公众号。

参考文献

1.       Panos A, Mavridis D. TableOne: an online web application and R package for summarising and visualising data. Evid Based Ment Health 2020;23(3):127-130.doi:10.1136/ebmental-2020-300162