nlme程序包是基于廣義最小二乘法和線性混合效應模型研發的、可通過R軟件實現廣義線性和非線性混合效應模型下的Meta分析。該程序包實現Meta分析時，需要對數據先行轉化為效應量的對數值才可進行。本文介紹了使用R軟件nlme程序包實現網狀Meta分析的過程，詳細呈現了如何轉化數據這一核心步驟。
引用本文： 張超, 牛玉明, 曾憲濤. R軟件nlme程序包在網狀Meta分析中的應用. 中國循證醫學雜志, 2014, 14(3): 355-360. doi: 10.7507/1672-2531.20140060 復制
隨著網狀Meta分析的發展及各種軟件的成功研發，運算模型的選擇也日益受到關注，廣義線性模型（generalized linear model，GLM）就是其中之一。GLM是線性模型的擴展，其特點是不強行改變數據的自然度量，數據可以具有非線性和非恒定方差結構；與線性模型相比，GLM模型中Y的分布可以是任何形式的指數分布（如高斯分布、泊松分布、二項式分布），聯結函數可以是任何單調可微函數（如對數函數logarithm或邏輯函數logit），這些優點使得GLM模型可處理多種變量，日益受到關注（欲詳細了解GLM的讀者建議閱讀參考文獻[1, 2]）。
nlme是一款基于S語言在S-PLUS軟件中實現混合效應模型（mixed-effects models）分析的軟件。因此，其可以同時進行線性和非線性混合效應模型（linear and nonlinear mixed effects models，NLME）分析[3, 4]。當前，混合效應模型也可通過nlme程序包在R語言中實現[5, 6]。nlme程序包亦可以實現網狀Meta分析[6, 7]，本文仍以《R軟件R2WinBUGS程序包在網狀Meta分析中的應用》一文中實例[8, 9]為例進行展示。
1 軟件及程序包的安裝/加載
R軟件的安裝已在前文講述，本文使用的是R-3.0.1 [8, 10, 11]。
此外，需要安裝nlme程序包，具體命令為：install.packages(“nlme”)。在彈出的對話框中選擇某個鏡像（CRAN）安裝，安裝完成后再由library(“nlme”)命令完成加載。
2 數據的加載與預處理
2.1 數據加載
R軟件nlme程序包實現網狀Meta分析的基本思路是將因變量y設置為效應量，自變量x設置為干預措施，結合相關的干預措施之間的比較來計算得出線性關系[7]。本文將y設為logOR，x設為13種藥物以及安慰劑，并對相關藥物之間的比較作一定標記。
因實例中含有3臂試驗，所以，首先將數據中所有的3臂試驗拆分成2臂試驗數據進行排列（見表 1）。拆分之后，使用如下命令加載數據：
?
linedata < -read.table(" D:\\Users\\Administrator\\Desktop\\Rwork\\linedata.txt" ,header=TRUE)。
?
表1
本文示例及預處理之后的數據
2.2 數據轉化
在進行線性模擬之前，需用到metafor程序包[10]來計算logOR及方差。數據轉化可分為三步。
第一步，使用install.packages(“metafor”)和library(“metafor”)命令安裝和加載Metafor程序包[10]。
第二步，計算因變量（即效應量）logOR。因實例為二分類數據，metafor程序包使用的是escalc命令：
?
metadata < -escalc(ai=r1, n1i=n1, ci=r2, n2i=n2, measure=“OR”, data=linedata)。
?
命令中：r1、n1、r2和n2分別表示試驗組的事件發生數與總數及對照組的事件發生數與總數；measure表示得出的效應量的對數［如本例選“OR”則結果自動轉換為log(OR)，但其結果是以變量名yi來儲存log(OR)數據集，而其方差集以變量名vi表示］。需要注意的是，得出結果數據集將直接寫入到原始數據中。
第三步，轉換變量名稱。為方便記憶，本例將Meta數據集中的yi及vi變量名稱改為logOR及var，具體代碼為：
?
names(metadata)[8:9] < -c(“logOR”, “var”)。
?
其中，[8:9]代表第8行和第9行對應的數據集。
2.3 創建啞變量矩陣
本例我們需要將88個研究組中的14個自變量設定成88×14的啞變量矩陣，具體命令如下：
?
matrixdata < -matrix(0, nrow=88, ncol=14, byrow=TRUE, dimnames=list(1:88, c(“A”, “B”, “C”, “D”, “E”, “F”, “G”, “H”, “I”, “J”, “K”, “L”, “M”, “N”)))
?
其中，0為矩陣中的數值（在此全部設置為0）；nrow為矩陣行數；ncol為矩陣列數；byrow為依據行來對矩陣進行賦值；dimnames為對矩陣的行和列進行命名；1:88為行的命名（從1至88），在此代表 88個研究組；c(“A”, “B”, “C”, “D”, “E”, “F”, “G”, “H”, “I”, “J”, “K”, “L”, “M”, “N”)為列的命名（從A到N），在此代表 13種藥物及安慰劑。
接著，設定配對變量矩陣，命令代碼如下：
?
for(aa in 1:88)
{ for(bb in 6:7)
??{
???matrixdata[aa, as.numeric(charToRaw(as.character(metadata$t1[aa])))-64] < -1
???matrixdata[aa, as.numeric(charToRaw(as.character(metadata$t2[aa])))-64] < --1
??}
}
?
上述代碼作用在于將兩種治療藥物相互配對，配對的方式為將metadata試驗組使用的藥物轉換成1放置于matrixdata矩陣所對應的藥物，相應的對照組設為-1。其中，as.character()和as.numeric()是指將括號里面的元素分別轉換成字符型和數字型，charToRaw()為將字符轉換成相應的ACSⅡ碼。
2.4 標識配對藥物
接著，需將metadata數據中的比較藥物做相關的配對標識，以便后面lme()函數的使用。實現的命令為：
?
mark < -paste(as.character(metadata$t1), as.character(metadata$t2))。
?
2.5 預處理數據匯總
在進行線性模型數據處理之前，我們需要搜集相關的預處理數據信息，具體命令為：totaldata < -cbind(metadata, matrixdata, mark)。
3 完成網狀Meta分析
在計算之前，可展示totaldata數據集（見表 1），命令為：totaldata。表 1中，變量logOR、A～N及mark為最終nlme程序包中lme()函數所需的有效數據。
接著，進行線性混合模型運算，命令如下：
?
lme < -lme(logOR~0+B+C+D+E+F+G+H+I+J+K+L+M+N, random=~1|mark, data=totaldata)。
?
命令中：logOR~0+B+C+D+E+F+G+H+I+J+K+L+M+N成為固定部分，其“～”左邊代表因變量（效應量），“～”右邊代表自變量（13種藥物及安慰劑）；在線性模型中，對于存在n個自變量時，通常需要設置n-1個啞變量，另一個則需要用數字“0”替代，此處，將“A”替換為“0”來進行運算。mark為配對比較藥物的標記，同時作為random～1的隨機效應范圍。lindata為儲存數據的變量名。
由于各種不同數據精準計算的需要，有時還需對以下3個方面作相應的設置，本處均為默認設置如下：①算法（method）默認：限制性極大對數似然法（REML）；②權重（weight）默認：NULL，組間同方差（值得注意的是，通常在設置的時候，采取方差的協變量）；③相關性（correlation）默認：NULL，組間不存在自相關性。
上述代碼運行之后，再運行匯總結果命令。命令為：summary(lme)。
命令執行后的結果見框1。
上述結果中，是以“A”藥物為基線藥物作參考；當然，還可以使用數字“0”依次替換“B”至“N”，替換后，即可得到相應的網狀計算結果。有興趣的讀者可自行嘗試。
4 繪制圖形
4.1 網狀關系圖繪制
同R2WinBUGS程序包一樣，nlme程序包當前尚無繪制網狀關系圖的功能，需要通過network程序包實現。有關使用此程序包及其他相關方法繪制網狀Meta分析關系圖的方法，詳見《網狀Meta分析中網狀關系圖的繪制》一文[12]。
4.2 相關圖形繪制
nlme程序包可繪制的圖形及命令分別如下（本文對這5種圖形不一一展示）：
（1）截距匯總圖：plot(ranef(lme))
（2）殘差分布圖：plot(lme)
（3）按mark分組，標準化殘差與擬合值關系圖：plot(lme, resid(., type=“p”)~fitted(.)|mark, abline=0)
（4）按mark分組，殘差的箱型圖：plot(lme, mark~resid(.))
（5）按mark分組，效應量與擬合值關系圖：plot(lme, logor~fitted(.)|mark, abline=c(0, 1))
5 結語
實際上，GLM在網狀Meta分析中的運用，雖較以前有所突破，但其在考慮間接比較方面仍較為欠缺，特別是在圖形的繪制上；此外，還需加入異質性及一致性的相關參數。再者，進行分析時，需先行對數據進行轉化才可繼續進行。而在這些方面，貝葉斯的線性混合模型較傳統線性混合模型具有更大優勢。
盡管如此，nlme程序包為網狀Meta分析提供了一種新的實現思路與方法。當前網狀Meta分析仍有許多需要研究與完善的地方，如傳遞性、異質性和一致性的檢測。nlme程序包是基于傳統線性模型，而非貝葉斯模型。該程序包在網狀Meta分析實現中的最大缺陷在于未考慮一致性、異質性及協方差參數的設定，亦不具備檢測功能；且有研究顯示使用非線性模型實現網狀Meta分析，方法學仍在突破。因此，介紹nlme程序包是必要的，該程序包為從貝葉斯模型之外的途徑完善網狀Meta分析提供新的可能途徑。隨著今后nlme程序包的不斷完善，其在網狀Meta分析中的運用將越發廣泛。
隨著網狀Meta分析的發展及各種軟件的成功研發，運算模型的選擇也日益受到關注，廣義線性模型（generalized linear model，GLM）就是其中之一。GLM是線性模型的擴展，其特點是不強行改變數據的自然度量，數據可以具有非線性和非恒定方差結構；與線性模型相比，GLM模型中Y的分布可以是任何形式的指數分布（如高斯分布、泊松分布、二項式分布），聯結函數可以是任何單調可微函數（如對數函數logarithm或邏輯函數logit），這些優點使得GLM模型可處理多種變量，日益受到關注（欲詳細了解GLM的讀者建議閱讀參考文獻[1, 2]）。
nlme是一款基于S語言在S-PLUS軟件中實現混合效應模型（mixed-effects models）分析的軟件。因此，其可以同時進行線性和非線性混合效應模型（linear and nonlinear mixed effects models，NLME）分析[3, 4]。當前，混合效應模型也可通過nlme程序包在R語言中實現[5, 6]。nlme程序包亦可以實現網狀Meta分析[6, 7]，本文仍以《R軟件R2WinBUGS程序包在網狀Meta分析中的應用》一文中實例[8, 9]為例進行展示。
1 軟件及程序包的安裝/加載
R軟件的安裝已在前文講述，本文使用的是R-3.0.1 [8, 10, 11]。
此外，需要安裝nlme程序包，具體命令為：install.packages(“nlme”)。在彈出的對話框中選擇某個鏡像（CRAN）安裝，安裝完成后再由library(“nlme”)命令完成加載。
2 數據的加載與預處理
2.1 數據加載
R軟件nlme程序包實現網狀Meta分析的基本思路是將因變量y設置為效應量，自變量x設置為干預措施，結合相關的干預措施之間的比較來計算得出線性關系[7]。本文將y設為logOR，x設為13種藥物以及安慰劑，并對相關藥物之間的比較作一定標記。
因實例中含有3臂試驗，所以，首先將數據中所有的3臂試驗拆分成2臂試驗數據進行排列（見表 1）。拆分之后，使用如下命令加載數據：
?
linedata < -read.table(" D:\\Users\\Administrator\\Desktop\\Rwork\\linedata.txt" ,header=TRUE)。
?
表1
本文示例及預處理之后的數據
2.2 數據轉化
在進行線性模擬之前，需用到metafor程序包[10]來計算logOR及方差。數據轉化可分為三步。
第一步，使用install.packages(“metafor”)和library(“metafor”)命令安裝和加載Metafor程序包[10]。
第二步，計算因變量（即效應量）logOR。因實例為二分類數據，metafor程序包使用的是escalc命令：
?
metadata < -escalc(ai=r1, n1i=n1, ci=r2, n2i=n2, measure=“OR”, data=linedata)。
?
命令中：r1、n1、r2和n2分別表示試驗組的事件發生數與總數及對照組的事件發生數與總數；measure表示得出的效應量的對數［如本例選“OR”則結果自動轉換為log(OR)，但其結果是以變量名yi來儲存log(OR)數據集，而其方差集以變量名vi表示］。需要注意的是，得出結果數據集將直接寫入到原始數據中。
第三步，轉換變量名稱。為方便記憶，本例將Meta數據集中的yi及vi變量名稱改為logOR及var，具體代碼為：
?
names(metadata)[8:9] < -c(“logOR”, “var”)。
?
其中，[8:9]代表第8行和第9行對應的數據集。
2.3 創建啞變量矩陣
本例我們需要將88個研究組中的14個自變量設定成88×14的啞變量矩陣，具體命令如下：
?
matrixdata < -matrix(0, nrow=88, ncol=14, byrow=TRUE, dimnames=list(1:88, c(“A”, “B”, “C”, “D”, “E”, “F”, “G”, “H”, “I”, “J”, “K”, “L”, “M”, “N”)))
?
其中，0為矩陣中的數值（在此全部設置為0）；nrow為矩陣行數；ncol為矩陣列數；byrow為依據行來對矩陣進行賦值；dimnames為對矩陣的行和列進行命名；1:88為行的命名（從1至88），在此代表 88個研究組；c(“A”, “B”, “C”, “D”, “E”, “F”, “G”, “H”, “I”, “J”, “K”, “L”, “M”, “N”)為列的命名（從A到N），在此代表 13種藥物及安慰劑。
接著，設定配對變量矩陣，命令代碼如下：
?
for(aa in 1:88)
{ for(bb in 6:7)
??{
???matrixdata[aa, as.numeric(charToRaw(as.character(metadata$t1[aa])))-64] < -1
???matrixdata[aa, as.numeric(charToRaw(as.character(metadata$t2[aa])))-64] < --1
??}
}
?
上述代碼作用在于將兩種治療藥物相互配對，配對的方式為將metadata試驗組使用的藥物轉換成1放置于matrixdata矩陣所對應的藥物，相應的對照組設為-1。其中，as.character()和as.numeric()是指將括號里面的元素分別轉換成字符型和數字型，charToRaw()為將字符轉換成相應的ACSⅡ碼。
2.4 標識配對藥物
接著，需將metadata數據中的比較藥物做相關的配對標識，以便后面lme()函數的使用。實現的命令為：
?
mark < -paste(as.character(metadata$t1), as.character(metadata$t2))。
?
2.5 預處理數據匯總
在進行線性模型數據處理之前，我們需要搜集相關的預處理數據信息，具體命令為：totaldata < -cbind(metadata, matrixdata, mark)。
3 完成網狀Meta分析
在計算之前，可展示totaldata數據集（見表 1），命令為：totaldata。表 1中，變量logOR、A～N及mark為最終nlme程序包中lme()函數所需的有效數據。
接著，進行線性混合模型運算，命令如下：
?
lme < -lme(logOR~0+B+C+D+E+F+G+H+I+J+K+L+M+N, random=~1|mark, data=totaldata)。
?
命令中：logOR~0+B+C+D+E+F+G+H+I+J+K+L+M+N成為固定部分，其“～”左邊代表因變量（效應量），“～”右邊代表自變量（13種藥物及安慰劑）；在線性模型中，對于存在n個自變量時，通常需要設置n-1個啞變量，另一個則需要用數字“0”替代，此處，將“A”替換為“0”來進行運算。mark為配對比較藥物的標記，同時作為random～1的隨機效應范圍。lindata為儲存數據的變量名。
由于各種不同數據精準計算的需要，有時還需對以下3個方面作相應的設置，本處均為默認設置如下：①算法（method）默認：限制性極大對數似然法（REML）；②權重（weight）默認：NULL，組間同方差（值得注意的是，通常在設置的時候，采取方差的協變量）；③相關性（correlation）默認：NULL，組間不存在自相關性。
上述代碼運行之后，再運行匯總結果命令。命令為：summary(lme)。
命令執行后的結果見框1。
上述結果中，是以“A”藥物為基線藥物作參考；當然，還可以使用數字“0”依次替換“B”至“N”，替換后，即可得到相應的網狀計算結果。有興趣的讀者可自行嘗試。
4 繪制圖形
4.1 網狀關系圖繪制
同R2WinBUGS程序包一樣，nlme程序包當前尚無繪制網狀關系圖的功能，需要通過network程序包實現。有關使用此程序包及其他相關方法繪制網狀Meta分析關系圖的方法，詳見《網狀Meta分析中網狀關系圖的繪制》一文[12]。
4.2 相關圖形繪制
nlme程序包可繪制的圖形及命令分別如下（本文對這5種圖形不一一展示）：
（1）截距匯總圖：plot(ranef(lme))
（2）殘差分布圖：plot(lme)
（3）按mark分組，標準化殘差與擬合值關系圖：plot(lme, resid(., type=“p”)~fitted(.)|mark, abline=0)
（4）按mark分組，殘差的箱型圖：plot(lme, mark~resid(.))
（5）按mark分組，效應量與擬合值關系圖：plot(lme, logor~fitted(.)|mark, abline=c(0, 1))
5 結語
實際上，GLM在網狀Meta分析中的運用，雖較以前有所突破，但其在考慮間接比較方面仍較為欠缺，特別是在圖形的繪制上；此外，還需加入異質性及一致性的相關參數。再者，進行分析時，需先行對數據進行轉化才可繼續進行。而在這些方面，貝葉斯的線性混合模型較傳統線性混合模型具有更大優勢。
盡管如此，nlme程序包為網狀Meta分析提供了一種新的實現思路與方法。當前網狀Meta分析仍有許多需要研究與完善的地方，如傳遞性、異質性和一致性的檢測。nlme程序包是基于傳統線性模型，而非貝葉斯模型。該程序包在網狀Meta分析實現中的最大缺陷在于未考慮一致性、異質性及協方差參數的設定，亦不具備檢測功能；且有研究顯示使用非線性模型實現網狀Meta分析，方法學仍在突破。因此，介紹nlme程序包是必要的，該程序包為從貝葉斯模型之外的途徑完善網狀Meta分析提供新的可能途徑。隨著今后nlme程序包的不斷完善，其在網狀Meta分析中的運用將越發廣泛。
