田 薇, 李曉毅, 付志慧

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

?

統(tǒng)計學

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

田 薇, 李曉毅, 付志慧

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

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

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

0 引 言

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

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

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

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

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

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

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

其中κ=a-b/2 。

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

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

(4)

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

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

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

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

2 模擬研究

在二項分布族中,logistic回歸模型是最重要的模型。對于響應變量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ù)進行估計,這也是logistic回歸模型屬于廣義線性模型的原因。

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

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

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

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

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

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

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

3 結 論

[ 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。

國家自然科學基金青年基金資助項目(11201313)。

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

1673-5862(2016)03-0321-04

O212.8

A

10.3969/ j.issn.1673-5862.2016.03.014