R语言|14. 森林图-3: 亚组分析三线表及森林图

解释：亚组分析公式：

cox1<-coxph(Surv(time,status==0)~RT, subset= (R14[,i]=='1'), data=R14)，time和ststus是数据的随访时间和结局，0为感兴趣事件；RT是要进行亚组分析的变量；subset=(x=='1')是亚组，即做关于x变量中第1个亚组的RT的单因素分析；数据来源是R14。 

for循环-1：for(i in 5:dim(R14)[2])意思是从R14文件中每次提取1个变量纳入亚组分析的subset=(R14[,i]=='1')，其中R14[,i]为循环变量，并注明从第5个变量即T开始纳入。这样我们就可以得到T/ER/PR/KI67/LVI/N/Grade变量中第1亚组的分析结果。

aa<-summary(cox1)：将结果放入aa中，等待提取 。

for循环-2：for( j in 1:dim(aa$coefficients)[1])即将aa中的各种系数（coefficients）提出，从第1个开始提。rbind():将aa的第2列coefficients提出（是HR）,第2-3列的conf.int（是95%CI），保留两位小数；coefficients的第5列（P）提出并保留3位小数。然后放入a表。每次提取都这么干，直至结束。关于为何coefficients第2列为HR，请看下，这是coxph()汇报结果的格式。

#   coxph(formula = Surv(time, status == 0) ~ RT, data = R14)
#   n= 350, number of events= 55 
# 
#          coef   exp(coef)  se(coef)  z    Pr(>|z|)
# RTYes 0.0003434 1.0003435 0.2714870 0.001  0.999
# 
#        exp(coef) exp(-coef)  lower .95  upper .95
# RTYes   1         0.9997      0.5876      1.703
# 
# Concordance= 0.513  (se = 0.036 )
# Likelihood ratio test= 0  on 1 df,   p=1
# Wald test            = 0  on 1 df,   p=1
# Score (logrank) test = 0  on 1 df,   p=1

#3-a表加入第一列名字，
#即变量T/ER/PR/KI67/LVI/N/Grade中亚变量名为1的实际含义
rownames(a)<-c("Characteristics",
               "T1",
               "Negative",
               "Negative",
               "Low",
               "Absent",
               "N0",
               "I");a <- data.frame(a)

#4- 将a表的"HR","5%CI","95%CI"组合成HR(5CI-95CI)格式
#上图X1列是HR，X2,X3是95%CI
a$HR.CI95<-paste0(a$X1," (",a$X2,"-",a$X3,")");a

经过以上4步，a表就做完了。其他表是同理。接下来做所有亚变量为2的亚组分析。

#1- 建立b表
b<-c("HR","5%CI","95%CI","P")

#2- 所有变量为2的亚组分析for循环
for(i in 5:dim(R14)[2]){ 
  cox1<-coxph(Surv(time,status==0)~RT,subset=(R14[,i]=='2'),data=R14)
  bb<-summary(cox1)
  for(j in 1:dim(bb$coefficients)[1]){
    b<-rbind(b,c(round(bb$coefficients[j,2],2),
                 round(bb$conf.int[j,3],2),
                 round(bb$conf.int[j,4],2),
                 round(bb$coefficients[j,5],3)))}}

#3- b表第一列名字，即变量T/ER/PR/KI67/LVI/N/Grade中亚变量名为2的实际含义
rownames(b)<-c("Characteristics",
              "T2",
              "Positive",
              "Positive",
              "High",
              "Present",
              "N1",
              "II")
b <- data.frame(b);b

#4- 加入HR(95% CI)格式
b$HR.CI95<-paste0(b$X1," (",b$X2,"-",b$X3,")");b

b表：所有亚变量为2的亚组分析完成。红色是需要根据自己数据改的的地方，蓝色为主要的地方，虽然与a表的代码是一样的，但命名要分为b和bb，不然会把数据做到a表上，导致错误。

注意，此时只剩N和Grade有第3的亚组，所以for循环要从第10列开始（N变量）循环。这就是在进行亚组分析前将数据按变量多少依次排列的原因。红色根据自己数据改的，蓝色是注意不需修改。

#1- 建立c表
c<-c("HR","5%CI","95%CI","P")

#2- 亚变量为3的亚组分析
for(i in 10:dim(R14)[2]){ 
  cox2<-coxph(Surv(time,status==0)~RT,subset=(R14[,i]=='3'),data=R14)
  cc<-summary(cox2)
  for(j in 1:dim(cc$coefficients)[1]){
    c<-rbind(c,c(round(cc$coefficients[j,2],2),
                 round(cc$conf.int[j,3],2),
                 round(cc$conf.int[j,4],2),
                 round(cc$coefficients[j,5],3)))}}

#3- c表的第一列名字，即变量N/Grade中亚变量名为3的实际含义
rownames(c)<-c("Characteristics",
               "N2~3",
               "III")
c <- data.frame(c);c
#4- 加入HR（95% CI）列
c$HR.CI95<-paste0(c$X1," (",c$X2,"-",c$X3,")");c

c表只有两条结果，因为只有两个变量为3分类变量。要是还有4分类，大家现在应该就会了，再复制一下上面的循环代码，稍微修改一下即可。

这次仅仅是一次单因素分析而已，做对照用。非常简单，循环函数都不需要。

#1- RT单因素分析
cox3<-coxph(Surv(time,status==0)~RT,data=R14)
dd<-summary(cox3)

#2- 提取信息
dd$coefficients 
dd$conf.int    
HR<- round(dd$coefficients[,2],2)     
PValue <- round(dd$coefficients[,5],3) 
CI5<- round(dd$conf.int[,3],2)
CI95<- round(dd$conf.int[,4],2)
HR.CI95<-paste0(HR," (",CI5,"-",CI95,")")

#3- d表是最后建的，就是把信息提取出来组合为d表
d<- data.frame("X1" = HR,
               "X2" = CI5,
               "X3" = CI95,
               "X4" = PValue,
               "HR.CI95" = HR.CI95)
#4- 行名，命名为所有人群
rownames(d)<-"All parents"           

这样，所需的HR、95%CI、P值都完成了，接下来是按照变量的分门别类把他们排好顺序。注：下一步要是嫌麻烦完全可以手动合成的。

#1- 查看数据变量名
names(R14)
#"status" "time""RT" "AGE""T" "ER" "PR" "KI67""LVI" "N" "Grade" 
#2- 指定需进行基线统计的变量，纳入的变量为上面进行的亚组分析的变量
myVars <- c("T","ER" ,"PR" ,"KI67","LVI","N","Grade")
#指定基线表中哪些变量是分类变量
catVars<-c("T","ER" ,"PR" ,"KI67","LVI","N","Grade")
#3- 构建table1函数，RT为分类标准
table<- print(CreateTableOne(vars=myVars,
                              data=R14,
                              strata="RT",
                              factorVars=catVars),
               catDigits = 2,contDigits = 2,showAllLevels=TRUE)     

此时，我们提取到了在放疗或未放疗两组中的患者人数，是一个17行5列的表。但我们不要p值和text列，level列要从数字变为英文格式。

#4- 去掉第4.5列（上图p和test列），新表命名为table1（17行3列）
table1<-table[,c(-4,-5)]

注：函数a[x  ,  y]意思是a表的行（x）和列（y）。例如：b<-a[1,]提取第一行给b，b<-a[-1,]a表删除第一行命名为b。b<-a[,1]提取第一列给b，b<-a[,-1]a表删除第一列并命名为b表。b<-a[,c（-1，-3）]a表删除第1.3列并命名为b表。

#1- 按照table1的格式将abcd表合成一个表
res1<-rbind(d[1,],#abcd表合为res1表
            a[2,],b[2,],#T的亚变量
            a[3,],b[3,],#ER
            a[4,],b[4,],#PR
            a[5,],b[5,],#KI67
            a[6,],b[6,],#LVI
            a[7,],b[7,],c[2,],#N
            a[8,],b[8,],c[3,])#G

使用的是rbind()和a[x,]函数。rbind()意思为上下行粘贴，a[x, ]是按行提取。abcd组合结果如下：res1表（一个17行5列表）

#2- res1表列名进入表格内，并命名为Characteristics
res1<-tibble::rownames_to_column(res1, var = "Characteristics")

此时的res1表变为17行6列表（多了亚变量名列）

#3- table与res1合成为res2
res2<-cbind(table1,res1)

前面介绍rbind()为上下的行粘贴，cbind()函数为左右的列粘贴。此时我们把17行3列的table1和17行6列的res2组合为res2（17行9列）

#4- 第一列不再需要，删除（res3,17行8列）
res3<-res2[,-1]

#5- 将合成的res3进行列排序，亚变量名放首位
res4<-data.frame(res3[c(3,2,1,4:6,8,7)])

至此，快要完成了。现在只有亚变量名缺变量名；没有列名；p值有=0的，应该为<0.001; N分期要放T分期后面。

#1- 每组亚变量上一行插入他们的变量名，根据自己数据来/或保存为Excel手写
N<-rbind(res4[1, ], #提取res4表第一行
         c("T stage",rep(NA, 7)),#插入一行第一个单元格为T stage，后7列为NA
         res4[c(2:3), ],#提取res4表第2，3行
         c("ER stage",rep(NA, 7)),#插入一行第一个单元格为ER stage，后7列为NA
         res4[c(4:5),],#以此类推
         c("PR stage",rep(NA, 7)),
         res4[c(6:7),],
         c("KI67",rep(NA, 7)), 
         res4[c(8:9),],
         c("LVI",rep(NA, 7)),
         res4[c(10:11),],
         c("N stage",rep(NA, 7)), 
         res4[c(12:14),],
         c("Grade",rep(NA, 7)),
         res4[c(15:17),])

#2-插入行名
for(i in 2:8) {N[, i] = as.character(N[, i])}#先让变量性质变为character类型
result <- rbind( c("Characteristics","PMRT","No_PMRT",NA,NA,NA,"HR (95%CI)","P Value"),
                 N[c(1:25),])

#3- p=0的变为<0.001
result$X4[result$X4==0]<-"<0.001"

#4- N分期放在T分期下面
result<-data.frame(result,row.names=NULL);result
result <- rbind(result [c(1:5,18:21,6:17,22:25),])

放疗(RT)与否在各亚组获益情况的亚组分析三线表完成

fig<-forestplot(result[,c(1,2,3,7,8)], #12378列显示为原数字格式
                 mean=result[,4],   
                 lower=result[,5],  
                 upper=result[,6], 
                 zero=1,            
                 boxsize=0.4,      
                 graph.pos= 4,#图放在第四列
                 hrzl_lines=list("1" = gpar(lty=1,lwd=2),
                                 "2" = gpar(lty=2),
                                 "26"= gpar(lwd=2,lty=1)),
                 graphwidth = unit(.25,"npc"),
                 xlab="HR<1,意味着放疗能改善局部区域复发风险",
                 xticks=c(0,1,3) ,
#----------------#字体
                is.summary=c(T,
                             F,
                             T,F,F,
                             T,F,F,F,
                             T,F,F,
                             T,F,F,
                             T,F,F,
                             T,F,F,
                             T,F,F,F
                             ),#T=粗体
                 txt_gp=fpTxtGp(label=gpar(cex=1),
                                ticks=gpar(cex=1.1), 
                                xlab=gpar(cex=1), 
                                title=gpar(cex=2)),
#----------------#线条粗细（x轴、置信区间）
                 lwd.zero=2,
                 lwd.ci=2,
                 lwd.xaxis=1, 
                 lty.ci=1,
                 ci.vertices =T,
                 ci.vertices.height=0.2, 
                 clip=c(0,3),
#----------------#行间距、字间距/box形状                 
                 ineheight=unit(8, 'mm'), 
                 line.margin=unit(8, 'mm'),
                 colgap=unit(6, 'mm'),
                 col=fpColors(zero = "#e22e2a",
                              box = '#048410', 
                              lines = 'black'),
                 fn.ci_norm="fpDrawCircleCI", 
                 title="亚组分析森林图")

总结：至此亚组分析三线表及其森林图以全部完成。难点其实就是表格组合，个人认为手动组合快点。 

若使用R，用到基础函数有rbind()数据上下粘贴组合；cbind()是数据左右组合；a[行，列]可以对a表行列进行挑选和删除等，虽然过程有些繁琐，但是完成之后感觉自己对R语言又增进了许多了解。

亚组分析是用以单因素Cox回归为核心的for循环来批量提取关键信息，并结合tableone包提取人群在个亚组的数目。

本期用的所有代码其实在前面早已介绍，换汤不换药。我相信只要掌握了代码含义和使用方法，大家都能完美解决临床研究中的问题，为论文的统计分析节省大量时间。

感谢阅读，如有错误请指正！

喜欢的小伙伴还请点赞分享一波，谢谢支持！！！