發生率是常用的效應指標。兩組發生率的比值稱為發生率比,可用于比較兩組人群中單位人時內某事件發生情況的差異。目前,RevMan軟件尚無法實現以發生率比為效應量的Meta分析。本文通過一組模擬數據展示了如何使用R Studio軟件的meta程序包實現以發生率比為效應量的Meta分析過程。
發生率(incidence rate,IR),也稱為人時率(person-time rate),是指一定觀察期內,某人群中新發生某事件的人數與觀察的總人時之比[1]。這個事件可以是疾病發生、死亡等。兩組IR的比值稱為發生率比(incidence rate ratio,IRR),用于比較兩組人群中單位人時內某事件發生情況的差異,常見于隊列研究、臨床試驗、篩查研究等。目前,RevMan軟件尚無法實現以IRR為效應量的Meta分析。為此,本研究使用一組模擬數據,詳細介紹基于R Studio軟件的meta程序包實現以IRR為效應量的Meta分析過程,并討論了以IRR為效應量進行Meta分析的應用場景。
1 模擬數據
本研究以“某急性傳染病的流行對乳腺癌篩查的影響”為例,假設:某急性傳染病的流行會對單位時間內的乳腺癌篩查人數產生影響,模擬生成從18項研究中提取的不同地區某急性傳染病流行前后乳腺癌篩查數據,詳見表1。在真實情況下,急性傳染病流行前后記錄的乳腺癌篩查時長可能不同。所以,本研究的模擬數據中,設定某急性傳染病流行前后乳腺癌篩查的時長不同。“Study”代表不同研究,“Year”代表出版年,“Country”代表研究地區所在國家,“Continent”代表國家所在洲,“Preevents”代表某地區某急性傳染病流行前的一段時間內進行乳腺癌篩查的人數,“Pretime”代表某地區某急性傳染病流行前記錄的乳腺癌篩查人數所對應的篩查時長,單位為天。“Duringevents”代表某地區某急性傳染病流行后的一段時間內進行乳腺癌篩查的人數,“Duringtime”代表某地區某急性傳染病流行后記錄的乳腺癌篩查人數所對應的篩查時長,單位為天。使用Microsoft Excel存儲模擬數據,命名為breast_cancer。

本研究中的發生率為單位時間(天)內的乳腺癌篩查率,傳染病流行前單位時間內的乳腺癌篩查率IR1=Preevents/Pretime;傳染病流行后單位時間內的乳腺癌篩查率IR2=Duringevents/Duringtime,IRR為傳染病流行后乳腺癌篩查率與傳染病流行前乳腺癌篩查率的比值,即:IRR=IR2/IR1。
2 R Studio軟件和meta程序包的安裝
R語言軟件是一款免費的開源軟件[2],可以從官方網站根據自身電腦運行環境選擇合適版本的R程序包。根據安裝提示,即可完成R軟件的安裝。R Studio軟件是一組有助于便捷高效利用R的專業集成工具[3],本研究分析過程所用R Studio軟件版本為4.2.1。安裝完成后,打開R Studio軟件,在左上方的命令窗口輸入相應代碼安裝meta程序包,代碼如下:
install.packages("meta") ##安裝meta程序包
library(meta) ##加載meta程序包
輸入代碼后可點擊命令窗口右上方“Run”,運行代碼。
3 以IRR為效應量的Meta分析操作過程
3.1 數據輸入
R Studio軟件可以采用多種方式導入數據。由于本研究的模擬數據存儲在Excel中,可以利用readxl程序包將文件名為breast_cancer的Excel數據導入R Studio軟件,代碼如下:
install.packages(“readxl”) ##安裝readxl程序包
library(readxl) ##加載readxl程序包
breast_cancer <- read_excel("~/Desktop/breast_cancer.xlsx") ##數據導入,括號內為模擬數據存儲位置,具體操作時應按數據存儲的實際位置進行修改。
點擊“Run”運行代碼。
3.2 Meta分析過程
metainc()函數是R Studio軟件的meta分析程序包中專門用于計算合并IRR的程序。其基本命令如下:
metainc(event.e,time.e,event.c,time.c,data,sm,studlab,method)
event.e為試驗組或暴露組的事件發生數,time.e為試驗組或暴露組的觀察時長,event.c為對照組或非暴露組的事件發生數,time.c為對照組或非暴露組的觀察時長,data為數據集,sm為合并效應量,studlab為標簽,method為合并效應量的算法。
本研究通過調用metainc()函數,實現以IRR為效應量的Meta分析,具體運行程序如下:
dataA<-metainc(Duringevents,Duringtime,Preevents,Pretime,data=breast_cancer,sm="IRR",studlab=paste(breast_cancer$Study,breast_cancer$Year,sep="-"),method="Inverse") ##metainc函數默認同時使用固定效應模型(common)和隨機效應模型(random)計算合并效應量。當納入的研究異質性較大時,研究者可以選擇隨機效應模型;默認的合并效應量算法為Mantel Haenszel(MH),其它算法還包括Cochran、廣義線性模型(generalised linear mixed model,GLMM)。
dataA
點擊“Run”運行代碼。運行結果如表2所示,I2>50%表明各研究之間的統計學異質性較大,故選擇隨機效應模型的Meta分析結果。合并效應量IRR=0.621 4[95%Cl(0.542 7,0.711 5),P<0.01],結果表明,急性傳染病流行后乳腺癌篩查率是流行前乳腺癌篩查率的62.14%,差異有統計學意義。

3.3 繪制森林圖
使用R Studio軟件meta程序包中的forest()命令繪制森林圖,代碼如下:
forest(dataA,label.e="During",label.c="Before",common=FALSE,col.square="orange",col.diamond.random="skyblue",test.overall=TRUE) ##通過label.e和label.c分別定義試驗組/暴露組和對照組/非暴露組的標簽;通過common=FALSE實現只顯示隨機效應模型的Meta分析結果,通過random=FALSE可以實現只顯示固定效應模型的Meta分析結果;通過col.square調整森林圖中代表單個研究效應量的方塊顏色,orange代表橘黃色;通過col.diamond.random調整森林圖中代表合并效應量的菱形顏色,skyblue代表天藍色;通過輸入test.overall=TRUE即可在森林圖中顯示顯著性檢驗結果,不輸入test.overall參數時則不會顯示顯著性檢驗結果。
點擊“Run”運行代碼,結果見圖1。

3.4 亞組分析
本研究以數據來源的國家所在洲(continent)為分組因素進行亞組分析,代碼如下:
dataB<-metainc(Duringevents,Duringtime,Preevents,Pretime,data=breast_cancer,sm="IRR", studlab=paste(breast_cancer$Study,breast_cancer$Year,sep="-"),common=FALSE,method="Inverse",subgroup=breast_cancer$Continent) ##通過subgroup進行亞組分析
dataB
點擊“Run”運行代碼,結果見表3。與流行前相比,四個洲的急性傳染病流行后乳腺癌篩查率均顯著降低。但是,亞組間比較的P<0.01,提示不同洲急性傳染病流行后乳腺癌篩查率的下降程度不同。

繪制亞組分析結果的森林圖,代碼如下:
forest(dataB,label.e="During",label.c="Before",col.square="orange",col.diamond.random="skyblue",test.overall=TRUE,test.effect.subgroup=TRUE) ##通過輸入test.effect.subgroup=TRUE即可在森林圖中顯示各亞組的顯著性檢驗結果,不輸入此參數則不會顯示各亞組的顯著性檢驗結果。
點擊“Run”運行代碼,結果見圖2。

3.5 敏感性分析
使用R Studio軟件meta程序包中的metainf()函數逐一剔除研究后計算合并效應量IRR,進行敏感性分析,代碼如下:
metainf(dataA,pooled="random") ##通過輸入pooled="random"選擇隨機效應模型進行Meta分析,不輸入該參數時默認固定效應模型。
點擊“Run”運行代碼,運行結果如表4。

繪制敏感性分析的森林圖,代碼如下:
forest(metainf(dataA,pooled = "random"))
點擊“Run”運行代碼,結果如圖3。

3.6 發表偏倚檢驗
在meta程序包中,可以通過funnel()函數繪制漏斗圖,定性判斷發表偏倚,代碼如下:
funnel(dataA,common=FALSE) ##通過輸入common=FALSE實現只顯示使用隨機效應模型的計算結果。
點擊“Run”運行代碼,結果如圖4。

研究者也可以通過meta程序包中metabias()函數定量檢驗漏斗圖的對稱性,代碼如下:
metabias(dataA,method.bias="Egger") ##metabia()函數中常用的算法包括“Begg”、“Egger”和“Thompson”。如果不對算法進行規定(缺少參數method.bias),除OR值作為合并效應量的Meta分析外,默認method.bias="Egger"。此外,默認在研究數量為10或更大的情況下才進行不對稱性檢驗[4]。
結果如表5所示,P=0.034,漏斗圖不對稱,提示可能存在發表偏倚。

4 討論
本文以癌癥篩查研究為例,詳細介紹了以IRR為效應量進行Meta分析的R Studio軟件實現過程。以某急性傳染病的流行作為暴露因素,將該急性傳染病流行前的人群作為非暴露組,將該急性傳染病流行后的人群作為暴露組。IRR為該急性傳染病流行后乳腺癌篩查率與該急性傳染病流行前乳腺癌篩查率的比值。合并后的IRR反映了該急性傳染病流行對乳腺癌篩查的影響程度。在比較兩組IR時,metainc()函數還可以計算合并發生率差異(Incidence Rate Difference, sm="IRD"),平方根變換的發生率差(Square root transformed Incidence Rate Difference, sm="IRSD"),疫苗效力或疫苗有效性(Vaccine efficacy or vaccine effectiveness, sm="VE")。除R Studio軟件外,也有研究使用Stata[5-6]和Commercial Meta-Analysis Software[7-8]進行以IRR為效應量的Meta分析。
除了篩查研究以外,IRR也常見于隊列研究[9-11]和隨機對照試驗[7,12-13]。但是,不同研究對IRR的定義可能不同[14-16]。例如,一項納入隊列研究的Meta分析[5]評估了炎癥性腸病患者發生帶狀皰疹的風險。暴露組的IR為炎癥性腸病患者發生帶狀皰疹的人數除以隨訪患者人年數。非暴露組的IR為非炎癥性腸病患者發生帶狀皰疹的人數除以隨訪患者人年數。IRR=暴露組的IR/非暴露組的IR。合并后的IRR反映了患有炎癥性腸病對發生帶狀皰疹的影響程度。Marilly等[17]開展了一項評估鈉-葡萄糖共轉運蛋白2抑制劑(SGLT2i)治療2型糖尿病風險和獲益的Meta分析。該研究納入了5個隨機對照試驗,以總體死亡率作為Meta分析的主要結局指標。試驗組的IR為試驗組死亡人數除以試驗組隨訪人年數。對照組的IR為對照組死亡人數除以對照組隨訪人年數。IRR=試驗組的IR/對照組的IR。合并后的IRR反映了SGLT2i與總體死亡率的關聯程度。
IRR、RR(risk ratio)和HR(hazard ratio)均是反映事件發生相對風險的指標[18-20]。RR的計算過程通常不涉及時間變量。當試驗/暴露組與對照/非暴露組的隨訪時間不同時,合并后的RR可能存在估計不精確的問題。HR的計算涉及終點事件和發生終點事件的時間。但是,一般需要通過Cox回歸估計HR[21]。當一些原始研究未提供具體的HR值時,可能導致這些研究的結果無法進行Meta分析。IR的分子是在觀察期內新發生終點事件的人數,分母是人群的總觀察人時,描述了終點事件在人群中的發生速度。常用于需要長期隨訪的隊列研究,特別是當不斷有研究對象進入、退出研究或失訪時,IR可以準確描述終點事件在人群中的發生概率[1]。因此,如果研究提供了事件發生數及篩查或隨訪時長,可以考慮計算和合并IRR。當終點事件重復發生時,在IRR的計算過程中,IR的分子可以是終點事件的發生頻次。例如:一項Meta分析[22]調查了運動員受傷率的性別差異。在隨訪期內,同一位運動員可能出現多次受傷的情況。IRR=(女運動員受傷次數/女運動員隨訪時間)/(男運動員受傷次數/男運動員隨訪時間)。但是,RR和HR不能用于終點事件重復發生的相對風險估計。
綜上所述,本文詳細介紹了使用R Studio軟件的meta程序包實現以IRR為效應量的Meta分析過程,可為開展同類研究提供重要的方法學參考。
發生率(incidence rate,IR),也稱為人時率(person-time rate),是指一定觀察期內,某人群中新發生某事件的人數與觀察的總人時之比[1]。這個事件可以是疾病發生、死亡等。兩組IR的比值稱為發生率比(incidence rate ratio,IRR),用于比較兩組人群中單位人時內某事件發生情況的差異,常見于隊列研究、臨床試驗、篩查研究等。目前,RevMan軟件尚無法實現以IRR為效應量的Meta分析。為此,本研究使用一組模擬數據,詳細介紹基于R Studio軟件的meta程序包實現以IRR為效應量的Meta分析過程,并討論了以IRR為效應量進行Meta分析的應用場景。
1 模擬數據
本研究以“某急性傳染病的流行對乳腺癌篩查的影響”為例,假設:某急性傳染病的流行會對單位時間內的乳腺癌篩查人數產生影響,模擬生成從18項研究中提取的不同地區某急性傳染病流行前后乳腺癌篩查數據,詳見表1。在真實情況下,急性傳染病流行前后記錄的乳腺癌篩查時長可能不同。所以,本研究的模擬數據中,設定某急性傳染病流行前后乳腺癌篩查的時長不同。“Study”代表不同研究,“Year”代表出版年,“Country”代表研究地區所在國家,“Continent”代表國家所在洲,“Preevents”代表某地區某急性傳染病流行前的一段時間內進行乳腺癌篩查的人數,“Pretime”代表某地區某急性傳染病流行前記錄的乳腺癌篩查人數所對應的篩查時長,單位為天。“Duringevents”代表某地區某急性傳染病流行后的一段時間內進行乳腺癌篩查的人數,“Duringtime”代表某地區某急性傳染病流行后記錄的乳腺癌篩查人數所對應的篩查時長,單位為天。使用Microsoft Excel存儲模擬數據,命名為breast_cancer。

本研究中的發生率為單位時間(天)內的乳腺癌篩查率,傳染病流行前單位時間內的乳腺癌篩查率IR1=Preevents/Pretime;傳染病流行后單位時間內的乳腺癌篩查率IR2=Duringevents/Duringtime,IRR為傳染病流行后乳腺癌篩查率與傳染病流行前乳腺癌篩查率的比值,即:IRR=IR2/IR1。
2 R Studio軟件和meta程序包的安裝
R語言軟件是一款免費的開源軟件[2],可以從官方網站根據自身電腦運行環境選擇合適版本的R程序包。根據安裝提示,即可完成R軟件的安裝。R Studio軟件是一組有助于便捷高效利用R的專業集成工具[3],本研究分析過程所用R Studio軟件版本為4.2.1。安裝完成后,打開R Studio軟件,在左上方的命令窗口輸入相應代碼安裝meta程序包,代碼如下:
install.packages("meta") ##安裝meta程序包
library(meta) ##加載meta程序包
輸入代碼后可點擊命令窗口右上方“Run”,運行代碼。
3 以IRR為效應量的Meta分析操作過程
3.1 數據輸入
R Studio軟件可以采用多種方式導入數據。由于本研究的模擬數據存儲在Excel中,可以利用readxl程序包將文件名為breast_cancer的Excel數據導入R Studio軟件,代碼如下:
install.packages(“readxl”) ##安裝readxl程序包
library(readxl) ##加載readxl程序包
breast_cancer <- read_excel("~/Desktop/breast_cancer.xlsx") ##數據導入,括號內為模擬數據存儲位置,具體操作時應按數據存儲的實際位置進行修改。
點擊“Run”運行代碼。
3.2 Meta分析過程
metainc()函數是R Studio軟件的meta分析程序包中專門用于計算合并IRR的程序。其基本命令如下:
metainc(event.e,time.e,event.c,time.c,data,sm,studlab,method)
event.e為試驗組或暴露組的事件發生數,time.e為試驗組或暴露組的觀察時長,event.c為對照組或非暴露組的事件發生數,time.c為對照組或非暴露組的觀察時長,data為數據集,sm為合并效應量,studlab為標簽,method為合并效應量的算法。
本研究通過調用metainc()函數,實現以IRR為效應量的Meta分析,具體運行程序如下:
dataA<-metainc(Duringevents,Duringtime,Preevents,Pretime,data=breast_cancer,sm="IRR",studlab=paste(breast_cancer$Study,breast_cancer$Year,sep="-"),method="Inverse") ##metainc函數默認同時使用固定效應模型(common)和隨機效應模型(random)計算合并效應量。當納入的研究異質性較大時,研究者可以選擇隨機效應模型;默認的合并效應量算法為Mantel Haenszel(MH),其它算法還包括Cochran、廣義線性模型(generalised linear mixed model,GLMM)。
dataA
點擊“Run”運行代碼。運行結果如表2所示,I2>50%表明各研究之間的統計學異質性較大,故選擇隨機效應模型的Meta分析結果。合并效應量IRR=0.621 4[95%Cl(0.542 7,0.711 5),P<0.01],結果表明,急性傳染病流行后乳腺癌篩查率是流行前乳腺癌篩查率的62.14%,差異有統計學意義。

3.3 繪制森林圖
使用R Studio軟件meta程序包中的forest()命令繪制森林圖,代碼如下:
forest(dataA,label.e="During",label.c="Before",common=FALSE,col.square="orange",col.diamond.random="skyblue",test.overall=TRUE) ##通過label.e和label.c分別定義試驗組/暴露組和對照組/非暴露組的標簽;通過common=FALSE實現只顯示隨機效應模型的Meta分析結果,通過random=FALSE可以實現只顯示固定效應模型的Meta分析結果;通過col.square調整森林圖中代表單個研究效應量的方塊顏色,orange代表橘黃色;通過col.diamond.random調整森林圖中代表合并效應量的菱形顏色,skyblue代表天藍色;通過輸入test.overall=TRUE即可在森林圖中顯示顯著性檢驗結果,不輸入test.overall參數時則不會顯示顯著性檢驗結果。
點擊“Run”運行代碼,結果見圖1。

3.4 亞組分析
本研究以數據來源的國家所在洲(continent)為分組因素進行亞組分析,代碼如下:
dataB<-metainc(Duringevents,Duringtime,Preevents,Pretime,data=breast_cancer,sm="IRR", studlab=paste(breast_cancer$Study,breast_cancer$Year,sep="-"),common=FALSE,method="Inverse",subgroup=breast_cancer$Continent) ##通過subgroup進行亞組分析
dataB
點擊“Run”運行代碼,結果見表3。與流行前相比,四個洲的急性傳染病流行后乳腺癌篩查率均顯著降低。但是,亞組間比較的P<0.01,提示不同洲急性傳染病流行后乳腺癌篩查率的下降程度不同。

繪制亞組分析結果的森林圖,代碼如下:
forest(dataB,label.e="During",label.c="Before",col.square="orange",col.diamond.random="skyblue",test.overall=TRUE,test.effect.subgroup=TRUE) ##通過輸入test.effect.subgroup=TRUE即可在森林圖中顯示各亞組的顯著性檢驗結果,不輸入此參數則不會顯示各亞組的顯著性檢驗結果。
點擊“Run”運行代碼,結果見圖2。

3.5 敏感性分析
使用R Studio軟件meta程序包中的metainf()函數逐一剔除研究后計算合并效應量IRR,進行敏感性分析,代碼如下:
metainf(dataA,pooled="random") ##通過輸入pooled="random"選擇隨機效應模型進行Meta分析,不輸入該參數時默認固定效應模型。
點擊“Run”運行代碼,運行結果如表4。

繪制敏感性分析的森林圖,代碼如下:
forest(metainf(dataA,pooled = "random"))
點擊“Run”運行代碼,結果如圖3。

3.6 發表偏倚檢驗
在meta程序包中,可以通過funnel()函數繪制漏斗圖,定性判斷發表偏倚,代碼如下:
funnel(dataA,common=FALSE) ##通過輸入common=FALSE實現只顯示使用隨機效應模型的計算結果。
點擊“Run”運行代碼,結果如圖4。

研究者也可以通過meta程序包中metabias()函數定量檢驗漏斗圖的對稱性,代碼如下:
metabias(dataA,method.bias="Egger") ##metabia()函數中常用的算法包括“Begg”、“Egger”和“Thompson”。如果不對算法進行規定(缺少參數method.bias),除OR值作為合并效應量的Meta分析外,默認method.bias="Egger"。此外,默認在研究數量為10或更大的情況下才進行不對稱性檢驗[4]。
結果如表5所示,P=0.034,漏斗圖不對稱,提示可能存在發表偏倚。

4 討論
本文以癌癥篩查研究為例,詳細介紹了以IRR為效應量進行Meta分析的R Studio軟件實現過程。以某急性傳染病的流行作為暴露因素,將該急性傳染病流行前的人群作為非暴露組,將該急性傳染病流行后的人群作為暴露組。IRR為該急性傳染病流行后乳腺癌篩查率與該急性傳染病流行前乳腺癌篩查率的比值。合并后的IRR反映了該急性傳染病流行對乳腺癌篩查的影響程度。在比較兩組IR時,metainc()函數還可以計算合并發生率差異(Incidence Rate Difference, sm="IRD"),平方根變換的發生率差(Square root transformed Incidence Rate Difference, sm="IRSD"),疫苗效力或疫苗有效性(Vaccine efficacy or vaccine effectiveness, sm="VE")。除R Studio軟件外,也有研究使用Stata[5-6]和Commercial Meta-Analysis Software[7-8]進行以IRR為效應量的Meta分析。
除了篩查研究以外,IRR也常見于隊列研究[9-11]和隨機對照試驗[7,12-13]。但是,不同研究對IRR的定義可能不同[14-16]。例如,一項納入隊列研究的Meta分析[5]評估了炎癥性腸病患者發生帶狀皰疹的風險。暴露組的IR為炎癥性腸病患者發生帶狀皰疹的人數除以隨訪患者人年數。非暴露組的IR為非炎癥性腸病患者發生帶狀皰疹的人數除以隨訪患者人年數。IRR=暴露組的IR/非暴露組的IR。合并后的IRR反映了患有炎癥性腸病對發生帶狀皰疹的影響程度。Marilly等[17]開展了一項評估鈉-葡萄糖共轉運蛋白2抑制劑(SGLT2i)治療2型糖尿病風險和獲益的Meta分析。該研究納入了5個隨機對照試驗,以總體死亡率作為Meta分析的主要結局指標。試驗組的IR為試驗組死亡人數除以試驗組隨訪人年數。對照組的IR為對照組死亡人數除以對照組隨訪人年數。IRR=試驗組的IR/對照組的IR。合并后的IRR反映了SGLT2i與總體死亡率的關聯程度。
IRR、RR(risk ratio)和HR(hazard ratio)均是反映事件發生相對風險的指標[18-20]。RR的計算過程通常不涉及時間變量。當試驗/暴露組與對照/非暴露組的隨訪時間不同時,合并后的RR可能存在估計不精確的問題。HR的計算涉及終點事件和發生終點事件的時間。但是,一般需要通過Cox回歸估計HR[21]。當一些原始研究未提供具體的HR值時,可能導致這些研究的結果無法進行Meta分析。IR的分子是在觀察期內新發生終點事件的人數,分母是人群的總觀察人時,描述了終點事件在人群中的發生速度。常用于需要長期隨訪的隊列研究,特別是當不斷有研究對象進入、退出研究或失訪時,IR可以準確描述終點事件在人群中的發生概率[1]。因此,如果研究提供了事件發生數及篩查或隨訪時長,可以考慮計算和合并IRR。當終點事件重復發生時,在IRR的計算過程中,IR的分子可以是終點事件的發生頻次。例如:一項Meta分析[22]調查了運動員受傷率的性別差異。在隨訪期內,同一位運動員可能出現多次受傷的情況。IRR=(女運動員受傷次數/女運動員隨訪時間)/(男運動員受傷次數/男運動員隨訪時間)。但是,RR和HR不能用于終點事件重復發生的相對風險估計。
綜上所述,本文詳細介紹了使用R Studio軟件的meta程序包實現以IRR為效應量的Meta分析過程,可為開展同類研究提供重要的方法學參考。