農(nóng)業(yè)經(jīng)濟調(diào)查缺失數(shù)據(jù)的多重插補及應(yīng)用

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

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

1 問題的提出

農(nóng)業(yè)在中國的國民經(jīng)濟和社會發(fā)展中發(fā)揮著基礎(chǔ)作用，所以在中國的社會科學(xué)研究中關(guān)于農(nóng)業(yè)經(jīng)濟以及農(nóng)業(yè)管理的研究占有很大的比重，而這些農(nóng)業(yè)經(jīng)濟管理研究很多都基于農(nóng)業(yè)經(jīng)濟調(diào)查以獲取數(shù)據(jù)。跟其他社會調(diào)查一樣，農(nóng)業(yè)經(jīng)濟調(diào)查數(shù)據(jù)也會產(chǎn)生一個幾乎無可避免的問題：缺失值[1]。缺失值的產(chǎn)生可能是由于農(nóng)戶無法回答調(diào)查問題或拒絕作出回答，比如當(dāng)問農(nóng)戶的收入時，農(nóng)戶一般很少精確統(tǒng)計過自己的收入所以無法作答，或者不想告訴調(diào)查者自己的收入而拒絕作答。導(dǎo)致缺失值還有一種情況是調(diào)查人員的遺漏，或由于疏忽造成的明顯錯誤記錄等；甚至，缺失值就是變量取值之一，即所觀測變量樣本空間的一個點[2]。

缺失值產(chǎn)生是無可避免的，但在大部分研究中，都將其作為無效數(shù)據(jù)刪除，即使不人為刪除，在諸如SPSS、SAS、Stata、R的很多統(tǒng)計軟件的數(shù)據(jù)處理中，也會默認刪除缺失值從而使整個數(shù)據(jù)處理過程能流暢進行[3]。因此盡管大量的農(nóng)業(yè)經(jīng)濟管理研究都基于調(diào)查，但幾乎無人對農(nóng)業(yè)經(jīng)濟調(diào)查數(shù)據(jù)的缺失值問題展開研究。成列刪除也叫完整個案刪除，即刪掉所有含缺失值的單元數(shù)據(jù)，是對缺失值刪除的一般方法，這種刪除在缺失比重非常小的時候是可取的[4]，特別是在在數(shù)據(jù)的缺失是完全隨機（MCAR）時甚至有較高的效率[5，6]。但當(dāng)缺失比重較高時成列刪除會導(dǎo)致過量信息被刪除，特別是當(dāng)數(shù)據(jù)有很多變量時單個變量哪怕是極小的缺失比重，成列刪除都會導(dǎo)致整體很大比例的數(shù)據(jù)遭刪除[7]。

當(dāng)然可以采取成對刪除來降低數(shù)據(jù)刪除[8]、用加權(quán)調(diào)整方法來修正成列刪除在數(shù)據(jù)非完全隨機缺失（MMAR）下的有偏估計[9]。但更好的方法是用插補來取代刪除，即根據(jù)數(shù)據(jù)的后驗分布，為缺失值插入合理的估計值，如均值、回歸值等，為了提高插補值的差異性，還可以為插補值加上隨機干擾項[10,11]，這就是單一插補。

但單一插補的問題往往低估了統(tǒng)計量的方差，使統(tǒng)計量的估計置信區(qū)間變窄，或者參數(shù)估計的檢驗顯著性下降。解決該問題的一個思路是采取重抽樣方法，如Jackknife法[12]、Bootstrap法[13]對單一插補后估計量的標(biāo)準誤差進行修正。對大部分農(nóng)業(yè)經(jīng)濟研究而言，并不需要完整的數(shù)據(jù)，只需要數(shù)據(jù)的模型估計和檢驗結(jié)果，所以一個更好的解決思路是采取多重插補，對每一個缺失值插入不止一個插補值，由于不同插補值之間存在差異，估計量的方差也會增加[14，15]。

圖1 多重插補思想

圖1顯示由于對同一缺失值插補了m次，這樣就產(chǎn)生m個“完整”數(shù)據(jù)，先根據(jù)每個“完整”數(shù)據(jù)單獨參數(shù)估計，然后將m個估計匯總成總估計，在匯總估計的標(biāo)準差計算中加入不同“完整”數(shù)據(jù)估計的差異信息。多重插補模型可根據(jù)含缺失值變量的數(shù)量分為一元模型和多元模型，多元模型在農(nóng)業(yè)經(jīng)濟調(diào)查的缺失數(shù)據(jù)中更為普遍。不過多元模型的多重插補并非單一插補的簡單疊加，需解決兩個問題：一是針對多變量的數(shù)據(jù)缺失如何同時產(chǎn)生插補值；二是很難產(chǎn)生穩(wěn)定的插補值，因為不同缺失變量的插補值會相互影響，而不同插補值又會影響模型參數(shù)。針對第一個問題，在多元正態(tài)線性假設(shè)下，Schafer提出了基于聯(lián)合分布的多重插補方法[16]，而Van Buuren等提出了基于全條件分布的多重插補方法[17]，都有較好的效果。針對第二個問題，不能用固定算法產(chǎn)生插補值而只能用基于馬爾科夫鏈蒙特卡羅（MCMC）的迭代算法，比如Gibbs抽樣[18]、數(shù)據(jù)擴增（DA）算法[19]。

本文結(jié)合農(nóng)業(yè)經(jīng)濟調(diào)查缺失數(shù)據(jù)特點，重點研究多元正態(tài)模型下的聯(lián)合分布多重插補方法。先提出該方法的假設(shè)和理論，然后模擬產(chǎn)生農(nóng)業(yè)經(jīng)濟調(diào)查缺失數(shù)據(jù)進行應(yīng)用分析。

2 假設(shè)條件、模型和方法

2.1 模型和方法的假設(shè)條件

本文所采用的模型和方法基于以下假設(shè)：

（1）Y=(Y1，Y1，…，Yp)來自一個 p 元正態(tài)總體：

其中 μ=(μ1，μ2，…，μp)為Y 的 p 維均值向量，而：

是Y的p×p維協(xié)方差矩陣，為一個正定矩陣，具有聯(lián)合分布密度函數(shù)：

（3）Yj(j=1，2，…，p) 的 缺 失 只 跟 Y-j(j=1，2，…，j-1，j+1，…，p)相關(guān)而跟自身不相關(guān)，即Y為隨機缺失（MAR）。

（4）數(shù)據(jù)Y=(Y1，Y1，…，Yp)為一般缺失，且其 p 個變量Yj(j=1，2，…，p)中，至少超過一個變量存在缺失值。

2.2 具體缺失模式

基于上述假設(shè)，農(nóng)業(yè)經(jīng)濟調(diào)查數(shù)據(jù)不同變量的缺失彼此獨立，而聯(lián)合分布法對缺失值進行插補需要在每一個具體的缺失模式中展開，需知道數(shù)據(jù)的具體缺失模式。比如一個維度為6×3農(nóng)業(yè)經(jīng)濟調(diào)查數(shù)據(jù)Y=(Y1，Y1，Y3)，如果缺失指示矩陣R為：

那么其具體的缺失模式如表1所示。表1說明有1個觀測無缺失為 (1，1，1)，2個觀測具體缺失模式為(1，1，0)，在該模式中Y1和Y2未缺失，而Y3缺失，那么就用Y1和 Y2為 Y3插補，即從后驗分布 P(Y3|Y1，Y2，?)中產(chǎn)生插補值；3個觀測具體缺失模式為(1，0，1)，在該模式中Y1和Y3未缺失，而Y2缺失，那么就用Y3和Y1為Y2插補。

表1 缺失模式具體信息

在P步中，根據(jù)I步的插補值結(jié)果重新估計條件分布參數(shù)：

將I步和P步重復(fù)迭代，直至插補值和參數(shù)最后收斂。

設(shè)農(nóng)業(yè)經(jīng)濟調(diào)查缺失數(shù)據(jù) Y={yij;i=1，2，…n;j=1，2，…p}，模型參數(shù) θ=(μ，Σ)的初始值 θ0=(μ0，Σ0)可由成列刪除或成對刪除后的均值向量和協(xié)方差矩陣組成。當(dāng)然參數(shù)的初始值也可由基于多元正態(tài)后驗分布最大似然估計產(chǎn)生。設(shè)迭代次數(shù)為t=1，2，…，T，T的取值取決于插補值和參數(shù)的收斂速度，收斂速度越快需要的T就越小，一般而言T取值為20，參數(shù)即可取得較好的收斂結(jié)果。

然后將數(shù)據(jù)按照具體缺失模式進行排序，根據(jù)其k個具體缺失模式分成k部分Y(k)(k=1，2，…，K)。在第t次迭代的第k部分Y(k)中，假設(shè)有ok個不含缺失值的變量記為O(k)，mk個變量含缺失值記為M(k)，則：

2.3 模型和方法

根據(jù)聯(lián)合分布法完成多重插補，需要知道數(shù)據(jù)不同的具體缺失模式，在每一模式中根據(jù)缺失變量基于未缺失變量的后驗聯(lián)合分布隨機產(chǎn)生插補值。由于不同變量插補值之間相互影響、而插補值又影響后驗分布參數(shù)，因此需要用迭代算法產(chǎn)生插補值而不能用固定算法。聯(lián)合分布法多重插補中引入的算法叫數(shù)據(jù)擴增（Data Augmentation,DA）法。

數(shù)據(jù)擴增法與期望最大法（Expectation-Maximum，EM）類似，均屬于馬爾科夫鏈蒙特卡洛法（Markov Chain Monte Carlo,MCMC）的一種貝葉斯統(tǒng)計方法。其算法主要由隨機抽取的I步和P步構(gòu)成，類似于EM算法中的E步和M步。在I步中隨機抽取產(chǎn)生各缺失變量的插補值，根據(jù)缺失變量基于未缺失變量的聯(lián)合條件分布：

將參數(shù)θ=(μ，Σ)更改為(p+1)×(p+1)階矩陣形式：

那么在Y(k)中，根據(jù)數(shù)據(jù)具體缺失情況，將參數(shù)矩陣進一步分解為：

接著清理參數(shù)矩陣中的未缺失變量部分子陣[16]，參數(shù)矩陣清理后為：

式（11）中:

為后驗回歸方程的截距項，而

為后驗回歸方程的斜率部分，而

為回歸模型的殘差項。清理后的參數(shù)矩陣（記為A）又可表示如下：

然后將 A中由 M(k)對應(yīng)的m(k)×m(k)階子方陣AM(k)進行Cholesky分解：

于是，在Y(k)可由以下公式產(chǎn)生插補值：

其中z˙為一m(k)維向量，由標(biāo)準正態(tài)分布N(0，1)的m(k)次抽取產(chǎn)生。

然后以此類推為每一個Y(k)(k=1，2，…，K)完成插補，數(shù)據(jù)擴增（DA）算法的I步就完成了，接下來根據(jù)I步的插補重新計算模型參數(shù)完成P步。根據(jù)I步插補后的數(shù)據(jù)，計算參數(shù)后驗分布P(θ |Yobs，Ymis)參數(shù)：均值向量 μ 和協(xié)方差陣Σ。其中協(xié)方差陣Σ的后驗分布服從Inv-Wishart分布[20]：

其中，S為k×k階對稱正定矩陣，v為自由度。其概率密度函數(shù)為：

在給定協(xié)方差陣Σ下，μ的后驗分布服從多元正態(tài)分布：

從式（18）和式（20）中各進行一次隨機抽取 θ˙=(μ˙，Σ˙)就可以完成數(shù)據(jù)擴增算法的P步。然后將I步和P步重復(fù)迭代t次，直至插補值和參數(shù)收斂，就產(chǎn)生了一次插補。再將前面的過程重復(fù)m次，就得到m個“完整”數(shù)據(jù)。

2.4 模擬分析方法

利用計算機的強大計算能力進行模擬分析，是現(xiàn)代科學(xué)研究的一種重要理論分析檢驗方法。本文的模擬分析基本思路是，先根據(jù)假設(shè)模擬產(chǎn)生缺失數(shù)據(jù)，然后利用聯(lián)合分布法進行多重插補，并對參數(shù)估計檢驗，觀測能否達到預(yù)期的處理效果，或者采用不同模型方法處理同一缺失值進行處理效果比較。

其實從理論上說，利用模擬數(shù)據(jù)分析比實際數(shù)據(jù)分析往往有更好的效果。這是因為模擬分析可以事先知道總體參數(shù)、缺失值的真值，可以讓缺失值處理模型方法的假設(shè)條件得到完全滿足，可以采取大樣本或抽取大量樣本分析，以避免單個樣本或者單次模擬分析結(jié)果的偶然性。而實際數(shù)據(jù)僅是樣本無數(shù)可能取值中的一次觀測，且無法斷定其是否符合假設(shè)條件，譬如從理論上假設(shè)數(shù)據(jù)是隨機產(chǎn)生，但在實際中調(diào)查樣本往往不是隨機的。因此利用實際缺失數(shù)據(jù)進行處理只能得到最終結(jié)果，但不能評估缺失值處理模型方法的效果，這就是本文采取模擬分析方法的原因。本文所有的數(shù)據(jù)模擬和分析都基于R語言程序。

3 模擬分析

3.1 基本步驟

模擬數(shù)據(jù)變量為6個：化肥成本（萬元）、農(nóng)藥成本（萬元）、柑橘銷售額（萬元）、柑橘產(chǎn)量（千斤）、柑橘種植面積（畝）、家庭收入（萬元）。模擬步驟為：①確定各變量的均值以及相關(guān)系數(shù)；②模擬產(chǎn)生樣本，從設(shè)定均值向量和相關(guān)系數(shù)矩陣的正態(tài)總體中隨機抽取產(chǎn)生；③在各變量中利用二項分布（失敗概率為10%）挖空觀測值，形成缺失數(shù)據(jù)；④用聯(lián)合分布方法對模擬缺失數(shù)據(jù)進行m次插補；⑤匯總估計多重插補結(jié)果，并跟未挖空之前的完整數(shù)據(jù)估計結(jié)果；⑥將模擬農(nóng)業(yè)經(jīng)濟調(diào)查缺失數(shù)據(jù)成列刪除，進行傳統(tǒng)回歸分析，并與聯(lián)合分布模型多重插補后的匯總估計結(jié)果對比。

3.2 原始數(shù)據(jù)模擬

設(shè) μ=(7.6，4.3，5.2，6.4，0.5，0.6)為6個農(nóng)業(yè)經(jīng)濟調(diào)查變量的均值向量，相關(guān)系數(shù)矩陣為：

以均值向量和相關(guān)系數(shù)矩陣構(gòu)建多元正態(tài)總體，并從中隨機抽取100個觀測作為樣本，繪制相關(guān)散點圖矩陣于圖2。

圖2 模擬數(shù)據(jù)的相關(guān)散點矩陣

從圖2相關(guān)散點圖來看，6個變量之間的相關(guān)關(guān)系明顯；主對角圖為6個變量的核密度曲線和軸須線，從中可見6個變量基本呈正態(tài)分布。圖中線性擬合曲線（細線）與平滑擬合曲線（粗線）非常接近，可見變量間的線性回歸關(guān)系也較為明顯，因為種植面積和產(chǎn)量、銷售額和產(chǎn)量、化肥成本和農(nóng)藥成本之間相關(guān)程度較高；種植面積和收入、產(chǎn)量和收入之間相關(guān)程度較低；其他變量之間相關(guān)程度居中。

各變量的均值為：

進一步計算樣本相關(guān)系數(shù)矩陣展示在表2中?？砂l(fā)現(xiàn)模擬效果較好，樣本均值基本接近總體均值，樣本相關(guān)系數(shù)矩陣與總體相關(guān)系數(shù)矩陣雖有差異，但可忽略。

表2 模擬樣本的相關(guān)矩陣

3.3 缺失數(shù)據(jù)模擬

對樣本隨機挖空產(chǎn)生缺失值，挖空的方法是利用二項分布，二項分布的成功概率為90%，這樣各變量的缺失比例大約都為10%。表3顯示了具體的缺失信息。

根據(jù)表3該缺失數(shù)據(jù)完整觀測數(shù)有57個，含缺失值的觀測數(shù)為43個，具體的缺失模式為12個。

表3 模擬缺失數(shù)據(jù)的具體缺失模式信息表

圖3進一步展示更多的缺失信息，從左圖可見各變量缺失比重并不完全為事先設(shè)定的10%，比如化肥成本，缺失比例為12%，缺失比例最低的農(nóng)藥成本只有3%。

表3的信息可進一步由圖3的右圖展示，從圖中可以更輕易地看出不同變量的缺失是彼此獨立而且完全隨機。

圖3 模擬缺失數(shù)據(jù)的缺失信息

3.4 聯(lián)合分布模型插補

接著模擬產(chǎn)生缺失數(shù)據(jù)，然后對其進行多重插補，利用聯(lián)合分布法，插補次數(shù)設(shè)定為4次。模擬數(shù)據(jù)的初始參數(shù)θ0=(μ0，Σ0)設(shè)定為成列刪除后的均值向量和協(xié)方差陣，計算結(jié)果如下：

每一次插補中都設(shè)定迭代次數(shù)為30次，在數(shù)據(jù)擴增的I步計算中，依照表3的具體缺失信息，根據(jù)缺失變量基于未缺失變量的條件聯(lián)合分布對缺失值進行插補；在P步計算中，根據(jù)I步插補值重修估算參數(shù)。其他3次插補同理產(chǎn)生。

模型估計參數(shù)為：其中α為回歸模型的截距項，β為回歸模型的斜率項，ε為回歸模型的殘項。根據(jù)前面聯(lián)合分布法多重插補后的4個“完整”數(shù)據(jù)，將參數(shù)?匯總估計在表4中。該回歸模型確定部分可以表示如下：

銷售收入=1.2802+0.5533柑橘銷售額+0.6941柑橘產(chǎn)量-0.984柑橘種植面積+5.6410農(nóng)藥成本+6.3488化肥成本

表4 缺失數(shù)據(jù)聯(lián)合分布法多重插補的匯總估計

通過表4可以發(fā)現(xiàn)，根據(jù)聯(lián)合分布法的多重插補結(jié)果計算的估計檢驗量有很好的效果，除了銷售額的回歸系數(shù)顯著性較低外，其他分量的回歸系數(shù)都非常顯著，雙側(cè)檢驗（原假設(shè)H0為：β=0）的P值都很小。

3.5 比較分析

接著比較分析聯(lián)合分布法多重插補后匯總估計的效果。先將多重插補結(jié)果與模擬缺失之前的完整數(shù)據(jù)分析結(jié)果進行比較，完整數(shù)據(jù)的分析結(jié)果顯示在表5中。通過對比發(fā)現(xiàn)，兩者的參數(shù)估計結(jié)果很接近，只是由于數(shù)據(jù)的缺失，根據(jù)多重插補結(jié)果估計的參數(shù)顯著性更低，P值更大。

表5 原始樣本的回歸參數(shù)估計

最后將模擬缺失數(shù)據(jù)成列刪除并進行回歸估計，其結(jié)果見表6。將其結(jié)果跟聯(lián)合分布法多重插補后的結(jié)果比較可發(fā)現(xiàn)，成列刪除的點估計結(jié)果跟多重插補的點估計結(jié)果以及未挖空的完整數(shù)據(jù)的點估計結(jié)果都很接近，這是因為數(shù)據(jù)的缺失是獨立的且完全隨機（MCAR）。但從參數(shù)估計的檢驗來看，較之多重插補，成列刪除的估計參數(shù)的顯著性普遍較低，P值更大，柑橘產(chǎn)量對應(yīng)的系數(shù)即使在0.1的顯著水平下也未拒絕原假設(shè)，這是由于成列刪除額外刪除了更多的有用信息造成的。

表6 缺失數(shù)據(jù)成對刪除后的回歸檢驗

4 結(jié)論

正如Allison所言，任何人在進行數(shù)據(jù)分析的時候，早晚（通常是早）都會遇到缺失數(shù)據(jù)問題[1]。農(nóng)業(yè)經(jīng)濟調(diào)查數(shù)據(jù)的缺失值也是難以避免的，遺憾的是大部分場合下都被研究人員和數(shù)據(jù)分析軟件作簡單的刪除處理。插補是一個更好的缺失值處理方法，比較分析發(fā)現(xiàn)成列刪除會增加額外信息的丟失，從而引致參數(shù)估計的顯著性下降。由于單一插補會低估估計量的標(biāo)準誤差，而多重插補可以用同一缺失值不同插補值之間的差異來彌補標(biāo)準誤差的低估。

在數(shù)據(jù)為一般缺失模式下，利用多元正態(tài)模型下的聯(lián)合分布法對農(nóng)業(yè)經(jīng)濟調(diào)查缺失數(shù)據(jù)進行多重插補，擁有很好的估計檢驗效果。模擬分析顯示，聯(lián)合分布模型多重插補后的估計量跟完整數(shù)據(jù)的估計量非常接近，只是由于數(shù)據(jù)缺失造成的誤差損失使檢驗顯著性下降；跟成列刪除后數(shù)據(jù)的估計檢驗結(jié)果相比，其準確性更高、檢驗顯著性更強。

不過需要指出的是插補方法雖然沒有像刪除方法那樣丟失信息，但是加入了額外信息，那么其加入額外信息的準確與否就非常關(guān)鍵。對缺失數(shù)據(jù)進行插補一般有良好的效果但是也存在風(fēng)險，它會讓人們以為獲得了完整數(shù)據(jù)，但其插補值的正確性很難界定，只能寄希望于插補值和真實值之間的差距不要太大[21]。所以在實際農(nóng)業(yè)經(jīng)濟調(diào)查缺失數(shù)據(jù)的分析時，對數(shù)據(jù)是否符合模型假設(shè)進行檢驗是非常必要的。此外多重插補還有一個問題值得注意，那就是其只產(chǎn)生估計檢驗結(jié)果，不產(chǎn)生插補后數(shù)據(jù)，如果農(nóng)業(yè)研究人員要展示數(shù)據(jù)就不能用多重插補方法。