前情提要

• 在上這節課之前，我們假設你已經有這些能力了：

– 搜索文獻&擷取Paper中需要的資訊

– 在R之內讀取檔案&安裝套件

– 簡單的合併分析&畫出森林圖

– 畫出漏斗圖&檢定它是否對稱

• 今天的課程我們將會依據步驟完成這些部份：

– 非2*2Table的合併

– 次組分析

– Meta-regression

• 由於課程時間的限制，下面這些都只有簡單的範例：

– Trial sequence analysis

– Network meta-analysis

Topic 1 非2*2Table的合併(連續變項)

• 有許多研究的Outcome並非事件的發生，而是一個性狀的改變(如血壓的改變)

– 請點這裡下載練習檔案

• 這個資料是出自這篇文章的範例：

Topic 1 連續變項的分析

• 下面這串是讀取檔案的指令

Example1=read.csv("d:/Example1.csv") #讀取檔案後存成Example1物件(資料表格式)

Example1 #看Example1的內容

##           study n.amlo mean.amlo var.amlo n.plac mean.plac var.plac
## 1  Protocol 154     46    0.2316   0.2254     48   -0.0027   0.0007
## 2  Protocol 156     30    0.2811   0.1441     26    0.0270   0.1139
## 3  Protocol 157     75    0.1894   0.1981     72    0.0443   0.4972
## 4 Protocol 162A     12    0.0930   0.1389     12    0.2277   0.0488
## 5  Protocol 163     32    0.1622   0.0961     34    0.0056   0.0955
## 6  Protocol 166     31    0.1837   0.1246     31    0.0943   0.1734
## 7 Protocol 303A     27    0.6612   0.7060     27   -0.0057   0.9891
## 8  Protocol 306     46    0.1366   0.1211     47   -0.0057   0.1291

• 下面是在meta-analysis中分析連續變項的程式碼，主要是依靠metacont()函數

– 注意裡面的prediction = TRUE，這是在說我希望呈現prediction interval。

Result1 = metacont(n.e = n.amlo, mean.e = mean.amlo, sd.e = sqrt(var.amlo),
n.c = n.plac, mean.c = mean.plac, sd.c = sqrt(var.plac),
data = Example1, studlab = study, prediction = TRUE)
Result1 #看分析結果

##                    MD            95%-CI %W(fixed) %W(random)
## Protocol 154   0.2343 [ 0.0969; 0.3717]     21.22      17.47
## Protocol 156   0.2541 [ 0.0663; 0.4419]     11.35      12.74
## Protocol 157   0.1451 [-0.0464; 0.3366]     10.92      12.45
## Protocol 162A -0.1347 [-0.3798; 0.1104]      6.67       9.04
## Protocol 163   0.1566 [ 0.0072; 0.3060]     17.94      16.21
## Protocol 166   0.0894 [-0.1028; 0.2816]     10.85      12.40
## Protocol 303A  0.6669 [ 0.1758; 1.1580]      1.66       2.90
## Protocol 306   0.1423 [-0.0015; 0.2861]     19.39      16.79
## 
## Number of studies combined: k=8
## 
##                          MD            95%-CI     z  p.value
## Fixed effect model   0.1619 [ 0.0986; 0.2252] 5.013 < 0.0001
## Random effects model 0.1589 [ 0.0710; 0.2467] 3.544   0.0004
## Prediction interval         [-0.0680; 0.3858]               
## 
## Quantifying heterogeneity:
## tau^2 = 0.0066; H = 1.33 [1; 2]; I^2 = 43.2% [0%; 74.9%]
## 
## Test of heterogeneity:
##      Q d.f.  p.value
##  12.33    7   0.0902
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

Topic 1 連續變項的分析(續)

• 接著當然是想辦法畫出森林圖&漏斗圖，聰明的妳一定想的到，這與分析類別變項的程式碼相同：

– 森林圖

forest(Result1)

– 漏斗圖&漏斗圖的檢定

metabias(Result1,k.min=3)

## 
##  Linear regression test of funnel plot asymmetry
## 
## data:  Result1
## t = 0.3378, df = 6, p-value = 0.747
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
##    bias se.bias   slope 
##  0.5885  1.7422  0.1104

funnel(Result1)

Topic 1 非2*2Table的合併(整合參數)

• 在大部分的情況下，我們恐怕都只能拿到各研究的整合指標，如Hazzard ratio(HR)之類的。

• 事實上，這樣的資訊就夠了，還記得上禮拜講到的公式嗎？Meta-analysis一樣是先計算出各研究的RR，在加權合併它。

– 如果你忘了，請點這裡複習一下上禮拜的部份(Page 11 to 14)

• 然而，由於RR/OR是log尺度的，因此我們必須給的是log轉換過後的點估計，但通常paper上給的是原始的點估計&信賴區間，因此我們需要自行幫他轉換。

– 下面這是轉換的公式

Topic 1 整合參數的分析

• 請點這裡下載檔案

– 這是上週我們使用的Folic acid的資料，只是現在只留下了RR&他的95% CI

Example2=read.csv("d:/RR_95%CI.csv") #讀取檔案後存成Example2物件(資料表格式)
Example2$logRR=log(Example2$RR) #計算logRR
Example2$SE1=(log(Example2$RR)-log(Example2$lower))/1.96 #使用下界計算SE
Example2$SE2=(log(Example2$upper)-log(Example2$RR))/1.96 #使用上界計算SE

• 接著，我們將使用metagen()來使用這些數字

Result2=metagen(TE=logRR,seTE=SE1,sm="RR",studlab=Study,data=Example2)
Result2 #看看分析結果

• 聰明的妳一定猜的到森林圖和漏斗圖怎樣畫…,這裡我就只放上森林圖了

forest(Result2)

練習1

• 請點這裡下載練習檔案，並利用每個研究的HR以及95% CI，算出其log(HR)以及其標準誤，並在Excel/R裡面增加兩個新變項，名稱自訂

• 接著，請你利用R讀取檔案，並且利用類似於剛剛的教學，並畫出森林圖&漏斗圖

• 最後，請你順便想一下，如果使用下界計算出來的標準誤，和使用上界計算出來的標準誤有很大的不同，那該怎麼辦?

Topic 2 - 次組分析

• 我們時常會發現研究之間存在很高的異質性，這時候該怎麼辦呢`?

– 也許在合併這些文獻之前，我們早就懷疑有些因子可能會影響治療效果了！(注意，按照Cochrane的Guildline，次組分析必須在事前有機轉/流行病學等研究提出合理的懷疑)

• 我們使用上週用的檔案作為範例，請點這裡下載檔案

– 不要忘記怎麼讀檔&分析，我們使用上週的程式碼

Folic=read.csv("d:/Folic acid.csv")
Folic
Result=metabin(event.e = ev.trt, n.e = n.trt, event.c = ev.ctrl, n.c =  n.ctrl, studlab = name,
data = Folic, sm = "RR", method.tau="DL")
Result

• 假定我們現在想要使用“food”這個變項進行次組分析，將研究分為3組，我們只需要這樣修改

Result1=metabin(event.e = ev.trt, n.e = n.trt, event.c = ev.ctrl, n.c =  n.ctrl, studlab = name,
data = Folic, sm = "RR", method.tau="DL", byvar=food)
Result1
forest(Result1, bylab="Food fortification")

Topic 2 - 次組分析(續)

• 假定我們想要以出版年代做區分，區分為2007年以前和2008年以後的文章，我們可以創造出一個新變項再使用它進行次組分析：

Folic$pubyear[Folic$year<=2007]="Early than 2007"
Folic$pubyear[Folic$year>=2008]="Later than 2008"
Result2=metabin(event.e = ev.trt, n.e = n.trt, event.c = ev.ctrl, n.c =  n.ctrl, studlab = name,
data = Folic, sm = "RR", method.tau="DL", byvar=pubyear)
forest(Result2, bylab="Year of publish")

練習2

• 請利用剛剛的方式，將這些研究分為糖尿病盛行率在30%以下/以上的兩個次組，進行分析並畫出森林圖。

– 你可以在Excel裡面修改你的原始檔，增加新變數後重新讀檔進R，之後再進行分析。

Topic 3 - Meta-regression

• 學會次組分析以後，我們還有個最大的問題，那就是難道我們不能使用連續變項來解釋異質性，一定得把研究分成幾個組別？

– Meta-regression基本上就是次組分析的延伸，而且它能夠容許同時存在兩個以上的因子。

• 這樣聽起來還是很難懂，我們用下面這個例子來學習：

• 這是你剛剛畫出來的圖，注意到兩個次組的RR各是0.89以及1.00。

– 因此，若我們將DM = < 30%做為參考組，那我們的交互作用項理論上應該要是0.99/0.89=1.11 (以Fixed model為例)

– 當然，在Meta-regression中是log尺度的，所以實際上的數字應該是log(1.11)=0.11

• 了解這點以後，我們迅速的使用Meta-regression測試看看數字是否就是0.11左右。

Topic 3 - Meta-regression(續)

• 下面這裡是分析的程式碼，我們必須使用另外一個套件(metafor)

library(metafor) #載入metafor套件，若尚未安裝則需先安裝
Folic$DM=c(0,1,0,0,1,0,0,1,0,0,1,1,1)
res1=rma(ai = ev.trt, ci = ev.ctrl, n1i = n.trt, n2i = n.ctrl, mods=~DM, data=Folic, method="FE",measure="RR")
res1

• 我們可以看到下面這張報表：

• 根據這張報表，我們將可以得到下列這些資訊：

– 當DM=0(也就是DM盛行率<30%)時，合併的結果為exp(-0.1173)=0.889

– 當DM=1(也就是DM盛行率>30%)時，合併的結果為exp(-0.1173+0.1092)=0.992

– DM的調節效果(或稱交互作用效果)為exp(0.1092)=1.115

Topic 3 - 更複雜的例子

• 當然，我們可以多增加幾個變項，假設我們同時想要看DM&food，其中DM改成連續變項，food是個三類別變項。

– 另外，這次我們用Mixed effects model：

library(metafor) #載入metafor套件，若尚未安裝則需先安裝
Folic$DMr=c(0.136,0.455,0.193,0.232,0.55,0,0.226,0.4,0.11,0.24,0.4,1,0.405) #產生一個連續變項
res2=rma(ai = ev.trt, ci = ev.ctrl, n1i = n.trt, n2i = n.ctrl, mods=~DMr+factor(food), data=Folic, method="DL",measure="RR")
res2

##Mixed-Effects Model (k = 13; tau^2 estimator: DL)
##
##tau^2 (estimated amount of residual heterogeneity):     0 (SE = 0.0079)
##tau (square root of estimated tau^2 value):             0
##I^2 (residual heterogeneity / unaccounted variability): 0.00%
##H^2 (unaccounted variability / sampling variability):   1.00
##R^2 (amount of heterogeneity accounted for):            100.00%
##
##Test for Residual Heterogeneity: 
##QE(df = 9) = 8.7522, p-val = 0.4605
##
##Test of Moderators (coefficient(s) 2,3,4): 
##QM(df = 3) = 4.3740, p-val = 0.2238
##
##Model Results:
##
##                                 se     zval    pval    ci.lb    ci.ub   
##intrcpt             -0.1621  0.0763  -2.1251  0.0336  -0.3117  -0.0126  *
##DMr                  0.4535  0.3464   1.3091  0.1905  -0.2255   1.1326   
##factor(food)partly  -0.0830  0.1370  -0.6056  0.5448  -0.3515   0.1856   
##factor(food)yes     -0.0728  0.1488  -0.4889  0.6249  -0.3644   0.2189   
##
##---
##Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

• 根據上面這張報表，我們可以寫出一個這樣的meta-regression

練習3

• 利用同一個資料表，請創造出兩個新的變項（可以在Excel之內作完再重新讀檔）

– 將dose(葉酸劑量)分為5mg/day以下&以上

– 將hcy(介入組血液中HCY的下降率)分為-20%以上&以下

• 請利用這兩個變項創造出一個Meta-regression，並試著去解釋他的內涵，你必須要告訴我…

– dose < 5mg/day 和 hcy > -20% 之研究的合併RR

– dose > 5mg/day 和 hcy > -20% 之研究的合併RR

– dose < 5mg/day 和 hcy < -20% 之研究的合併RR

– dose > 5mg/day 和 hcy < -20% 之研究的合併RR

– 如果你很快的做完了，試著想一下Ecological fallacy，他是否會發生在meta-analysis之中？

其他meta-analysis related method的介紹 - Trial sequence analysis(1)

• 現在有越來越多人在做meta-analysis之前，會像做其他研究一樣，先進行樣本數估計，而在meta-analysis中，由於樣本的增加是隨著Ttial的增加而變多的，所以我們稱作Trial sequence analysis

– 另外，由於是Ttial的增加是需要較長時間的，因此可以將其視為有非常多個個期中分析的臨床試驗研究。在一般的臨床試驗中，對期中分析的顯著水準會做一定程度的懲罰(避免多重比較問題)，因此Trial sequence analysis同樣會懲罰樣本不足的meta-analysis。

• 下面這張圖是使用上週下載的檔案分析的結果。

其他meta-analysis related method的介紹 - Trial sequence analysis(2)

• R裡面目前沒有相關的套件，目前最常用的是TSA軟體。

– 下面這是這個軟體的官方網站：

其他meta-analysis related method的介紹 - Network meta analysis

• 有的時候針對同一種族群，我們可能同時有多個藥物可以治療，這時候我們可能需要一次整合這些資訊。

– R裡面的netmeta套件可以幫助我們進行Network meta analysis。(下面是範例)

library(netmeta)
data(Senn2013)
net1 <- netmeta(TE, seTE, treat1, treat2, studlab,data=Senn2013, sm="MD", reference="plac")

summary(net1)

## Number of studies: k=26
## Number of treatments: n=10
## Number of pairwise comparisons: m=28
## 
## Fixed effect model
## 
## Treatment estimate (sm='MD', reference.group='plac'):
##          MD            95%-CI
## acar -0.827  [-1.040; -0.615]
## benf -0.905  [-1.154; -0.656]
## metf -1.114  [-1.231; -0.997]
## migl -0.944  [-1.193; -0.695]
## piog -1.066  [-1.215; -0.918]
## plac  0.000  [ 0.000;  0.000]
## rosi -1.202  [-1.295; -1.108]
## sita -0.570  [-0.823; -0.317]
## sulf -0.439  [-0.619; -0.260]
## vild -0.700  [-0.950; -0.450]
## 
## Quantifying heterogeneity/inconsistency:
## tau^2 = 0.1087; I^2 = 81.4%
## 
## Test of heterogeneity/inconsistency:
##      Q d.f.  p.value
##  96.99   18 < 0.0001

netgraph(net1)

小結

• 這兩堂課之中，除了讀取檔案之外，我們總共學了下列這些程式碼：

– metabin(), metacont(), metagen()

– forest(), funnel(), metabias()

– rma()

• 你可以訪問下列三個網站以取得更詳細的資訊

– Cochrane Handbook for Systematic Reviews)

– meta package說明

– metafor package說明