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

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

1 問題的提出

農(nóng)業(yè)在中國的國民經(jīng)濟和社會發(fā)展中發(fā)揮著基礎作用，所以在中國的社會科學研究中關于農(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)查問題或拒絕作出回答，比如當問農(nóng)戶的收入時，農(nóng)戶一般很少精確統(tǒng)計過自己的收入所以無法作答，或者不想告訴調(dià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]。但當缺失比重較高時成列刪除會導致過量信息被刪除，特別是當數(shù)據(jù)有很多變量時單個變量哪怕是極小的缺失比重，成列刪除都會導致整體很大比例的數(shù)據(jù)遭刪除[7]。

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

但單一插補的問題往往低估了統(tǒng)計量的方差，使統(tǒng)計量的估計置信區(qū)間變窄，或者參數(shù)估計的檢驗顯著性下降。解決該問題的一個思路是采取重抽樣方法，如Jackknife法[12]、Bootstrap法[13]對單一插補后估計量的標準誤差進行修正。對大部分農(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個估計匯總成總估計，在匯總估計的標準差計算中加入不同“完整”數(shù)據(jù)估計的差異信息。多重插補模型可根據(jù)含缺失值變量的數(shù)量分為一元模型和多元模型，多元模型在農(nóng)業(yè)經(jīng)濟調(diào)查的缺失數(shù)據(jù)中更為普遍。不過多元模型的多重插補并非單一插補的簡單疊加，需解決兩個問題：一是針對多變量的數(shù)據(jù)缺失如何同時產(chǎn)生插補值；二是很難產(chǎn)生穩(wěn)定的插補值，因為不同缺失變量的插補值會相互影響，而不同插補值又會影響模型參數(shù)。針對第一個問題，在多元正態(tài)線性假設下，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)合分布多重插補方法。先提出該方法的假設和理論，然后模擬產(chǎn)生農(nóng)業(yè)經(jīng)濟調(diào)查缺失數(shù)據(jù)進行應用分析。

2 假設條件、模型和方法

2.1 模型和方法的假設條件

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

（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)相關而跟自身不相關，即Y為隨機缺失（MAR）。

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

2.2 具體缺失模式

基于上述假設，農(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步重復迭代，直至插補值和參數(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é)方差矩陣組成。當然參數(shù)的初始值也可由基于多元正態(tài)后驗分布最大似然估計產(chǎn)生。設迭代次數(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)中，假設有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步構成，類似于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)對應的m(k)×m(k)階子方陣AM(k)進行Cholesky分解：

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

其中z˙為一m(k)維向量，由標準正態(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步重復迭代t次，直至插補值和參數(shù)收斂，就產(chǎn)生了一次插補。再將前面的過程重復m次，就得到m個“完整”數(shù)據(jù)。

2.4 模擬分析方法

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

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

3 模擬分析

3.1 基本步驟

模擬數(shù)據(jù)變量為6個：化肥成本（萬元）、農(nóng)藥成本（萬元）、柑橘銷售額（萬元）、柑橘產(chǎn)量（千斤）、柑橘種植面積（畝）、家庭收入（萬元）。模擬步驟為：①確定各變量的均值以及相關系數(shù)；②模擬產(chǎ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ù)模擬

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

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

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

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

各變量的均值為：

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

表2 模擬樣本的相關矩陣

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

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

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

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

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

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

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

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

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

每一次插補中都設定迭代次數(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è)檢驗（原假設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)量對應的系數(shù)即使在0.1的顯著水平下也未拒絕原假設，這是由于成列刪除額外刪除了更多的有用信息造成的。

表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ù)估計的顯著性下降。由于單一插補會低估估計量的標準誤差，而多重插補可以用同一缺失值不同插補值之間的差異來彌補標準誤差的低估。

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

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