徐州醫學院流行病與衛生統計教研室(221002) 王 婷 曾 平 何 鵬
Probit回歸最先由Bliss在1934年提出〔1〕,在進行殺蟲劑的毒理實驗時,Bliss發現了一個有趣的現象:無論他配制的殺蟲劑濃度有多高,在用藥之后總會有一兩只昆蟲還活著;無論他怎么稀釋殺蟲劑,即便只是用裝過殺蟲劑的容器,試驗結果也總會有幾只昆蟲死掉。Bliss原創性地采用概率分布這種新的數學思想來解決殺蟲劑實驗時所遭遇的困境。Probit回歸建立了“劑量”與“使用該劑量時一只動物會死掉的概率”這兩者間的關系,因此稱為Probit(單位概率)模型,現在已經越來越多地應用到二分類數據的回歸分析中。本文將主要在貝葉斯統計框架內討論Probit回歸和參數后驗分布的潛變量Gibbs抽樣。
1.Probit回歸和后驗分布
設解釋變量為X,回歸系數向量為β,根據廣義線性模型原理〔2〕,可建立 Probit回歸:

這里假定事件發生服從參數為p的Bernoulli分布,n表示樣本量,Φ表示標準正態的累積概率分布函數,如Φ(1.96)=0.975,Φ-1為累積分布函數的逆函數,如Φ-1(0.975)=1.96,這樣通過Probit連接函數Φ-1將取值為0~1之間的p映射到了整個實數空間。似然。設 g(β)為回歸參數β的先驗分布(prior distribution),貝葉斯統計和頻率統計最大的區別之一就在于假定參數為隨機變量,當有關于未知參數的歷史知識、主觀認識或者專家意見時,可以選擇有信息先驗。當對未知參數的信息一無所知、或先驗分布有太多參數需要指定時,認為參數在其空間內有等可能的取值概率而不特別偏向某些取值,則選取均勻分布作為先驗分布,又稱貝葉斯假定。根 據 貝 葉 斯 原 理, 后 驗 分 布〔3〕(posterior distribution)p(β|Y,X)為:


3.潛變量Probit回歸的Gibbs抽樣
設存在一個潛在的連續變量 y*,y*稱為特征(trait)或傾向得分(propensity score)〔5〕。考慮以下的模型:均數(-1.5897)和中位數(-1.5018)來看,男孩比女 孩更易發胖。


圖1 腰圍后驗樣本的直方圖、核密度圖、軌跡圖和自相關圖

表1 參數的后驗樣本描述
本文討論了醫學領域中二分類數據分析的貝葉斯Probit回歸,像絕大部分貝葉斯模型一樣,貝葉斯Probit回歸參數的后驗分布異常復雜,需要采用MCMC模擬抽樣,Gibbs抽樣是眾多MCMC算法中最常用的模擬方法。Gibbs抽樣中需要已知參數的滿條件分布,在此基礎上迭代抽樣生成參數的Markov鏈,但從Probit回歸的后驗分布卻不能得到簡單并且抽樣方便的滿條件分布,因此執行Gibbs抽樣也就不大可行,本文通過增加潛在變量解決了這個問題。
潛在變量并不能被觀察到,因此在貝葉斯統計中當作未知量看待,則此時Probit回歸的后驗分布為g(β,Y*|Y),是回歸參數和潛變量的聯合密度函數。在這個后驗分布中,β,Y*各自的滿條件分布分別是g(β|Y*,Y)和 g(Y*|β,Y),在 Probit回歸中,前者為多元正態分布,后者為截尾正態分布,這兩個滿條件分布都比較簡單而且容易直接進行模擬抽樣,因此執行Gibbs抽樣也就沒有困難。在給定回歸參數初始值后在這兩個滿條件分布之間反復迭代生成參數和潛變量的Markov鏈,在算法收斂后則可認為生成的參數序列來自Probit回歸后驗分布。構造潛變量的Gibbs抽樣可以看作是一種數據擴增技術,通過在模型的后驗分布中增加輔助變量使得Gibbs抽樣更加容易。另一個好處是,生成的潛在變量向量可以進一步作為模型診斷的信息加以利用。
1.Salsburg D.The lady tasting tea:how statistics revolutionized science in the twentieth century.Holt Paperbacks,2002.
2.Dobson AJ,Barnett A.An introduction to generalized linear models,third edition.Chapman & Hall,2009.
3.Gelman A,Carlin JB,Stern HS,et al.Bayesian data analysis,2nd ed.London:Chapman & Hall,2004.
4.Albert J.Bayesian computation with R,2nd Ed.New York:Springer,2009.
5.Lynch SM.Introduction to applied bayesian.Statistics and Estimation for Social Scientists,New York:Springer,2009.
6.http://www.r-project.org/