基于Maximin方法的異質(zhì)性數(shù)據(jù)分位數(shù)回歸模型研究

蔡 超 黃聰聰

（山東工商學院 統(tǒng)計學院，山東 煙臺 264005）

一、引言

Fan J和Han F等（2014）[1]認為大數(shù)據(jù)除了具有海量化這一重要特征之外，異質(zhì)性也是大數(shù)據(jù)的重要特征，即大數(shù)據(jù)通常由具有多個來源的數(shù)據(jù)組合而成，不同來源的數(shù)據(jù)呈現(xiàn)不同的影響模式。譬如在空氣質(zhì)量研究中，影響空氣質(zhì)量的氣象因素影響模式在不同季節(jié)存在差異。如果在建立空氣質(zhì)量影響因素模型時忽略上述異質(zhì)性影響，可能會導致模型估計出現(xiàn)錯誤。特別是當多個來源的數(shù)據(jù)混合在一起無法分離時，建模和估計都具有一定難度。因此，研究具有多個來源的異質(zhì)性數(shù)據(jù)的建模問題有重要的理論意義和實踐價值。

針對具有多個來源的異質(zhì)性數(shù)據(jù)的建模研究，PinheiroJ C和Bates DM（2001）[2]提出了混合效應模型，假定變量的回歸系數(shù)為常系數(shù)，通過引入可加形式的固定和隨機效應來刻畫數(shù)據(jù)的異質(zhì)性。但有研究表明上述假定過于嚴格，往往與現(xiàn)實不符。因此，Hastie T和Tibshirani R（1993）[3]、Zhang W 和 Fan J（1999）[4]、Cai Z和Fan J等（2000）[5]等提出了變系數(shù)模型，假定模型中的參數(shù)與變量之間具有函數(shù)關(guān)系，數(shù)據(jù)的異質(zhì)性通過這種函數(shù)關(guān)系來刻畫。但變系數(shù)模型過多的待估參數(shù)會損失自由度。鑒于上述問題，Lin CC和 NgS（2014）[6]、Su L和Shi Z等（2016）[7]提出了組異質(zhì)性回歸模型，通過對不同來源的數(shù)據(jù)設(shè)定不同的回歸系數(shù)刻畫數(shù)據(jù)的異質(zhì)性。馬雙鴿和王小燕等（2015）[8]、斯介生和李揚等（2017）[9]、方匡南和趙夢巒（2018）[10]提出了整合分析模型，通過融合不同來源數(shù)據(jù)，同時求解多個不同來源數(shù)據(jù)的回歸模型。然而，上述模型復雜性高，計算量很大。Meinshausen N 和 Bü hlmnn P（2015）[11]提出了均值回歸的Maximin估計，即最大化所有來源數(shù)據(jù)最小的可解釋方差。這個模型本質(zhì)上是尋找一個簡單模型來提取不同來源數(shù)據(jù)的共同屬性，保證整個數(shù)據(jù)集上都有很好的預測精度。Rothenhausler D和Meinshausen N等（2016）[12]理論證明了均值回歸Maximin估計的漸近性質(zhì)，并給出了其置信區(qū)間。秦磊和夏傳信等（2018）[13]將Maximin估計推廣到廣義線性模型，并且利用模擬數(shù)據(jù)和應用研究說明了這種方法具有較好的預測效果。

Meinshausen N 和 Bü hlmnn P（2015）[11]的研究工作是建立在均值框架下，僅能預測響應變量的條件均值，無法揭示響應變量條件分布的變化規(guī)律。Koenker R和Bassett G（1978）[14]提出的分位數(shù)回歸模型，不僅能夠獲得響應變量關(guān)于給定解釋變量的多個條件分位數(shù)，而且能夠不受異方差的限制，得到更加穩(wěn)健的結(jié)果。因此，本文將Maximin估計推廣到分位數(shù)回歸模型，構(gòu)建異質(zhì)性數(shù)據(jù)的分位數(shù)回歸估計方法：分位數(shù)回歸的Maximin估計（Maximin Estimator of Quantile Regression，Maximin QR），給出了其數(shù)學表示、參數(shù)估計、模型檢驗與預測方法，并通過數(shù)值模擬檢驗Maximin QR方法的預測效果，最后將其應用于北京地區(qū)PM2.5的條件密度預測研究。

二、模型與方法

（一）模型表示

（二）參數(shù)估計

（三）系數(shù)檢驗

（四）模型預測

三、數(shù)值模擬

本節(jié)檢驗Maximin QR方法的實際表現(xiàn)，并將其與傳統(tǒng)分位數(shù)回歸方法（以下簡稱Classical QR）、傳統(tǒng)的均值回歸方法（以下簡稱Classical MR）和Meinshausen N和 Bü hlmnn P（2015）[11]提出的 Maximin MR 方法進行對比，數(shù)值模擬在R中編程實現(xiàn)。

（一）數(shù)據(jù)生成

考慮誤差項為獨立同分布的模型：

和誤差項為非獨立同分布的模型

式中，ng=104，xig∈ i5由標準正態(tài)分布生成，xig，1為 xig的第一個元素。設(shè)置隨機誤差項εig三種不同的分布：N（0，1）、t（3）和χ2（3）。

在誤差項的不同分布情形下，數(shù)據(jù)有8種來源，即G=8，回歸系數(shù)的取值如表1所示。由表1可知，異質(zhì)性通過對不同來源數(shù)據(jù)設(shè)置不同的回歸系數(shù)來刻畫，而且每組來源數(shù)據(jù)具有共同結(jié)構(gòu) β0=（1，0，0，0，0）T。為檢驗Maximin QR方法的預測效果，選取前5組來源數(shù)據(jù)作為樣本內(nèi)數(shù)據(jù)進行估計，后3組來源數(shù)據(jù)作為樣本外數(shù)據(jù)進行預測。這樣，樣本外數(shù)據(jù)會含有樣本內(nèi)數(shù)據(jù)中的某些結(jié)構(gòu)，同時含有樣本內(nèi)數(shù)據(jù)中沒有的結(jié)構(gòu)。

表1 回歸系數(shù)取值

（二）評價指標

選取平均絕對誤差（Mean Absolute Deviations，MAD）和分位數(shù)平均絕對誤差（Quantile Average Absolute Error，QAAE）來評價模型不同分位點處的預測能力，其定義如下：

（三）結(jié)果比較

圖1報告了在不同誤差分布情形下Classical QR和Maximin QR方法預測誤差MAD的箱線圖（由于各分位點的結(jié)果類似，因此只報告了τ=0.5時的結(jié)果）。為了方便對比，Classical MR和Maximin MR方法預測誤差的箱線圖也列在圖1中。圖2報告了在不同誤差分布情形下Classical QR和Maximin QR方法預測誤差QAAE的箱線圖（由于QAAE依賴于分位點τ，因此Classical MR和Maximin MR方法無法計算QAAE的值）。需要說明的是，預測誤差均是基于100次重復生成的數(shù)據(jù)集計算獲得的，圖1和圖2中的散點是100次重復試驗的預測誤差。

圖1 預測誤差MAD的箱線圖

圖2 預測誤差QAAE的箱線圖

由圖1和圖2可知：第一，在不同誤差分布情形下，4種方法的預測結(jié)果表現(xiàn)基本相同，這表明數(shù)值模擬的結(jié)果是穩(wěn)定的。第二，在各分位點處，相對于Classical QR方法，Maximin QR具有較好的預測結(jié)果，表現(xiàn)為更低的MAD和QAAE值，因此可以認為Maximin QR方法在預測上更具優(yōu)勢。第三，與Maximin MR方法相比，中位點處的Maximin QR方法預測的MAD值都低于Maximin MR方法的預測值，意味著Maximin QR方法比Maximin MR方法更能得到準確的預測結(jié)果。此外，當誤差項為非對稱分布或異方差時，Maximin QR方法與Maximin MR方法預測誤差的差異更為明顯，這一結(jié)論表明，Maximin QR方法在誤差項為非對稱分布或異方差時的預測能力更具優(yōu)勢。

四、應用研究

（一）數(shù)據(jù)選取與描述

PM2.5（細懸浮顆粒物）是霧霾天氣的主要成因。它不僅能夠形成霧霾，影響大氣能見度，而且因其顆粒小能攜帶有毒有害物質(zhì)，進入人體會影響身體健康。PM2.5一般受到污染物本地排放和氣象條件的影響，深入研究各種因素與PM2.5之間的關(guān)系有助于理解PM2.5分布的特點和揭示PM2.5聚集和消散過程背后的規(guī)律。由于秋冬季節(jié)是霧霾高發(fā)期，而且冬季是北方地區(qū)的供暖時期，燃煤污染更易引發(fā)霧霾天氣。因此，可以認為各因素對PM2.5的影響模式在每個月份都不一致，即不同月份的數(shù)據(jù)具有不同的來源，整個數(shù)據(jù)集是異質(zhì)性的。如果將所有數(shù)據(jù)進行合并分析，會忽略各月份數(shù)據(jù)的差異性；如果分月份單獨建模分析，則會忽略各月份數(shù)據(jù)的關(guān)聯(lián)性。基于此，本文考慮數(shù)據(jù)集間的關(guān)聯(lián)性和異質(zhì)性，運用Maximin QR方法進行估計，并進行條件密度預測。

以加州大學歐文分校機器學習資源庫（UCI Machine Learning Repository）中的北京PM2.5數(shù)據(jù)集為研究對象（關(guān)于數(shù)據(jù)集更詳細的信息請見Liang X和Zou T等（2015）[16]），尋求PM2.5與影響因素之間的關(guān)系。剔除缺失值后該數(shù)據(jù)集包含41757個觀測值，以日為單位記錄了2010年1月1日至2014年12月31日的PM2.5、露點溫度、溫度、氣壓、風速、每小時降雪量、每小時降水量、是否西北風、是否東南風、是否東北風等10個變量。數(shù)據(jù)的描述統(tǒng)計結(jié)果如表2所示。由表2知，PM2.5的偏度系數(shù)大于0且峰度系數(shù)大于0，說明PM2.5為右偏尖峰分布，表明運用均值回歸模型分析各因素對PM2.5的影響并不奏效，需要利用分位數(shù)回歸模型研究各因素與PM2.5之間的關(guān)系。

表2 北京PM2.5數(shù)據(jù)集描述統(tǒng)計

（二）模型預測比較

為了驗證Maximin QR方法的預測性能，在北京PM2.5數(shù)據(jù)集上隨機抽取9個月的數(shù)據(jù)作為樣本內(nèi)數(shù)據(jù)進行估計，其余3個月的數(shù)據(jù)作為樣本外數(shù)據(jù)進行預測，以上過程重復試驗100次。表3報告了100次重復試驗的平均預測誤差，圖3報告了100次重復試驗預測誤差的箱線圖。由表3和圖3可知，在各分位點處，Maximin QR方法的預測誤差最小，且在中分點處小于Maximin MR方法的預測誤差。綜合來看，Maximin QR方法的預測效果優(yōu)于傳統(tǒng)的分位數(shù)回歸方法和Maximin MR方法。

表3 北京PM2.5數(shù)據(jù)集的平均預測誤差

圖3 北京PM2.5數(shù)據(jù)集預測誤差的箱線圖

（三）模型估計結(jié)果

對全部數(shù)據(jù)進行Maximin QR估計，估計結(jié)果如表4所示。為方便對比，Maximin MR估計的結(jié)果也列在表4中。由表4的估計結(jié)果可得到如下結(jié)論：

1.露點溫度對PM2.5的增加有顯著影響，且隨著分位點的提高，系數(shù)值不斷增大。這說明露點溫度對PM2.5具有聚集作用，而且PM2.5越高時，露點溫度使PM2.5聚集的作用越強。這主要由于PM2.5越高越能吸附更多的水汽和顆粒物，而露點溫度增加更使PM2.5聚集，從而形成惡性循環(huán)。

表4 系數(shù)估計結(jié)果

2.溫度、氣壓、風速、每小時降雪量和每小時降雨量對PM2.5的減少有顯著影響，且隨著分位點的提高，系數(shù)值不斷減小。這表明溫度、氣壓、風速、每小時降雪量和每小時降雨量等氣象因素增加時，會使PM2.5消散，且PM2.5越高時，這些氣象因素使PM2.5消散的越多。這主要由于PM2.5越高越有較多的顆粒物聚集在一起，這些氣象因素就消散得越多。

3.東北風和西北風對PM2.5具有負向影響，且隨著分位點的提高，系數(shù)值不斷減小。而東南風對PM2.5具有正向影響，且隨著分位點的提高，系數(shù)值不斷減小。這意味著北風使PM2.5有下降趨勢，而南風使PM2.5有上升趨勢。這主要是由于北京以北是太行山脈和燕山山脈，污染工業(yè)較少，北風帶來的是相對潔凈的空氣；而北京的東南方向，廣泛分布著消耗大量煤炭的重工業(yè)企業(yè)，東南風會把北京以南的污染物傳送到北京。

（四）條件密度預測

首先，分位點在0.01和0.99之間每隔0.01連續(xù)取值，共設(shè)置99個分位點；其次，選取全部數(shù)據(jù)利用Maximin QR方法估計99個分位點處的回歸系數(shù)；最后，分別選取一個解釋變量的較低水平（低）、中等水平（中）和較高水平（高）（如果解釋變量是離散變量，則選取0（否）和1（是）），其他解釋變量取值不變，連續(xù)變量取其均值，離散變量取其眾數(shù)，計算PM2.5在各分位點處的條件分位數(shù)預測值，進而討論PM2.5的條件密度變化。各個解釋變量不同水平的具體取值見表5。圖4給出了9個解釋變量分別變化時PM2.5的條件密度預測與條件均值預測，其中基于Maximin QR方法的條件密度預測用曲線表示，基于Maximin MR方法的條件均值預測用垂線表示。表5報告了PM2.5條件密度預測的描述統(tǒng)計量：均值、標準差、偏度和峰度。為便于比較，條件均值預測的結(jié)果也列于表5中。由表5可知，第一，基于Maximin MR方法的預測只能給出響應變量的一個條件均值水平預測結(jié)果，而基于Maximin QR方法的預測能夠得到響應變量整個條件分布情況，能夠獲取更為全面的信息；第二，PM2.5條件密度預測的偏度都大于0，即都呈現(xiàn)右偏狀態(tài)，表明預測的PM2.5存在非對稱性；第三，PM2.5條件密度預測的峰度都大于0，即預測的PM2.5呈尖峰分布，表明PM2.5的預測值較為密集地分布在眾數(shù)的周圍，預測其眾數(shù)可能更為準確。

表5 條件密度預測的描述統(tǒng)計

圖4 PM2.5的條件密度預測

由圖4可知，當露點溫度增加時，條件密度曲線向右移動，散布逐漸變大，這表明露點溫度對PM2.5具有正向影響，即露點溫度越高，PM2.5越大；當溫度、氣壓、風速、每小時降雪量和每小時降雨量分別增加時，條件密度曲線向左移動，散布逐漸變小，這表明溫度、氣壓、風速、每小時降雪量和每小時降雨量對PM2.5具有負向影響且逐漸向其中心值集中，即這些天氣因素增大時，PM2.5會降低。當東北風和西北風從無到有時，條件密度曲線向左移動，散布逐漸變小，而當東南風從無到有時，條件密度曲線向右移動，散布逐漸變大，這意味著北風使北京的PM2.5有下降趨勢，而南風使PM2.5有上升趨勢。條件均值預測結(jié)果也印證了各因素對PM2.5的上述影響，但遠沒有條件密度預測提供的信息豐富。

五、結(jié)論與啟示

針對具有多個來源的異質(zhì)性數(shù)據(jù)，本文提出了分位數(shù)回歸的Maximin估計方法，并給出了其數(shù)學表示、參數(shù)估計、模型檢驗與預測方法。它的基本思想是最大化所有來源數(shù)據(jù)的最小可解釋殘差，構(gòu)建一個簡單的共性模型，以減少數(shù)據(jù)來源較多而呈現(xiàn)的復雜性。數(shù)值模擬的結(jié)果顯示：Maximin QR方法與傳統(tǒng)的分位數(shù)回歸方法和Maximin MR方法相比，更能獲得精確的預測結(jié)果，證明了該估計方法的正確性和有效性。最后，將Maximin QR方法應用于北京PM2.5數(shù)據(jù)集，研究各因素對PM2.5的影響，并在此基礎(chǔ)上給出PM2.5的條件密度預測，結(jié)果表明：Maximin QR方法不僅能夠比傳統(tǒng)分位數(shù)回歸方法和Maximin MR方法更好地預測PM2.5，而且通過條件密度預測曲線可知，PM2.5的條件分布呈右偏尖峰分布，預測其眾數(shù)可能比預測其均值更為有效。針對PM2.5條件密度預測所呈現(xiàn)的分布特征，地方政府可以根據(jù)氣象因素的不同水平準確預測PM2.5的未來走向，制定和調(diào)整本地空氣質(zhì)量治理政策。

本文提出的方法適用于數(shù)據(jù)來源已知的情況，但當數(shù)據(jù)來源有多個且未知時，本文提出的方法將不再適用。在未來的研究過程中，可以考慮將本文的估計方法延伸到數(shù)據(jù)來源未知的異質(zhì)性數(shù)據(jù)。例如，首先用交叉驗證等準則確定數(shù)據(jù)來源的組數(shù)，然后，按照無放回抽樣將數(shù)據(jù)平分為G組，最后，分組后的數(shù)據(jù)按照數(shù)據(jù)來源已知的方式進行處理。