R语言使用forestploter包绘制单组及双组森林图（2）（附有超详细备注）

森林图在论文中很常见，多用于表示多因素分析中的变量与结果变量的比值效应，可以用图示的方法比较直观的绘制出来。既往我们在文章《R语言快速绘制多因素回归分析森林图（1）》已经介绍了怎么绘制森林图，但是绘图比较简单，不够美观，不能绘制相对复杂的森林图。今天我们来介绍一下forestploter包，它等于是在forestplot包的基础上进一步强化功能，制作方法也相对简单一点，而且加强了对图形的精细控制，而且可以绘制单组和多组森林图。

我们先导入R包和数据

library(grid)
library(forestploter)
dt<-read.csv("E:/r/test/forest2.csv",sep=',',header=TRUE)
1
2
3

这是一个森林图数据（公众号回复：森林图数据2，可以获得数据）这个数据，嫌麻烦的可以在这里下载：https://download.csdn.net/download/dege857/86945654?spm=1001.2014.3001.5501
前几个变量我解释一下Subgroup：分组，就是的组别和亚组，Treatment：治疗组例数，Placebo：使用安慰剂，也就是对照组的例数。Est：效应值，可以是OR或者HR，low：效应值的最低值，可以认为是可信区间的下限，hi：效应值的最高值，可以认为是可信区间的上限。
我们先来画一个最基础的森林图，后面的几个数据暂时不需要，我们先精简一下数据

dt <- dt[,1:6]
View(dt)
1
2

数据变成了一个精简的数据（上图），我们先把数据格式整理一下，先让它缩进一格

dt$Subgroup <- ifelse(is.na(dt$Placebo), 
                      dt$Subgroup,
                      paste0("   ", dt$Subgroup))
1
2
3

接下来我们把治疗组和对照组NA（有缺失）的地方变成一个空格

dt$Treatment <- ifelse(is.na(dt$Treatment), "", dt$Treatment)
dt$Placebo <- ifelse(is.na(dt$Placebo), "", dt$Placebo)
1
2

生成一个变量se，它在绘图的时候表示正方形的大小

dt$se <- (log(dt$hi) - log(dt$est))/1.96
1

生成一个绘图区间，等下用来绘图

dt$` ` <- paste(rep(" ", 20), collapse = " ")
1

生成HR和可信区间

dt$`HR (95% CI)` <- ifelse(is.na(dt$se), "",
                           sprintf("%.2f (%.2f to %.2f)",
                                   dt$est, dt$low, dt$hi))#sprintF返回字符和可变量组合
1
2
3

最后得到绘图数据

整理好数据我们就可以绘图了

p <- forest(dt[,c(1:3, 8:9)],
            est = dt$est,       #效应值
            lower = dt$low,     #可信区间下限
            upper = dt$hi,      #可信区间上限
            sizes = dt$se,     #黑框的大小
            ci_column = 4,   #在那一列画森林图，要选空的那一列
            ref_line = 1,
            arrow_lab = c("Placebo Better", "Treatment Better"),
            xlim = c(0, 4),
            ticks_at = c(0.5, 1, 2, 3),
            footnote = "This is the demo data. Please feel free to change\nanything you want.")
p
1
2
3
4
5
6
7
8
9
10
11
12

有些人可能想问，有些森林图是有P值的，这个没有P值怎么办，我们可以给它加上去

dt$p <- paste(rep("<0.05", 22))
1

重新绘图

p <- forest(dt[,c(1:3, 8:10)],
            est = dt$est,       #效应值
            lower = dt$low,     #可信区间下限
            upper = dt$hi,      #可信区间上限
            sizes = dt$se,     #黑框的大小
            ci_column = 4,   #在那一列画森林图，要选空的那一列
            ref_line = 1,
            arrow_lab = c("Placebo Better", "Treatment Better"),
            xlim = c(0, 4),
            ticks_at = c(0.5, 1, 2, 3),
            footnote = "This is the demo data. Please feel free to change\nanything you want.")
p
1
2
3
4
5
6
7
8
9
10
11
12

这样P值就出来了

生成图形后我们还可以对图形进行细节调整，我们把第一行总病例数放到最后改成汇总，然后把它的名字改成Overall，

dt_tmp <- rbind(dt[-1, ], dt[1, ])
dt_tmp[nrow(dt_tmp), 1] <- "Overall"
1
2

forest_theme函数可以对森林图图片的细节进行调整，我们可以在forest_theme函数上设定预定的森林图样式模块，到时直接绘图就好了

tm <- forest_theme(base_size = 10,  #文本的大小
                   # Confidence interval point shape, line type/color/width
                   ci_pch = 15,   #可信区间点的形状
                   ci_col = "#762a83",    #CI的颜色
                   ci_fill = "blue",     #ci颜色填充
                   ci_alpha = 0.8,        #ci透明度
                   ci_lty = 1,            #CI的线型
                   ci_lwd = 1.5,          #CI的线宽
                   ci_Theight = 0.2, # Set an T end at the end of CI  ci的高度，默认是NULL
                   # Reference line width/type/color   参考线默认的参数，中间的竖的虚线
                   refline_lwd = 1,       #中间的竖的虚线
                   refline_lty = "dashed",
                   refline_col = "grey20",
                   # Vertical line width/type/color  垂直线宽/类型/颜色   可以添加一条额外的垂直线，如果没有就不显示
                   vertline_lwd = 1,              #可以添加一条额外的垂直线，如果没有就不显示
                   vertline_lty = "dashed",
                   vertline_col = "grey20",
                   # Change summary color for filling and borders   更改填充和边框的摘要颜色
                   summary_fill = "yellow",       #汇总部分大菱形的颜色
                   summary_col = "#4575b4",
                   # Footnote font size/face/color  脚注字体大小/字体/颜色
                   footnote_cex = 0.6,
                   footnote_fontface = "italic",
                   footnote_col = "red")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

设定好模块后直接绘图就可以了

pt <- forest(dt_tmp[,c(1:3, 8:9)],
             est = dt_tmp$est,
             lower = dt_tmp$low, 
             upper = dt_tmp$hi,
             sizes = dt_tmp$se,
             is_summary = c(rep(FALSE, nrow(dt_tmp)-1), TRUE),
             ci_column = 4,
             ref_line = 1,
             arrow_lab = c("Placebo Better", "Treatment Better"),
             xlim = c(0, 4),
             ticks_at = c(0.5, 1, 2, 3),
             footnote = "This is the demo data. Please feel free to change\nanything you want.",
             theme = tm)

plot(pt)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

我们还可以对图片细节进行修改，比如我们想把第三行变成红色

g <- edit_plot(p, row = 3, gp = gpar(col = "red", fontface = "italic"))
g
1
2

更改3,6,10,11行方块及长条的颜色为绿色

g <- edit_plot(g,
               row = c(3, 6, 11, 13),
               col = 4,
               which = "ci",
               gp = gpar(col = "green"))
g
1
2
3
4
5
6

把2, 5, 10, 13, 17, 20行的文本变成粗体

g <- edit_plot(g,
               row = c(2, 5, 10, 13, 17, 20),
               gp = gpar(fontface = "bold"))
g
1
2
3
4

把第五行的背景改成绿色

g <- edit_plot(g, row = 5, which = "background",
               gp = gpar(fill = "darkolivegreen1"))
g
1
2
3

在顶部插入文本

g <- insert_text(g,
                 text = "Treatment group",
                 col = 2:3,
                 part = "header",
                 gp = gpar(fontface = "bold"))
g
1
2
3
4
5
6

在标题下方加一条下划线

g <- add_underline(g, part = "header")
g
1
2

在第十行位置插入文本

g <- insert_text(g,
                 text = "This is a long text. Age and gender summarised above.\nBMI is next",
                 row = 10,
                 just = "left",
                 gp = gpar(cex = 0.6, col = "green", fontface = "italic"))

g
1
2
3
4
5
6
7

下面来介绍一下怎么绘制多组的森林图，绘制对组森林图就要涉及多组的数据，我们把数据重新导入一下

dt<-read.csv("E:/r/test/forest2.csv",sep=',',header=TRUE)
1

设置一下缩进，这样好看一点，这步和前面是一样的

dt$Subgroup <- ifelse(is.na(dt$Placebo), 
                      dt$Subgroup,
                      paste0("   ", dt$Subgroup))#######如果变量没有缺失，就缩进一格，也就是前进一格
1
2
3

因为是双组变量，所以要设置2个n，这步和前面基本一样

dt$n1 <- ifelse(is.na(dt$Treatment), "", dt$Treatment)###将缺失的部分变为空格
dt$n2 <- ifelse(is.na(dt$Placebo), "", dt$Placebo)
1
2

因为是要画两个森林图，所以要增加两个空地方来画图

dt$`CVD outcome` <- paste(rep(" ", 20), collapse = " ")
dt$`COPD outcome` <- paste(rep(" ", 20), collapse = " ")
1
2

设置一些森林图的基本参数，这步和前面是一样的

tm <- forest_theme(base_size = 10,
                   refline_lty = "solid",   #参考线类型
                   ci_pch = c(15, 18),
                   ci_col = c("#377eb8", "#4daf4a"),
                   footnote_col = "blue",
                   legend_name = "Group",   #设置标题名字
                   legend_value = c("Trt 1", "Trt 2"),   #设置分组名字
                   vertline_lty = c("dashed", "dotted"),
                   vertline_col = c("#d6604d", "#bababa"))
1
2
3
4
5
6
7
8
9

最后绘图，这里我要说一下，ci_column = c(3, 5)是指在第3和5列绘图，est_gp1和est_gp2为一组，est_gp3和est_gp4为一组，其他的依次类推

p <- forest(dt[,c(1, 19, 21, 20, 22)],
            est = list(dt$est_gp1,
                       dt$est_gp2,
                       dt$est_gp3,
                       dt$est_gp4),
            lower = list(dt$low_gp1,
                         dt$low_gp2,
                         dt$low_gp3,
                         dt$low_gp4), 
            upper = list(dt$hi_gp1,
                         dt$hi_gp2,
                         dt$hi_gp3,
                         dt$hi_gp4),
            ci_column = c(3, 5),
            ref_line = 1,
            vert_line = c(0.5, 2),
            nudge_y = 0.2,
            theme = tm)

p
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

最后双组森林图也绘制好了，也可以像单组一样调整。这应该是我写得最详细的1次了，整整搞了15页。