R语言| 亚组分析森林图（基于HR结果）
欢迎大家关注我的公众号：一只勤奋的科研喵
亚组分析森林图
本期介绍使用R语言制作亚组分析三线表及其森林图。
注：由于未能找到专门做亚组分析的包或相关代码，本期代码是小编根据以往所学而写的，算是抛砖引玉了，之后有更好更简洁的代码时会及时更新本文。
更多代码解释及示例数据请阅读原文：R语言|14. 森林图-3: 亚组分析三线表及森林图
载入R包和数据
#1.载入包
library(survival)
library(plyr)
##制作table1
library(tableone)
##制作森林图
library(forestplot)
#2.清理工作环境
rm(list = ls())
#3.数据放入工作路径，加载
R14<- read.csv('R14.csv')
#4.查看数据前6行
head(R14)
#5.查看数据数据性质
str(R14)
#6.查看结局:0为局部区域复发（55人）
R14$status<-factor(R14$status)
summary(R14$status)
3.png
- 相关变量定义：
RT：放疗
T：1=T1，2=T2；
ER: 1=阴性，2=阳性；
PR: 1=阴性，2=阳性；
KI67: 1=低表达（<=30%）,2=高表达；
LVI(脉管癌栓)：1=无，2=有；
Grade: 1=I级，2=2级，3=级。
- 本文数据的特殊处理：
1、分类变量要全变为数字，并统一变为1.2.3.格式（稍后解释）；
2、变量2分类在前，然后是3分类，...（稍后解释）。一定要这么处理！！！
4.png
所有需要修改的代码均用红色标出（详见原文链接,简书无法正常显示）
一、亚组分析：HR、95%CI及P值批量获取
1、所有亚变量都为“1”的亚组分析
#1-建一个名为"HR","5%CI","95%CI","P"的4列表：a表。
#用于放批量获取的结果。
a<-c("HR","5%CI","95%CI","P")
##2建两个for循环，第一个是批量Cox亚组分析，
##第二个是将结果批量提取放入a表
for(i in 5:dim(R14)[2]){
cox1<-coxph(Surv(time,status==0)~RT,subset=(R14[,i]=='1'),data=R14)
aa<-summary(cox1)
for(j in 1:dim(aa$coefficients)[1]){
a<-rbind(a,c(round(aa$coefficients[j,2],2),
round(aa$conf.int[j,3],2),
round(aa$conf.int[j,4],2),
round(aa$coefficients[j,5],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")
5.png
#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
6.png
经过以上4步，a表就做完了。其他表是同理。接下来做所有亚变量为2的亚组分析。
2、所有亚变量都为“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
7.png
b表：所有亚变量为2的亚组分析完成。虽然与a表的代码是一样的，但命名要分为b和bb，不然会把数据做到a表上，导致错误。
3、所有亚变量都为“3”的亚组分析
注意，此时只剩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
8.png
4、在所有患者中做RT的单因素分析
#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"
9.png
二、放疗患者在各亚组的数目
这里代码意义参考：R语言|4. 轻松绘制临床基线表Table 1
#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)
10.png
此时，我们提取到了在放疗或未放疗两组中的患者人数，是一个17行5列的表。但我们不要p值和text列，level列要从数字变为英文格式。
#4- 去掉第4.5列（上图p和test列），新表命名为table1（17行3列）
table1<-table[,c(-4,-5)]
11.png
三、4个亚组分析表合成一个表并与table1表合并
12.png
#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列表）
13.png
#2- res1表列名进入表格内，并命名为Characteristics
>此时的res1表变为17行6列表（多了亚变量名列）
res1<-tibble::rownames_to_column(res1, var = "Characteristics")
14.png
#3- table与res1合成为res2
res2<-cbind(table1,res1)
前面介绍rbind()为上下的行粘贴，cbind()函数为左右的列粘贴。此时我们把17行3列的table1和17行6列的res2组合为res2（17行9列）
15.png
