農業經濟調查缺失數據的多重插補及應用

潘傳快，韓京芳，熊 巍，祁春節

（1.武漢紡織大學 經濟學院，武漢 430200；2.華中農業大學 經濟管理學院，武漢 430070）

1 問題的提出

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

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

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

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

圖1 多重插補思想

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

本文結合農業經濟調查缺失數據特點，重點研究多元正態模型下的聯合分布多重插補方法。先提出該方法的假設和理論，然后模擬產生農業經濟調查缺失數據進行應用分析。

2 假設條件、模型和方法

2.1 模型和方法的假設條件

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

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

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

是Y的p×p維協方差矩陣，為一個正定矩陣，具有聯合分布密度函數：

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

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

2.2 具體缺失模式

基于上述假設，農業經濟調查數據不同變量的缺失彼此獨立，而聯合分布法對缺失值進行插補需要在每一個具體的缺失模式中展開，需知道數據的具體缺失模式。比如一個維度為6×3農業經濟調查數據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，?)中產生插補值；3個觀測具體缺失模式為(1，0，1)，在該模式中Y1和Y3未缺失，而Y2缺失，那么就用Y3和Y1為Y2插補。

表1 缺失模式具體信息

在P步中，根據I步的插補值結果重新估計條件分布參數：

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

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

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

2.3 模型和方法

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

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

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

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

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

式（11）中:

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

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

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

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

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

其中z˙為一m(k)維向量，由標準正態分布N(0，1)的m(k)次抽取產生。

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

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

在給定協方差陣Σ下，μ的后驗分布服從多元正態分布：

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

2.4 模擬分析方法

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

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

3 模擬分析

3.1 基本步驟

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

3.2 原始數據模擬

設 μ=(7.6，4.3，5.2，6.4，0.5，0.6)為6個農業經濟調查變量的均值向量，相關系數矩陣為：

以均值向量和相關系數矩陣構建多元正態總體，并從中隨機抽取100個觀測作為樣本，繪制相關散點圖矩陣于圖2。

圖2 模擬數據的相關散點矩陣

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

各變量的均值為：

進一步計算樣本相關系數矩陣展示在表2中。可發現模擬效果較好，樣本均值基本接近總體均值，樣本相關系數矩陣與總體相關系數矩陣雖有差異，但可忽略。

表2 模擬樣本的相關矩陣

3.3 缺失數據模擬

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

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

表3 模擬缺失數據的具體缺失模式信息表

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

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

圖3 模擬缺失數據的缺失信息

3.4 聯合分布模型插補

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

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

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

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

表4 缺失數據聯合分布法多重插補的匯總估計

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

3.5 比較分析

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

表5 原始樣本的回歸參數估計

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

表6 缺失數據成對刪除后的回歸檢驗

4 結論

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

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

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