农业经济调查缺失数据的多重插补及应用

潘传快，韩京芳，熊 巍，祁春节

（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]。所以在实际农业经济调查缺失数据的分析时，对数据是否符合模型假设进行检验是非常必要的。此外多重插补还有一个问题值得注意，那就是其只产生估计检验结果，不产生插补后数据，如果农业研究人员要展示数据就不能用多重插补方法。