基于R語言的互補(bǔ)雙對(duì)數(shù)模型分析

戴建國(guó) 楊劍紅

摘 要： 針對(duì)分類數(shù)據(jù)模型的誤差不滿足正態(tài)分布，Logistic分布等常用分布，而是滿足極值分布這類特殊情況時(shí)，采用互補(bǔ)雙對(duì)數(shù)模型進(jìn)行分析。以R語言作為分析工具，并將其與幾種常用模型進(jìn)行對(duì)比分析，結(jié)果表明互補(bǔ)雙對(duì)數(shù)模型優(yōu)于其他模型。

關(guān)鍵詞： 互補(bǔ)雙對(duì)數(shù)模型； R語言； 有序分類

中圖分類號(hào)： O 212.1 文獻(xiàn)標(biāo)志碼： A 文章編號(hào)： 1671-2153（2017）04-0087-03

1 問題提出

大數(shù)據(jù)時(shí)代，數(shù)據(jù)挖掘與數(shù)據(jù)分析過程中遇到的數(shù)據(jù)類型主要有連續(xù)型，離散型，以及混合型，其涉及的方法眾多，如機(jī)器學(xué)習(xí)，深度學(xué)習(xí)等。當(dāng)然，對(duì)于離散型數(shù)據(jù)，回歸模型也是比較流行的統(tǒng)計(jì)分析方法之一，而且在多數(shù)情況下，因變量有可能是有序二分類或多分類的離散型變量，在醫(yī)療衛(wèi)生統(tǒng)計(jì)或社會(huì)統(tǒng)計(jì)尤為常見，如對(duì)某事的評(píng)價(jià)可能分為不好、好、非常好這樣三個(gè)等級(jí)。對(duì)于這類數(shù)據(jù)常用的模型有Logistic回歸模型（包括累積Logistic，鄰近Logistic等模型），Probit回歸模型，但它們分別是假設(shè)模型誤差分布是Logistic分布和正態(tài)分布的，而有些情況誤差并不滿足這樣的假設(shè)，其誤差分布是非對(duì)稱的極值分布或稱Gumbel 分布[1]。如某年某河面最高水位這樣一個(gè)極值，這時(shí)需要尋求一個(gè)更為有效的模型，此時(shí)采用互補(bǔ)雙對(duì)數(shù)線性模型與其他幾種模型進(jìn)行對(duì)比分析。并且近年來R語言已成為當(dāng)下最流行做數(shù)據(jù)分析和挖掘的工具之一，它不僅有豐富的函數(shù)包，而且能與Hadoop，Python等軟件結(jié)合，使得成為重要的數(shù)據(jù)分析工具。下面基于R語言用互補(bǔ)雙對(duì)數(shù)模型對(duì)一組實(shí)際數(shù)據(jù)進(jìn)行對(duì)比分析。

2 互補(bǔ)雙對(duì)數(shù)線性模型[1]

該模型的一般形式為

log{-log[1-P（Y≤j）]}=αj+β'X， （1）

從而可得求概率的表達(dá)式為

式中：αj為解釋變量不能解釋的部分；β'為參數(shù)構(gòu)成的向量。當(dāng)只有一個(gè)自變量，并且當(dāng)X=xi（i=1，2）時(shí)，有性質(zhì)：

P（Y>j|x1）=P（Y>j|x2）exp[β（x1-x2）]。

對(duì)于各自變量的回歸系數(shù)通常采用極大似然法來估計(jì)。由于互補(bǔ)雙對(duì)數(shù)函數(shù)是嚴(yán)格凹的，海森矩陣便是負(fù)的，因此通過牛頓—拉夫森方法[2]得到參數(shù)估計(jì)，對(duì)參數(shù)進(jìn)行檢驗(yàn)，模型擬合優(yōu)度檢驗(yàn)[3]，這些都可以由R語言的summary（）函數(shù)得到。

3 實(shí)例分析與R程序

下面實(shí)例是研究膚色和性別對(duì)壽命的影響，其中壽命作為因變量分為5個(gè)類別（0～20歲、20～40歲、40～50歲、50～60歲、>60歲），分別記為1，2，3，4，5；性別和膚色作為兩個(gè)自變量，記為X，Z。對(duì)于自變量X包含男性、女性兩個(gè)類別，分別用0，1表示；自變量Z包含白人、黑人兩個(gè)類別，同樣用0，1表示。且每行的總數(shù)均為1000，數(shù)據(jù)見表1。由表1可以看出，在壽命的給類別上分布是很不平衡的，如比較壽命分別對(duì)應(yīng)0～20和>60的兩列，下面給出幾種模型的對(duì)比分析。

3.1 建立數(shù)據(jù)文件

首先將表1中數(shù)據(jù)的每一列構(gòu)成因子，然后用cbind（）函數(shù)將其按列合并，最后用as.data.frame（）將合并后的數(shù)據(jù)轉(zhuǎn)化成數(shù)據(jù)框，并命名為w，程序如下：

>gender<-c（0， 0， 1， 1）； race<-c（0， 1， 0， 1）；y1<-c（13， 26， 9， 18）；y2<-c（28， 49， 13， 24）；y3<-c（32， 56， 19， 37）；>y4<-c（122， 201， 80， 129）； y5<-c（805， 668， 879， 792）

>w<-as.data.frame（cbind（gender，race，y1，y2，y3，y4，y5））

3.2 安裝和載入VGAM包

在數(shù)據(jù)文件建立后，需要安裝和載入所需要的函數(shù)包，以下用一種簡(jiǎn)單有效的方法載入和安裝所需的VGAM包。

>install.package（“VGAM”）； library（VGAM）

3.3 模型的構(gòu)建及預(yù)測(cè)

在完成前面的準(zhǔn)備工作后，使用VGAM包中的vglm（）函數(shù)對(duì)表中的數(shù)據(jù)用上述模型和其他兩種常用模型進(jìn)行分析及預(yù)測(cè)。首先構(gòu)建互補(bǔ)雙對(duì)數(shù)模型，模型的表達(dá)式已在式（1）給出。以下代碼中設(shè)置link = cloglog表示互補(bǔ)雙對(duì)數(shù)，將構(gòu)建的模型記為CLL_fit，可用summary（CLL_fit）得到模型參數(shù)的估計(jì)和檢驗(yàn)，并用predict（）所構(gòu)建模型對(duì)數(shù)據(jù)進(jìn)行預(yù)測(cè)，用round（）結(jié)果保留一位小數(shù)，預(yù)測(cè)結(jié)果見表1中括號(hào)內(nèi)的數(shù)據(jù)。相關(guān)代碼如下：

>CLL_fit<-vglm（cbind（y1，y2，y3，y4，y5）～gender+race，family = cumulative（link = cloglog，parallel = T），data = w）

>summary（CLL_fit）

>fitted<-predict（CLL_fit，type = "response"）；round（fitted*1000，1）

3.3.1 鄰近logit回歸模型

該模型的一般形式為

式中：c為因變量的類別數(shù)；X為解釋變量構(gòu)成的向量；πj為類別j在給定解釋變量下對(duì)應(yīng)的概率；βj'為對(duì)應(yīng)的回歸系數(shù)，該模型的假設(shè)是模型誤差服從Logistic分布。同樣使用VGAM包中的vglm（）函數(shù)對(duì)表中的數(shù)據(jù)進(jìn)行分析及預(yù)測(cè)，此時(shí)需設(shè)置family=acat（），將構(gòu)建后的模型記為AL_fit，由predict（）命令得到預(yù)測(cè)結(jié)果見表2。代碼如下：

>AL_fit<-vglm（cbind（y1，y2，y3，y4，y5）～gender+race，family=acat（reverse=TRUE，parallel=FALSE）， data=w）

>fitted<-predict（AL_fit，type ="response"）；round（fitted*1000，1）

3.3.2 累積Probit模型

模型的一般形式為

P（Y≤j）=？椎（αj+β'X） （j=1，…，c-1）， （3）

式中：c為因變量的類別數(shù)；X為解釋變量構(gòu)成的向量；β'為對(duì)應(yīng)的回歸系數(shù)；？椎為正態(tài)分布的分布函數(shù)；P（Y≤j）為累積概率，該模型的假設(shè)是模型誤差服從正態(tài)分布?？赏パa(bǔ)雙對(duì)數(shù)模型一樣用代碼構(gòu)建模型，只需將連接函數(shù)設(shè)為probit，即link=probit，同樣，得到的擬合值如表2所示，代碼如下：

>fit.cprobit<-vglm（cbind（y1，y2，y3，y4，y5）～gender+race，family=cumulative（link=probit，parallel=TRUE）， data=w）

>fitted<-predict（fit.cprobit，type = "response"）；round（fitted*1000，1）

3.4 模型擬合評(píng)價(jià)

本文分別用BIC和AIC準(zhǔn)則[4-5]對(duì)上述模型進(jìn)行評(píng)價(jià)，對(duì)于具體的公式及證明在此不討論，有興趣者可參考相關(guān)文獻(xiàn)。當(dāng)模型的BIC和AIC值越小，說明模型擬合的越好。下面分別給出代碼計(jì)算各模型的AIC和BIC值，結(jié)果如表3所示，代碼如下：

>AIC（CLL_fit）；BIC（CLL_fit）；AIC（AL_fit）；BIC（AL_fit）；AIC（fit.cprobit）；BIC（fit.cprobit）

綜合上面代碼得到表2和表3的結(jié)果，比較發(fā)現(xiàn)：表1中模型（1）的擬合數(shù)據(jù)比表2中另外兩種模型的擬合值更接近真實(shí)值，尤其是在年齡大于60的那一列，而且比較表3中三個(gè)模型的AIC和BIC值可知模型（1）的值明顯小于另外兩種模型，因此，表明對(duì)于這樣一類特殊數(shù)據(jù)，互補(bǔ)雙對(duì)數(shù)模型具有明顯的優(yōu)勢(shì)。

由互補(bǔ)雙對(duì)數(shù)模型擬合的輸出結(jié)果，可以得到4個(gè)回歸方程分別為

log{-log[1-P（Y≤1）]}=-4.21-0.54X+0.61Z，

log{-log[1-P（Y≤2）]}=-3.19-0.54X+0.61Z，

log{-log[1-P（Y≤3）]}=-2.58-0.54X+0.61Z，

log{-log[1-P（Y≤4）]}=-1.52-0.54X+0.61Z。

由互補(bǔ)雙對(duì)數(shù)模型給出的性質(zhì)有P（Y>j|X=1，Z）=P（Y>j|X=0，Z）exp（0.54），即在給定膚色的情況下，男性壽命長(zhǎng)于固定時(shí)間的比例是女性比例的exp（0.54）=1.71次方；在給定性別時(shí)，黑人壽命長(zhǎng)于某一固定時(shí)間的比例是白人比例的exp（0.61）=1.84次方，并且可以得到對(duì)于壽命長(zhǎng)于某一固定時(shí)間，黑種男人的比例是白種女人比例的3.55次方。

4 結(jié)束語

對(duì)于互補(bǔ)雙對(duì)數(shù)線性模型主要用于誤差分布是極值分布的情況，如上例所示，其可以對(duì)自變量進(jìn)行選擇，也可以估計(jì)概率的大小，從而得到擬合值，當(dāng)因變量是二分變量時(shí)，可以用于預(yù)測(cè)分類，在文獻(xiàn)[6-7]中給出了相關(guān)實(shí)例分析，說明了在誤差分布非對(duì)稱時(shí)，互補(bǔ)雙對(duì)數(shù)模型體現(xiàn)一定的優(yōu)勢(shì)，而且R語言能夠很有效，簡(jiǎn)潔完成分析過程，從而也說明R語言是一款很好的統(tǒng)計(jì)分析工具。

參考文獻(xiàn)：

[1] AGRESTI，ALAN. Analysis Of Ordinal Categorical Data，Second Edition[M]. Canada：Wiley，2010.

[2] BOLARINWA，OLUBUSOYE O E. On the sensitivity of probit and complementary log-log models to violated tolerance[J]. 2013，12：351-360.

[3] AGRESTI，ALAN. An Introduction to Categorical Data Analysis[M]. Wiley，New York，Duxbury Press，1996.

[4] BURNHAM K P，ANDERSON D R. Multimodel Inference understanding AIC and BIC in model selection[J]. Sociological Methods & Research，2004，33（33）：261-304.

[5] KUHA J. AIC and BIC comparisons of assumptions and performance[J]. Sociological Methods & Research，2004，33（2）：188-229.

[6] 肖珉，周宗放. 基于互補(bǔ)雙對(duì)數(shù)模型的企業(yè)集團(tuán)財(cái)務(wù)困境預(yù)測(cè)研究[C]. 中國(guó)災(zāi)害防御協(xié)會(huì)風(fēng)險(xiǎn)分析專業(yè)委員會(huì)第四屆年會(huì)論文集，2010：976-981.

[7] 史芳芳，楊海珍，徐昭，等. 基于互補(bǔ)雙對(duì)數(shù)模型的國(guó)際資本流入風(fēng)險(xiǎn)影響因素研究[J].數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，2016，46（10）：91-96.