【原】逻辑回归亚组分析森林图绘制

💡专注R语言在🩺生物医学中的使用

---

设为“星标”，精彩不错过

---

我之前以COX回归为例演示了亚组分析和森林图的绘制：

• R语言亚组分析及森林图绘制
• R语言亚组分析1行代码实现！

但是并没有演示逻辑回归亚组分析和森林图的绘制，因为道理都是相通的，而且逻辑回归比COX回归还要简单，但是没想到有许多朋友完全不会逻辑回归的亚组分析和森林图绘制。。。

本着帮人帮到底，送佛送到西的原则，我还是再写一篇逻辑回归亚组分析森林图绘制的推文吧。。

因为有非常好用的jstable可以一行代码实现，所以我就不演示手动实现的过程了。

安装

install.packages("jstable")

## From github: latest version
remotes::install_github('jinseob2kim/jstable')
library(jstable)

准备数据

还是使用上次演示的数据。

使用survival包中的colon数据集用于演示，这是一份关于结肠癌患者的生存数据，共有1858行，16列，共分为3个组，1个观察组+2个治疗组，观察他们发生终点事件的差异。

各变量的解释如下：

• id：患者id
• study：没啥用，所有患者都是1
• rx：治疗方法，共3种，Obs(观察组), Lev(左旋咪唑), Lev+5FU(左旋咪唑+5-FU)
• sex：性别，1是男性 age：年龄
• obstruct：肠梗阻，1是有 perfor：肠穿孔，1是有
• adhere：和附近器官粘连，1是有
• nodes：转移的淋巴结数量
• status：生存状态，0代表删失，1代表发生终点事件
• differ：肿瘤分化程度，1-well,2-moderate,3-poor
• extent：局部扩散情况，1-submucosa，2-muscle，3-serosa，4-contiguous-structures
• surg：手术后多久了，1-long,2-short
• node4：是否有超过4个阳性淋巴结，1代表是
• time：生存时间
• etype：终点事件类型，1-复发，2-死亡

rm(list = ls())
library(survival)

str(colon)

'data.frame':   1858 obs. of  16 variables:
 $ id      : num  1 1 2 2 3 3 4 4 5 5 ...
 $ study   : num  1 1 1 1 1 1 1 1 1 1 ...
 $ rx      : Factor w/ 3 levels "Obs","Lev","Lev+5FU": 3 3 3 3 1 1 3 3 1 1 ...
 $ sex     : num  1 1 1 1 0 0 0 0 1 1 ...
 $ age     : num  43 43 63 63 71 71 66 66 69 69 ...
 $ obstruct: num  0 0 0 0 0 0 1 1 0 0 ...
 $ perfor  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ adhere  : num  0 0 0 0 1 1 0 0 0 0 ...
 $ nodes   : num  5 5 1 1 7 7 6 6 22 22 ...
 $ status  : num  1 1 0 0 1 1 1 1 1 1 ...
 $ differ  : num  2 2 2 2 2 2 2 2 2 2 ...
 $ extent  : num  3 3 3 3 2 2 3 3 3 3 ...
 $ surg    : num  0 0 0 0 0 0 1 1 1 1 ...
 $ node4   : num  1 1 0 0 1 1 1 1 1 1 ...
 $ time    : num  1521 968 3087 3087 963 ...
 $ etype   : num  2 1 2 1 2 1 2 1 2 1 ...

分类变量需要变为因子型，这样在进行回归时会自动进行哑变量设置。

为了演示，我们只选择Obs组和Lev+5FU组的患者，所有的分类变量都变为factor，把年龄也变为分类变量并变成factor。

而且由于要要是逻辑回归，所以我们不用其中的生存时间这个变量。

suppressMessages(library(tidyverse))

df <- colon %>% 
  mutate(rx=as.numeric(rx)) %>% 
  filter(etype == 1, !rx == 2) %>%  #rx %in% c("Obs","Lev+5FU"), 
  select(time, status,rx, sex, age,obstruct,perfor,adhere,differ,extent,surg,node4) %>% 
  mutate(sex=factor(sex, levels=c(0,1),labels=c("female","male")),
         age=ifelse(age >65,">65","<=65"),
         age=factor(age, levels=c(">65","<=65")),
         obstruct=factor(obstruct, levels=c(0,1),labels=c("No","Yes")),
         perfor=factor(perfor, levels=c(0,1),labels=c("No","Yes")),
         adhere=factor(adhere, levels=c(0,1),labels=c("No","Yes")),
         differ=factor(differ, levels=c(1,2,3),labels=c("well","moderate","poor")),
         extent=factor(extent, levels=c(1,2,3,4),
                       labels=c("submucosa","muscle","serosa","contiguous")),
         surg=factor(surg, levels=c(0,1),labels=c("short","long")),
         node4=factor(node4, levels=c(0,1),labels=c("No","Yes")),
         rx=ifelse(rx==3,0,1),
         rx=factor(rx,levels=c(0,1))
         )

str(df)

'data.frame':   619 obs. of  12 variables:
 $ time    : num  968 3087 542 245 523 ...
 $ status  : num  1 0 1 1 1 1 0 0 0 1 ...
 $ rx      : Factor w/ 2 levels "0","1": 1 1 2 1 2 1 2 1 1 2 ...
 $ sex     : Factor w/ 2 levels "female","male": 2 2 1 1 2 1 2 1 2 2 ...
 $ age     : Factor w/ 2 levels ">65","<=65": 2 2 1 1 1 2 2 1 2 2 ...
 $ obstruct: Factor w/ 2 levels "No","Yes": 1 1 1 2 1 1 1 1 1 1 ...
 $ perfor  : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
 $ adhere  : Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 1 1 1 ...
 $ differ  : Factor w/ 3 levels "well","moderate",..: 2 2 2 2 2 2 2 2 3 2 ...
 $ extent  : Factor w/ 4 levels "submucosa","muscle",..: 3 3 2 3 3 3 3 3 3 3 ...
 $ surg    : Factor w/ 2 levels "short","long": 1 1 1 2 2 1 1 2 2 1 ...
 $ node4   : Factor w/ 2 levels "No","Yes": 2 1 2 2 2 2 1 1 1 1 ...

亚组分析

使用jstable，只要1行代码即可！！！

library(jstable)

res <- TableSubgroupMultiGLM(
  
  # 指定公式，不要乱写！
  formula = status ~ rx, 
  
  # 指定哪些变量有亚组
  var_subgroups = c("sex","age","obstruct","perfor","adhere",
                    "differ","extent","surg","node4"), 
  data = df #指定你的数据
  )

res

               Variable Count Percent   OR Lower                 Upper P value
rx1             Overall   619     100 1.99  1.45                  2.75  <0.001
1...1               sex  <NA>    <NA> <NA>  <NA>                  <NA>    <NA>
init...2         female   312    50.4 1.43  0.92                  2.24   0.114
...3               male   307    49.6 2.93  1.84                  4.71  <0.001
1...4               age  <NA>    <NA> <NA>  <NA>                  <NA>    <NA>
init...5            >65   224    36.2 2.48  1.45                  4.28   0.001
...6               <=65   395    63.8 1.76  1.18                  2.63   0.005
1...7          obstruct  <NA>    <NA> <NA>  <NA>                  <NA>    <NA>
init...8             No   502    81.1 1.94  1.36                  2.77  <0.001
...9                Yes   117    18.9 2.24  1.07                  4.75   0.034
1...10           perfor  <NA>    <NA> <NA>  <NA>                  <NA>    <NA>
init...11            No   602    97.3 1.94   1.4                  2.69  <0.001
...12               Yes    17     2.7 5.83  0.78                 61.72   0.104
1...13           adhere  <NA>    <NA> <NA>  <NA>                  <NA>    <NA>
init...14            No   533    86.1 2.03  1.44                  2.88  <0.001
...15               Yes    86    13.9  1.7  0.72                  4.05   0.229
...16            differ  <NA>    <NA> <NA>  <NA>                  <NA>    <NA>
...17              well    56     9.2 3.78  1.28                 11.92   0.019
...18          moderate   444    73.3 1.95  1.34                  2.86   0.001
...19              poor   106    17.5 1.61  0.74                  3.54   0.228
...20            extent  <NA>    <NA> <NA>  <NA>                  <NA>    <NA>
...21         submucosa    18     2.9    0  <NA> 1.23555810375665e+213   0.996
...22            muscle    70    11.3 2.83  0.97                  9.05   0.064
...23            serosa   500    80.8 2.07  1.45                  2.96  <0.001
...24        contiguous    31       5 1.94  0.42                  9.23   0.392
1...25             surg  <NA>    <NA> <NA>  <NA>                  <NA>    <NA>
init...26         short   452      73 2.17  1.49                  3.17  <0.001
...27              long   167      27 1.54  0.83                  2.86   0.169
1...28            node4  <NA>    <NA> <NA>  <NA>                  <NA>    <NA>
init...29            No   453    73.2 2.21  1.51                  3.26  <0.001
...30               Yes   166    26.8 1.61  0.84                  3.11   0.155
          P for interaction
rx1                    <NA>
1...1                 0.031
init...2               <NA>
...3                   <NA>
1...4                 0.317
init...5               <NA>
...6                   <NA>
1...7                 0.734
init...8               <NA>
...9                   <NA>
1...10                0.316
init...11              <NA>
...12                  <NA>
1...13                  0.7
init...14              <NA>
...15                  <NA>
...16                 0.447
...17                  <NA>
...18                  <NA>
...19                  <NA>
...20                 0.087
...21                  <NA>
...22                  <NA>
...23                  <NA>
...24                  <NA>
1...25                0.351
init...26              <NA>
...27                  <NA>
1...28                0.407
init...29              <NA>
...30                  <NA>

直接就得出了结果！除了亚组分析的各种结果，还给出了交互作用的P值！

很多人问能不能在公式的右边添加多个协变量，答案是不能啊，我想你可能是搞混了亚组分析和多因素回归的概念，毕竟二者的森林图可以说是长得一模一样！赶紧去看这篇推文：亚组分析和多因素回归的森林图比较

画森林图

这个结果不需要另存为csv也能直接使用（除非你是细节控，需要修改各种大小写等信息），当然如果你需要HR(95%CI)这种信息，还是需要自己添加一下的。

我们选择自己想要展示的信息，比如我这里选择了：变量名字、人数、OR值、OR值的可信区间、P值、交互作用的P值：

plot_df <- res[,c("Variable","Count","OR","Lower","Upper","P value",
                  "P for interaction")]
plot_df

               Variable Count   OR Lower                 Upper P value
rx1             Overall   619 1.99  1.45                  2.75  <0.001
1...1               sex  <NA> <NA>  <NA>                  <NA>    <NA>
init...2         female   312 1.43  0.92                  2.24   0.114
...3               male   307 2.93  1.84                  4.71  <0.001
1...4               age  <NA> <NA>  <NA>                  <NA>    <NA>
init...5            >65   224 2.48  1.45                  4.28   0.001
...6               <=65   395 1.76  1.18                  2.63   0.005
1...7          obstruct  <NA> <NA>  <NA>                  <NA>    <NA>
init...8             No   502 1.94  1.36                  2.77  <0.001
...9                Yes   117 2.24  1.07                  4.75   0.034
1...10           perfor  <NA> <NA>  <NA>                  <NA>    <NA>
init...11            No   602 1.94   1.4                  2.69  <0.001
...12               Yes    17 5.83  0.78                 61.72   0.104
1...13           adhere  <NA> <NA>  <NA>                  <NA>    <NA>
init...14            No   533 2.03  1.44                  2.88  <0.001
...15               Yes    86  1.7  0.72                  4.05   0.229
...16            differ  <NA> <NA>  <NA>                  <NA>    <NA>
...17              well    56 3.78  1.28                 11.92   0.019
...18          moderate   444 1.95  1.34                  2.86   0.001
...19              poor   106 1.61  0.74                  3.54   0.228
...20            extent  <NA> <NA>  <NA>                  <NA>    <NA>
...21         submucosa    18    0  <NA> 1.23555810375665e+213   0.996
...22            muscle    70 2.83  0.97                  9.05   0.064
...23            serosa   500 2.07  1.45                  2.96  <0.001
...24        contiguous    31 1.94  0.42                  9.23   0.392
1...25             surg  <NA> <NA>  <NA>                  <NA>    <NA>
init...26         short   452 2.17  1.49                  3.17  <0.001
...27              long   167 1.54  0.83                  2.86   0.169
1...28            node4  <NA> <NA>  <NA>                  <NA>    <NA>
init...29            No   453 2.21  1.51                  3.26  <0.001
...30               Yes   166 1.61  0.84                  3.11   0.155
          P for interaction
rx1                    <NA>
1...1                 0.031
init...2               <NA>
...3                   <NA>
1...4                 0.317
init...5               <NA>
...6                   <NA>
1...7                 0.734
init...8               <NA>
...9                   <NA>
1...10                0.316
init...11              <NA>
...12                  <NA>
1...13                  0.7
init...14              <NA>
...15                  <NA>
...16                 0.447
...17                  <NA>
...18                  <NA>
...19                  <NA>
...20                 0.087
...21                  <NA>
...22                  <NA>
...23                  <NA>
...24                  <NA>
1...25                0.351
init...26              <NA>
...27                  <NA>
1...28                0.407
init...29              <NA>
...30                  <NA>

我们添加个空列用于显示可信区间，并把不想显示的NA去掉即可，还需要把P值，可信区间这些列变为数值型。

plot_df$` ` <- paste(rep(" ", nrow(plot_df)), collapse = " ")
plot_df[,2:7] <- apply(plot_df[,2:7],2,as.numeric)

Warning in apply(plot_df[, 2:7], 2, as.numeric): NAs introduced by coercion

plot_df[,c(2,6,7)][is.na(plot_df[,c(2,6,7)])] <- " "

画图就非常简单！

library(forestploter)
library(grid)

p <- forest(
  data = plot_df[,c(1,2,8,6,7)],
  lower = plot_df$Lower,
  upper = plot_df$Upper,
  est = plot_df$OR,
  ci_column = 3,
  sizes = (plot_df$OR+0.001)*0.3, 
  ref_line = 1, 
  xlim = c(0,4)
  )

Warning in forest(data = plot_df[, c(1, 2, 8, 6, 7)], lower = plot_df$Lower, :
Missing lower and/or upper limit on column3 row 22

print(p)

这样就搞定了，easy。

---