基于Maximin方法的異質性數據分位數回歸模型研究

蔡 超 黃聰聰

（山東工商學院 統計學院，山東 煙臺 264005）

一、引言

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

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

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

二、模型與方法

（一）模型表示

（二）參數估計

（三）系數檢驗

（四）模型預測

三、數值模擬

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

（一）數據生成

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

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

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

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

表1 回歸系數取值

（二）評價指標

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

（三）結果比較

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

圖1 預測誤差MAD的箱線圖

圖2 預測誤差QAAE的箱線圖

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

四、應用研究

（一）數據選取與描述

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

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

表2 北京PM2.5數據集描述統計

（二）模型預測比較

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

表3 北京PM2.5數據集的平均預測誤差

圖3 北京PM2.5數據集預測誤差的箱線圖

（三）模型估計結果

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

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

表4 系數估計結果

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

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

（四）條件密度預測

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

表5 條件密度預測的描述統計

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

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

五、結論與啟示

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

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