实用技巧——绘制森林图

目录

1. 数据格式整理

（1）单组数据格式

a.cox回归数据提取

b. COX回归数据的精修

（2）多组数据格式

2.绘制单组森林图

3.绘制多组森林图

4.双变量森林图

5.常用的细节调整

a.改变某行颜色

b.改变置信区间CI条的颜色

c.改变文本字体，字号

d.改变某行背景颜色

e.顶部插入文本

f.列名添加横线

g.某行插入文本

---

1. 数据格式整理

（1）单组数据格式

a.cox回归数据提取

以R语言自带的数据集lung为例： 这是一份关于肺癌患者的生存数据。time是生存时间，以天为单位，status是生存状态，1代表删失，2代表死亡。

#加载系统自带的生存分析数据集
library(survival)
str(lung)

 inst：机构代码
time：以天为单位的生存时间
status：删失状态1 = 删失，2 = 出现失效事件
age：岁
sex：性别，男= 1女= 2
ph.ecog：ECOG评分（0 =好，5 =死）
ph.karno：医师进行的Karnofsky评分（0 = 差，100 = 好）
pat.karno：患者自行进行的Karnofsky评分（0 = 差，100 = 好）
meal.cal：用餐时消耗的卡路里
wt.loss：最近六个月的体重减轻

使用cox回归探索危险因素。分类变量需要变为因子型，这样在进行回归时会自动进行哑变量设置。

lung$sex <- factor(lung$sex, labels = c("male","female"))
lung$ph.ecog <- factor(lung$ph.ecog, labels = c("asymptomatic", "symptomatic",'in bed <50%','in bed >50%'))
str(lung)

拟合多因素Cox回归模型，这里我们只用sex/age/ph.ecog/ph.karno4个变量做演示： 

（1）常规的cox分析如下：

fit.cox<-coxph(Surv(time,status)~sex+age+ph.karno+ph.ecog,data=lung)
#查看结果
summary(fit.cox)

 (2) 现状需要对上述数据进行提取形成一个可以绘制森林图的表格(重点)

常见的格式：变量+HR+95%CI_low+95%CI_upper+P value

方法一：

思想：分别提取p value，HR，low CI，high CI

x <- summary(fit.cox)
pvalue <- signif(as.matrix(x$coefficients)[,5],2)#把x$coefficients作为矩阵，提取第五列
HR <- signif(as.matrix(x$coefficients)[,2],2)
low <- signif(x$conf.int[,3],2)
high <- signif(x$conf.int[,4],2)

x$coefficients；x$conf.int 分别是cox回归结果的第一个表和第二个表

 将上述内容整合成为一个数据库格式

result_mul_cox <- data.frame(p.value=pvalue,HR=HR,low=low,high=high,stringsAsFactors = F)
result_mul_cox

 给数据框增加一列为HR（95%CI）

result_mul_cox$'HR(95%CI)' <- paste0(HR," (",low,"-",high,")",sep=" ")
result_mul_cox

方法二（TableOne函数）：

第一步：提取变量HR+95%CI+95%CI

x <- summary(fit.cox)
colnames(x$conf.int)
multi1 <- as.data.frame(round(x$conf.int[,c(1,3,4)],2)) 
#提取出HR，95%CI，95%CI

第二步：提取HR（95%CI）和pvalue

library(tableone)
multi2 <- ShowRegTable(fit.cox,#fit.cox是你的cox回归模型
                       exp=T,  #TRUE表示在表格中显示指数（HR）
                       digits = 2,  #小数点后保留的位数
                       pDigits = 3,  #P值的小数点后显示的位数
                       printToggle = T, #TRUE表示打印结果表格
                       quote = F,   #False表示不对列名添加引号
                       ciFun = confint) #表示使用confint函数计算置信区间

第三步：将两次提取的结果合并

result <- cbind(multi1,multi2)
result

第四步：将行名转为第一列******

#将行名转为表格第一列
library(tibble)
result <- rownames_to_column(result,var="Characteristics")
result

b. COX回归数据的精修

（1）对第一列的数据名称进行修改

从表格中可以看到分类变量的名称是有点问题，由原变量的名称+分类变量名称组合而成，所以需要对原变量名称的字符串进行去除。

#修改变量名称
library(stringr)
result$Characteristics <- str_remove(result$Characteristics,"ph.ecog|sex")
result$Characteristics

（2）对数据的格式进行修改

#将第五行和第六行HR（95%CI）和P值改为字符型
str(result)
result[,5:6] <- lapply(result[,5:6],as.character)
str(result)

#将2，3，4列转化为数值型，便于绘图
result[,2:4] <- lapply(result[,2:4],as.numeric)

（3）创建新的result，补充空行，为了留给分类变量的名称

result <- rbind(
  c("characteristic",NA,NA,NA,"HR (95%CI)","p"),
  result[2,],
  c("Sex",NA,NA,NA,NA,NA),
  c("female",1.00,NA,NA,"1.00 (ref)",NA),
  result[1,],
  result[3,],
  c("ph.ecog",NA,NA,NA,NA,NA),
  c("asymptomatic",1.00,NA,NA,"1.00 (ref)",NA),
  result[4,],
  result[5,],
  result[6,],
  rep(NA,ncol(result)) #插入全为NA的空行
)
result

因为中间三列是数值型，所以不能给第一列标注名称！！！！！

C.总结每个变量的病例数

library(tableone)
str(lung)
myVars <- c("sex","age","ph.karno","ph.ecog") #所有需要统计的数据
catVars <- c("sex","ph.ecog") #所有的分类变量
table1 <- print(CreateTableOne(
  vars=myVars,
  data = lung,
  factorVars = catVars),
  showAllLevels=T) #设置showAlllevels参数可以显示所有的级别

插入基线表table1中的空行，使行数与变量result的行数一致

N <- rbind(c(NA,NA),
            table1[4,],
           c(NA,NA),
           table1[2:3,],
           table1[5,],
           c(NA,NA),
           table1[6:9,],
           c(NA,NA));N
N <- N[,-1]
N <- as.data.frame(N)#注意对应相应的变量
head(N)

将result和N合并

result1 <- cbind(result[,1],N,result[,2:6])
#设置第一行名
result1[1,] <- c("characteristic","number (%)",NA,NA,NA,"HR (95%CI)","p") #因为中间三行为数值型不能命名

---

（2）多组数据格式

格式要求：

方法一：excel直接录入

方法二：使用上面的单组数据的方法，在列其他的HR

2.绘制单组森林图

（1）forestplot包

以前面的格式部分产生的result1为例，可以没有结果，可以把前面代码直接复制到R中跑一遍，生成result1

result1
library(forestplot)
fig <- forestplot(
  result1[,c(1,2,6,7)],#需要在森林图中展示的列
  mean=result1[,3],    #均值列（HR）
  lower=result1[,4],  #95%置信区间的下限数据列
  upper=result1[,5],  #95%置信区间的上限数据列
  zero=1,              #均值HR为1为参考线
  boxsize=0.4,        #方框的大小
  graph.pos="right",   #森林图在右侧
  hrzl_lines=list(      #水平线样式设置（三线表的线）
    "1"=gpar(lty=1,lwd=2),  #均值线
    "2"=gpar(lty=2),        #下限和上限之间的虚线
    "13"=gpar(lwd=2,lty=1,columns=c(1:4)) #下限和上限线
  ),
  graphwidth=unit(.25,"npc"),   #设置森林图宽度，整个图形的宽度
  xlab="HR",       #X轴标签
  xticks=c(-0.8,1,3,6,10),   #x轴刻度
  is.summary=c(T,F,T,F,F,F,T,F,F,F,F), #判断是否为汇总行，不参与绘图的行设置为T，是否整行加粗
  txt_gp=fpTxtGp(        #设置文本样式
    label = gpar(cex=0.8),   #设置标签大小（正文内容字体大小）
    ticks = gpar(cex=1),     #设置刻度大小
    xlab = gpar(cex=0.9),    #x轴标签大小
    title = gpar(cex=1.2)    #标题大小
  ),
  lwd.zero=1,                #HR=1时线条宽度，零刻度线得宽度
  lwd.ci=2,                #置信区间线得宽度
  lwd.xaxis=2,               #x轴得宽度
  lty.ci=1.5,                #置信区间得样式
  ci.vertices=T,             #是否显示置信区间得顶点
  ci.vertices.height=0.2,    #置信区间顶点的高度
  clip=c(0.1,8),             #X轴的截断范围
  ineheight=unit(8,"mm"),    #森林图内线的高度
  line.margin=unit(8,"mm"),  #森林图线的边距
  colgap=unit(6,"mm"),       #森林图中方框的间距
  fn.ci_norm="fpDrawDiamondCI",#置信区间的形状
  title="多因素COX回归森林图",#森林图的标题
  col=fpColors(               #设置颜色
    box = "blue4",           #方框的颜色
    lines = "blue4",         #线条颜色
    zero = "black"           #均值为0的颜色
  )
)
fig

（2）forestploter包（重点）

a.格式要求

1.基本数据格式

重点：forestploter包格式不能包括第一行的变量名称

 对result1进行修改，将第一行设置为列名

library(forestploter)
result2 <- result1 
colnames(result2) <- result2[1,] #将第一行赋值为列名
#colnames(result2) <- c("characteristic","number","HR","low","high","HR(95%CI)","p"  )
result2 <- result2[-1,]
result2 <- result2[-nrow(result2),]

 标准格式如下：

 前面可以通过加粗的方式修改分类变量的名称，这里可以通过分变量进行缩进体现

result2$characteristic <- ifelse(is.na(result2$number),
                                 result2$characteristic,
                                 paste0("    ",result2$characteristic))

 2.把NA转为空字符串

result2$number <- ifelse(is.na(result2$number)," ",result2$number)
result2$`HR(95%CI)` <- ifelse(is.na(result2$`HR(95%CI)`)," ",result2$`HR(95%CI)`)
result2$p <- ifelse(is.na(result2$p)," ",result2$p)

HR，low，high三行的NA不能随意修改为字符串，不然数据格式就被改变不是数值型，不能计算绘图。

3.计算标准误差SE，绘图时表现为正方形的大小

result2[,3:5] <- lapply(result2[,3:5],as.numeric)
str(result2)
result2$se <- -(log(result2$high)-log(result2$low))/1.96

4.为森林图添加空白列，为了产生一个绘图区间

result2$" " <- paste(rep(" ",12),collapse = " ")

b.设置主题

tm <- forest_theme(base_size = 10,
                   refline_col = "red4",
                   arrow_type = "closed",
                   footnote_col = "blue4")

主题参数介绍：

所有主题参数都通过上述forest_theme( )函数调整

base_size=10                  设置文本的基础大小

设置可信区间的外观：

ci_pch= 15                         设置可信区间点的形状

ci_col="blue4"                    设置可信区间的边框颜色

ci_fill="blue4"                     设置可信区间的填充颜色

ci_alpha=0.8                      设置可信区间的透明度

ci_lty=1,                              设置可信区间的线条类型

ci_lwd=1.5                          设置可信区间的线条宽度

ci_Theight=0.2                  设置T字在可信区间末端的高度，默认NULL

设置x轴的粗细：

xaxis_lwd=0.6                    设置x轴的粗细

xaxis_cex=1                       设置x轴文本大小

设置参考线的外观：

refline_lwd=1                      设置宽度

refline_lty="dashed"            设置线的类型

refline_col="grey20"            设置线的颜色

设置垂直线的外观：

vertline_lwd=1

vertline_lty="dashed"

vertline_col="grey20"

设置脚注的字体大小，字体样式和颜色：

footnote_cex=0.6                设置脚注字体大小

footnote_fontface="italic"    设置字体样式

footnote_col="red4"             设置脚注文本颜色

设置箭头格式：

arrow_type="closed"            设置箭头类型为闭合箭头

arrow_label_just="end"         设置箭头标签对齐方式为末尾

设置图例信息：

legend_name= "group"           设置图例标题

legend_value=c("trt1","trt2")    设置图例值

c.绘制森林图

p <- forest(result2[,c(1,2,9,6,7)],         #选择要在森林图中显示的数据列，包括变量名称，样本数，绘图的空白列，HR（95%ci），p值
            est = result2$HR,
            lower = result2$low,
            upper = result2$high,
            sizes = result2$se,
            ci_column = 3,
            ref_line = 1,
            arrow_lab = c("Low risk","High risk"),
            xlim = c(-1,10),
            ticks_at = c(-0.5,1,3,6,10),
            theme = tm,
            footnote = "XXXXXXX")
p

3.绘制多组森林图

以格式中提到的格式为例：

 1.读取数据，设置缩进，设置绘图的空白区域

library(readxl)
library(forestploter)
slt <- read_excel("GRS森林图数据.xlsx",sheet=1)
str(slt)
slt$Var <- ifelse(is.na(slt$HR1),slt$Var,paste0("    ",slt$Var))
slt$HR1 <- ifelse(is.na(slt$HR1)," ",slt$HR1)
slt$HR2 <- ifelse(is.na(slt$HR2)," ",slt$HR2)
slt$HR3 <- ifelse(is.na(slt$HR3)," ",slt$HR3)
slt$plot1 <- paste0(rep(" ",20),collapse = " ")
slt$plot2 <- paste0(rep(" ",20),collapse = " ")
slt$plot3 <- paste0(rep(" ",20),collapse = " ")

 设置主题

tm <- forest_theme(
  base_size = 10,
  ci_pch = 18,
  ci_col = "black",
  ci_lty = 1,
  ci_lwd =1.8,
  ci_Theight = 0.2,
  refline_lwd = 1,
  refline_lty = "dashed",
  refline_col = "grey20",
  footnote_cex = 0.8,
  footnote_fontface = "italic",
  footnote_col = "red4"
)

 绘图函数，用list列表来描述每次需要描绘的结局：

p <- forest(slt[,c(1,2,17,7,18,12,19)],
           est=list(slt$HR1_mean,
                    slt$HR2_mean,
                    slt$HR3_mean),
           lower=list(slt$HR1_low,
                      slt$HR2_low,
                      slt$HR3_low),
           upper=list(slt$HR1_high,
                      slt$HR2_high,
                      slt$HR3_high),
           sizes = 0.6,
           ci_column=c(3,5,7),
           ref_line=c(1,1,1),
       arrow_lab = c("Low risk","High risk"),
       xlim = list(c(0.6,2),
                   c(0.6,2.2),
                   c(0.6,1.7)),
       ticks_at = list(c(0.8,1,1.3,1.6),
                       c(0.8,1,1.3,1.6,2),
                       c(0.8,1,1.2,1.4)),
       xlab = "HR",
       theme=tm,
       footnote = "
       
       SSBs, sugar-sweetened beverages; 
       ASBs, artificial sweetened beverages; 
       PJs, pure fruit/vegetable juices;
       HR, hazard ratio; CI, confidence interval."
       )
p

 对细节进行修饰，参考第五节部分内容细节调整

PP <- edit_plot(p,
                part = "header",
                gp=gpar(fontsize=12,
                        fontface="bold")) #标题栏放大
PP <- edit_plot(PP,
                row=c(1,6,11),
                gp=gpar(fontface="bold")) #变量名称加粗
PP <- add_border(PP,part = "header",where = "bottom") #增加边际线
PP <- add_border(PP,part = "header",where = "top")
PP <- edit_plot(PP,
                row = c(5,10),
                col=2,
                gp=gpar(fontface="bold"))#统计学意义的HR加粗
PP <- edit_plot(PP,
                row = c(5,10),
                col=3,
                which = "ci",
                gp=gpar(col="green4",
                        lwd=2))
PP <- edit_plot(PP,
                row = c(4,5,9,10),
                col=4,
                gp=gpar(fontface="bold"))
PP <- edit_plot(PP,
                row = c(4,5,9,10),
                col=5,
                which = "ci",
                gp=gpar(col="green4",
                        lwd=2))
PP <- edit_plot(PP,
                row = c(3,8,13,15),
                col=6,
                gp=gpar(fontface="bold"))
PP <- edit_plot(PP,
                row = c(3,8,13,15),
                col=7,
                which = "ci",
                gp=gpar(col="green4",
                        lwd=2))
PP <- edit_plot(PP,
                which = "background",
                gp=gpar(fill="white"));PP#设置背景色

4.双变量森林图

5.常用的细节调整

以下内容只能修改forestplot包产生的森林图：

基本介绍：

使用edit_plot函数进行修改，基本修改格式为

edit_plot（p，

                 row=，

                 which=，

                  gp=gpar（））

p表示需要编辑的图片名称

row需要修改的行

which需要修改的内容：包括text文本，background（背景），ci(置信区间)三个内容

gpar（）中可以修改颜色col=，字体fontface=，填充颜色fill=

a.改变某行颜色

PP <- edit_plot(p,row=5,gp=gpar(col="red4",fontface="italic"))
PP

b.改变置信区间CI条的颜色

将7到10行的第三列置信区间改为绿色

PP <- edit_plot(PP,row=c(7:10),col=3,
                which="ci",
                gp=gpar(col="green4"));PP

c.改变文本字体，字号

PP <- edit_plot(PP,
                row=c(2,6),
                gp=gpar(fontface="bold"));PP

d.改变某行背景颜色

PP <- edit_plot(PP,
                row=5,
                which = "background",
                gp=gpar(fill="yellow"));PP

e.顶部插入文本（森林图标题）

col设置标题在位置，在哪几列之间

PP <- insert_text(PP,
                text="COX回归分析森林图",
                col=1:5,
                part="header",
                gp=gpar(fontface="bold"));PP

 

f.列名添加横线（三线表的横线）

PP <- add_border(PP,
                 part = "header",
                 where="bottom");PP

PP <- add_border(PP,
                 row = 10,
                 where="bottom");PP

 

g.某行插入文本

just表示文本的位置，可以设置为center，left，right

PP <- insert_text(PP,
                  text="Smoking",
                  row=5,
                  just = "left",
                  gp=gpar(col="blue4",
                          cex=1,
                          fontface="italic"));PP
PP <- insert_text(PP,
                  text="80(58.56)",
                  row=5,
                  col=2,
                  just = "center",
                  gp=gpar(col="blue4",
                          cex=1,
                          fontface="italic"));PP

h.添加文本

PP <- add_text(p,"HR (95%CI)",col=2,part = "header",just = "left");PP

I.修改列名****（重要）

因为列名重复造成的命名错误

plot <- edit_plot(
  plot,
  row = 1,
  col = c(3,5), #修改第一行，第三列和第五列
  part = "header",#列名，如果修改其他部分part设置为body
  which="text",#设置为修改文本
  label="Placebo\r\n(N=2841)"#填写修改的内容 \r\n表示换行
)