不同缺失机制并存时偏倚校正的模拟研究*

赵俊康 王 彤 荣惠英 孟繁龙

弱势人群的医疗救助问题一直以来备受世界各国政府关注[1]。这部分弱势人群的特点是收入偏低，极易陷入因贫致病和因病致贫的恶性循环中。根据第四次国家卫生服务调查结果[2]，我国约有38%的居民生病不去看医生，经医生诊断该住院治疗而未住院的达21%，其中70.3%的人未住院的主要原因仍然是“经济困难”。这些潜在患者不选择就医使得从医院收集数据仍然很难估计出这部分非从业人群的全部医疗费用需求。低收入非从业人群更可能由于贫困等原因得到高于普通人群的致病机会，因此若用一个总的平均水平来估计弱势人群医疗费用的实际需求将明显低估这种需求。

这种由于个体自我行为(因经济困难自主选择不就医)所导致的样本选择偏倚，单靠好的抽样设计是无法消除的。需要注意的是患病但自我选择未就医者应答表现出的0消费与真正未患病而不就医者的消费真值是不同的，即自我选择未就医者的医疗消费真值未知，应视为缺失数据[3]。将这类真值未知的0消费数据删除或者直接取因变量为0来应用多元线性回归等常规的统计学分析方法就忽视了这种无应答偏倚；同时，像这类较大规模的社会学或流行病学调查中无应答偏倚也是常态而不是偶然[4]，故而针对不同缺失机制下的无应答偏倚探讨其校正方法成为国内外学者长期以来关注的问题。

Rubin等人于1976年提出的缺失机制主要包括完全随机缺失MCAR(missing completely at random)、随机缺失MAR (missing at random)和非随机缺失NMAR(not missing at random)三类[5]。在MCAR假定下，对完全观测个体使用的任何分析方法仍然有效；在MAR假定下,主流观点是采用多重填补MI(multiple imputation)对随机缺失进行填补继而得出无偏估计[6]；而对于由于自主选择不就医而导致的NMAR，本研究选用适合于该类型数据的受限因变量(limited dependent variable)统计模型来进行校正[7-9]。

随机缺失机制下的多重填补方法

1.多重填补的具体步骤

多重填补(MI)主要由三个独立的步骤组成：填补阶段、分析阶段和合并阶段。MI其实是包含了一组方法的一个广义的术语，在其框架内的所有的方法中都含有这三步过程。图1描述了整个过程。

填补阶段 分析阶段 合并阶段

图1 MI的三个步骤

(1)填补阶段为每个缺失值抽取m个估计值进行填补，从而构成m个完整数据集，这m个数据集中只有观测数据是相同的，填补值一般不等。(2)分析阶段：分析步的主要分析对象就是填补好的数据集，这一步将应用数据原本完整时所用到的相同的方法来分析。唯一的区别在于要对每个完整数据集分别使用该方法处理，因此将分析m次。对于本研究中包含自主选择性缺失的医疗费用数据，这一步可使用填补后的多个数据集与选择性偏倚导致的因变量为虚假0(视为缺失)的数据合并进行样本选择模型分析以解决选择性偏倚导致的那部分缺失，继而得出m个样本选择模型拟合结果。(3)合并阶段：综合这m个拟合结果，根据Rubin(1987)提出的针对参数估计值与标准误的合并准则[10]，最终得到对目标变量的统计推断。

2.填补模型

(1)预测均数匹配法(PMM)

预测均数匹配法(PMM)是处理单调缺失模式中定量变量缺失的多重填补方法之一，PMM法的具体填补步骤如下：

令Yj为有缺失值的定量变量，用Yj及其协变量X1,X2,…,Xk均被观测到的个体观察值建立回归模型：

Yj=β0+β1X1+β2X2+…+βkXk

(1)

每次填补产生填补值的步骤如下：

(2)

其中nj是变量Yj未缺失的观测个体数。

然后，从标准正态分布N(0,1)中抽取k+1个独立的变量，组成一个有k+1个元素的向量Z，得到新的回归系数：

(3)

②对于每个缺失值，其预测值为

(4)

④最后从这k0个观测值中随机抽取一个值填补缺失值。

(2) 倾向性得分法

倾向性得分(PS)是指对给定的观察到的协变量条件下，每个观察值被分配到某特定处理组的条件概率。PS法的具体填补步骤如下：

①产生一个指示变量Rj，当Rj为0时，表示变量Yj中有缺失值的个体；当Rj为1时，则表示变量Yj中被观测到的个体。

②拟合logistic回归模型

logit(pj)=β0+β1X1+β2X2+…+βkXk

(5)

其中X1,X2,…,Xk是Yj的协变量，

pj=Pr(Rj=0|X1,X2,…,Xk)

logit(pi)=log(pi/(1-pi))

③根据模型计算变量Yj上每个个体数据缺失的倾向性得分logit(pj)，并根据该得分将所有的观测分组，一般为5组，如果观测数量较多，可分为更多的组。

⑤重复以上步骤，直到每一个缺失变量都得到填补。

(3)基于Bootstrap的EM算法

基于Bootstrap的EM算法(EMB)是通过bootstrap算法从参数的后验密度中抽取新的参数，从而代替其他方法中复杂的抽取过程。该方法不必估计参数的方差矩阵，也不用像期望最大化重要性抽样EMis(expectation maximization importance sampling)算法那样进行重要性抽样，甚至不需要像数据增广DA (data augmentation ) 算法一样推导Markov链并检查收敛性，而且还可以应用于非常大型的数据[11]。

EM算法是一种迭代算法，广泛运用于寻找参数的最大似然估计值，尤其是在缺失数据的问题中非常有用[12]。它的每一次迭代都由两步组成：E步(求期望)和M步(极大化)。其中E步是在给定已观测到的数据和当前参数下，求缺失数据的条件期望，然后用这些条件期望值去填补缺失值。M步是当缺失数据被填补之后就像没有缺失一样进行的极大似然估计。重复以上两步,直至前后两次计算结果达到规定的收敛标准。

将经典的EM算法和Bootstrap算法相结合来从后验分布中进行抽取。具体填补步骤如下：

①从含量为n的完整数据集Yobs中有放回地抽取m个大小为n的样本。

②对每一个样本运行稳定而快速的EM算法，得到一组参数的点估计值μ和∑，共m组。

③用Ycom中的观测数据分别结合每组参数估计值得出缺失数据的条件分布，并从中抽取填补值，继而得到m个经填补后的完整数据集。

(4)Markov Chain Monte Carlo方法

Markov Chain Monte Carlo(MCMC)方法在应用于缺失值领域时称为数据增广DA算法[13]，同EM算法一样，DA算法也是依次填补缺失值和推断未知参数的一种迭代方法。区别就在于DA是以随机的方式对缺失值和参数进行抽取，而EM算法只是缺失值和参数的点估计。

该算法通过填充(imputation)及后验(posterior)两步迭代来实现：

①填补步(I-step)

填补的数据是从给定观测数据、均数和协方差矩阵后的缺失数据的条件分布中随机抽取得到，从贝叶斯角度来看，该分布又称为后验预测分布。

在每次迭代过程中，从事先给定的均值向量μ和协方差阵Σ的初始估计值开始，在给定Yobs下的条件分布P(Ymis|Yobs)中抽取Ymis。

在某个观测个体具有类似的缺失模式时，令Yobs=y1，就得到均值向量及条件协方差矩阵分别为

(10)

(11)

的多元正态分布P(Ymis|Yobs=y1)，也就是Ymis的条件分布。整个I-step可以表述成：

(12)

②后验步(P-step)

由于在多重填补过程中需要产生多个完整的数据集，因此在每个I-step需要不同的均数向量和协方差矩阵，因此P-step的目的就是轮流产生参数估计值。

(13)

(14)

(15)

这样从后验分布中抽取新的参数后，接下来的I-step使用这些新的参数值来产生新的填补值；然后新的填补数据继续用于下一个P-step，继而再抽取另一组新的参数估计值。如此循环往复重复这两个步骤一定的次数，产生一个足够长的随机序列：

(16)

该随机序列是就一条马尔科夫链，并且在一定的正则条件下会收敛到一个稳定分布[14]。当该链收敛到一个稳定的分布P(Ymis,θ|Yobs)时，就可以近似独立地从该分布中为缺失值抽取填补值。

两种缺失机制并存时的两阶段校正模拟研究

1.模拟设计

首先需要构建出含结果等式和选择等式的样本选择模型，如下所示：

y0=x1+ε

(17)

d0=x2+v

(18)

d=1(d0≥c)，d=0(d0

(19)

y1=y0·d

(20)

取10000例观测值(n=10000)，根据样本选择模型的结果等式(17)和选择等式(18)模拟六个变量x1、x2、y0、d0、ε、v，x1和x2分别取自均值为0，标准差为1，相关系数为0的双变量正态分布，而ε和v取自均值为0，标准差为1，相关系数为0.75的双变量正态分布。y0和d0分别通过公式y0=e1+x1+ε和d0=1+x2+v求出。第一步首先对全部10000例观测值的因变量分别进行轻度、中度和重度截取，即以5%、30%和70%的比例向下截取产生对应于调查中生病但自主选择不就医者发生的虚假0消费因变量。可以通过以上截取比例给y0定义一个相应的界值c，当d0>c时，令y1=y0且d=1；当d0≤c时，令y1为缺失且d=0,最后对y1进行对数转换，对应于使医疗费用值近似服从正态分布。通过调整c的值就可以获得针对样本选择模型不同程度缺失率的非随机缺失数据。第二步分别在上述三个不同缺失率下，令d=1的个体以5%、30%和70%的比例随机产生缺失(对应于调查中可能发生的随机缺失)。这样就产生9种不同的组合数据。在上述不同组合下，首先对d=1的个体(仅存在随机缺失)分别应用PMM、PS、MCMC和EMB法进行多重填补，然后把填补后的数据与d=0的数据合并，应用样本选择模型的两步似然估计来获得各自的回归系数估计量来校正虚假0消费产生的选择性偏倚。最后，重复抽样100次,计算9种组合下两阶段校正方法所获得结果等式中自变量x1的回归系数和标准误。本次模拟分析中，多重填补技术的EMB算法选用了R软件，PMM、PS和MCMC法选用了SAS软件，样本选择模型分析也选用了SAS软件。

2．评价标准

在比较两阶段校正方法下样本选择模型结果等式回归系数的优劣时，选用以下三个评价标准[15-16]。

(1) 标准偏倚(Standardized bias)

当标准偏倚落在±0.4区间之外时，偏倚就会对功效、可信区间覆盖率和误差率产生明显的负面影响。因此，将±0.4作为评价标准偏倚的上下界值，即若某方法的标准偏倚绝对值超出0.4，此方法便无法接受。标准偏倚做为评价准确度的指标是方法评价指标中的首要观测指标。

(2)可信区间平均长度(length)

如果一个方法与另一个方法相比，有相同的或更高的准确度，但得出的可信区间平均却更短，那么此方法的精确度就更高。

(3)均方误差的平方根(RMSE)

3．模拟分析结果

表1 四种填补方法下样本选择模型结果等式的回归系数估计值的各项评价标准比较(一)

从表1到表3可知，四种填补方法的标准偏倚绝对值均不等，其中PS法超过了所规定的界值，故该法效果相对不理想；其余三个方法中，均方误差的平方根和可信区间平均长度均相差不大，因此根据标准偏倚绝对值大小便可判断出不同缺失机制组合下的填补方法优劣。

综上，各种情况下不同方法的推荐结果如表4。

表2 四种填补方法下样本选择模型结果等式的回归系数估计值的各项评价标准比较(二)

表3 四种填补方法下样本选择模型结果等式的回归系数估计值的各项评价标准比较(三)

表4 不同缺失机制组合下的填补方法选择

讨 论

数据缺失现象在调查研究中非常普遍，它不仅会降低参数估计的效率，同时也给统计分析带来很大偏倚。根据数据缺失机制，可将数据缺失分为三类：完全随机缺失(MCAR)、随机缺失(MAR)和非随机缺失(NMAR)。针对MAR机制，统计学家们提出了多种方法来校正这种缺失带来的偏倚，MI就是被广为推崇的方法之一；针对NMAR机制，由于该机制的复杂性，当前还没有一种统一的方法来校正这种偏倚，不过当回归模型中的应变量为非随机缺失时，某些情况下可以应用样本选择模型来纠正这种NMAR带来的偏倚；但当两种缺失机制并存时的偏倚纠正方法尚未见有介绍。对此，本研究提出了一个两阶段策略以纠正不同缺失机制造成的偏倚。第一阶段首先利用只包含随机缺失数据的个体对单纯无应答缺失按照MAR机制进行多重填补，在第二阶段中使用填补后的多个数据集与选择性偏倚导致的因变量为虚假0(视为缺失)数据合并进行样本选择模型分析以校正由于非随机缺失所造成的偏倚，最后对多个样本选择模型拟合结果进行合并。

模拟研究结果表明：当非随机缺失为轻度时，PS法由于标准偏倚绝对值远远超过了规定的界值，所以该法的结果相对不理想；而MCMC、EMB和PMM法均得出较好的结果。不同程度随机缺失情况下的填补方法选择为：随机缺失也为轻度时，MCMC法最好；随机缺失为中度时，EMB法最好；在随机缺失为重度时，PMM法最好。

当非随机缺失为中度时，PS法由于标准偏倚绝对值远远超过了规定的界值，所以仍不可取，而MCMC、EMB和PMM法均得出较好的结果。此时，无论随机缺失程度如何，MCMC法都是最好的方法。

当非随机缺失为重度时，PS法由于标准偏倚绝对值远远超过了规定的界值，所以仍不可取，而MCMC、EMB和PMM法均得出较好的结果。此时，无论随机缺失程度如何，PMM法都是最好的方法。

本文以医疗费用调查研究中可能出现的两种缺失为假设背景，探索性地提出两阶段策略纠正这两种偏倚，希望能为以后在缺失值处理方面的应用提供一些方法学依据。

参 考 文 献

1.Fisher ES,Bynum JP,Skinner JS.Slowing the growth of health care costs—lessons from regional variation.New England Journal of Medicine,2009,360(9):849-852.

2.卫生部统计信息中心.第四次国家卫生服务调查主要结果.[cited 2010年3月16日];http://www.moh.gov.cn/publicfiles/business/htmlfiles/mohbgt/s3582/200902/39201.htm].

3.Baer OCJ.Bradley JC,et al.Testing and correcting for non-random selection bias due to censoring:an application to medical costs.Health Services and Outcomes Research Methodology,2003,4(2):93-107.

4.Peytchev A,Baxter RK,Carley-Baxter LR.Not All Survey Effort is Equal.Public Opinion Quarterly,2009,73(4):785-806.

5.Rubin DB.Inference and missing data.Biometrika,1976,63(3):581-592.

6.Little RJA.Rubin DB.Statistical analysis with missing data.2nd.Vol.2.2002:Wiley New York:2002.

7.薛小平,史东平，王彤.受限因变量模型及其半参数估计.中国卫生统计,2007,24(2):211-213.

8.张磊.样本选择模型的似然估计与两步估计.现代预防医学,2007,34(9):1607-1609.

9.张磊.样本选择模型及其在医疗费用研究中的应用.[硕士学位论文].山西：山西医科大学,2007.

10.Donald B.Multiple Imputation for Nonresponse in Surveys.American Journal of Sociology,1987,76:346.

11.Honaker J,King G.What to Do about Missing Values in Time‐Series Cross‐Section Data.American Journal of Political Science,2010,54(2):561-581.

12.Dempster AP,Laird NM,Rubin DB.Maximum likelihood from incomplete data via the EM algorithm.Journal of the Royal Statistical Society.Series B (Methodological),1977:1-38.

13.Tanner MA,Wong WH.The calculation of posterior distributions by data augmentation.Journal of the American Statistical Association,1987,82(398):528-540.

14.Schunk D.A Markov chain Monte Carlo multiple imputation procedure for dealing with item nonresponse in the German save survey.2007.

15.Co1lins LM，Schafer JL，Kam CM.A Comparison of Inclusive and Restrictive Strategies in Modem Missing Data Procedures.Psychological methods,2001,6(4):330.

16.Burton A,et al.The design of simulation studies in medical statistics.Statistics in Medicine,2006,25(24):4279-4292.