農(nóng)業(yè)經(jīng)濟(jì)調(diào)查缺失數(shù)據(jù)的貝葉斯和Bootstrap多重插補的比較

熊 巍，潘傳快，祁春節(jié)

（1.華中農(nóng)業(yè)大學(xué) 經(jīng)濟(jì)管理學(xué)院，武漢430070；2武漢紡織大學(xué) 經(jīng)濟(jì)學(xué)院，武漢430200）

1 問題的提出

對農(nóng)業(yè)、農(nóng)民、農(nóng)村問題的研究很多都需要開展農(nóng)業(yè)經(jīng)濟(jì)調(diào)查以獲取數(shù)據(jù)。理想的情況是這些數(shù)據(jù)都是完整無缺失的，但其實是很難做到的[1]。通常情況下，農(nóng)業(yè)經(jīng)濟(jì)研究人員遇到缺失值都會采取成列刪除的方法，即作為無效數(shù)據(jù)刪除所有含缺失值的個案，留下完整數(shù)據(jù)作為分析[2]。即使農(nóng)業(yè)研究人員忽略了這個問題，但在大部分分析軟件中，為了數(shù)據(jù)分析得以流暢進(jìn)行，都會自動成列刪除缺失數(shù)據(jù)[3]。如果數(shù)據(jù)的缺失比重很低，并且是完全隨機缺失（MCAR），成列刪除相當(dāng)于在損失很小的情況下得到一個完整數(shù)據(jù)，而且這個完整數(shù)據(jù)還是原樣本的簡單隨機子樣，就意味著用該數(shù)據(jù)估計是無偏的[4,5]。但遺憾的是農(nóng)業(yè)經(jīng)濟(jì)調(diào)查的缺失比重一般較高，而且很多時候即使缺失比重不高，但變量太多，采取成列刪除也會引起很大的數(shù)據(jù)損失[4]。此外，數(shù)據(jù)缺失往往不是完全隨機的，比如產(chǎn)量越高的種植戶越不傾向于透露自己的產(chǎn)量，那么成列刪除后的數(shù)據(jù)就是有偏的。所以在處理農(nóng)業(yè)經(jīng)濟(jì)調(diào)查缺失數(shù)據(jù)時，刪除一般都不是最好的處理方法[6]。更值得推薦的方法是不刪除任何數(shù)據(jù)，而結(jié)合數(shù)據(jù)的經(jīng)驗分布，對其缺失部分進(jìn)行插補[7]，或?qū)笔е档臄?shù)據(jù)進(jìn)行最大似然估計[8]。插補特別是多重插補由于能更精準(zhǔn)地估計模型參數(shù)而更多地被采用[9]。因為多重插補可以利用同一缺失值不同填補值之間的差異來抵消估計和檢驗誤差的減少，使估計檢驗更準(zhǔn)確[10]。特別是當(dāng)農(nóng)業(yè)經(jīng)濟(jì)研究人員不需要數(shù)據(jù)本身，只需要根據(jù)數(shù)據(jù)建立農(nóng)業(yè)經(jīng)濟(jì)計量模型以完成研究。一個基本的多重插補思想是將根據(jù)完整未缺失數(shù)據(jù)的均值或回歸值加上均值為0的“白噪聲”干擾，但其插補值的離散程度仍然偏低，因為利用該方法的模型參數(shù)是固定的[11]。可以利用貝葉斯法和Bootstrap法讓回歸參數(shù)隨機化，貝葉斯法的基本思想是讓插補模型參數(shù)來自其后驗分布的一個隨機抽取[12]，而Bootstrap法的基本思想是讓線形模型的參數(shù)來自原完整樣本的多個Bootstrap樣本[13]。Bootstrap法跟貝葉斯法相比，可以不需要兩個步驟：Choleski分解和從χ2分布中抽取以完成殘差估計，因此效率上會更高一點。

本文主要對一元線性模型下貝葉斯多重插補和Bootstrap多重插補進(jìn)行比較。作為實證分析案例，本文選取在中國柑橘主產(chǎn)區(qū)對種植戶所做的調(diào)查數(shù)據(jù)，分別利用兩種方法對缺失值進(jìn)行多重插補，再匯總計算農(nóng)業(yè)經(jīng)濟(jì)計量模型所需要參數(shù)的估計和檢驗結(jié)果并加以比較。本文所有的數(shù)據(jù)分析均使用R程序語言軟件。

2 模型和方法

2.1 假設(shè)

為了方便兩種方法的比較，本文的插補模型和方法基于以下假設(shè)：

（1）目標(biāo)缺失變量為單一變量Y且呈正態(tài)分布，而用以輔助插補的完整變量X可以是多元的，但仍成正態(tài)分布。如果數(shù)據(jù)中存在多個變量缺失，那就分解成多個插補模型，即每一個缺失變量構(gòu)建一個插補模型。

（2）目標(biāo)變量Y為隨機缺失（MAR），即Y的缺失跟本身未觀測值無關(guān)，但跟輔助變量X有關(guān)。這意味著只要借助輔助變量而不需要借助缺失變量來完成對缺失值的插補。

（3）Y跟X之間顯著存在正態(tài)線性相關(guān)關(guān)系，即：

式（1）中，β為回歸系數(shù)，是一個q維向量，而σ2為殘差標(biāo)準(zhǔn)誤差，也即插補模型隨機干擾項的平方。

2.2 貝葉斯多重插補

農(nóng)業(yè)經(jīng)濟(jì)調(diào)查數(shù)據(jù)由缺失變量Y和未缺失變量X構(gòu)成，如果未缺失變量不唯一，則X為矩陣。Y有n1個未缺失的觀測值記為Yo，其對應(yīng)的X部分記為Xo；其n0缺失值記為Ym，其對應(yīng)的X部分記為Xm。根據(jù)式（1），利用數(shù)據(jù)的未缺失部分Yo和Xo，可以獲得模型的最小平方估計和：

其中為模型回歸參數(shù)向量，而殘差方差估計值為：

根據(jù)貝葉斯思想，β和的后驗分布可由構(gòu)建，其中的后驗分布為χ2分布：

而基于的β的條件后驗分布為正態(tài)分布：

缺失變量的每個缺失值Ym可以采取以下模型產(chǎn)生任意多個插補值：

式（6）中，為根據(jù)式（4）確定的后驗分布的任意個隨機抽取，而來自式（5）給定的后驗分布的任意個隨機抽取，Z1為n0個標(biāo)準(zhǔn)正態(tài)分布值。即：

式（7）中，g是來自χ2(n1-q)分布的一個隨機值；式（8）中Z2為q個標(biāo)準(zhǔn)正態(tài)分布隨機值為的Cholesky分解。

2.3 Bootstrap多重插補

Bootstrap多重插補的一個核心思想是，讓插補模型的回歸參數(shù)來自Bootstrap樣本的估計，由于Bootstrap樣本不唯一，所以根據(jù)每一個Bootstrap樣本估計的參數(shù)也不唯一。設(shè)(Yb，Xb)為來自觀測部分(Yo，Xo)一個的Bootstrap樣本,其容量為n1。根據(jù)該Bootstrap樣本，模型參數(shù)的最小平方估計為：

而殘差方差的估計為：

根據(jù)式（9）和式（10），可以利用以下模型產(chǎn)生缺失部分Ym的任意插補值插補：

式（11）中，Z1為n0個標(biāo)準(zhǔn)正態(tài)分布值。

2.4 匯總估計方法

（1）插補次數(shù)的選擇

多重插補的思想是利用對同一缺失值的多次插補來彌補估計標(biāo)準(zhǔn)誤差的損失，那么需要多少次的插補才能彌補這種損失呢？從理論上來說，插補次數(shù)越多估計越準(zhǔn)確，如果插補次數(shù)趨于無窮大，估計標(biāo)準(zhǔn)誤差的估計幾乎是完全準(zhǔn)確的。但是在實際上，這意味著巨大的計算量，縱使對計算機而言，如果插補次數(shù)設(shè)定很大計算量也是驚人的。很多研究顯示，當(dāng)插補次數(shù)超過10次時，對估計標(biāo)準(zhǔn)誤差的改進(jìn)已經(jīng)非常有限，因此在實際操作中，只要將插補次數(shù)設(shè)定為5～10次即可。

（2）多重插補的點估計

假設(shè)在一次農(nóng)業(yè)經(jīng)濟(jì)調(diào)查中建立農(nóng)業(yè)經(jīng)濟(jì)計量模型所需要的總體參數(shù)為θ，θ可以是單一的標(biāo)量或者是向量，取決于感興趣參數(shù)的個數(shù)，比如想知道農(nóng)戶的平均家庭收入θ=μ。多重插補后θ的匯總估計為：

式（12）中，m為插補次數(shù)對，而參數(shù)的內(nèi)部方差：

其中，Uj為插補后數(shù)據(jù)的的協(xié)方差矩陣。而同一缺失值不同插補值之間的方差為：

式（14）中，m-1為自由度。因此估計參數(shù)θ的總方差為：

式（15）為有限次多重插補的參數(shù)方差計算公式，并不精確，可以用下式進(jìn)行調(diào)整：

（3）多重插補的區(qū)間估計和檢驗統(tǒng)計量

統(tǒng)計量ˉ符合一個自由度為νnc的Student-t分布：

式（17）中，自由度νnc的計算方法為：

根據(jù)式（17），θ的一個置信水平為(1-α)的置信區(qū)間為：

式（19）中，tα2(νnc)為Student-t分布在自由度νnc下的第100(1-α/2)分位數(shù)。

在檢驗時，θ的檢驗t值如下：

其對應(yīng)的P值為：

如果利用方差分析檢驗?zāi)Ｐ蛥?shù)的顯著性，則F檢驗統(tǒng)計量構(gòu)建如下：

其相應(yīng)的P值為：

3 實證分析

本文的實證數(shù)據(jù)為國家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系（柑橘）團(tuán)隊在中國柑橘主產(chǎn)區(qū)開展的一項調(diào)查。調(diào)查內(nèi)容為柑橘種植戶的柑橘生產(chǎn)銷售成本收益，調(diào)查對象為主產(chǎn)區(qū)的柑橘種植戶，調(diào)查方式為入戶訪問調(diào)查，調(diào)查問卷以開放性問題為主，調(diào)查樣本為調(diào)查工作人員篩選后的典型種植戶。

3.1 缺失信息

本次調(diào)查搜集的數(shù)據(jù)較為龐大，由于篇幅所限和分析需要，本文選取其中的部分典型變量來進(jìn)行實證分析，選取的這些變量中既有完整數(shù)據(jù)也有缺失數(shù)據(jù)，而且變量之間普遍存在相關(guān)關(guān)系，這些變量的部分?jǐn)?shù)據(jù)顯示在表1中。

表1 選區(qū)的變量的部分單元數(shù)據(jù)

為了在R中運行，在表1中將缺失值表示為“NA”即為缺失值（Not Available）。表2進(jìn)一步描述了這些變量的缺失信息。

表2 每個變量的缺失信息

從表2可以發(fā)現(xiàn)，除了總產(chǎn)量以外，所有的變量都有缺失，但是有些變量缺失比重很高，比如銷售費用和人工成本，而有些變量缺失比重非常低，比如農(nóng)藥費和化肥費用，其他變量的缺失比重居中。

3.2 輔助變量選擇

本文選取總產(chǎn)值為目標(biāo)缺失變量，按照模型假設(shè)輔助變量必須是完整的，但該數(shù)據(jù)只有總產(chǎn)量一個變量達(dá)到要求，為了提高插補的準(zhǔn)確性，再將缺失比例較低的化肥費用和農(nóng)藥費用加進(jìn)來。這兩個變量為數(shù)不多的幾個缺失值用各自的簡單隨機抽取進(jìn)行單一插補，插補效果如圖1所示，圖中“×”表示插補值。

圖1 農(nóng)藥費和化肥費的簡單隨機插補

為進(jìn)一步考察目標(biāo)缺失變量和輔助變量之間的相關(guān)關(guān)系，將這4個變量的完整部分?jǐn)?shù)據(jù)繪制成相關(guān)散點圖,如圖2所示。從圖中可以發(fā)現(xiàn)，4個變量均成正態(tài)分布，基本符合正態(tài)假設(shè)，只是由于數(shù)據(jù)來自典型種植戶調(diào)查，而這些典型種植戶有很多是種植大戶，因此分布明顯右偏。雖然極端值存在且破壞了相關(guān)關(guān)系，但從圖中仍可以看出，各個變量之間存在線性相關(guān)關(guān)系，也基本符合假設(shè)需要。

圖2 變量之間的相關(guān)散點圖矩陣

3.3 多重插補

現(xiàn)利用R編程，采用貝葉斯法和Bootstrap法對缺失數(shù)據(jù)進(jìn)行多重插補，然后估計經(jīng)濟(jì)計量模型的參數(shù)，事先假定該經(jīng)濟(jì)計量模型為總產(chǎn)值關(guān)于化肥費用、農(nóng)藥費用以及總產(chǎn)量的線性回歸模型。表3為采取貝葉斯多重插補后估計的模型參數(shù)以及相應(yīng)檢驗。

表3 貝葉斯多重插補后數(shù)據(jù)的估計檢驗結(jié)果

根據(jù)貝葉斯多重插補的分析結(jié)果，可以得到如下農(nóng)業(yè)經(jīng)濟(jì)計量模型：

總產(chǎn)值=2074.37+1.92283×總產(chǎn)量-4.5269×化肥費+2.95065×農(nóng)藥費

表3不僅給出了該模型參數(shù)的點估計和區(qū)間估計，還給出了標(biāo)準(zhǔn)誤差、檢驗t值和P值。表中的lambda值和fmi值表示因缺失比重的增加而增加的誤差。作為比較，再將該缺失數(shù)據(jù)利用Bootstrap法進(jìn)行多重插補，其模型的估計檢驗結(jié)果展示在表4中。

表4 Bootstrap多重插補后數(shù)據(jù)的估計檢驗結(jié)果

根據(jù)Bootstrap多重插補的分析結(jié)果，可以得到如下農(nóng)業(yè)經(jīng)濟(jì)計量模型：

總產(chǎn)值=2066.00+1.93053×總產(chǎn)量-4.3052×化肥費+2.85846×農(nóng)藥費

3.4 比較分析

通過表3和表4的比較分析，采取貝葉斯法和Bootstrap法多重插補后的參數(shù)點估計區(qū)別非常小，幾乎可以忽略，估計量的標(biāo)準(zhǔn)差也非常接近，檢驗的t統(tǒng)計量也差別很小，從這點上來看兩者似乎沒有太大的差別。但可以發(fā)現(xiàn)兩者對自由度的估計差別還是很大的，采用貝葉斯法估計的參數(shù)自由度似乎比Bootstrap法估計的自由度大（農(nóng)藥費用對應(yīng)的自由度除外），說明貝葉斯法的估計更可靠些（有更大的樣本支持），從檢驗P值上也可以看出，在相同的t值水平下，采用貝葉斯法由于更大的自由度，P值更偏小些，從區(qū)間估計上看也具有更窄的置信區(qū)間，說明同等置信水平下估計的準(zhǔn)確性更高。

兩者的fmi值以及λ值區(qū)別不大，可以看出由于目標(biāo)變量缺失的比例不高，因此由于數(shù)據(jù)缺失而引致的回歸模型額外方差比重并不算高。

4 結(jié)論

在缺失模式為單一缺失、缺失機制為隨機缺失（MAR）的條件下，假設(shè)缺失變量和輔助變量的后驗分布為正態(tài)分布，彼此關(guān)系為線性關(guān)系，那么就可以構(gòu)建一元正態(tài)線性多重插補模型。在該模型下，通過兩個方面來增加插補值的離散型，一是為插補值加入隨機干擾項，二是讓模型參數(shù)隨機產(chǎn)生。模型參數(shù)的隨機產(chǎn)生方法有兩種：貝葉斯法和Bootstrap方法，兩者產(chǎn)生隨機參數(shù)的思想截然不同。

根據(jù)對中國柑橘主產(chǎn)區(qū)的種植戶的入戶調(diào)查數(shù)據(jù)的實證分析發(fā)現(xiàn)，雖然貝葉斯法和Bootstrap法的思路不同，但最終根據(jù)兩者構(gòu)建的農(nóng)業(yè)經(jīng)濟(jì)計量模型的點估計結(jié)果卻差異很小。不過值得注意的是，雖然點估計結(jié)果差別很小，但是兩者對參數(shù)自由度的估計卻差別較大，采用貝葉斯法估計似乎有更大的自由度，也就意味者有更準(zhǔn)確的區(qū)間估計和更顯著的檢驗效果。但是從效率上說，Bootstrap方法會更有優(yōu)勢一點，因為其比貝葉斯法省了兩個步驟，使計算機編程的難度降低。此外，跟成列刪除相比，插補雖然沒有損失信息，但加入了額外信息，這需要充分驗證缺失值插補模型是否符合假設(shè)，否則其插補效果甚至不如刪除。