基于潛變量Gibbs抽樣的Probit回歸*

徐州醫(yī)學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)教研室(221002) 王 婷 曾 平 何 鵬

Probit回歸最先由Bliss在1934年提出〔1〕，在進(jìn)行殺蟲(chóng)劑的毒理實(shí)驗(yàn)時(shí)，Bliss發(fā)現(xiàn)了一個(gè)有趣的現(xiàn)象:無(wú)論他配制的殺蟲(chóng)劑濃度有多高，在用藥之后總會(huì)有一兩只昆蟲(chóng)還活著;無(wú)論他怎么稀釋殺蟲(chóng)劑，即便只是用裝過(guò)殺蟲(chóng)劑的容器，試驗(yàn)結(jié)果也總會(huì)有幾只昆蟲(chóng)死掉。Bliss原創(chuàng)性地采用概率分布這種新的數(shù)學(xué)思想來(lái)解決殺蟲(chóng)劑實(shí)驗(yàn)時(shí)所遭遇的困境。Probit回歸建立了“劑量”與“使用該劑量時(shí)一只動(dòng)物會(huì)死掉的概率”這兩者間的關(guān)系，因此稱(chēng)為Probit(單位概率)模型，現(xiàn)在已經(jīng)越來(lái)越多地應(yīng)用到二分類(lèi)數(shù)據(jù)的回歸分析中。本文將主要在貝葉斯統(tǒng)計(jì)框架內(nèi)討論P(yáng)robit回歸和參數(shù)后驗(yàn)分布的潛變量Gibbs抽樣。

方法和原理

1．Probit回歸和后驗(yàn)分布

設(shè)解釋變量為X，回歸系數(shù)向量為β，根據(jù)廣義線(xiàn)性模型原理〔2〕，可建立 Probit回歸:

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

3．潛變量Probit回歸的Gibbs抽樣

設(shè)存在一個(gè)潛在的連續(xù)變量 y*，y*稱(chēng)為特征(trait)或傾向得分(propensity score)〔5〕?？紤]以下的模型:均數(shù)(－1.5897)和中位數(shù)(－1.5018)來(lái)看，男孩比女 孩更易發(fā)胖。

圖1 腰圍后驗(yàn)樣本的直方圖、核密度圖、軌跡圖和自相關(guān)圖

表1 參數(shù)的后驗(yàn)樣本描述

討 論

本文討論了醫(yī)學(xué)領(lǐng)域中二分類(lèi)數(shù)據(jù)分析的貝葉斯Probit回歸，像絕大部分貝葉斯模型一樣，貝葉斯Probit回歸參數(shù)的后驗(yàn)分布異常復(fù)雜，需要采用MCMC模擬抽樣，Gibbs抽樣是眾多MCMC算法中最常用的模擬方法。Gibbs抽樣中需要已知參數(shù)的滿(mǎn)條件分布，在此基礎(chǔ)上迭代抽樣生成參數(shù)的Markov鏈，但從Probit回歸的后驗(yàn)分布卻不能得到簡(jiǎn)單并且抽樣方便的滿(mǎn)條件分布，因此執(zhí)行Gibbs抽樣也就不大可行，本文通過(guò)增加潛在變量解決了這個(gè)問(wèn)題。

潛在變量并不能被觀察到，因此在貝葉斯統(tǒng)計(jì)中當(dāng)作未知量看待，則此時(shí)Probit回歸的后驗(yàn)分布為g(β，Y*|Y)，是回歸參數(shù)和潛變量的聯(lián)合密度函數(shù)。在這個(gè)后驗(yàn)分布中，β，Y*各自的滿(mǎn)條件分布分別是g(β|Y*，Y)和 g(Y*|β，Y)，在 Probit回歸中，前者為多元正態(tài)分布，后者為截尾正態(tài)分布，這兩個(gè)滿(mǎn)條件分布都比較簡(jiǎn)單而且容易直接進(jìn)行模擬抽樣，因此執(zhí)行Gibbs抽樣也就沒(méi)有困難。在給定回歸參數(shù)初始值后在這兩個(gè)滿(mǎn)條件分布之間反復(fù)迭代生成參數(shù)和潛變量的Markov鏈，在算法收斂后則可認(rèn)為生成的參數(shù)序列來(lái)自Probit回歸后驗(yàn)分布。構(gòu)造潛變量的Gibbs抽樣可以看作是一種數(shù)據(jù)擴(kuò)增技術(shù)，通過(guò)在模型的后驗(yàn)分布中增加輔助變量使得Gibbs抽樣更加容易。另一個(gè)好處是，生成的潛在變量向量可以進(jìn)一步作為模型診斷的信息加以利用。

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/