Topic 9. SCI 文章第二张表—单因素回归分析表

在SCI论文中，我们经常可以看见一些这样的表格，大多数命名表格 2.，主要用来表示原因和结果的单因素分析的关系或者是分组变量的关系，也就是单因素分析，那这样一张表格该怎么完成呢？

01 单因素回归分析方法

单因素回归分析就是在构建回归模型时，只纳入一个因素进入到回归模型中进行拟合，其理论上也应该属于单因素分析的范畴，只不过是用到了稍微高级一点的回归方法而已。

单因素回归分析有三种方法，那么三种方法该怎么选择，以及什么样的数据适合哪种回归分析方法，这其中都存在很复杂的过程，但是经我细细的道来，自然就明了许多，单因素回归分析三种方法包括如下：  

1. 单因素 Cox 回归

2. 单因素 Logistic 回归

3. 多线性回归

   先看第一种 Cox 回归分析，这种方法最常用，常常与 Logisitc 回归分析方法混淆，那么什么样的数据应当选择 Cox 回归分析方法呢？如下数据中的有 time-to-event, 对应着 sex 二分类方法，如下：

library("survival")
library("survminer")
#导入数据
data("lung")
head(lung)
  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      2  74   1       1       90       100     1175      NA
2    3  455      2  68   1       0       90        90     1225      15
3    3 1010      1  56   1       0       90        90       NA      15
4    5  210      2  57   1       1       90        60     1150      11
5    1  883      2  60   1       0      100        90       NA       0
6   12 1022      1  74   1       1       50        80      513       0

#COX分析需要的生存时间time和生存状态status
res.cox<- coxph(Surv(time, status)~sex,data = lung)
summary(res.cox)
Call:
coxph(formula = Surv(time, status) ~ sex, data = lung)

  n= 228, number of events= 165 

       coef exp(coef) se(coef)      z Pr(>|z|)   
sex -0.5310    0.5880   0.1672 -3.176  0.00149 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    exp(coef) exp(-coef) lower .95 upper .95
sex     0.588      1.701    0.4237     0.816

Concordance= 0.579  (se = 0.021 )
Likelihood ratio test= 10.63  on 1 df,   p=0.001
Wald test            = 10.09  on 1 df,   p=0.001
Score (logrank) test = 10.33  on 1 df,   p=0.001

02 批量单因素 Cox 回归

一般我们的关注的特征都比较多，用上面的代码一个一个来做单因素cox回归分析效率太低了，下面我们来看看如何批量做单因素cox回归分析。对于上面的肺癌数据，我们对除了time-to-event 之外的8特征做单因素cox回归分析，然后对数据进行整理，因为数据有缺失，去除缺失的数据，如下：

lung<-na.omit(lung)

对单因素数据结果进行整理，之后达到一张表格，如下：

#对除了time-to-event 之外的8特征做单因素cox回归分析
covariates <- c("inst","age", "sex",  "ph.karno", "ph.ecog", "wt.loss","meal.cal","pat.karno")
#分别对每一个变量，构建生存分析的公式
univ_formulas <- sapply(covariates,
                        function(x) as.formula(paste('Surv(time, status)~', x)))

#循环对每一个特征做cox回归分析
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})
univ_models
#提取HR，95%置信区间和p值
univ_results <- lapply(univ_models,
                       function(x){ 
                         x <- summary(x)
                         #获取p值
                         p.value<-signif(x$wald["pvalue"], digits=2)
                         #获取HR
                         HR <-signif(x$coef[2], digits=2);
                         #获取95%置信区间
                         HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
                         HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
                         HRa <- paste0(HR, " (", 
                                       HR.confint.lower, "-", HR.confint.upper, ")")
                         res<-c(p.value,HRa,HR,HR.confint.lower,HR.confint.upper)
                         # names(res)<-c("P-value","HR (95% CI for HR)","","","")
                         return(res)
                       })
#转换成数据框，并转置
res <- t(as.data.frame(univ_results, check.names = FALSE))
res<-as.data.frame(res)
names(res)=c("Variants","P-value","Hazard Ratio (95%CI)","","")
res
         Variants         P-value Hazard Ratio (95%CI)         
inst          0.18    1 (0.99-1.1)                    1 0.99 1.1
age          0.044   0.97 (0.94-1)                 0.97 0.94   1
sex          0.089    1.7 (0.93-3)                  1.7 0.93   3
ph.karno     0.063       1 (1-1.1)                    1    1 1.1
ph.ecog       0.18 0.73 (0.46-1.2)                 0.73 0.46 1.2
wt.loss       0.63   0.99 (0.97-1)                 0.99 0.97   1
meal.cal      0.73         1 (1-1)                    1    1   1
pat.karno     0.38      1 (0.99-1)                    1 0.99   1

write.table(file="univariate_cox_result.xls",res,quote=F,sep=",",row.names = T)

03 单因素回归结果可视化

对于这个结果一般来说文章分析中都会存在，有时是单独一张表格，有时是与多因素放在一起，不过我们这里先说的就是批量整理单因素回归分析的结果，可视化结果我们利用R包forestplot，如下：

#install.packages("forestplot")
library(forestplot)
rs_forest <- read.csv('univariate_cox_result.xls',header = FALSE)
rs_forest
# 读入数据的时候大家一定要把header设置成FALSE，保证第一行不被当作列名称〿
 tiff('Figure 2.tiff',height = 1800,width = 3000,res= 600)
forestplot(labeltext = as.matrix(rs_forest[,c(1,3,2)]),
           #设置用于文本展示的列，此处我们用数据的前四列作为文本，在图中展示
           mean = rs_forest$V4, #设置均倿
           lower = rs_forest$V5, #设置均值的lowlimits陿
           upper = rs_forest$V6, #设置均值的uplimits陿
           #is.summary = c(T,T,T,T,T,T,T,T,T),
           #该参数接受一个逻辑向量，用于定义数据中的每一行是否是汇总值，若是，则在对应位置设置为TRUE，若否，则设置为FALSE；设置为TRUE的行则以粗体出现
           zero = 1.5, #设置参照值，此处我们展示的是HR值，故参照值是1，而不昿0
           boxsize = 0.3, #设置点估计的方形大小
           lineheight = unit(8,'mm'),#设置图形中的行距
           colgap = unit(2,'mm'),#设置图形中的列间跿
           lwd.zero = 2,#设置参考线的粗绿
           lwd.ci = 2,#设置区间估计线的粗细
           col=fpColors(box='#458B00', summary= "#8B008B",lines = 'black',zero = '#7AC5CD'),
           #使用fpColors()函数定义图形元素的颜色，从左至右分别对应点估计方形，汇总值，区间估计线，参考线
           xlab="The estimates",#设置x轴标筿
           lwd.xaxis=2,#设置X轴线的粗绿
           lty.ci = "solid",
           graph.pos = 3)#设置森林图的位置，此处设置为4，则出现在第四列
dev.off()

怎么样，学会了吗？当然后面还会有多因素表格，但是我们一般会将单因素与多因素的表格合并为一张表格，下期将说说怎么自动化合并一张这么复杂的表并进行可视化展示，敬请期待！

关注公众号，每天有更新，记得扫码进群，每天有讨论，免费交流哦！

Reference:

1. Liu Y, Zhao P, Cheng M, Yu L, Cheng Z, Fan L, Chen C. AST to ALT ratio and arterial stiffness in non-fatty liver Japanese population:a secondary analysis based on a cross-sectional study. Lipids Health Dis. 2018 Dec 3;17(1):275.

2. Hosmer, D. & Lemeshow, S. (2000). Applied Logistic Regression (Second Edition). New York: John Wiley & Sons, Inc.

3. Long, J. Scott (1997). Regression Models for Categorical and Limited Dependent Variables. Thousand Oaks, CA: Sage Publications.

桓峰基因

生物信息分析，SCI文章撰写及生物信息基础知识学习：R语言学习，perl基础编程，linux系统命令，Python遇见更好的你

30篇原创内容

公众号