基于bootstrap法的我國糧食產量回歸分析

劉 薇，常振海，張德生

(1.天水師范學院 數學與統計學院，甘肅 天水 741001；2.西安理工大學 理學院,西安 710054)

0 引言

糧食產量和糧食安全是人們一直關注的熱點,學者們從不同角度討論了糧食問題。因為能夠用于統計分析的糧食年度數據較少,導致建立起的模型統計性質缺乏可信度,bootstrap方法給我們提供了增加樣本量的便捷途徑,關于bootstrap方法的應用研究近年來也熱度不減[1～4]。將bootstrap方法用于研究糧食產量的文獻所見不多,考慮到糧食產量前后幾十年增加較快,可能存在距離均值較遠的數據,本文結合較為穩健的M-估計[5],建立起關于我國糧食產量的bootstrap法回歸模型。

1 基本知識與方法

考慮最簡單的回歸模型

在c=1.345時，它對正態數據的效率為95%。

(2)計算bootstrap法估計T*=β(S*)。

(3)重復(1)和(2)R次,本文取R=1999,這樣得到的R個T*r(r=1，…，R)。

(4)利用R個T*r得到T的近似分布,從而可以得到待估計參數β的點估計,區間估計等統計量。

很顯然,這種方法省去了對分布的假定及基于假定分布的統計推斷,利用重復抽樣來擴大樣本量,建立起基于樣本的統計分布。

2 實證分析

樣本范圍是1980～2012年共33個樣本,初步選用變量為：

x1：農業機械總動力(萬千瓦),

x2：受災面積(千公頃),

x3：國家財政用于農業的支出(億元),

x4：有效灌溉面積(千公頃),

x5：農用化肥施用折純量(萬噸),

x6：農村用電量(億千瓦小時),

x7：糧食作物播種面積(千公頃),

Y：糧食產量(萬噸),

數據來源于國家統計局網站。計算主要基于R語言。

自變量個數相對于樣本個數來說較多,都選入模型顯然不合適,為了選取對糧食產量有重要影響的因素,我們首先對這些自變量進行了降維分析,結果見圖1。

圖1 因子分析碎石圖

圖2 最大方差法因子旋轉圖

從圖1能看出,提取兩個因子,方差累計貢獻率就達到了88.95%,這樣,在這7個自變量中抽取出兩個公因子,僅損失約11%的信息,效果還是很不錯的。為了能更加清晰地顯示提取的兩個因子及其實際代表意義,我們做了最大方差法因子旋轉,結果見圖2。

從圖2中能看出,變量x1,x3,x4,x5,x6代表了第一因子,x2,x7代表了第二因子。我們結合這些變量的實際意義,很容易看出第一因子與糧食單位面積產量有關,第二因子于糧食作物播種面積有關。進一步分析,在第一因子中的變量,對糧食產量影響最大的因素是x5,在第二因子的變量中則是x7。因此,我們分別以這兩個變量為自變量,建立回歸模型,同時為了便于比較自變量對因變量的貢獻大小,這里首先對數據進行了中心化處理,處理后的變量分別稱為zx5,zx7,zy,結果見表1(前4列)。

表1 M-估計結果

從表1中的t值能看出,系數是顯著的,不過這個推斷是基于大樣本去推斷的,很顯然,穩妥起見,需要擴大樣本量,采用bootstrap方法,結果見表1最后三列,這里的偏差和標準誤差分別是

和表1中M-估計方法相比,bootstrap M-估計方法的系數幾乎沒什么改變,偏差很小,標準誤略有增大,這個可能和離群值有關系,后面再討論離群值的問題。這些bootstrap方法的估計系數分布見圖3。

圖3 bootstrap法的M-估計系數分布

從圖3能看出,zx5系數接近于正態分布,尾部稍有點厚,zx7系數明顯右偏,和正態分布相去甚遠。在預測中這兩個系數的關系可用圖4來顯示。

圖4 兩個系數的聯合分布散點圖

從圖4能看出,這兩個系數近似于不相關,說明在對zy預測時它們相對獨立,這和開始我們采用降維的方法選取的變量結論是一致的。有了它們各自的分布和聯合分布,我們就可以計算這些系數的M-估計區間,結果見表2。

表2 bootstrap法的M-估計區間

從表2能看出,系數zx5的三種區間結果相差不是很大,和該系數的分布近似于正態分布有關系,但系數zx7的則相差較大,圖3中顯示了其分布左尾較長,這個結果和圖3中呈現的分布相一致,比較來看,顯然BCa區間的結果更加可信。

這樣,建立起的基于bootstrap法和M-估計方法的回歸模型為

比較起見,我們將基于普通最小二乘方法的回歸結果列出,見式(5)。

比較式(4)和(5),M-估計方法調低了農用化肥施用折純量對糧食產量的影響,同時增加了糧食作物播種面積的影響。式(4)中系數的t值更大,標準誤更小,說明M-估計更可靠。從回歸結果看,在固定一個自變量的情形下,農用化肥施用折純量和糧食作物播種面積每增加一個單位,糧食產量將分別增加1.0645和0.3612個單位,說明化肥在糧食種植中的作用較大,在我國經濟快速發展的形勢下,播種面積增加較為困難,重視化肥的使用將有助于糧食產量的增加。模型擬合的殘差見圖5。

圖5 回歸擬合的殘差

從圖5能看出,1980年、1981年和2003年的數據明顯偏離樣本均值,說明原始數據中可能存在異常值,為了直觀顯示每個觀測點對系數的影響,采用Jackknife-after-Bootstrap方法,結果見（圖6）。

圖6 Jackknife-after-Bootstrap方法散點圖

從圖6能看出,觀測點1(1980年),2(1981年)增大了系數zx5,減小了系數zx7,而點33(2012年)則使兩個系數均增大。對照原始數據去看,因為我國糧食產量在近30年來增加較快,所以從數據上就使得1980年和1980年的糧食產量數據顯得異常,而第24個數據(2003年)是近些年來糧食產量唯一比上一年下降的年份,所以也顯得異常,再有就是第33個數據(2012年)也顯示異常,說明去年的糧食產量數據較大,或者說距離均值較遠,也說明我國糧食產量持續增高。

3 結論

(1)在對糧食產量有眾多影響的因素中,經過降維分析,最終我們選取出了兩個主要的因素,即x5：農用化肥施用折純量(萬噸)和x7：糧食作物播種面積(千公頃),同時為了能夠比較這兩個因素對糧食產量的影響大小,我們對數據進行了中心化處理。

(2)因為樣本觀測值只有33個,實際中又不可能增加,為了得到待估計系數的優良統計性質,我們采用了近些年來較為流行的bootstrap方法來增加樣本量,同時采用較為穩健的M-估計方法,得到了回歸結果。

(3)從回歸結果看,在固定一個自變量的情形下,農用化肥施用折純量和糧食作物播種面積每增加一個單位,糧食產量將分別增加1.0645和0.3612個單位,說明化肥在糧食種植中的作用較大,在我國經濟快速發展的形勢下,播種面積增加較為困難,重視化肥的使用將有助于糧食產量的增加。

(4)最后,基于Jackknife-after-Bootstrap方法,并結合殘差對糧食產量進行了異常值分析,結果顯示觀測點1(1980年),2(1981年)增大了系數 zx5,減小了系數 zx7,而點33(2012年)則使兩個系數均增大。對照原始數據去看,因為我國糧食產量在近30年來增加較快,所以從數據上就使得1980年和1980年的糧食產量數據顯得異常,而第24個數據(2003年)是近些年來糧食產量唯一比上一年下降的年份,所以也顯得異常,再有就是第33個數據(2012年)也顯示異常,說明去年的糧食產量數據較大,或者說距離均值較遠,也說明我國糧食產量持續增高。

[1]Peter H,Joel H.A Simple Bootstrap Method for Constructing Nonparametric Confidence Bands for Functions[J].The Annals of Statistics,2013,41(4).

[2]Michael P.F,Erica H B,Michael A.P.Pointwise Confidence Intervals for A Survival Distribution With Small Samples or Heavy Censoring Biostat[J].Biostatistics,2013,14(4).

[3]Hoai T T,France M,Nicholas H.G.H.A Comparison of Bootstrap Approaches for Estimating Uncertainty of Parameters in Linear Mixed-Effects Models[J].Pharmaceutical Statistics,2013,12(3).

[4]劉薇，常振海.基于自助法的統計泛函估計比較研究[J].統計與決策,2013,(2).

[5]Huber P J.Robust Regression:Asymptotic Conjectures and Monte Carlo[J].Ann.Statist,1973,1(5).