生存率是基于KM曲线函数结果，关于某变量的KM曲线及结果已在 R语言|8. Kaplan-Meier曲线及美化和 R语言|10. KM曲线代码进阶：各种亚组分析助你深入分析数据详细介绍。

但是，有些小伙伴想做某变量在所有亚组下的生存率，那就需要我们做很多组KM曲线，例如上图8个变量，若做KM曲线就要做17幅，会占用大量文章面积，得不偿失。

而亚组分析森林图能清晰、明了的总结研究变量在各个亚组的n年生存率差异。结合上期介绍的HR结果将会使文章增色不少。

设计思路 (此处可以先跳过)

1、大框架：同基于HR的亚组分析一样，基于生存率的亚组分析三线表也是3大类：变量名、各亚组人数及OS% (95%CI)。[本例中亚组人群包括了实际人数和死亡人数。OS也是包括了放疗组和非放疗组 .] 大框架是一致。

2、核心函数代码：

• survfit( )：构建亚组分析生存函数模型，来源于survival包；
• survdiff( ): 提取log-rank test 的P值，来源于survival包；
• tbl_survfit( ): 提取亚组分析生存率和死亡人数，gtsummary包【神包，五星推荐】点击看介绍：R语言| 日常笔记1
• CreateTableOne( ):提取各亚组患者数，来源于tableone包；
• forestplot( ):绘制森林图，来源于forestplot包

3、辅助函数代码：

• rbind( ), 数据上下组合；
• cbind( ), 数据左右组合；
• str_remove_all( )，删除单元格中符合条件的函数；
• separate( ), 将一列数据按条件分为几列，它的反面是paste0(),将几列数据组合为1列【即HR（95%CI）的组合与拆分】
• a[行，列]，对行列进行提取或删除，b<- a[,-1]，a表删除第一列并命名为b表；b<- a[-1,]，a表删除第一行并命名为b表。

4、注意同第14期一样，亚变量要处理为1234的数字格式，变量2分类在前，3分类在后，以此类推。

5、步骤：构建for循环，用上述函数提取信息并组合成表。先收集受放疗组在亚变量为1的信息，然后是2,3...。再收集不接受放疗组在亚变量为1的信息，然后是2,3...。也就是说要建6-8个for循环。但代码是一摸一样的，只需修改接受提取信息表的名字、亚组名称和是否放疗即可。 

难点还是表的组合，可以手动组合，更快更简单。

#0.安装包
install.packages("x")#x为未安装的包
#1.载入包
library(survival)
library(plyr)
library(gtsummary)
library(stringr)
library(tidyverse)
##制作table1
library(tableone)
##制作森林图
library(forestplot)
#2.清理工作环境
rm(list = ls()) 
#3.数据放入工作路径，加载
R15<- read.csv('R15.csv')
#4.查看数据前6行
head(R15)
#5.查看数据数据性质
str(R15)
#6.查看结局:0为死亡（652人）
R15$status<-factor(R15$status)
summary(R15$status)

数据处理格式，亚变量处理为数字（研究变量rt可以不变），二、三分类变量排序

相关变量定义

rt：放疗 

t：1=T1，2=T2；3=T3

er: 1=阴性，2=阳性；

pr: 1=阴性，2=阳性；

her2: 1=阴性，2=阳性； 

sur(手术): 1=根治术 ,2=保乳术； 

che（化疗）：1=无，2=有； 

g（分级）: 1=I级，2=2级，3=级。

红色标记是根据自己数据修改的代码

#1.构建容纳死亡人数，OS% (5%CI,95%CI)及P值的表（a表）
a<-c("Event N","OS% (5%CI, 95%CI)","P")
#2for循环，提取信息至a表
for(i in 5:dim(R15)[2]){ 
  fit1<-survfit(Surv(time,status==0)~rt,subset=(R15[,i]==1),data=R15)
  p<-survdiff(Surv(time,status==0)~rt,rho=0,subset=(R15[,i]==1),data=R15)
  tbl1<-tbl_survfit(fit1,times=60)%>%add_nevent()
  a<-rbind(a,c(
      inline_text(tbl1,column = nevent),
      inline_text(tbl1,time =60, level = "Yes"),
      1-pchisq(p$chisq,length(p$n)-1)))}

红色是需要修改，蓝色需要注意。

• survfit( )：构建亚组分析生存函数模型
• survdiff( ): 计算log-rank test 的P值
• tbl_survfit( ): 提取亚组分析生存率和死亡人数
  %>%add_nevent()，增加计算死亡人数
• 1-pchisq(p$chisq,length(p$n)-1)，提取p值
• subset=(R15[,i]==1)：R15中亚变量为1的组
  
  
   a<-rbind( ) 把提取的3个数据组合成一行，结合for循环从R15数据第5个（er）变量开始纳入至R15[,i]，每纳入一个变量取它为1的亚变量进入上述3个函数计算OS，死亡人数，p值。每结束一次循环，结果放入a表，指导R15数据随后一个变量
• time、status是自己数据的时间和结局变量名
  
• ==0是结局事件（）死亡
• time=60，计算60月的生存率
• rt，本次研究变量（rt=放疗，本次目的是研究放疗在个亚组的生存时间）
• #
• level = "Yes"：先计算放疗=Yes在各亚组的生存时间，最后几个for循环后面改为“No”

#3.给a表每行命名，即亚变量为1的含义（7个亚变量）
a<-data.frame(cbind(c("Characteristics",
                      " Negative",
                      " Negative",
                      " Negative",
                      " No",
                      " No",
                      " I", 
                      " T1"),a))

#4. 删除X3列中的括号和逗号                  
a <- data.frame(a,str_remove_all(a$X3, "[(,)%]"))

#5. 以空格为条件，将OS（95%CI），分为OS,5%CI和95%CI格式，完成yes1表
yes1<- separate(data=a, col=5,into =c("OS%", "5%CI", "95%CI"), sep = " ");yes1
#col=5：分裂a表第5列为OS%,5%ci,95%CI,分裂条件为空格，注意双引号的空格不要删除

之后的所有的for循环结构一模一样，不再细述

与a表相比，只需建立b表，1变为2，写入=2的变量真是含义，生成yes2表。

若使用自己数据，可以做完yes1表后，直接复制自己修改完的yes1代码，然后,按照紫色和红色修改。

#1、
b<-c("Event N","OS% (5%CI, 95%CI)","P")
#2、
for(i in 5:dim(R15)[2]){ 
  fit1<-survfit(Surv(time,status==0)~rt,subset=(R15[,i]==2),data=R15)
  p<-survdiff(Surv(time,status==0)~rt,rho=0,subset=(R15[,i]==2),data=R15)
  tbl1<-tbl_survfit(fit1,times=60)%>%add_nevent()
  b<-rbind(b,c(
    inline_text(tbl1,column = nevent),
    inline_text(tbl1,time =60, level = "Yes"),
    1-pchisq(p$chisq,length(p$n)-1)))}
#3、7变量为2的真实含义
b<-data.frame(cbind(c("Characteristics",
                      " Positive",
                      " Positive",
                      " Positive",
                      " Yes",
                      " Yes",
                      " II", 
                      " T2"),b))
#4、
b <- data.frame(b,str_remove_all(b$X3, "[(,)%]"))
#5、
yes2<- separate(data=b, col=5,into =c("OS%", "5%CI", "95%CI"), sep = " ");yes2

注意，此时只有两个变量有3分类。所以结果只有两条。

仍然可以按自己数据修改好yes1表的代码后，复制一次，然后修改下面紫色红色处。

#1
c<-c("Event N","OS% (5%CI, 95%CI)","P")
#2、因为之前为2分类，所以for循环从第10（g）个变量开始取
for(i in 10:dim(R15)[2]){ 
  fit1<-survfit(Surv(time,status==0)~rt,subset=(R15[,i]==3),data=R15)
  p<-survdiff(Surv(time,status==0)~rt,rho=0,subset=(R15[,i]==3),data=R15)
  tbl1<-tbl_survfit(fit1,times=60)%>%add_nevent()
  c<-rbind(c,c(
    inline_text(tbl1,column = nevent),
    inline_text(tbl1,time =60, level = "Yes"),
    1-pchisq(p$chisq,length(p$n)-1)))}
#3、只有两个变量
c<-data.frame(cbind(c("Characteristics",
                      " III", 
                      " T3"),c))
#4
c <- data.frame(c,str_remove_all(c$X3, "[(,)%]"))
#5
yes3<- separate(data=c, col =5,into =c("OS%", "5%CI", "95%CI"), sep = " ");yes3

因为只有一条记录，无需for循环，无需亚组。

此步骤建议用下面的代码，然后修改红色处。

d<-c("Event N","OS% (5%CI, 95%CI)","P")

fit1<-survfit(Surv(time,status==0)~rt,data=R15)
tbl1<-tbl_survfit(fit1,times=60)%>%add_nevent()
p<-survdiff(Surv(time,status==0)~rt,rho=0,data=R15)

d<-rbind(d,c(inline_text(tbl1,column = nevent),
         inline_text(tbl1,time =60, level = "Yes"),
         1-pchisq(p$chisq,length(p$n)-1)))
         
d<-data.frame(cbind(c("Characteristics",
                      "All parents"),d))
                      
d<- data.frame(d,str_remove_all(d$X3, "[(,)%]"))

yes4<- separate(data=d, col =5,into =c("OS%", "5%CI", "95%CI"), sep = " ");yes4

res1<-rbind(yes4[c(1,2),],#所以人
            yes1[2,],yes2[2,],#er的亚变量1和2，yes1的第1行，yes2的第1行
            yes1[3,],yes2[3,],#pr
            yes1[4,],yes2[4,],#her2
            yes1[5,],yes2[5,],#surg
            yes1[6,],yes2[6,],#che
            yes1[7,],yes2[7,],yes3[2,],#g
            yes1[8,],yes2[8,],yes3[3,])#t

下面的4个表就更简单了，把上面修改好的4个代码复制一遍。按照说明删除无用代码

例no1（黑色背景），红色是需要修改的。

修改流程见白色背景，与上4个表不同的是，后四个的表无需再提取死亡人数和p值，所以相关的代码删除。

然后，将yes改为no。yes1改为no1。之后的表也是如此，yes变no，表名为no2，no3.。。。

a<-c("OS% (5%CI, 95%CI)")
for(i in 5:dim(R15)[2]){ 
  fit1<-survfit(Surv(time,status==0)~rt,subset=(R15[,i]==1),data=R15)
  tbl1<-tbl_survfit(fit1,times=60)
  a<-rbind(a,inline_text(tbl1,time =60, level = "No"))}

a<-data.frame(cbind(c("Characteristics",
                      " Negative",
                      " Negative",
                      " Negative",
                      " No",
                      " No",
                      " I", 
                      " T1"),a))
a <- data.frame(a,str_remove_all(a$X2, "[(,)%]"))
no1<- separate(data=a, col=3,into =c("OS%", "5%CI", "95%CI"), sep = " ");no1

a<-c("Event N","OS% (5%CI, 95%CI)","P value (Log-rank)")
for(i in 5:dim(R15)[2]){ 
  fit1<-survfit(Surv(time,status==0)~rt,subset=(R15[,i]==1),data=R15)
  p<-survdiff(Surv(time,status==0)~rt,rho=0,subset=(R15[,i]==1),data=R15)
  tbl1<-tbl_survfit(fit1,times=60)%>%add_nevent()
   a<-rbind(a,c(
      inline_text(tbl1,column = nevent),
      inline_text(tbl1,time =60, level = "Yes"),
      1-pchisq(p$chisq,length(p$n)-1)))}

a<-data.frame(cbind(c("Characteristics",
                      " Negative",
                      " Negative",
                      " Negative",
                      " No",
                      " No",
                      " I", 
                      " T1"),a))
a <- data.frame(a,str_remove_all(a$X3, "[(,)%]"))
yes1<- separate(data=a, col=5,into =c("OS%", "5%CI", "95%CI"), sep = " ");yes1

#no2表按照yes2改
b<-c("OS% (5%CI, 95%CI)")

for(i in 5:dim(R15)[2]){ 
  fit1<-survfit(Surv(time,status==0)~rt,subset=(R15[,i]==2),data=R15)
  tbl1<-tbl_survfit(fit1,times=60)
  b<-rbind(b,inline_text(tbl1,time =60, level = "No"))}

b<-data.frame(cbind(c("Characteristics",
                      " Positive",
                      " Positive",
                      " Positive",
                      " Yes",
                      " Yes",
                      " II", 
                      " T2"),b))
                      
b <- data.frame(b,str_remove_all(b$X2, "[(,)%]"))

no2<- separate(data=b, col=3,into =c("OS%", "5%CI", "95%CI"), sep = " ");no2

#no3表按照yes3改
c<-c("OS% (5%CI, 95%CI)")
for(i in 10:dim(R15)[2]){ 
  fit1<-survfit(Surv(time,status==0)~rt,subset=(R15[,i]==3),data=R15)
  tbl1<-tbl_survfit(fit1,times=60)
  c<-rbind(c,inline_text(tbl1,time =60, level = "No"))}

c<-data.frame(cbind(c("Characteristics",
                      "III", 
                      "T3"),c))
c <- data.frame(c,str_remove_all(c$X2, "[(,)%]"))
no3<- separate(data=c, col=3,into =c("OS%", "5%CI", "95%CI"), sep = " ");no3

d<-c("OS% (5%CI, 95%CI)")
fit1<-survfit(Surv(time,status==0)~rt,data=R15)
tbl1<-tbl_survfit(fit1,times=60)
d<-rbind(d,inline_text(tbl1,time =60, level = "No"))
d<-data.frame(cbind(c("Characteristics","No"),d))
d<- data.frame(d,str_remove_all(d$X2, "[(,)%]"))
no4<- separate(data=d, col=3,into =c("OS%", "5%CI","95%CI"), sep= " ");no4

#4表合1
res2<-rbind(no4[c(1,2),],#所有患者生存率
            no1[2,],no2[2,],#er的亚变量
            no1[3,],no2[3,],#pr
            no1[4,],no2[4,],#her2
            no1[5,],no2[5,],#surg
            no1[6,],no2[6,],#che
            no1[7,],no2[7,],no3[2,],#g
            no1[8,],no2[8,],no3[3,])#t

#cbind()数据左右粘贴
res3<-cbind(res1,res2)

• 多了1列亚变量列，删除。列名不规范，删了重写，即删除第1行和第8列

res4<- res3[-1,-8]

•  p值保留3位小数。p值=0.000的改为p<0.001。

data.frame(res4)
res4$X4 <- as.numeric(as.character(res4$X4))#因子变数值
res4$X4 <- round(res4$X4,3)
res4$X4[res4$X4==0]<-"<0.001"

• 提取亚变量人群数目，用tableone包的CreateTableOne（）函数

#提取变量名和患者数目
names(R15)
#[1]"status" "time" "age","er","pr" ,"her2", "surg","che","g","t" 
myVars <- c("er","pr" ,"her2", "surg","che","g","t"     )
catVars<-c("er","pr" ,"her2", "surg","che","g","t" )
table<- print(CreateTableOne(vars=myVars,
                             data=R15,
                             factorVars=catVars),
              catDigits = 2,contDigits = 2,showAllLevels=TRUE)
#去掉第1列，新表命名为table1
table1 <- data.frame(table[,-1])

•  上述提取的亚变量人群总数加入到亚组分析三线表中

res5<-cbind(table1,res4)

• 调整列的顺序，把亚变量名放首位。

res5<-data.frame(res5[c(2,1,3:12)])

• 把变量名写上，第一行写上列名，完成三线表（result表）

每组亚变量上一行插入他们的变量名，
for(i in 1:12) {res5[, i] = as.character(res5[, i])}
result<-rbind(c("Characteristics","Number (%)" ,"Event", "5Y-OS (RT)","P" , 
                NA,NA,NA,"5Y-OS (No_RT)",NA,NA,NA),
              res5[1, ],
              c("ER stage",rep(NA, 11)),
              res5[c(2:3), ],
              c("PR stage",rep(NA, 11)),
              res5[c(4:5),],
              c("HER2",rep(NA, 11)), 
              res5[c(6:7),],
              c("Surgery",rep(NA, 11)),
              res5[c(8:9),],
              c("Chemotherapy",rep(NA, 11)), 
              res5[c(10:11),],
              c("Grade",rep(NA, 11)),
              res5[c(12:14),],
              c("T stage",rep(NA, 11)),
              res5[c(15:17),]);result

for (i in names(result)[c(6:8,10:12)]){result[,i] <- as.numeric(result[,i])}#批量数值转因子
fig<-forestplot(result[,c(1,2,3,4,9,5)], 
                mean=cbind(result[,6],  result[,10]),   
                lower=cbind(result[,7],  result[,11]),  
                upper=cbind(result[,8],  result[,12]),  
                zero=90,            
                boxsize=0.2,  
                graph.pos= 4,
                graphwidth = unit(.2,"npc"),
                xlab="是否放疗在各亚组的5年生存率 (%)",
                xticks=c(75,80,85,90,95) ,
#----------------#字体
                is.summary=c(T,
                              F,
                              T,F,F,
                              T,F,F,
                              T,F,F,
                              T,F,F,
                              T,F,F,
                              T,F,F,F,
                              T,F,F,F),#T=粗体
                txt_gp=fpTxtGp(label=gpar(cex=0.8),
                                ticks=gpar(cex=0.8), 
                                xlab=gpar(cex=1), 
                                title=gpar(cex=2)),
#----------------#线条粗细（x轴、置信区间）
                hrzl_lines=list("1" = gpar(lty=1,lwd=2),
                                "2" = gpar(lty=2),
                                "26"= gpar(lwd=2,lty=1)),
                lwd.zero=1,
                lty.zero =2,
                lwd.ci=1,
                lwd.xaxis=1, 
                lty.ci =1,
                ci.vertices =T,
                ci.vertices.height=0.1, 
                clip=c(75,95),
                grid=structure(c(75,80,85,95), 
                               gp=gpar(col="grey",lty=2,lwd=1)),
#----------------#行间距、字间距/box形状                 
                ineheight=unit(8, 'mm'), 
                line.margin=unit(8, 'mm'),
                colgap=unit(6, 'mm'),
                col=fpColors(zero = "grey",
                              box = c("#f31017", "blue"),
                              lines =c("#f31017", "blue")),
                fn.ci_norm =fpDrawCircleCI,
                legend = c("PMRT", "No_PMRT"),
                title="亚组分析森林图-2 (基于生存率)")

• 关于临床回顾性研究的森林图基本就结束了。
• 大部分代码是自制，虽然有的地方比较麻烦，但是都有解释，亲测可用。之后如有更简洁高效的代码会实时更新！请勿将代码用于商业！
• 下期开始更新临床回顾性研究可视化预测模型：Nomogram（列线图）。首先是通过几篇文献介绍列线图的基本构成，然后是变量选择，模型绘制，验证与比较方法，网页版动态Nomogram（https://impcofmxd.shinyapps.io/IMPC2/）等。
  
  
  感谢阅读，如有错误请指正！