當醫學研究需要驗證兩個變量之間關系時，研究者通過改變和調試自變量來觀測因變量的變化；如果因變量隨著自變量的變化而變化，說明兩者是相關的，至于如何相關，則需要進一步利用不同的統計模型來做出判斷。因此構建合理的統計學模型探索研究變量間的關系，是研究者目前需要解決的問題。以往進行臨床醫學研究時，研究者通常使用最小二乘回歸和一般分位數回歸進行統計分析。最小二乘回歸關注的是均值，使各樣本殘差平方和最小，分位數回歸采用加權最小絕對離差和法進行估計。最小二乘回歸置信區間較大，若對數據進行不同分位數下的研究，便可使其方差盡可能小，參數估計精準度更高。分位數回歸在非正態分布臨床醫學數據中有應用條件更為寬松、挖掘的信息更豐富、估計結果可信性更高等特點。醫學研究中分位數回歸多用于費用數據[1-5]、年齡數據[6]、生命質量數據[7-9]中。Yu等[10]最早提出使用貝葉斯方法進行分位數回歸的統計分析方法。研究表明[11]貝葉斯分位數回歸參數估計結果精確程度和可信程度更高，在小樣本情形下也可以得到可靠參數信息[12]。目前該方法在臨床醫學研究領域應用較少，為介紹貝葉斯分位數回歸在臨床醫學數據分析中的應用，本文結合實例，使用R Studio軟件實踐臨床醫學研究數據的貝葉斯分位數回歸，展現其實際應用價值，為貝葉斯分位數回歸在臨床醫學數據研究中的應用提供參考。
1 資料與方法
1.1 數據來源與加載
R語言常用于統計計算和圖像繪制，廣泛應用于統計學和數據科學領域。R是世界上最強大的統計計算、機器學習和圖形編程語言，是數據科學領域最流行的語言，并且R語言是完全面向數據的，比其他的編程語言更加注重從數據的角度思考問題[13]。R是一套開源的統計計算機語言軟件（操作環境），其特點是免費、源代碼開放，適用于統計分析和繪圖的語言和操作環境。R Studio是R語言最常用的集成開發環境。
本文研究所用數據來自首都臨床特色應用研究專項250例膝骨關節炎患者臨床資料（表1）。評估膝骨關節炎患者血清IgG值與年齡之間的關系。數據需錄入Excel軟件中并保存為csv格式文件，存入相應電腦硬盤分區（此例筆者將文件列入D盤）建立名為IgG~Age.csv的Excel表格，采用英文及英文縮寫便于代碼書寫及程序運行，“年齡”定義為“Age”，錄入對應數據，實現數據清洗及數據整理。
表1
250例膝骨關節炎患者年齡及IgG值
在安裝R Studio之前須安裝最新版本R語言軟件。使用加載程序函數依次加載多個數據包（dplyr、ggplot2、quantreg等）。
1.2 原始數據讀入與格式轉換
建立R系統臨時數據框dataframe進行數據加載。
執行命令：
> fichier<-"D:/IgG~Age.csv"
> mydata <-read.table（fichier, header=TRUE, sep=","）
> mydata
讀取D盤目標文件IgG~Age.csv數據，header=TRUE為邏輯值，其表示文件第一行已包含變量名稱，即數值讀取時不會將其納入而出現運行錯誤。sep=","即用逗號來分隔行內數據值[14]。
1.3 數據提取與賦值
執行命令：
>IgG<-mydata$IgG
>Age<-mydata$Age
$表示從R Studio臨時構建的dataframe中讀取出對應列數據。
2 結果
2.1 貝葉斯分位數回歸估計
分別建立貝葉斯分位數雙變量一次項、二次項回歸方程式：
IgG[i]==beta[0]+beta[1]*age[i]
IgG[i]==beta[0]+beta[1]*Age[i]+beta[2]*Age[i]^2
執行命令：
>fit=Brq（IgG~age, tau=0.5, data=mydata, runs=10000）
>summary（fit）
>main=expression（IgG[i]==beta[0]+beta[1]*age[i]）
>plot（age, IgG, main=main, cex=1.5, axes=F, pch=20, col='#cc1e1e55'）
>axis（side=1, lwd=2）
>axis（side=2, lwd=2）
>legend（"topleft", y=max（IgG）+0.8, legend=c（.05,0.25,0.50,0.75,0.95）,
lty=c（1,1,1,1,1）, lwd=c（1,1,1,1,1）, col=c（1:5）, title="Quantile"）
>estimate=data.frame（）
>for（i in 1:5） {
taus=c（0.05, 0.25, 0.5, 0.75, 0.95）
fit=Brq（IgG~age, tau=taus[i], runs=10000, burn=1000）
est=summary（fit）
estimate_i <- data.frame（est$coefficients）
estimate_i['variables']=rownames（estimate_i）
estimate_i['quantile']=fit$tau
estimate=rbind（estimate, estimate_i）
abline（fit, col=i, lwd=2）
print（summary（fit））
}
>p2 <- ggplot（data=estimate, aes（x=quantile, y=Estimate））+
geom_point（）+
geom_smooth（method='loess'）+
geom_ribbon（aes（ymin=L.CredIntv, ymax = U.CredIntv）, fill='gray', alpha=0.5）+
facet_wrap（'variables', nrow=1, scales='free'）+
theme_bw（base_size=16）
>plot（fit, plottype="tracehist", Coefficients=c（1:2）, D=c（1:2）,col1=2,col2=2,col3=3,col4=4）
重復上述代碼進行二次函數貝葉斯分位數回歸。
我們使用R語言的quantreg包進行分位數回歸，quantreg即Quantile Regression。命令表達的是IgG與Age之間的函數關系。tau表示計算不同分位點的參數，改變tau值，可以得到不同分位數下的參數估計值。此處選取tau=0.5進行運算。main函數為主函數，是程序運行的入口。expression為載入的表達式，本例我們通過分別載入一次項和二次項方程式來展現直線擬合與曲線擬合的區別。為方便進行兩圖對比，本文將兩函數圖放入一頁。使用plot（）函數進行畫圖，使用par（）函數進行繪圖參數特定，使用mfrow參數設定繪圖布局，定義為一頁兩圖。調整cex改變點符號大小，axes=F表示禁止生成坐標軸，col=2表示使用紅色，pch=20代表實心點，使用axis（）函數添加軸線，線條和刻度線的顏色定義為黃色。設置標簽，將其置于右上方，定義不同分位數下線條顏色、類型及寬度，標簽圖例標題為“分位數”。將估計值輸出至數據框中。使用for in語句，在0.05～0.95分位數下依次進行代碼執行，遍歷完不同分位數即循環結束。每一次對抽樣過程的重復稱為一次“迭代”，每一次迭代所得的結果將會作為下一次迭代的初始值。構造馬爾科夫鏈進行重復抽樣時，在收斂前的一段時間，分布還未達到穩定分布，為消除其對抽樣結果的影響，需將初始n個迭代值去掉。本文使用run語句對樣本迭代抽樣10 000次，使用burn語句拋去前1 000次抽樣，實現對模型的退火（burn in），利用第1 001次到10 000次的抽樣數據進行參數估計。本文使用abline函數繪制一次函數圖像中的直線，使用curve函數繪制二次函數對應的曲線，將led設置線條寬度為默認寬度的二倍。使用print語句輸出參數估計。退火次數越高，估計值越準確。輸出結果見圖1。
圖1
曲線擬合圖
將結果輸出至p2，aes函數配合ggplot（）使用，設定自變量X為quantile，因變量Y為estimate。使用geom_point（）繪制圖中散點，使用geom_smooth（）函數向圖中添加擬合曲線，默認平滑算法，使用geom_ribbon（）函數添加灰色置信區間。將結果輸出至p2，aes（）函數配合ggplot（）使用，設定自變量X為quantile，因變量Y為estimate。使用geom_point（）繪制圖中散點，使用geom_smooth（）函數向圖中添加擬合曲線，默認平滑算法，使用geom_ribbon（）函數添加灰色置信區間。輸出結果見圖2。
圖2
系數變化情況
將type參數設定為矩陣直方圖形式，進行吉布斯抽樣收斂性的路徑圖檢查及后驗直方圖評估。輸出結果見圖3、圖4。
圖3
一次函數迭代運行軌跡圖及密度函數圖
圖4
二次函數迭代運行軌跡圖及密度函數圖
圖1可以看出曲線擬合圖中在不同分位下得到的值差距很大，而直線擬合圖中不同分位下得到的值差距較小，故曲線擬合更能體現出圖像變化趨勢，參數估計更加準確。
計算后驗分布期望的傳統數值計算方法（頻率法）是數值積分、拉普萊斯近似計算和蒙特卡洛重要抽樣。MCMC方法，即馬爾可夫鏈—蒙特卡羅（Markov chain Monte Carlo）方法目前是最流行的貝葉斯計算方法，一方面是由于它處理非常復雜問題的效率，另一方面是因為它的編程方法相對容易，是能夠解決復雜概率統計問題的有效工具[15]。貝葉斯估計尤其是采用MCMC方法，在小樣本性質、假設檢驗以及預測方面具有極大的優勢[16]。MCMC方法目前常的主要有Gibbs抽樣和Metropolis-Hastings算法。Gibbs抽樣是應用最廣泛的MCMC方法，本文使用Gibbs抽樣進行后驗分布。后驗Gibbs抽樣軌跡圖判斷參數迭代過程是否收斂。通過各參數的后驗密度函數圖，可進一步判斷抽樣送代的收斂性。如圖所示，由迭代軌跡圖可以看出抽樣序列比較集中，在小范圍內波動，沒有明顯的周期性和趨勢性，密度函數圖近似于正態，說明馬爾科夫鏈是收斂的，MCMC模擬過程處于平穩狀態，利用10 000次迭代，1 000次退火進行參數估計是有效的。
使用貝葉斯分位數回歸進行參數估計，將所得系數代入回歸公式，得到不同分位數下回歸公式：Y1=?6.022 063 47+2.026 913 73X?0.015 077 69X2……Y5=24.610 542 414?0.395 059 497X+0.004 205 064X2，據此可以發現膝骨關節炎患者血清IgG值有隨著年齡升高的趨勢。
3 討論
3.1 貝葉斯分位數回歸估計
貝葉斯分位數回歸是將目前應用廣泛的貝葉斯方法與分位數回歸結合，構建貝葉斯分位數回歸模型，采用分位數回歸分析各分位點參數的變化，通過貝葉斯理論讓參數值估計更精確。充分發揮二者的優勢，是一種強有力的參數估計方法。Yu等[10]給出貝葉斯分位數回歸基本原理，利用非對稱Laplace分布，將對損失函數的求解轉化為對似然函數的求解，再根據貝葉斯推斷得到估計參數的后驗分布，引出了貝葉斯分位數回歸模型。其核心思想是將不符合正態分布的數據，經過反復抽樣，將其抽成真實世界情況。研究時需要考慮臨床醫學數據中非正態分布等特點，因此貝葉斯分位數回歸是研究者很好的選擇。以往因其統計學方法太過復雜，對計算機技術要求較高，一直以來未能廣泛使用。
在臨床醫學研究實踐中，常常會遇到兩個變量間的關系不呈直線而是呈某種曲線關系又或者兩變量在某局部范圍內呈現直線關系，但整個范圍又顯現出曲線趨勢。此時若用直線強行擬合則會產生較多離散點，使得結果與真實情況相差較大。因此我們常會根據樣本散點分布的形狀選擇合適的曲線函數。曲線擬合的目的就是找到一個合適的曲線回歸方程來更合理得描述變量間的關系。本文我們通過一次函數與最基本的二次函數曲線圖像的對比，展現了曲線擬合相對于直線擬合的優勢。
3.2 分位數回歸優勢
普通最小二乘回歸分析關注的是根據給定自變量平均值來預測因變量的平均值。當回歸殘差是正態分布時，那么關于回歸系數的推論就是有效的、無偏的。然而，在進行臨床醫學數據統計分析時通常會出現嚴重偏離的數據，數據存在離群值，在因變量處于尖峰、厚尾分布等情況下，均值回歸結果可信度就會降低，同時可能會遺漏一些重要信息[17]。因此我們需要一種能夠得到更佳估計效果的回歸方法。
Koenker和Bassett（1978）首次提出了分位數回歸思想。分位數回歸是估計自變量與因變量不同分位數之間線性關系的一種建模方法。將因變量拆分為多個分位數點（如25%、50%、75%），研究不同分位點下回歸影響的關系。其中，50%分位點的回歸也叫中位數回歸。分位數回歸采用的是加權的最小絕對離差（weighted least squares absolute deviation，WLAD）進行估計。使用分位數回歸處理臨床醫學樣本數據時無需符合正態分布或近似正態分布假定，因此其應用條件較最小二乘回歸分析更為寬松。對于符合正態分布或近似正態分布的臨床醫學數據可以使用最小二乘回歸，對于不符合正態分布的數據可選擇分位數回歸進行描述分析。普通最小二乘回歸只能提供平均數，假設自變量對因變量的影響在全局分布上是恒定的。分位數回歸能提供不同分位數下的估計結果，能更加全面地刻畫整個條件分布全貌和各個分位處的局部特征，能夠展現自變量與因變量之間關系的細節，檢測出可能不會被注意到的關聯，若只是檢查均值，會遺漏這些細節。因此，分位數回歸可以更精確地展現變量間關系的變化。此外，分位數回歸本身的特性決定了其受醫學樣本數據中異常值的影響較小，因而其估計更為穩健。簡言之，分位數回歸在非正態分布臨床醫學數據中有應用條件更為寬松、挖掘的信息更豐富、估計結果可信性更高等特點。隨著分位數回歸算法不斷發展，分位數越來越廣泛地被應用到醫學領域[1-9,18-23]。
3.3 貝葉斯分位數回歸優勢
分位數回歸常用參數估計方法有頻率派估計和貝葉斯估計兩種。貝葉斯相較于頻率派具有不可比擬的優勢[24]。頻率學派將模型參數視為固定常數，而貝葉斯將模型參數視作隨機變量，充分考慮到參數的不確定性。貝葉斯估計綜合利用先驗信息和樣本信息通過貝葉斯公式得到后驗信息，再根據后驗信息去推斷未知參數。基于抽樣的貝葉斯方法，在小樣本條件下也能得到較優的參數估計。以往缺乏技術支持，分位數回歸研究多使用頻率法，貝葉斯法在眾多研究中不便實現。隨著計算機技術的發展，計算機運算速度越來越快，貝葉斯法可以輕松實現，因此貝葉斯流派在逐漸興起，大有超過頻率學派的趨勢。貝葉斯估計可以利用依據貝葉斯定理得到的新信息迭代處理得到進一步的信息，數據處理過程更加科學。貝葉斯法加入了參數的先驗信息，合適的先驗分布可提高臨床醫學數據參數估計精確度。當無先驗信息時，亦不會影響參數估計的準確性。貝葉斯估計能夠充分利用樣本信息以及參數的先驗信息，是一種能夠有效地獲得參數后驗分布的方法。
樣本數據較小是眾多臨床醫學數據研究中常遇到的問題，在眾多臨床研究當中很難實現大樣本數，難以直接代表整體。不同于頻率法對樣本量的敏感程度，貝葉斯法定了先驗，然后進行反復抽樣，更具有細膩性，對小數據具有更好的推斷性[25]。貝葉斯方法在醫學研究小樣本條件下所得結論更加可靠，提高了估計準確度和預測的精度[26-27]。近年來大量的研究表明貝葉斯方法是一種合理有效的統計推斷方法。貝葉斯原理在大數據時代有獨特的價值，是目前研究者在沒有充分或準確信息時最優的推理結構[28]。現如今基于貝葉斯原理發展而來的貝葉斯網絡不確定性知識表達和推理方法，在醫學領域已得到廣泛應用[29-33]。
3.4 貝葉斯分位數回歸應用
貝葉斯分位數回歸參數估計已廣泛應用于金融、互聯網等各大領域，但在臨床醫學研究領域未得到足夠重視。張穎等[12]將基于Gibbs抽樣的MCMC算法引入GJR-CAViaR模型中，利用上證綜指日收益率數據進行實證分析，展現了該方法的優勢。謝婷婷等[11]運用貝葉斯分位數法對科技創新、金融發展對產業結構升級的關系進行實證檢驗，分析其影響程度。羅玉波等[34]分別使用簡單最小二乘法、普通的分位數回歸、貝葉斯類等6種方法對模擬數據進行估計，發現在絕大多數情形下，貝葉斯類方法優于其他的估計方法。并分析了在貝葉斯框架下，運用Gibbs抽樣法對模型進行條件分位數回歸的可行性。結果發現模擬研究和實證研究中，其估計效果均穩健合理。陸平等[35]使用貝葉斯分位數回歸分析互聯網產品眾籌行為與績效之間的關系。方麗婷等[36]提出了空間滯后分位數回歸模型的貝葉斯估計方法，結果表明其提出的貝葉斯估計方法在小樣本條件下仍有良好的估計效果，參數估計精確度較高，應用實例展示了貝葉斯分位數回歸方法的實際應用價值。姚萍等[37]建立分位數SV模型，進行貝葉斯框架下的統計分析，以中國金融市場數據進行實證。Kuhudzai等[38]通過對比一般和貝葉斯分位數回歸在血壓風險因素對血壓分布的不同分位數下的影響中的應用，展現了貝葉斯法估計更加精確。Alampi等[39]利用貝葉斯分位數回歸分析妊娠期有毒物質暴露與自閉癥行為間的關系，同時展現了分位數回歸可以揭示變量間更加細微的關系。HyungBok等[40]使用貝葉斯分位數回歸機器學習法研究分析冠狀動脈疾病不同階段與多個心血管危險因素的關系。Suleiman等[41]使用貝葉斯有序分位數回歸研究強迫癥患者對氟伏沙明的反應。
目前，分位數回歸在醫學領域常用于以下數據情況：① 數據有極端值；② 數據呈現尖峰或者厚尾的分布情況；③ 數據存在嚴重異方差；④ 數據呈偏態分布以及方差不齊，數據方差齊性及正態性不滿足。分位數回歸的使用在連續數據[20]、分類數據[2,21]、生存數據[22]方面均有涉及，因其應用條件更為寬松，涉及數據不需要滿足正態分布，故連續數據多為非正態分布數據，也用于正態分布數據；分類數據中二分類數據、有序多分類數據、無序多分類數據均有所涉及；離散數據也有涉及，也可將離散變量轉化為連續變量進行分析，本文所選實例中年齡變量即為離散變量。分位數回歸解決的問題多以線性相關關系為主，較少涉及非線性相關關系[23]。
當臨床醫學研究樣本量足夠大可以代表整體時，研究者對于該樣本數據可直接進行描述性統計，數據為正態分布時，可以直接使用方差分析進行一般線性回歸。樣本數量越少，統計模型的設計就越困難，對醫學統計方法要求就越高，研究者進行假設檢驗就越難，越不符合真實世界，這種情況下研究者可以使用貝葉斯方法進行反復抽樣將其抽成真實世界，以提高研究準確性。然而即使貝葉斯法適合于小樣本醫學臨床數據。若樣本量過小，即使抽樣迭代很多次也可能出現馬爾科夫鏈不收斂、抽樣序列不集中的情況，這樣得到的參數的估計值就會較偏離真實值，估計的不確定性較大。
綜上所述，本文實例分析展現了貝葉斯分位數回歸方法參數估計結果精確，可信程度較高，在小樣本條件下也可以得到可靠參數信息，并且展現了曲線擬合在非線性關系中相對較于直線擬合的優勢。貝葉斯分位數回歸在臨床醫學中的應用較少，更多體現在理論層面，還需要更多高質量的案例對貝葉斯分位數回歸法的優勢進行應用、分析和總結。本文通過介紹貝葉斯分位數回歸方法并使用R Studio軟件進行實例分析，展現了貝葉斯分位數回歸在臨床醫學研究的應用價值，并提供具體代碼編寫，為臨床醫學數據研究者統計方法的使用提供參考。
當醫學研究需要驗證兩個變量之間關系時，研究者通過改變和調試自變量來觀測因變量的變化；如果因變量隨著自變量的變化而變化，說明兩者是相關的，至于如何相關，則需要進一步利用不同的統計模型來做出判斷。因此構建合理的統計學模型探索研究變量間的關系，是研究者目前需要解決的問題。以往進行臨床醫學研究時，研究者通常使用最小二乘回歸和一般分位數回歸進行統計分析。最小二乘回歸關注的是均值，使各樣本殘差平方和最小，分位數回歸采用加權最小絕對離差和法進行估計。最小二乘回歸置信區間較大，若對數據進行不同分位數下的研究，便可使其方差盡可能小，參數估計精準度更高。分位數回歸在非正態分布臨床醫學數據中有應用條件更為寬松、挖掘的信息更豐富、估計結果可信性更高等特點。醫學研究中分位數回歸多用于費用數據[1-5]、年齡數據[6]、生命質量數據[7-9]中。Yu等[10]最早提出使用貝葉斯方法進行分位數回歸的統計分析方法。研究表明[11]貝葉斯分位數回歸參數估計結果精確程度和可信程度更高，在小樣本情形下也可以得到可靠參數信息[12]。目前該方法在臨床醫學研究領域應用較少，為介紹貝葉斯分位數回歸在臨床醫學數據分析中的應用，本文結合實例，使用R Studio軟件實踐臨床醫學研究數據的貝葉斯分位數回歸，展現其實際應用價值，為貝葉斯分位數回歸在臨床醫學數據研究中的應用提供參考。
1 資料與方法
1.1 數據來源與加載
R語言常用于統計計算和圖像繪制，廣泛應用于統計學和數據科學領域。R是世界上最強大的統計計算、機器學習和圖形編程語言，是數據科學領域最流行的語言，并且R語言是完全面向數據的，比其他的編程語言更加注重從數據的角度思考問題[13]。R是一套開源的統計計算機語言軟件（操作環境），其特點是免費、源代碼開放，適用于統計分析和繪圖的語言和操作環境。R Studio是R語言最常用的集成開發環境。
本文研究所用數據來自首都臨床特色應用研究專項250例膝骨關節炎患者臨床資料（表1）。評估膝骨關節炎患者血清IgG值與年齡之間的關系。數據需錄入Excel軟件中并保存為csv格式文件，存入相應電腦硬盤分區（此例筆者將文件列入D盤）建立名為IgG~Age.csv的Excel表格，采用英文及英文縮寫便于代碼書寫及程序運行，“年齡”定義為“Age”，錄入對應數據，實現數據清洗及數據整理。
表1
250例膝骨關節炎患者年齡及IgG值
在安裝R Studio之前須安裝最新版本R語言軟件。使用加載程序函數依次加載多個數據包（dplyr、ggplot2、quantreg等）。
1.2 原始數據讀入與格式轉換
建立R系統臨時數據框dataframe進行數據加載。
執行命令：
> fichier<-"D:/IgG~Age.csv"
> mydata <-read.table（fichier, header=TRUE, sep=","）
> mydata
讀取D盤目標文件IgG~Age.csv數據，header=TRUE為邏輯值，其表示文件第一行已包含變量名稱，即數值讀取時不會將其納入而出現運行錯誤。sep=","即用逗號來分隔行內數據值[14]。
1.3 數據提取與賦值
執行命令：
>IgG<-mydata$IgG
>Age<-mydata$Age
$表示從R Studio臨時構建的dataframe中讀取出對應列數據。
2 結果
2.1 貝葉斯分位數回歸估計
分別建立貝葉斯分位數雙變量一次項、二次項回歸方程式：
IgG[i]==beta[0]+beta[1]*age[i]
IgG[i]==beta[0]+beta[1]*Age[i]+beta[2]*Age[i]^2
執行命令：
>fit=Brq（IgG~age, tau=0.5, data=mydata, runs=10000）
>summary（fit）
>main=expression（IgG[i]==beta[0]+beta[1]*age[i]）
>plot（age, IgG, main=main, cex=1.5, axes=F, pch=20, col='#cc1e1e55'）
>axis（side=1, lwd=2）
>axis（side=2, lwd=2）
>legend（"topleft", y=max（IgG）+0.8, legend=c（.05,0.25,0.50,0.75,0.95）,
lty=c（1,1,1,1,1）, lwd=c（1,1,1,1,1）, col=c（1:5）, title="Quantile"）
>estimate=data.frame（）
>for（i in 1:5） {
taus=c（0.05, 0.25, 0.5, 0.75, 0.95）
fit=Brq（IgG~age, tau=taus[i], runs=10000, burn=1000）
est=summary（fit）
estimate_i <- data.frame（est$coefficients）
estimate_i['variables']=rownames（estimate_i）
estimate_i['quantile']=fit$tau
estimate=rbind（estimate, estimate_i）
abline（fit, col=i, lwd=2）
print（summary（fit））
}
>p2 <- ggplot（data=estimate, aes（x=quantile, y=Estimate））+
geom_point（）+
geom_smooth（method='loess'）+
geom_ribbon（aes（ymin=L.CredIntv, ymax = U.CredIntv）, fill='gray', alpha=0.5）+
facet_wrap（'variables', nrow=1, scales='free'）+
theme_bw（base_size=16）
>plot（fit, plottype="tracehist", Coefficients=c（1:2）, D=c（1:2）,col1=2,col2=2,col3=3,col4=4）
重復上述代碼進行二次函數貝葉斯分位數回歸。
我們使用R語言的quantreg包進行分位數回歸，quantreg即Quantile Regression。命令表達的是IgG與Age之間的函數關系。tau表示計算不同分位點的參數，改變tau值，可以得到不同分位數下的參數估計值。此處選取tau=0.5進行運算。main函數為主函數，是程序運行的入口。expression為載入的表達式，本例我們通過分別載入一次項和二次項方程式來展現直線擬合與曲線擬合的區別。為方便進行兩圖對比，本文將兩函數圖放入一頁。使用plot（）函數進行畫圖，使用par（）函數進行繪圖參數特定，使用mfrow參數設定繪圖布局，定義為一頁兩圖。調整cex改變點符號大小，axes=F表示禁止生成坐標軸，col=2表示使用紅色，pch=20代表實心點，使用axis（）函數添加軸線，線條和刻度線的顏色定義為黃色。設置標簽，將其置于右上方，定義不同分位數下線條顏色、類型及寬度，標簽圖例標題為“分位數”。將估計值輸出至數據框中。使用for in語句，在0.05～0.95分位數下依次進行代碼執行，遍歷完不同分位數即循環結束。每一次對抽樣過程的重復稱為一次“迭代”，每一次迭代所得的結果將會作為下一次迭代的初始值。構造馬爾科夫鏈進行重復抽樣時，在收斂前的一段時間，分布還未達到穩定分布，為消除其對抽樣結果的影響，需將初始n個迭代值去掉。本文使用run語句對樣本迭代抽樣10 000次，使用burn語句拋去前1 000次抽樣，實現對模型的退火（burn in），利用第1 001次到10 000次的抽樣數據進行參數估計。本文使用abline函數繪制一次函數圖像中的直線，使用curve函數繪制二次函數對應的曲線，將led設置線條寬度為默認寬度的二倍。使用print語句輸出參數估計。退火次數越高，估計值越準確。輸出結果見圖1。
圖1
曲線擬合圖
將結果輸出至p2，aes函數配合ggplot（）使用，設定自變量X為quantile，因變量Y為estimate。使用geom_point（）繪制圖中散點，使用geom_smooth（）函數向圖中添加擬合曲線，默認平滑算法，使用geom_ribbon（）函數添加灰色置信區間。將結果輸出至p2，aes（）函數配合ggplot（）使用，設定自變量X為quantile，因變量Y為estimate。使用geom_point（）繪制圖中散點，使用geom_smooth（）函數向圖中添加擬合曲線，默認平滑算法，使用geom_ribbon（）函數添加灰色置信區間。輸出結果見圖2。
圖2
系數變化情況
將type參數設定為矩陣直方圖形式，進行吉布斯抽樣收斂性的路徑圖檢查及后驗直方圖評估。輸出結果見圖3、圖4。
圖3
一次函數迭代運行軌跡圖及密度函數圖
圖4
二次函數迭代運行軌跡圖及密度函數圖
圖1可以看出曲線擬合圖中在不同分位下得到的值差距很大，而直線擬合圖中不同分位下得到的值差距較小，故曲線擬合更能體現出圖像變化趨勢，參數估計更加準確。
計算后驗分布期望的傳統數值計算方法（頻率法）是數值積分、拉普萊斯近似計算和蒙特卡洛重要抽樣。MCMC方法，即馬爾可夫鏈—蒙特卡羅（Markov chain Monte Carlo）方法目前是最流行的貝葉斯計算方法，一方面是由于它處理非常復雜問題的效率，另一方面是因為它的編程方法相對容易，是能夠解決復雜概率統計問題的有效工具[15]。貝葉斯估計尤其是采用MCMC方法，在小樣本性質、假設檢驗以及預測方面具有極大的優勢[16]。MCMC方法目前常的主要有Gibbs抽樣和Metropolis-Hastings算法。Gibbs抽樣是應用最廣泛的MCMC方法，本文使用Gibbs抽樣進行后驗分布。后驗Gibbs抽樣軌跡圖判斷參數迭代過程是否收斂。通過各參數的后驗密度函數圖，可進一步判斷抽樣送代的收斂性。如圖所示，由迭代軌跡圖可以看出抽樣序列比較集中，在小范圍內波動，沒有明顯的周期性和趨勢性，密度函數圖近似于正態，說明馬爾科夫鏈是收斂的，MCMC模擬過程處于平穩狀態，利用10 000次迭代，1 000次退火進行參數估計是有效的。
使用貝葉斯分位數回歸進行參數估計，將所得系數代入回歸公式，得到不同分位數下回歸公式：Y1=?6.022 063 47+2.026 913 73X?0.015 077 69X2……Y5=24.610 542 414?0.395 059 497X+0.004 205 064X2，據此可以發現膝骨關節炎患者血清IgG值有隨著年齡升高的趨勢。
3 討論
3.1 貝葉斯分位數回歸估計
貝葉斯分位數回歸是將目前應用廣泛的貝葉斯方法與分位數回歸結合，構建貝葉斯分位數回歸模型，采用分位數回歸分析各分位點參數的變化，通過貝葉斯理論讓參數值估計更精確。充分發揮二者的優勢，是一種強有力的參數估計方法。Yu等[10]給出貝葉斯分位數回歸基本原理，利用非對稱Laplace分布，將對損失函數的求解轉化為對似然函數的求解，再根據貝葉斯推斷得到估計參數的后驗分布，引出了貝葉斯分位數回歸模型。其核心思想是將不符合正態分布的數據，經過反復抽樣，將其抽成真實世界情況。研究時需要考慮臨床醫學數據中非正態分布等特點，因此貝葉斯分位數回歸是研究者很好的選擇。以往因其統計學方法太過復雜，對計算機技術要求較高，一直以來未能廣泛使用。
在臨床醫學研究實踐中，常常會遇到兩個變量間的關系不呈直線而是呈某種曲線關系又或者兩變量在某局部范圍內呈現直線關系，但整個范圍又顯現出曲線趨勢。此時若用直線強行擬合則會產生較多離散點，使得結果與真實情況相差較大。因此我們常會根據樣本散點分布的形狀選擇合適的曲線函數。曲線擬合的目的就是找到一個合適的曲線回歸方程來更合理得描述變量間的關系。本文我們通過一次函數與最基本的二次函數曲線圖像的對比，展現了曲線擬合相對于直線擬合的優勢。
3.2 分位數回歸優勢
普通最小二乘回歸分析關注的是根據給定自變量平均值來預測因變量的平均值。當回歸殘差是正態分布時，那么關于回歸系數的推論就是有效的、無偏的。然而，在進行臨床醫學數據統計分析時通常會出現嚴重偏離的數據，數據存在離群值，在因變量處于尖峰、厚尾分布等情況下，均值回歸結果可信度就會降低，同時可能會遺漏一些重要信息[17]。因此我們需要一種能夠得到更佳估計效果的回歸方法。
Koenker和Bassett（1978）首次提出了分位數回歸思想。分位數回歸是估計自變量與因變量不同分位數之間線性關系的一種建模方法。將因變量拆分為多個分位數點（如25%、50%、75%），研究不同分位點下回歸影響的關系。其中，50%分位點的回歸也叫中位數回歸。分位數回歸采用的是加權的最小絕對離差（weighted least squares absolute deviation，WLAD）進行估計。使用分位數回歸處理臨床醫學樣本數據時無需符合正態分布或近似正態分布假定，因此其應用條件較最小二乘回歸分析更為寬松。對于符合正態分布或近似正態分布的臨床醫學數據可以使用最小二乘回歸，對于不符合正態分布的數據可選擇分位數回歸進行描述分析。普通最小二乘回歸只能提供平均數，假設自變量對因變量的影響在全局分布上是恒定的。分位數回歸能提供不同分位數下的估計結果，能更加全面地刻畫整個條件分布全貌和各個分位處的局部特征，能夠展現自變量與因變量之間關系的細節，檢測出可能不會被注意到的關聯，若只是檢查均值，會遺漏這些細節。因此，分位數回歸可以更精確地展現變量間關系的變化。此外，分位數回歸本身的特性決定了其受醫學樣本數據中異常值的影響較小，因而其估計更為穩健。簡言之，分位數回歸在非正態分布臨床醫學數據中有應用條件更為寬松、挖掘的信息更豐富、估計結果可信性更高等特點。隨著分位數回歸算法不斷發展，分位數越來越廣泛地被應用到醫學領域[1-9,18-23]。
3.3 貝葉斯分位數回歸優勢
分位數回歸常用參數估計方法有頻率派估計和貝葉斯估計兩種。貝葉斯相較于頻率派具有不可比擬的優勢[24]。頻率學派將模型參數視為固定常數，而貝葉斯將模型參數視作隨機變量，充分考慮到參數的不確定性。貝葉斯估計綜合利用先驗信息和樣本信息通過貝葉斯公式得到后驗信息，再根據后驗信息去推斷未知參數。基于抽樣的貝葉斯方法，在小樣本條件下也能得到較優的參數估計。以往缺乏技術支持，分位數回歸研究多使用頻率法，貝葉斯法在眾多研究中不便實現。隨著計算機技術的發展，計算機運算速度越來越快，貝葉斯法可以輕松實現，因此貝葉斯流派在逐漸興起，大有超過頻率學派的趨勢。貝葉斯估計可以利用依據貝葉斯定理得到的新信息迭代處理得到進一步的信息，數據處理過程更加科學。貝葉斯法加入了參數的先驗信息，合適的先驗分布可提高臨床醫學數據參數估計精確度。當無先驗信息時，亦不會影響參數估計的準確性。貝葉斯估計能夠充分利用樣本信息以及參數的先驗信息，是一種能夠有效地獲得參數后驗分布的方法。
樣本數據較小是眾多臨床醫學數據研究中常遇到的問題，在眾多臨床研究當中很難實現大樣本數，難以直接代表整體。不同于頻率法對樣本量的敏感程度，貝葉斯法定了先驗，然后進行反復抽樣，更具有細膩性，對小數據具有更好的推斷性[25]。貝葉斯方法在醫學研究小樣本條件下所得結論更加可靠，提高了估計準確度和預測的精度[26-27]。近年來大量的研究表明貝葉斯方法是一種合理有效的統計推斷方法。貝葉斯原理在大數據時代有獨特的價值，是目前研究者在沒有充分或準確信息時最優的推理結構[28]。現如今基于貝葉斯原理發展而來的貝葉斯網絡不確定性知識表達和推理方法，在醫學領域已得到廣泛應用[29-33]。
3.4 貝葉斯分位數回歸應用
貝葉斯分位數回歸參數估計已廣泛應用于金融、互聯網等各大領域，但在臨床醫學研究領域未得到足夠重視。張穎等[12]將基于Gibbs抽樣的MCMC算法引入GJR-CAViaR模型中，利用上證綜指日收益率數據進行實證分析，展現了該方法的優勢。謝婷婷等[11]運用貝葉斯分位數法對科技創新、金融發展對產業結構升級的關系進行實證檢驗，分析其影響程度。羅玉波等[34]分別使用簡單最小二乘法、普通的分位數回歸、貝葉斯類等6種方法對模擬數據進行估計，發現在絕大多數情形下，貝葉斯類方法優于其他的估計方法。并分析了在貝葉斯框架下，運用Gibbs抽樣法對模型進行條件分位數回歸的可行性。結果發現模擬研究和實證研究中，其估計效果均穩健合理。陸平等[35]使用貝葉斯分位數回歸分析互聯網產品眾籌行為與績效之間的關系。方麗婷等[36]提出了空間滯后分位數回歸模型的貝葉斯估計方法，結果表明其提出的貝葉斯估計方法在小樣本條件下仍有良好的估計效果，參數估計精確度較高，應用實例展示了貝葉斯分位數回歸方法的實際應用價值。姚萍等[37]建立分位數SV模型，進行貝葉斯框架下的統計分析，以中國金融市場數據進行實證。Kuhudzai等[38]通過對比一般和貝葉斯分位數回歸在血壓風險因素對血壓分布的不同分位數下的影響中的應用，展現了貝葉斯法估計更加精確。Alampi等[39]利用貝葉斯分位數回歸分析妊娠期有毒物質暴露與自閉癥行為間的關系，同時展現了分位數回歸可以揭示變量間更加細微的關系。HyungBok等[40]使用貝葉斯分位數回歸機器學習法研究分析冠狀動脈疾病不同階段與多個心血管危險因素的關系。Suleiman等[41]使用貝葉斯有序分位數回歸研究強迫癥患者對氟伏沙明的反應。
目前，分位數回歸在醫學領域常用于以下數據情況：① 數據有極端值；② 數據呈現尖峰或者厚尾的分布情況；③ 數據存在嚴重異方差；④ 數據呈偏態分布以及方差不齊，數據方差齊性及正態性不滿足。分位數回歸的使用在連續數據[20]、分類數據[2,21]、生存數據[22]方面均有涉及，因其應用條件更為寬松，涉及數據不需要滿足正態分布，故連續數據多為非正態分布數據，也用于正態分布數據；分類數據中二分類數據、有序多分類數據、無序多分類數據均有所涉及；離散數據也有涉及，也可將離散變量轉化為連續變量進行分析，本文所選實例中年齡變量即為離散變量。分位數回歸解決的問題多以線性相關關系為主，較少涉及非線性相關關系[23]。
當臨床醫學研究樣本量足夠大可以代表整體時，研究者對于該樣本數據可直接進行描述性統計，數據為正態分布時，可以直接使用方差分析進行一般線性回歸。樣本數量越少，統計模型的設計就越困難，對醫學統計方法要求就越高，研究者進行假設檢驗就越難，越不符合真實世界，這種情況下研究者可以使用貝葉斯方法進行反復抽樣將其抽成真實世界，以提高研究準確性。然而即使貝葉斯法適合于小樣本醫學臨床數據。若樣本量過小，即使抽樣迭代很多次也可能出現馬爾科夫鏈不收斂、抽樣序列不集中的情況，這樣得到的參數的估計值就會較偏離真實值，估計的不確定性較大。
綜上所述，本文實例分析展現了貝葉斯分位數回歸方法參數估計結果精確，可信程度較高，在小樣本條件下也可以得到可靠參數信息，并且展現了曲線擬合在非線性關系中相對較于直線擬合的優勢。貝葉斯分位數回歸在臨床醫學中的應用較少，更多體現在理論層面，還需要更多高質量的案例對貝葉斯分位數回歸法的優勢進行應用、分析和總結。本文通過介紹貝葉斯分位數回歸方法并使用R Studio軟件進行實例分析，展現了貝葉斯分位數回歸在臨床醫學研究的應用價值，并提供具體代碼編寫，為臨床醫學數據研究者統計方法的使用提供參考。
