遺傳模型的選擇在基因-疾病關聯研究 Meta 分析中是較為關鍵的難點。本文對 Minelli 等提出的無遺傳模型約束的貝葉斯 Meta 分析方法進行深入討論，利用 JAGS 和 R 語言對其進行實現，并在實例中比較了不同參數先驗信息對合并效應量的影響，特別是外部臨床先驗的影響。實例研究結果顯示，無論是無信息先驗還是外部臨床先驗，在納入研究個數較多的情況下，對合并效應量的影響均較小。
引用本文： 宋嘉麒, 金志超. 無遺傳模型約束的基因-疾病關聯研究的貝葉斯 Meta 分析方法及應用. 中國循證醫學雜志, 2017, 17(9): 1106-1110. doi: 10.7507/1672-2531.201703116 復制
在對慢性復雜疾病進行病因學的研究過程中，除考慮常見的環境危險因素外，往往需要考慮遺傳因素對疾病發病的影響。目前基因-疾病關聯研究常利用單核苷酸多態性（single nucleotide polymorphism，SNP）篩查與疾病相關的基因位點。由于測序技術的進步和成本的下降，基因-疾病關聯研究已大量開展，其文獻數量也呈爆發式增長[1]。由于單個基因多態性對于疾病的效應值往往較小[2]，因此需要更大樣本量來提高統計檢驗效能。因此，利用 Meta 分析方法整合各個獨立的基因-疾病關聯研究就顯得非常必要，目前已有研究者開展了此類研究的 Meta 分析方法學及應用研究[3, 4]。
在進行基于單核苷酸的基因-疾病關聯研究 Meta 分析中，比較突出的問題是如何選擇遺傳模型。研究者常假定一種遺傳模型（顯性或隱性遺傳模型）進行效應量的合并，但是這種方法缺乏生物學支持且較為武斷，如遺傳模型選用錯誤將導致基于基因-疾病關聯研究的 Meta 分析結論錯誤。Minelli 等[5, 6]提出通過計算不同基因型 OR 值的比值，并在此基礎上計算某基因型相對于參照基因型的 OR 值來避免選擇遺傳模型。Minelli 等[5, 6]分別利用極大似然估計和貝葉斯方法進行了估計，貝葉斯方法由于能夠利用外部信息作為先驗，估計結果更為精確，但這些研究未就外部先驗信息進行分析，只進行了多種模糊先驗的敏感性分析。
本研究在 Minelli 等研究的基礎上，利用各類模糊先驗、外部先驗信息，進行基因-疾病關聯研究 Meta 分析，利用 R 和 JAGS 軟件編寫 Meta 分析程序，采用 ACE I/D 基因多態性對于血管緊張素轉化酶抑制劑類降壓藥所致的咳嗽是否有影響的實例進行說明。
1 方法
1.1 基本原理
基因-疾病關聯研究 Meta 分析的特點之一是至少有 3 種可能的基因型分組，并且分析結果與所選擇的遺傳模型有關。在最簡單的雙等位基因情況下（野生基因型 A，突變基因型 G），一般通過收集 3 種不同基因型（AA，GA，GG）人群的目標疾病發生頻率來估計基因-疾病關聯性大小。因此，將會得到 2 個 OR 值，分別為 GG vs. AA 和 GA vs. AA，這 2 個 OR 值在遺傳模型下相互獨立。傳統的基因-疾病關聯研究 Meta 分析方法通過假設遺傳模型，將 3 組（AA，GA，GG）減少為 2 組，如：假設為隱性模型，則合并 GA 與 AA（GG vs. GA+AA）；假設為顯性模型，則合并 GG 與 GA（GG+GA vs. AA）；假設為超顯性模型，則合并 GG 與 AA（GG+AA vs. GA）等。
然而，遺傳模型通常是未知的，而錯誤的遺傳模型得出的 Meta 分析結果通常有誤。本研究提供的方法可避免選擇特定的遺傳模型，但需要計算 2 個 OR 值對數的比 λ：
 |
假設 λ 在研究間是固定值，即 λ 取值 0、0.5、1 時則分別對應隱性模型、共顯性模型和顯性模型。但實際上 λ 可取 0 到 1 間的任何值，因此不受特定遺傳模型限制。
基因-疾病關聯性研究的設計主要是病例-對照研究，本文利用回顧性似然構建相應的似然函數，同時假定病例組和對照組的基因型服從多項式分布（multinomial distribution）：
y0i~Multinomial (n0, p0i), y1i~Multinomial (n1, p1i)
其中，i=1, 2, 3，分別對應 AA，GA，GG 三種基因型，y0i 為對照組某基因型頻率，y1i 為病例組基因型頻率，n0 為對照組總樣本例數，n1 為病例組總樣本例數。
 |
其中，d 取值為 1 表示病例，取值為 0 表示對照，即 p0i 和 p1i 分別為某基因型在對照組和病例組中的概率。可將對照組進行以下的參數化過程：α1=1，α2=p02/p01，α3=p03/p01，進一步可構建 p01=1/（1+α2+α3），p02=α2/（1+α2+α3），p03=α3/（1+α2+α3）。同理，病例組參數化過程如下：β1=1，β2=p12/p11，β3=p13/p11。上述過程可得，GA vs. AA 的 OR 值為 OR2=β2/α2，GG vs. AA 的 OR 值為 OR3=β3/α3，其對數形式為 δ2=log（OR2），δ3=log（OR3），則 λ=δ2/δ3。如前所述，可以根據 λ 的值決定相應的遺傳模型（取值 0、0.5 或 1），或忽略遺傳模型直接進行效應量的合并（取值 0 到 1 間任何值）。病例組中各基因型概率可通過下式構建：
 |
則每個研究的似然函數可寫為：
 |
其中，假設研究間互相獨立，且第 j 個研究 GG vs. AA 的 OR 值的對數 δ3j 服從均數為 θ，方差為 τ2 的隨機效應正態分布：
 |
GA vs. AA 的 OR 值的對數 δ2j 可由 δ3j·λ 計算得出。因此，需要明確未知參數 θ、τ 和 λ 的先驗分布。
其中，均數 θ 在所有模型中服從正態分布：
 |
對于研究間的標準差 τ，本研究采用 3 種先驗分布：
分布 ①
分布 ②
分布 ③
其中，分布 ① 中，方差的倒數服從 Gamma 分布，雖然這是最常用的模糊先驗，但最近有學者推薦直接采用標準差的先驗[7]。分布 ② 中，研究間標準差 τ 服從標準半正態分布，這種先驗分布使得 τ>2 的概率很低。分布 ③ 為，0 到 2 范圍內的均勻分布，這種先驗分布排除了 τ>2 的可能性。
對于參數 λ，本研究采用 4 種先驗分布：
分布 ①
分布 ②
分布 ③
分布 ④
其中，分布 ① 為 2 個參數都為 1 的 Beta 分布，這種分布在 0 到 1 間是均勻的。但是，當參數取值接近 0 或 1 并且數據分布稀疏時，這種先驗分布傾向于將后驗估計值取為 0.5。例如，當近似于隱性模型時（λ 的真實值約等于 0），這種先驗分布更容易扭曲后驗估計值，因為有 90% 的概率得到大于 0.1 的值。分布 ② 為 2 個參數都為 0.5 的 Beta 分布，這種先驗分布在 λ 取接近 0 或 1 的值時，能有更大的先驗概率。但當模型近似于共顯性模型（λ=0.5）且數據分布稀疏時，這種分布可能會增加 λ 附近的不確定性。分布 ③ 為 0 到 1 范圍內的均勻分布，是一種常用的先驗分布。分布 ④ 為根據 Fujisawa 等[8]及 Minelli 等[5]的研究確定的外部先驗為均數 0.98，方差為 0.1 的正態分布。
1.2 實例
實例數據來自本文作者先前發表的 Cochrane 系統評價[9]，主要從基因多態性角度研究血管緊張素轉化酶抑制劑（ACEI）類降壓藥與咳嗽的關聯，共納入 17 個研究，其基本信息如表 1 所示。ACE 基因是高度多態的，其 16 內含子存在 287 bp 的 DNA 序列存在（插入，I）或沒有（缺失，D）的基因多態性，即研究最多的 ACE 基因插入/缺失（I/D）多態性。研究發現，D 等位基因與血漿 ACE 的水平呈正相關，DD 基因型者平均 ACE 水平大約是 II 基因型者的 2 倍。
表1
ACE I/D 基因多態性對血管緊張素轉化酶抑制劑類降壓藥所致咳嗽研究的基本特征
1.3 軟件實現
本文采用 R 和 JAGS 軟件進行編程，JAGS 軟件可在 http://mcmc-jags.sourceforge.net/ 下載，貝葉斯部分的 JAGS 代碼如下，以 txt 格式存儲在文件夾中即可：
model
{
for(i in 1:17 ) {
p0[i, 1]<-1/(1+a[i, 1]+a[i, 2]) #a 為公式中的 α
p0[i, 2]<-a[i, 1]/(1+a[i, 1]+a[i, 2])
p0[i, 3]<-a[i, 2]/(1+a[i, 1]+a[i, 2])
p1[i, 1]<-1/(1+a[i, 1]*exp(lambda*d[i])+a[i, 2]*exp(d[i])) #d 為 δ3
p1[i, 2]<-a[i, 1]*exp(lambda*d[i])/(1+a[i, 1]*exp(lambda*d[i])+a[i, 2]*exp(d[i]))
p1[i, 3]<-a[i, 2]*exp(d[i])/(1+a[i, 1]*exp(lambda*d[i])+a[i, 2]*exp(d[i]))
d[i]~dnorm(theta, prec)
ncont[i, 1:3]~dmulti(p0[i, 1:3], tcont[i])
ncase[i, 1:3]~dmulti(p1[i, 1:3], tcase[i])
a[i, 1]~dnorm(0, 0.0001)
a[i, 2]~dnorm(0, 0.0001)
}
lambda~dbeta(1, 1) #Beta 分布
#lambda~dbeta(0.5, 0.5)T(0.01, 0.99) #Beta 分布
#lambda~dunif(0, 1) #均勻分布
#lambda~dnorm(0.98, 0.1) #外部先驗
theta~dnorm(0, 0.000001)
tau~dgamma(0.001, 0.001) #Gamma 分布
#tau2~dunif(0, 1) #均勻分布
#tau2~dnorm(0, 1)I(0, ) #半正態分布
prec=1/(tau*tau)
or3<-exp(theta)
or2<-exp(theta*lambda)
}
在 R 中利用 R2jags 包調用 JAGS 程序，R 程序及相應的初始值如下：
library(R2jags)
library(dplyr)
ACEI <-read_excel("C:/work/modelfree/ACEI.xls") #讀取數據
attach(ACEI)
tcase=case_II+case_ID+case_DD
tcont=control_II+control_ID+control_DD
n<-17
ncase<-select(ACEI, contains("case"))
ncont<-select(ACEI, contains("control"))
data <-list("ncase", "ncont", "tcase", "tcont") #構建分析數據集
params <-c("lambda", "theta", "tau", "a", "or2", "or3") #構建參數集
#初始值設置，使用兩條鏈
a=matrix(c(rep(0.5, time=34)), ncol=2)
d<-c(rep(0, time=17))
inits1 <-list("lambda"=0.98, "theta"=0.5, "tau"=0.5, "a"=a, "d"=d)
inits2 <-list("lambda"=0.50, "theta"=1, "tau"=0.5, "a"=a, "d"=d)
inits <-list(inits1, inits2) #模型擬合
set.seed(20130801)
fit <-jags(data=data, inits=inits, parameters.to.save=params, n.chains=2, n.iter=50000, n.burnin=5000, model.file="C:/work/modelfree/ modelfree.txt")
print(fit)
plot(fit)
traceplot(fit)
本部分代碼主要涉及分析數據列表 data 的建立和初始值的設定，利用 R2jags 包中的 jags 函數進行抽樣和迭代，并利用 traceplot 函數進行收斂性評價，dplyr 包主要用作數據整理。在初始值的設置中，如果 λ 先驗信息為 dbeta（0.5，0.5）T（0.01，0.99），由于設置了抽樣的最大值和最小值，這時 λ 的初始值不能大于 0.99 或小于 0.01。
2 結果
不同的先驗信息下，ACE I/D 基因多態性對 ACEI 類降壓藥所致的咳嗽影響如表 2 所示。
表2
ACE I/D 基因多態性對 ACEI 類降壓藥所致咳嗽在不同先驗信息下的效應量合并結果
從表 2 中可以得出，相較于基因型Ⅱ，基因型 DD 和 DI 的攜帶者服用 ACEI 后引起干咳的可能性較小，與原文中采用顯性遺傳模型的結果相同[OR=0.70，95%CI（0.50，0.97）]。該結果在不同的先驗分布下非常穩定，但是研究間的異質性較大，可能引起的原因為納入研究的國家及人種不同引起，進一步亞組分析能夠減少研究間異質性。
3 討論
基因-疾病關聯 Meta 分析的無遺傳模型約束方法，同時考慮了遺傳效應的大小和遺傳模型的信息。雖然可用最大似然估計來實現，但貝葉斯方法是一種更具理論和實踐優勢的替代方法[5]。從理論角度看，貝葉斯方法能夠給予遺傳效應和遺傳模型明確的先驗信息。先驗信息可能來自于未被 Meta 分析納入的研究或專家意見，前者較為簡便而后者較難操作[10, 11]。從實際應用的角度看，采用馬爾可夫鏈蒙特卡羅（MCMC）算法估計的貝葉斯模型，更符合復雜遺傳模型下的基因疾病關聯。
Minelli 等[6]指出基因-疾病關聯 Meta 分析需要多個大樣本量的研究，使得 λ 先驗分布的不同選擇不會對結果造成過多的影響。同理，研究間方差 τ 的先驗分布也會影響結果。本研究結果顯示在研究數和樣本量較大的情況下，幾種無信息先驗和外部先驗對于結果的影響較小，也證實了 Minelli 等[6]的結論。另外，有研究指出，Gamma 分布并不是一個較好的研究間方差的先驗分布，而半正態分布和均勻分布可能更好[7]。從本文結果來看，3 種先驗分布下的結果差異不大，問題可能出現在數據分布稀疏或樣本量較小的情況下。因此，當納入研究較少的時，先驗信息對合并效應量的影響可能會較大，故在先驗信息的選擇上需謹慎，特別是在選取外部臨床先驗時，需要結合文獻和專家的意見并進行敏感性分析[12]。
本研究中還有兩個地方需要在接下來的研究中進行完善：一是異質性問題，在下一步研究中，將對原始研究中存在的協變量利用 Logistic 回歸進行檢驗，找出影響異質性的主要原因；二是基因型合并報告的問題，在部分原始研究中，研究者在報告結局的時候并不列出三種基因型的具體頻數，而是提供了在某種遺傳模型下的 OR 值，因此，無法直接利用文中方法進行合并，需要進一步利用外部信息，如人群某基因型的頻率等，利用貝葉斯方法進行合并。
綜上所述，本文對無遺傳模型約束的基因-疾病關聯研究貝葉斯 Meta 分析方法研究進行了實現，提供了 JAGS 和 R 代碼，并且在此基礎上通過實例比較了各種無信息先驗和外部先驗對于結果的影響。
在對慢性復雜疾病進行病因學的研究過程中，除考慮常見的環境危險因素外，往往需要考慮遺傳因素對疾病發病的影響。目前基因-疾病關聯研究常利用單核苷酸多態性（single nucleotide polymorphism，SNP）篩查與疾病相關的基因位點。由于測序技術的進步和成本的下降，基因-疾病關聯研究已大量開展，其文獻數量也呈爆發式增長[1]。由于單個基因多態性對于疾病的效應值往往較小[2]，因此需要更大樣本量來提高統計檢驗效能。因此，利用 Meta 分析方法整合各個獨立的基因-疾病關聯研究就顯得非常必要，目前已有研究者開展了此類研究的 Meta 分析方法學及應用研究[3, 4]。
在進行基于單核苷酸的基因-疾病關聯研究 Meta 分析中，比較突出的問題是如何選擇遺傳模型。研究者常假定一種遺傳模型（顯性或隱性遺傳模型）進行效應量的合并，但是這種方法缺乏生物學支持且較為武斷，如遺傳模型選用錯誤將導致基于基因-疾病關聯研究的 Meta 分析結論錯誤。Minelli 等[5, 6]提出通過計算不同基因型 OR 值的比值，并在此基礎上計算某基因型相對于參照基因型的 OR 值來避免選擇遺傳模型。Minelli 等[5, 6]分別利用極大似然估計和貝葉斯方法進行了估計，貝葉斯方法由于能夠利用外部信息作為先驗，估計結果更為精確，但這些研究未就外部先驗信息進行分析，只進行了多種模糊先驗的敏感性分析。
本研究在 Minelli 等研究的基礎上，利用各類模糊先驗、外部先驗信息，進行基因-疾病關聯研究 Meta 分析，利用 R 和 JAGS 軟件編寫 Meta 分析程序，采用 ACE I/D 基因多態性對于血管緊張素轉化酶抑制劑類降壓藥所致的咳嗽是否有影響的實例進行說明。
1 方法
1.1 基本原理
基因-疾病關聯研究 Meta 分析的特點之一是至少有 3 種可能的基因型分組，并且分析結果與所選擇的遺傳模型有關。在最簡單的雙等位基因情況下（野生基因型 A，突變基因型 G），一般通過收集 3 種不同基因型（AA，GA，GG）人群的目標疾病發生頻率來估計基因-疾病關聯性大小。因此，將會得到 2 個 OR 值，分別為 GG vs. AA 和 GA vs. AA，這 2 個 OR 值在遺傳模型下相互獨立。傳統的基因-疾病關聯研究 Meta 分析方法通過假設遺傳模型，將 3 組（AA，GA，GG）減少為 2 組，如：假設為隱性模型，則合并 GA 與 AA（GG vs. GA+AA）；假設為顯性模型，則合并 GG 與 GA（GG+GA vs. AA）；假設為超顯性模型，則合并 GG 與 AA（GG+AA vs. GA）等。
然而，遺傳模型通常是未知的，而錯誤的遺傳模型得出的 Meta 分析結果通常有誤。本研究提供的方法可避免選擇特定的遺傳模型，但需要計算 2 個 OR 值對數的比 λ：
|
假設 λ 在研究間是固定值，即 λ 取值 0、0.5、1 時則分別對應隱性模型、共顯性模型和顯性模型。但實際上 λ 可取 0 到 1 間的任何值，因此不受特定遺傳模型限制。
基因-疾病關聯性研究的設計主要是病例-對照研究，本文利用回顧性似然構建相應的似然函數，同時假定病例組和對照組的基因型服從多項式分布（multinomial distribution）：
y0i~Multinomial (n0, p0i), y1i~Multinomial (n1, p1i)
其中，i=1, 2, 3，分別對應 AA，GA，GG 三種基因型，y0i 為對照組某基因型頻率，y1i 為病例組基因型頻率，n0 為對照組總樣本例數，n1 為病例組總樣本例數。
|
其中，d 取值為 1 表示病例，取值為 0 表示對照，即 p0i 和 p1i 分別為某基因型在對照組和病例組中的概率。可將對照組進行以下的參數化過程：α1=1，α2=p02/p01，α3=p03/p01，進一步可構建 p01=1/（1+α2+α3），p02=α2/（1+α2+α3），p03=α3/（1+α2+α3）。同理，病例組參數化過程如下：β1=1，β2=p12/p11，β3=p13/p11。上述過程可得，GA vs. AA 的 OR 值為 OR2=β2/α2，GG vs. AA 的 OR 值為 OR3=β3/α3，其對數形式為 δ2=log（OR2），δ3=log（OR3），則 λ=δ2/δ3。如前所述，可以根據 λ 的值決定相應的遺傳模型（取值 0、0.5 或 1），或忽略遺傳模型直接進行效應量的合并（取值 0 到 1 間任何值）。病例組中各基因型概率可通過下式構建：
|
則每個研究的似然函數可寫為：
|
其中，假設研究間互相獨立，且第 j 個研究 GG vs. AA 的 OR 值的對數 δ3j 服從均數為 θ，方差為 τ2 的隨機效應正態分布：
|
GA vs. AA 的 OR 值的對數 δ2j 可由 δ3j·λ 計算得出。因此，需要明確未知參數 θ、τ 和 λ 的先驗分布。
其中，均數 θ 在所有模型中服從正態分布：
|
對于研究間的標準差 τ，本研究采用 3 種先驗分布：
分布 ①
分布 ②
分布 ③
其中，分布 ① 中，方差的倒數服從 Gamma 分布，雖然這是最常用的模糊先驗，但最近有學者推薦直接采用標準差的先驗[7]。分布 ② 中，研究間標準差 τ 服從標準半正態分布，這種先驗分布使得 τ>2 的概率很低。分布 ③ 為，0 到 2 范圍內的均勻分布，這種先驗分布排除了 τ>2 的可能性。
對于參數 λ，本研究采用 4 種先驗分布：
分布 ①
分布 ②
分布 ③
分布 ④
其中，分布 ① 為 2 個參數都為 1 的 Beta 分布，這種分布在 0 到 1 間是均勻的。但是，當參數取值接近 0 或 1 并且數據分布稀疏時，這種先驗分布傾向于將后驗估計值取為 0.5。例如，當近似于隱性模型時（λ 的真實值約等于 0），這種先驗分布更容易扭曲后驗估計值，因為有 90% 的概率得到大于 0.1 的值。分布 ② 為 2 個參數都為 0.5 的 Beta 分布，這種先驗分布在 λ 取接近 0 或 1 的值時，能有更大的先驗概率。但當模型近似于共顯性模型（λ=0.5）且數據分布稀疏時，這種分布可能會增加 λ 附近的不確定性。分布 ③ 為 0 到 1 范圍內的均勻分布，是一種常用的先驗分布。分布 ④ 為根據 Fujisawa 等[8]及 Minelli 等[5]的研究確定的外部先驗為均數 0.98，方差為 0.1 的正態分布。
1.2 實例
實例數據來自本文作者先前發表的 Cochrane 系統評價[9]，主要從基因多態性角度研究血管緊張素轉化酶抑制劑（ACEI）類降壓藥與咳嗽的關聯，共納入 17 個研究，其基本信息如表 1 所示。ACE 基因是高度多態的，其 16 內含子存在 287 bp 的 DNA 序列存在（插入，I）或沒有（缺失，D）的基因多態性，即研究最多的 ACE 基因插入/缺失（I/D）多態性。研究發現，D 等位基因與血漿 ACE 的水平呈正相關，DD 基因型者平均 ACE 水平大約是 II 基因型者的 2 倍。
表1
ACE I/D 基因多態性對血管緊張素轉化酶抑制劑類降壓藥所致咳嗽研究的基本特征
1.3 軟件實現
本文采用 R 和 JAGS 軟件進行編程，JAGS 軟件可在 http://mcmc-jags.sourceforge.net/ 下載，貝葉斯部分的 JAGS 代碼如下，以 txt 格式存儲在文件夾中即可：
model
{
for(i in 1:17 ) {
p0[i, 1]<-1/(1+a[i, 1]+a[i, 2]) #a 為公式中的 α
p0[i, 2]<-a[i, 1]/(1+a[i, 1]+a[i, 2])
p0[i, 3]<-a[i, 2]/(1+a[i, 1]+a[i, 2])
p1[i, 1]<-1/(1+a[i, 1]*exp(lambda*d[i])+a[i, 2]*exp(d[i])) #d 為 δ3
p1[i, 2]<-a[i, 1]*exp(lambda*d[i])/(1+a[i, 1]*exp(lambda*d[i])+a[i, 2]*exp(d[i]))
p1[i, 3]<-a[i, 2]*exp(d[i])/(1+a[i, 1]*exp(lambda*d[i])+a[i, 2]*exp(d[i]))
d[i]~dnorm(theta, prec)
ncont[i, 1:3]~dmulti(p0[i, 1:3], tcont[i])
ncase[i, 1:3]~dmulti(p1[i, 1:3], tcase[i])
a[i, 1]~dnorm(0, 0.0001)
a[i, 2]~dnorm(0, 0.0001)
}
lambda~dbeta(1, 1) #Beta 分布
#lambda~dbeta(0.5, 0.5)T(0.01, 0.99) #Beta 分布
#lambda~dunif(0, 1) #均勻分布
#lambda~dnorm(0.98, 0.1) #外部先驗
theta~dnorm(0, 0.000001)
tau~dgamma(0.001, 0.001) #Gamma 分布
#tau2~dunif(0, 1) #均勻分布
#tau2~dnorm(0, 1)I(0, ) #半正態分布
prec=1/(tau*tau)
or3<-exp(theta)
or2<-exp(theta*lambda)
}
在 R 中利用 R2jags 包調用 JAGS 程序，R 程序及相應的初始值如下：
library(R2jags)
library(dplyr)
ACEI <-read_excel("C:/work/modelfree/ACEI.xls") #讀取數據
attach(ACEI)
tcase=case_II+case_ID+case_DD
tcont=control_II+control_ID+control_DD
n<-17
ncase<-select(ACEI, contains("case"))
ncont<-select(ACEI, contains("control"))
data <-list("ncase", "ncont", "tcase", "tcont") #構建分析數據集
params <-c("lambda", "theta", "tau", "a", "or2", "or3") #構建參數集
#初始值設置，使用兩條鏈
a=matrix(c(rep(0.5, time=34)), ncol=2)
d<-c(rep(0, time=17))
inits1 <-list("lambda"=0.98, "theta"=0.5, "tau"=0.5, "a"=a, "d"=d)
inits2 <-list("lambda"=0.50, "theta"=1, "tau"=0.5, "a"=a, "d"=d)
inits <-list(inits1, inits2) #模型擬合
set.seed(20130801)
fit <-jags(data=data, inits=inits, parameters.to.save=params, n.chains=2, n.iter=50000, n.burnin=5000, model.file="C:/work/modelfree/ modelfree.txt")
print(fit)
plot(fit)
traceplot(fit)
本部分代碼主要涉及分析數據列表 data 的建立和初始值的設定，利用 R2jags 包中的 jags 函數進行抽樣和迭代，并利用 traceplot 函數進行收斂性評價，dplyr 包主要用作數據整理。在初始值的設置中，如果 λ 先驗信息為 dbeta（0.5，0.5）T（0.01，0.99），由于設置了抽樣的最大值和最小值，這時 λ 的初始值不能大于 0.99 或小于 0.01。
2 結果
不同的先驗信息下，ACE I/D 基因多態性對 ACEI 類降壓藥所致的咳嗽影響如表 2 所示。
表2
ACE I/D 基因多態性對 ACEI 類降壓藥所致咳嗽在不同先驗信息下的效應量合并結果
從表 2 中可以得出，相較于基因型Ⅱ，基因型 DD 和 DI 的攜帶者服用 ACEI 后引起干咳的可能性較小，與原文中采用顯性遺傳模型的結果相同[OR=0.70，95%CI（0.50，0.97）]。該結果在不同的先驗分布下非常穩定，但是研究間的異質性較大，可能引起的原因為納入研究的國家及人種不同引起，進一步亞組分析能夠減少研究間異質性。
3 討論
基因-疾病關聯 Meta 分析的無遺傳模型約束方法，同時考慮了遺傳效應的大小和遺傳模型的信息。雖然可用最大似然估計來實現，但貝葉斯方法是一種更具理論和實踐優勢的替代方法[5]。從理論角度看，貝葉斯方法能夠給予遺傳效應和遺傳模型明確的先驗信息。先驗信息可能來自于未被 Meta 分析納入的研究或專家意見，前者較為簡便而后者較難操作[10, 11]。從實際應用的角度看，采用馬爾可夫鏈蒙特卡羅（MCMC）算法估計的貝葉斯模型，更符合復雜遺傳模型下的基因疾病關聯。
Minelli 等[6]指出基因-疾病關聯 Meta 分析需要多個大樣本量的研究，使得 λ 先驗分布的不同選擇不會對結果造成過多的影響。同理，研究間方差 τ 的先驗分布也會影響結果。本研究結果顯示在研究數和樣本量較大的情況下，幾種無信息先驗和外部先驗對于結果的影響較小，也證實了 Minelli 等[6]的結論。另外，有研究指出，Gamma 分布并不是一個較好的研究間方差的先驗分布，而半正態分布和均勻分布可能更好[7]。從本文結果來看，3 種先驗分布下的結果差異不大，問題可能出現在數據分布稀疏或樣本量較小的情況下。因此，當納入研究較少的時，先驗信息對合并效應量的影響可能會較大，故在先驗信息的選擇上需謹慎，特別是在選取外部臨床先驗時，需要結合文獻和專家的意見并進行敏感性分析[12]。
本研究中還有兩個地方需要在接下來的研究中進行完善：一是異質性問題，在下一步研究中，將對原始研究中存在的協變量利用 Logistic 回歸進行檢驗，找出影響異質性的主要原因；二是基因型合并報告的問題，在部分原始研究中，研究者在報告結局的時候并不列出三種基因型的具體頻數，而是提供了在某種遺傳模型下的 OR 值，因此，無法直接利用文中方法進行合并，需要進一步利用外部信息，如人群某基因型的頻率等，利用貝葉斯方法進行合并。
綜上所述，本文對無遺傳模型約束的基因-疾病關聯研究貝葉斯 Meta 分析方法研究進行了實現，提供了 JAGS 和 R 代碼，并且在此基礎上通過實例比較了各種無信息先驗和外部先驗對于結果的影響。
