[統合分析] 如何生成single arm meta-analysis 的圖表

在統合分析裡面, 通常會有兩個arm, 不過有的時候會需要對症狀的發生率, 疾病的發生率進行統合, 此時就需要做所謂的single arm meta-analysis, R 裡面有非常好用的套件, 以下就由我來示範!

首先假設假想例子:

Study 1: 30個病人有5個中風, 15男有2個中風, 15女有3個中風

Study 2: 40個病人有7個中風, 20男有3個中風, 20女有4個中風

Study 3: 20個病人有9個中風, 10男有4個中風, 10女有5個中風

Study 4: 70個病人有10個中風, 30男有5個中風, 40女有5個中風

則可以使用R 的 metaprop 函式來生成森林圖(forest plot), 其程式碼如下:

library(meta)
library(metafor)
library("grid")

studies <- c("Study 1","Study 2","Study 3","Study 4")
obs <- c(5,7,9,10)
denom <- c(30,40,20,70)


m1 <- metaprop(obs, denom, studies, comb.fixed=T) #com.fixed 表示fixed effect model


forest(m1, leftlabs=c("Study","Stroke event","Total"), 
       ref=0.2, colgap=unit(3,"mm"), col.by="black",fontsize=10,
       rowgap=unit(3,"mm"))

生成圖片如下:

                                                                                                                                           

其中圖裡面的 I^2 就是統合分析中拿來評估異質性(heterogenity)的指標, 值越低表示異質性越低, 意指為統合分析中研究間的差異性越小

另外就是統合分析裡面, 經常要看所謂的publication bias, 程式碼也是非常簡單, 如下

funnel(m1)

metabias(m1) #Egger's test

如圖可以看到, 明顯兩邊不對稱, 雖然這麼少的文獻要看publication bias不太夠(依照目前cochrane handbook的建議, 至少要十篇論文來生成funnel plot會比較有意義), 不過一般來說不對稱就代表可能有publication bias!

subgroup analysis, 一般中文翻作次群體分析, 以這個例子為例, 上面有提供男女資訊, 那我們就會想知道男女是否會有差異, 一樣用程式碼做示範,

studies2 <- c("Study 1","Study 2","Study 3","Study 4",
              "Study 1","Study 2","Study 3","Study 4")

obs2 <- c(2,3,4,5,3,4,5,5)

denom2 <- c(15,20,10,30,15,20,10,40)

setting<-c("M", "M", "M","M",
           "F", "F", "F", "F")

m2 <- metaprop(obs2, denom2, studies2, comb.fixed=F,byvar=setting,
) #com.fixed 表示fixed effect model

forest(m2, leftlabs=c("Study","Stroke event","Total"), 
       ref=0.1, colgap=unit(5,"mm"), col.by="black",fontsize=10,
       rowgap=unit(3,"mm"))

之後再利用meta-regression, 替男女兩組做比較, 看是否有差異!

df1 <- structure(list(study_nr = c(1,2,3,4,1,2,3,4),
                     group = c("A", "A", "A", "A",
                               "B", "B", "B", "B"), 
                     n_tot = c(15,20,10,30,15,20,10,40), 
                     n_event = c(2,3,4,5,3,4,5,5)), 
                class = c("tbl_df", "tbl", "data.frame"),
                row.names = c(NA, -8L))
df1

# 計算 log(odds) 以及 corresponding sampling variances
df1 <- escalc(measure="PLO", xi=n_event, ni=n_tot, data=df1)
df1

# fit meta-regression model with 'group' as moderator
res1 <- rma(yi, vi, mods = ~ group, data=df1, method="DL")
res1

# 得到 odd ratio 和 95% CI
predict(res1, newmods=1, intercept=FALSE, transf=exp, digits=2)

Test for Residual Heterogeneity:

QE(df = 6) = 9.3173, p-val = 0.1565

p value >=0.05, 沒有顯著意義!

reference: 

1. metafor package

2. stackoverflow

3. https://rpubs.com/pekong/532068