在R語言中實(shí)現(xiàn)Bayes方法對logistic的回歸分析

田 薇, 李曉毅, 付志慧

(沈陽師范大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院, 沈陽 110034)

?

統(tǒng)計(jì)學(xué)

在R語言中實(shí)現(xiàn)Bayes方法對logistic的回歸分析

田 薇, 李曉毅, 付志慧

(沈陽師范大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院, 沈陽 110034)

對于logistic回歸分析的處理辦法,一直采用的都是極大似然估計(jì)的EM算法,由于計(jì)算方法的固定及計(jì)算過程的復(fù)雜性,例如,該算法對于初值的選取要求很高,否則收斂速度很慢。Gibbs抽樣法作為一種高效靈活的估計(jì)方法廣泛應(yīng)用于廣義線性回歸模型,其中Probit回歸模型由于聯(lián)系函數(shù)為正態(tài)分布,使得回歸系數(shù)的后驗(yàn)分布為共軛正態(tài),從而抽樣簡單快捷。而Logit模型的后驗(yàn)分布比較復(fù)雜,無法直接抽取。本文基于增加數(shù)據(jù)的Gibbs抽樣方法,通過引入Plya-Gamma分布族的潛在變量,使得模型中的回歸系數(shù)參數(shù)的滿條件分布為共軛正態(tài)分布,從而回歸系數(shù)的馬氏鏈很容易構(gòu)造,回歸系數(shù)的估計(jì)為后驗(yàn)均值估計(jì)。通過一組實(shí)際數(shù)據(jù),分別調(diào)用R語言Glm包和 BayesLogit包,并對比2種方法的估計(jì)結(jié)果,二者差別不大,表明Plya-Gamma潛變量Bayes估計(jì)法在處理logistic回歸模型時(shí)的可用性、準(zhǔn)確性。

logistic回歸模型; Plya-Gamma分布; MCMC; R語言

0 引 言

logistic回歸模型常用在尋找危險(xiǎn)因素、預(yù)測、判別3個方面。由于線性回歸模型只能分析連續(xù)型數(shù)據(jù),具有很強(qiáng)的局限性,而實(shí)際中的logistic回歸用途是極為廣泛的,logistic回歸幾乎已經(jīng)成了流行病學(xué)和醫(yī)學(xué)中最常用的分析方法,logit模型在處理屬性數(shù)據(jù)或分類數(shù)據(jù)方面極具優(yōu)勢。在估計(jì)logistic回歸模型時(shí),一般采用極大似然法。若引入合理的先驗(yàn)分布函數(shù),Bayes方法對于許多模型的參數(shù)估計(jì)問題一直特別有效,對提高統(tǒng)計(jì)推斷質(zhì)量具有實(shí)際意義。另外,忽略參數(shù)的先驗(yàn)信息,有時(shí)是一種浪費(fèi),甚至還會導(dǎo)致不合理的結(jié)論。logit回歸模型聯(lián)系函數(shù)為logistic分布函數(shù),回歸系數(shù)的后驗(yàn)分布無具體形式,需要采用Metropolis-Hastings抽樣法,該方法在應(yīng)用過程中若建議分布選取的不合理,會導(dǎo)致接受概率很小,因此馬氏鏈?zhǔn)諗亢苈?/p>

本文嘗試使用一種由Albert提出的數(shù)據(jù)添加新方法,該方法在一定程度上區(qū)別于常見的數(shù)據(jù)添加方式。引入一個來自Plya-Gamma分布的隨機(jī)變量(具體的構(gòu)造方法將在文中加以說明)。為了提高估計(jì)的質(zhì)量,除了當(dāng)前樣本數(shù)據(jù),還可以利用客觀信息和經(jīng)驗(yàn)累積的信息,先驗(yàn)信息的加入,參數(shù)估計(jì)更加穩(wěn)定,也更合理和符合實(shí)際?；貧w系數(shù)的共軛分布仍為正態(tài)分布,使得抽樣及后驗(yàn)估計(jì)很容易得到,更加方便計(jì)算。在R語言的BayesLogit程序包中,MCMC抽樣及數(shù)據(jù)處理均可得以實(shí)現(xiàn)。最后通過一組實(shí)際數(shù)據(jù),利用Glm包實(shí)現(xiàn)傳統(tǒng)似然估計(jì)方法,并與MCMC估計(jì)方法相比較。

定義1 隨機(jī)變量X是帶有參數(shù)b(b>0),和c∈R的Plya-Gamma分布,其中變量Χ分布記為X～PG(b,c), 形式為

這里gk～Ga(b,1),是獨(dú)立的伽瑪隨機(jī)變量。

其中:Vw=(XTΩX+B-1)-1;mω=Vω(XTκ+B-1b);κ=(y1-n1/2,…,yN-nN/2);Ω為ωi的對角線矩陣。

接下來開始對上面所得到的抽樣方法進(jìn)行演繹證明,先從一些定理和積分公式開始。

定理p(ω)為隨機(jī)變量ω的密度函數(shù),且ω～PG(b,0),b>0。對于所有的a∈R,有下列恒等式:

其中κ=a-b/2 。

對式子(2)的非正態(tài)化聯(lián)合密度處理,可得ω的條件分布

通過定理和積分,可以得到第i個觀測值的似然函數(shù)為

(4)

其中p(ωi|ni,0)為帶參數(shù)隨機(jī)變量(ni,0),服從Plya-Gamma分布的密度函數(shù)。

在n組數(shù)據(jù)下,β的后驗(yàn)條件分布為

從整理得到的分布形式知,β的后驗(yàn)條件分布服從正態(tài)分布,即P(β|ω,y)∝N(mω,Vω)。

其中:mω=Vω(XTκ+B-1b);Vw=(XTΩX+B-1)-1。其理論依據(jù)為正態(tài)分布(方差已知)的共軛先驗(yàn)還是正態(tài)分布。此處:z=(κ1/ω1,…,κN/ωN);Ω=diag(ω1,…,ωN)。由式(5)可知,β的條件分布為高斯似然,且先驗(yàn)p(β)也服從高斯分布,因此該線性模型得以簡單計(jì)算。

2 模擬研究

在二項(xiàng)分布族中,logistic回歸模型是最重要的模型。對于響應(yīng)變量Y有p個自變量(或稱為解釋變量),記為X1,X2,…,Xp。在p個自變量作用下出現(xiàn)成功的條件概率為P=P{Y=1|X1,X2,…,Xp},那么其logistic回歸模型可表示為

其中:稱β0為截距;稱β1,β2,…,βp為logistic回歸模型系數(shù)。

對上式作logit變換,logistic回歸模型可以變成下列線性形式:

可以使用線性回歸模型對參數(shù)進(jìn)行估計(jì),這也是logistic回歸模型屬于廣義線性模型的原因。

當(dāng)logistic回歸模型的分布函數(shù)為

農(nóng)村集體土地上不動產(chǎn)登記工作中，不動產(chǎn)登記權(quán)利人主體不一致的情況很常見，突出表現(xiàn)在以下3個方面：①規(guī)劃審批手續(xù)的建房人與土地審批手續(xù)的使用人不同；②土地使用權(quán)人與登記簿中房屋所有人不同；③房屋所有權(quán)與土地所有權(quán)人的主體不一致。

其聯(lián)合分布函數(shù)為

通常利用極大似然法,對該線性回歸模型中的參數(shù)進(jìn)行估計(jì)。

通過一組實(shí)際數(shù)據(jù),用R語言中BayesLogit包對數(shù)據(jù)進(jìn)行處理,得到spambase數(shù)據(jù)包的其他特征數(shù),詳見表1。并與原始方法得到的模型系數(shù)估計(jì)值進(jìn)行比照,詳見表2(logistic回歸模型系數(shù)估計(jì)值表)。其結(jié)果表明利用BayesLogit方法處理logistic回歸分析問題的可行性。

表1 spambase數(shù)據(jù)包其他特征數(shù)Tab.1 spambase packet number of other features

表2 logistic回歸模型系數(shù)估計(jì)值表Tab.2 logistic regression coefficient estimates table

3 結(jié) 論

[ 1 ]GAMERMAN D. Sampling from the posterior distribution in generalized linear mixed models[J]. Statistics and Computing, 1997(7):57-68.

[ 2 ]HOLMAN R, GLAS C A W. Modeling non-ignorable missing data mechanisms with item response theory models[J]. BRIT J MATH STAT PSY, 2005,58(1):1-17.

[ 3 ]HAMBLETON R K. Fundamentals of item response theory[M]. NewYork:Sage Publication, 1991.

[ 4 ]RUBIN D B. Inference and missing data[J]. Biometrika, 1976,63(3):581-592.

[ 5 ]LITTLE R J A, RUBIN D B. Statistical analysis with missing data[M]. Manhattan:John Wiley&Sons, 2014.

[ 6 ]MASTERS G N.ARasch model for partial credit scoring[J]. Psychometrika, 1982,47(2):149-174.

[ 7 ]ALBERT J H. Bayesian estimation of normal ogive item response curves using Gibbs sampling[J]. J EDUCBEHAV STAT, 1992,17(3):251-269.

[ 8 ]JONES D H, NEDIAK M S. Item parameter calibration of LSAT items using MCMC approximation of Bayes posterior distribution[M]. Newtown:Law School Admission Coucil, 2005.

[ 9 ]GELMAN A, RUBIN D B. Inference from iterative simulation using multiple sequences[J]. STAT SCI, 1992:457-472.

[10]MARIS G,BECHGER T M. An introduction to the DAT Gibbs sampler for the two-parameter logistic(2PL) model and beyond[J]. International Journal of Methodology and Experimental Psychology, 2005,26(2):327-352.

[11]LUDLOW L H, O’LEARY M. Scoring omitted and not-reached items: practical data analysis implications[J]. EDUC PSYCHOL MEAS, 1999,59(4):615-630.

[12]HUISMAN M. Imputation of missing itemresponses:Some simple techniques[J]. QUAL QUANT, 2000,34(4):331-351.

[13]MURAKI E, BOCK R D. PARSCALE:IRT based test scoring and item analysis for graded open-ended exercises and performance tasks[M]. Scientific Software International, 1993.

[14]LORD F M. Maximum likelihood and Bayesian parameter estimation in item response theory[J]. J EDUC MEAS, 1986,23(2):157-162.

[15]MOUSTAKI I, KNOTT M. Weighting for item non-response in attitude scales by using latent variable models with covariates[J]. J R STAT SOC B, 2000,163(3):445-459.

Bayesian inference for logistic models in R Language

TIAN Wei, LI Xiaoyi, FU Zhihui

(College of Mathemetics and Systems Science, Shenyang Normal University, Shenyang 110034, China)

For the approach to logistic regression analysis, using a maximum likelihood estimation are the EM, due to the complexity and fixity of calculation, for example, the initial value of the algorithm is demanding, otherwise the convergence rate is slow. Gibbs sampling as an efficient and flexible estimation is widely used for generalized linear regression models, due to the contact function is normal in Probit model, so that the posterior distribution of the regression coefficients is Conjugated Normality and sampling is easier.The posterior of Logit model is complex, unable to directly extract, based on Gibbs to increase data by introducing latent variables Plya-Gamma distribution families, making the regression coefficient parameters of full conditional distribution Conjugated Normality, thereby Markov chains regression coefficient is easy to construct the estimated regression coefficients for the posterior mean estimate. Through a set of actual data, respectively, calling R language package of BayesLogit and Glm, and comparing the results of the two methods, the difference is small, indicating Plya-Gamma latent variable Bayesian estimation in dealing with the accuracy of logistic regression model.

logistic regression model; Plya-Gamma distribution; MCMC; R language

2016-04-16。

國家自然科學(xué)基金青年基金資助項(xiàng)目(11201313)。

田 薇(1990-),女,遼寧葫蘆島人,沈陽師范大學(xué)碩士研究生; 通信作者: 李曉毅(1956-),女,遼寧葫蘆島人,沈陽師范大學(xué)教授。

1673-5862(2016)03-0321-04

O212.8

A

10.3969/ j.issn.1673-5862.2016.03.014