在系統評價或 Meta 分析中，風險比（Hazard ratio，HR）常作為效應量用于時間-事件結局分析，但在實踐中，常因忽視數據的對數轉換從而導致錯誤結果。本文以實例說明如何正確使用 Stata 和 R 軟件對獲得 HR 及其 95% 可信區間的生存資料進行 Meta 分析。
引用本文： 張天嵩. Stata 和 R 軟件在生存資料 Meta 分析中的正確使用. 中國循證醫學雜志, 2017, 17(10): 1237-1240. doi: 10.7507/1672-2531.201607004 復制
對于生存資料的 Meta 分析，如能從原始研究獲得風險比（Hazard ratio，HR）及其 95% 可信區間（confidence interval，CI），則可運用 Stata 和 R 等軟件輕松地獲得總的合并 HR[1, 2]。但在實踐中，常因忽視 HR 及其 95%CI 的對數轉換，從而導致錯誤結果[3]。為了避免系統評價研究者錯誤使用相關軟件，本文以兩個具體數據為例，說明正確使用 Stata 和 R 軟件進行生存資料 Meta 分析的方法。
1 資料和方法
1.1 數據來源
數據來源于兩個 Meta 分析研究。數據一來自于文獻[3]所引用 Xie 等[4]的 Meta 分析，主要觀察手術對老年人結直腸癌肝轉移患者是否有效，直接從原文提取數據，見表 1。數據二來源于 Le Chevalier 等[5]的 Meta 分析，觀察吉西他濱聯合鉑類藥物治療晚期非小細胞癌的療效，選擇其中第三代含鉑類藥物亞組數據，見表 2。每個數據均含有四個變量：study、hr、ll、ul，分別表示研究名稱、HR 及其 95%CI 的下限和上限。按表 1、表 2 所示格式輸入 Stata 數據管理器中，分別命名為 xie.dta 和 le.dta（以 Stata 12.0 版本以下數據格式）存儲于 C 盤根目錄下備用。
表1
結直腸癌 Meta 分析數據
表2
非小細胞癌 Meta 分析數據
1.2 Meta 分析方法和模型
Cochrane 干預措施系統評價手冊明確指出，對于時間-事件（生存數據），如果從單個原始研究 Cox 比例風險回歸模型的結果中獲得了 HR 的對數值（log HR）及其標準誤（se log HR），則可使用倒方差法進行合并，固定效應模型和隨機效應模型均可選用[6]。同理，如果單個研究獲得 HR 及其 95%CI，因 HR 常不服從正態分布，需要進行對數轉換，可通過以下公式計算出 log HR 及其標準誤，然后使用倒方差法進行合并。假設納入 Meta 分析的第 i 個研究的效應量及 95%CI 上、下限分別為 HRi、UCIi、LCIi，則單個研究的效應量為 log HRi，其標準誤為 se log HRi，方差為 var log HRi，則總的合并效應量為[7]：
，其中
，
，若要獲得合并的 HR，則只需將合并 log HR 指數化即可。
1.3 軟件命令或軟件包
為便于與文獻[3]中的結果相對比，本文仍采用 Stata 軟件的 metan 命令、metabias 命令；以及 R 軟件 meta 擴展包的 metagen（）函數、metabias（）函數。關于命令或軟件包的安裝，眾多文獻及書籍均有說明，可以參考相關文獻[1, 3]，此處不再贅述。
1.4 Stata 和 R 軟件在生存資料 Meta 分析中的應用
分別使用 Stata（13.1 版）和 R 軟件（3.2.2 版）合并 HR，均采用隨機效應模型，并比較對 HR 及其 95%CI 進行對數轉換后合并 HR 結果的異同。因為數據一僅有 8 個研究，所以對數據一不進行繪制漏斗圖和漏斗圖不對稱檢驗[6]。
1.4.1 Stata 軟件在生存資料 Meta 分析中的應用 為方便計算，筆者編寫了.do 文件，具體操作步驟和主要解釋見框 1，請注意在 Stata 軟件中把“*”放在命令或語句前面用于解釋和說明。有興趣的讀者，可在 Stata 軟件的菜單欄中按“Windows→Do-file Editor→New Do-file Editor”步驟操作，打開“Do-file”編輯器，將框 1 中代碼復制到已打開的“Do-file”編輯器中，并在編輯器的“Tools”下拉菜單中選擇“Excute（do）”執行即可[8]。請注意，代碼中“metan hr ll ul，randomi nograph”是本文為比較結果而使用的，請讀者在實際 Meta 分析操作中切勿使用該命令。
1.4.2 R 軟件在生存資料 Meta 分析中的應用 具體操作步驟代碼和主要解釋見框 2，請注意在 R 軟件中把“#”放在命令或語句之前用于解釋說明。
2 結果
針對不同的數據，通過對 HR 及其 95%CI 對數轉換后合并 HR，或通過合并 HR 的對數及對數標準誤，采用隨機效應模型，Stata 軟件和 R 軟件合并獲得的結果完全一致，而且與原文結果完全相同。但是未進行對數轉換直接合并 HR 及其 95%CI 結果與對數轉換后合并結果不同，特別是數據二，對數據進行對數轉換后，結果發生改變，差異由無統計學意義（P=0.077）變為有統計學意義（P<0.001），見表 3。Egger 法漏斗圖不對稱檢驗 t=?1.61，相應 P=0.142。
表3
是否進行對數轉換對合并結果的影響
3 討論
隨訪研究獲得的資料稱為生存資料，常包括隨訪事件、生存時間、分組變量以及其他協變量。生存資料的 Meta 分析在預后研究的系統評價中應用廣泛，常用的合并指標有生存率、中位生存期和 HR 等，應盡可能采用 HR 作為效應量進行合并[1]，眾多統計軟件均可用于生存資料的 Meta 分析，如 Stata、R 等。
對于生存資料 Meta 分析，應盡量使用個體參與者數據；如果不能獲得，再考慮使用聚合數據[6]。若獲得 HR 及其 95%CI 等聚合數據，一般可通過兩種途徑進行合并：
第一，將 HR 及其 95%CI 上限和下限分別取自然對數，然后可以用 Stata 的 metan 命令后跟三變量進行合并，基本命令格式如“metan loghr logll logul，[選擇項]”。如果數據未進行對數轉換就直接合并[3]，則結果與原文[4]有所差異；而本文通過對 HR 及其 95%CI 上限和下限進行對數轉換后合并，結果則與原文完全一致。特別是數據二，對數據是否進行對數轉換，其合并結果完全不一致，結論也相應發生改變。請注意，如果采用 Stata 第一個用于 Meta 分析的 meta 命令，則不需要進行對數轉換，如針對本文中的數據，使用命令“meta hr ll ul”得到的結果與表 3 中“進行對數轉換后合并結果”相同。
第二，根據 HR 及其 95%CI 獲得 HR 的自然對數及對數標準誤，然后使用 Stata 軟件的 metan 命令后跟兩變量格式進行合并，基本命令格式如“metan loghr seloghr，[選擇項]”。也可以采用 R 軟件 meta 擴展包的 metagen（）函數進行合并，如本文和文獻[3]均采用 metagen（）函數對數據一進行分析，但文獻[3]中描述直接合并 HR 及其標準誤有誤，在 meta 擴展包自帶的說明文件中明確的示例是合并 log HR 和 se log HR，從文獻[3]的圖 8 也可以看出，森林圖右側顯示的“HR1”的 95%CI 與單個的原始研究結果不同，而本文采用的方法計算結果與原文完全一致。
需要提醒的是，metabias 命令由 Steichen 編寫，最初發布于 1997 年，已數次更新；在 2009 年由 Harbord 作了重要修訂，并將以前 Steichen 編寫的版本命名為 metabias6，兩者命令行操作格式不同，如新版本的命令只能后跟效應量及其標準誤、二分類數據格式，已不支持效應量及其 95%CI 上、下限三變量格式，在應用時請讀者注意，詳細內容參看相關文獻[1, 2, 9]。
總之，對于生存資料的 Meta 分析，如果能獲得 HR 及其 95%CI，則采用 Stata、R 等軟件計算合并 HR 較為簡單，但應注意每種軟件的命令或包的正確使用方法，以免獲得錯誤的合并結果。
對于生存資料的 Meta 分析，如能從原始研究獲得風險比（Hazard ratio，HR）及其 95% 可信區間（confidence interval，CI），則可運用 Stata 和 R 等軟件輕松地獲得總的合并 HR[1, 2]。但在實踐中，常因忽視 HR 及其 95%CI 的對數轉換，從而導致錯誤結果[3]。為了避免系統評價研究者錯誤使用相關軟件，本文以兩個具體數據為例，說明正確使用 Stata 和 R 軟件進行生存資料 Meta 分析的方法。
1 資料和方法
1.1 數據來源
數據來源于兩個 Meta 分析研究。數據一來自于文獻[3]所引用 Xie 等[4]的 Meta 分析，主要觀察手術對老年人結直腸癌肝轉移患者是否有效，直接從原文提取數據，見表 1。數據二來源于 Le Chevalier 等[5]的 Meta 分析，觀察吉西他濱聯合鉑類藥物治療晚期非小細胞癌的療效，選擇其中第三代含鉑類藥物亞組數據，見表 2。每個數據均含有四個變量：study、hr、ll、ul，分別表示研究名稱、HR 及其 95%CI 的下限和上限。按表 1、表 2 所示格式輸入 Stata 數據管理器中，分別命名為 xie.dta 和 le.dta（以 Stata 12.0 版本以下數據格式）存儲于 C 盤根目錄下備用。
表1
結直腸癌 Meta 分析數據
表2
非小細胞癌 Meta 分析數據
1.2 Meta 分析方法和模型
Cochrane 干預措施系統評價手冊明確指出，對于時間-事件（生存數據），如果從單個原始研究 Cox 比例風險回歸模型的結果中獲得了 HR 的對數值（log HR）及其標準誤（se log HR），則可使用倒方差法進行合并，固定效應模型和隨機效應模型均可選用[6]。同理，如果單個研究獲得 HR 及其 95%CI，因 HR 常不服從正態分布，需要進行對數轉換，可通過以下公式計算出 log HR 及其標準誤，然后使用倒方差法進行合并。假設納入 Meta 分析的第 i 個研究的效應量及 95%CI 上、下限分別為 HRi、UCIi、LCIi，則單個研究的效應量為 log HRi，其標準誤為 se log HRi，方差為 var log HRi，則總的合并效應量為[7]：
，其中
，
，若要獲得合并的 HR，則只需將合并 log HR 指數化即可。
1.3 軟件命令或軟件包
為便于與文獻[3]中的結果相對比，本文仍采用 Stata 軟件的 metan 命令、metabias 命令；以及 R 軟件 meta 擴展包的 metagen（）函數、metabias（）函數。關于命令或軟件包的安裝，眾多文獻及書籍均有說明，可以參考相關文獻[1, 3]，此處不再贅述。
1.4 Stata 和 R 軟件在生存資料 Meta 分析中的應用
分別使用 Stata（13.1 版）和 R 軟件（3.2.2 版）合并 HR，均采用隨機效應模型，并比較對 HR 及其 95%CI 進行對數轉換后合并 HR 結果的異同。因為數據一僅有 8 個研究，所以對數據一不進行繪制漏斗圖和漏斗圖不對稱檢驗[6]。
1.4.1 Stata 軟件在生存資料 Meta 分析中的應用 為方便計算，筆者編寫了.do 文件，具體操作步驟和主要解釋見框 1，請注意在 Stata 軟件中把“*”放在命令或語句前面用于解釋和說明。有興趣的讀者，可在 Stata 軟件的菜單欄中按“Windows→Do-file Editor→New Do-file Editor”步驟操作，打開“Do-file”編輯器，將框 1 中代碼復制到已打開的“Do-file”編輯器中，并在編輯器的“Tools”下拉菜單中選擇“Excute（do）”執行即可[8]。請注意，代碼中“metan hr ll ul，randomi nograph”是本文為比較結果而使用的，請讀者在實際 Meta 分析操作中切勿使用該命令。
1.4.2 R 軟件在生存資料 Meta 分析中的應用 具體操作步驟代碼和主要解釋見框 2，請注意在 R 軟件中把“#”放在命令或語句之前用于解釋說明。
2 結果
針對不同的數據，通過對 HR 及其 95%CI 對數轉換后合并 HR，或通過合并 HR 的對數及對數標準誤，采用隨機效應模型，Stata 軟件和 R 軟件合并獲得的結果完全一致，而且與原文結果完全相同。但是未進行對數轉換直接合并 HR 及其 95%CI 結果與對數轉換后合并結果不同，特別是數據二，對數據進行對數轉換后，結果發生改變，差異由無統計學意義（P=0.077）變為有統計學意義（P<0.001），見表 3。Egger 法漏斗圖不對稱檢驗 t=?1.61，相應 P=0.142。
表3
是否進行對數轉換對合并結果的影響
3 討論
隨訪研究獲得的資料稱為生存資料，常包括隨訪事件、生存時間、分組變量以及其他協變量。生存資料的 Meta 分析在預后研究的系統評價中應用廣泛，常用的合并指標有生存率、中位生存期和 HR 等，應盡可能采用 HR 作為效應量進行合并[1]，眾多統計軟件均可用于生存資料的 Meta 分析，如 Stata、R 等。
對于生存資料 Meta 分析，應盡量使用個體參與者數據；如果不能獲得，再考慮使用聚合數據[6]。若獲得 HR 及其 95%CI 等聚合數據，一般可通過兩種途徑進行合并：
第一，將 HR 及其 95%CI 上限和下限分別取自然對數，然后可以用 Stata 的 metan 命令后跟三變量進行合并，基本命令格式如“metan loghr logll logul，[選擇項]”。如果數據未進行對數轉換就直接合并[3]，則結果與原文[4]有所差異；而本文通過對 HR 及其 95%CI 上限和下限進行對數轉換后合并，結果則與原文完全一致。特別是數據二，對數據是否進行對數轉換，其合并結果完全不一致，結論也相應發生改變。請注意，如果采用 Stata 第一個用于 Meta 分析的 meta 命令，則不需要進行對數轉換，如針對本文中的數據，使用命令“meta hr ll ul”得到的結果與表 3 中“進行對數轉換后合并結果”相同。
第二，根據 HR 及其 95%CI 獲得 HR 的自然對數及對數標準誤，然后使用 Stata 軟件的 metan 命令后跟兩變量格式進行合并，基本命令格式如“metan loghr seloghr，[選擇項]”。也可以采用 R 軟件 meta 擴展包的 metagen（）函數進行合并，如本文和文獻[3]均采用 metagen（）函數對數據一進行分析，但文獻[3]中描述直接合并 HR 及其標準誤有誤，在 meta 擴展包自帶的說明文件中明確的示例是合并 log HR 和 se log HR，從文獻[3]的圖 8 也可以看出，森林圖右側顯示的“HR1”的 95%CI 與單個的原始研究結果不同，而本文采用的方法計算結果與原文完全一致。
需要提醒的是，metabias 命令由 Steichen 編寫，最初發布于 1997 年，已數次更新；在 2009 年由 Harbord 作了重要修訂，并將以前 Steichen 編寫的版本命名為 metabias6，兩者命令行操作格式不同，如新版本的命令只能后跟效應量及其標準誤、二分類數據格式，已不支持效應量及其 95%CI 上、下限三變量格式，在應用時請讀者注意，詳細內容參看相關文獻[1, 2, 9]。
總之，對于生存資料的 Meta 分析，如果能獲得 HR 及其 95%CI，則采用 Stata、R 等軟件計算合并 HR 較為簡單，但應注意每種軟件的命令或包的正確使用方法，以免獲得錯誤的合并結果。
