缺失值处理：虚拟变量调整法靠谱吗？

作者：黄海嫦 (中山大学)
邮箱：[e-mail]

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：

 

编者按：本文部分整理自 PAUL ALLISON 的博客 Is Dummy Variable Adjustment Ever Good for Missing Data?，特此致谢！

---

目录

• 1. 引言

• 2. 简单了解一下DVA

• 3. DVA 适用的的两种情况

• 4. 对原文附录内容的展示

• 5. 笔者补充佐证

• 6. 相关推文

---

1. 引言

社会调查数据难免会有缺失值，常见情形之一是问卷题目跳转。比如，未婚群体受访时会跳过「婚姻幸福度」等类似题目。一般而言，缺失值有实际值，但由于受访者拒答、隐瞒或谎报造成答案不合理、追踪调查数据后期未能成功追访等原因无法观测。目前，数据缺失大致分为以下三种类型：

• 完全随机缺失 (Missing Completely at Random, MCAR)：表现为出现的数据缺失与已观测到的数据无关，并且与未观测到的数据也无关。例如，由于测量设备出故障，导致某些值缺失。
• 随机缺失 (Missing at Random, MAR)：此时出现的数据缺失的现象与已观测到的数据有关，而与未观测到的数据无关。例如，人们是否透露收入可能与性别、教育程度、职业等因素有关系。如果除了收入数据缺失，其他因素都能被观测到，并且收入缺失情况在不同性别、教育程度、职业的人群内存在差异，而与收入本身的值无关——那么收入就是随机缺失的。
• 非随机缺失 (Not Missing at Random, NMAR)：出现的数据缺失现象与未观测的数据有关。例如，在控制了性别、教育程度、职业等已观测因素之后，如果收入是否缺失还依赖于收入本身的值，像是高收入人群倾向于隐瞒自己的收入，那么收入就是非随机缺失的。

在数据缺失类型中，完全随机缺失和随机缺失属于可忽略的缺失，而非随机缺失属于不可忽略的缺失。总的来说，由于数据缺失原因复杂、影响难测，对缺失值进行合理处理就成为一项重要工作。目前，缺失数据的处理方法常被分为三类：加权法、删除法和插补法。

• 插补法又包括统计学插补法和机器学习插补法。
• 如果缺失数据是由单元无回答 (被调查者不愿意或者不能够回答整张问卷) 造成的，那么常用的方法是“加权法” (通过增加回答者的权重来弥补无回答，以减小无回答带来的偏倚)。
• 如果缺失数据是由项目无回答 (被调查者拒绝回答个别的调查项目) 造成的，那么常用的方法是“删除法”或者“插补法”。

在处理数据缺失问题上，虚拟变量调整法 (DVA, Dummy variable adjustment) 曾在许多年间大行其道。它主要针对回归分析中不完全的自变量，适用于任何类型的回归——线性、logistic、Cox 等，易于理解与实现。

然而， DVA 并不正确。Michael Jones (1996) 证明了：即使数据的缺失是完全随机的，DVA 也常常会产生有偏的系数估计。由此，目前在回顾缺失值处理方法时，DVA 常被弃置一旁。

本文将重新介绍 DVA。尽管 Jones 一再强调 DVA 的缺陷，但 Allison 本人认为，仍然存在两种情况，可以使 DVA 成为处理缺失值的有效方法。这种方法，甚至会优于 多重插补法和 (完全信息下的) 最大似然法。

2. 简单了解一下DVA

假设现有因变量 对两个自变量 和 进行线性回归。其中， 是一个不完全变量，有大量的缺失值。DVA 主要是指在进行回归之前，针对不完全自变量 的预处理：

• 首先，创建一个虚拟变量 。当该样本的 未缺失时，取值为 1，若缺失则为 0；
• 其次，将 所有缺失值替换为某个常数 C；
• 最后，将三个变量 (， 和 )，作为自变量进入回归模型。

从本质来看，常数 C 的选择不影响回归模型本身。准确地说，无论 C 是多少， 和 的系数都不会改变。

需要注意的是，Allison 本人在发表于 2010 年的一篇文章 中曾较正式地介绍了虚拟变量调整法中将缺失值替换为 0 这一特殊情形，即把上文中的常数 C 设为 0。

实际上，当有足够理由说明缺失数据的实际值为 0 时，可以将缺失值替换为 0。当理由不充分时，这一做法不可行。详见 When should missing data, in numerical variables, be replaced by zeros?。

在本文中，Allison 进一步指出， 的系数确实依赖于 C。在更一般的情况下，为了更好地解释 的系数，一个较好的 C 的设置是取 X 非缺失b部分的均值。回归模型如下：

模型给出了方程中系数的估计值。其中， 是控制了自变量 和关于变量 缺失情况的虚拟变量 后，变量 对 的影响。 是控制了变量 后，在可被观测到 的样本中，变量 对 的影响。

如果将 C 设为 非缺失部分的均值， 就是：在控制了其他变量的情况下，“具 非缺失部分均值的个体” 其 减去 “缺失 数据的个体” 其 的差值。也就是说，表示缺失情况的虚拟变量系数衡量着：“具备完全数据的样本”所对应的因变量的预测值，与“具备不完全数据的样本”所对应的因变量的预测值之间差异的大小。

如果变量 上也有缺失值呢？同上，创建另一个虚拟变量以表示 的缺失情况，然后用 非缺失部分的均值 (也可以是任何其他常数) 来替换 的所有缺失值。新的虚拟变量也进入回归模型。虚拟变量调整法还适用于存在更多 “不完全” 变量的情况。

虚拟变量调整法优点众多：

• 比较简单；
• 适用于任何类型的回归；
• 没有样本被忽略或舍弃；
• 所有可用的信息都被纳入了回归模型；
• 如果满足误差项的一般假设，虚拟变量调整法所估计的标准误将是“准确”的。(这一点很重要。许多传统方法，如成对删除和单一插补，在估计标准误方面不尽如人意。)

但正如 Jones 所指出的，DVA 所得系数往往是有偏的。这不难理解。如果用任意常数值替换掉 的缺失值，就会降低该变量的变异真实性。如果 和 相关， 的系数估计就有偏差，且进一步导致 的系数估计有偏。而且，即使数据的缺失是完全随机的，这种情况也会发生。从这个角度看，DVA 不能被称之为一个处理缺失值问题的好方法。

不过，Allison 在 Blog 文章中提出了两种例外情况。在这两种情况下，Jones 的结论不足以完全否定 DVA。

3. DVA 适用的的两种情况

第一种情况： 上的数据缺失是因为该变量对样本的这个子集“不适用”或没有意义。

Jones 的证明有一个隐含假设：存在未被观察到的 的真实值，且这些缺失值对因变量有影响。但是，如果这个假设不成立呢？

假设一个回归模型，其中， 是抑郁症的衡量标准， 是婚姻满意度的衡量标准。样本中有 30% 的人未婚，他们的婚姻满意度是缺失的。此时，用多重插补法来补全这 30% 样本的婚姻满意度没有意义。更好的方法是：假设数据生成机制由两个方程组成，一个用于已婚人士，另一个用于未婚人士：

Married：

Unmarried：

这两个方程之间的明显区别是 (婚姻满意度) 没有出现在未婚人群的方程中。此外， 变量包括任何其他潜在的影响抑郁程度的因素，比如教育。

分别估计已婚和未婚人士这两个方程是可行的。但 (或任何其他协变量) 对已婚和未婚人群的影响理论上应该是相同的。相比较于联立两个方程进行估计，分别估计两个方程会得到更大的标准误。另外，如果存在一个以上的“不适用” 变量，就很有可能需要估计几个不同的方程。

为了解决这个问题，本文直接假设 的系数对于已婚和未婚的人是相同的，即 。正如前文所述，当 未缺失 (即对应样本为已婚人士)，有 ，否则为 0。结合两个方程，有：

这个方程包含了 的“主效应”和 与 的相互作用。 与 相乘，相当于当 缺失时取 0。因而，这也只是 时的 DVA 模型。此时，如果 满足线性回归的一般假设，应用于该模型的 OLS 将产生所有系数的 最佳无偏估计。

还要注意，对于 ， 的系数是原始的两个方程 (已婚和未婚) 的截距之差。如果 被设为 非缺失 b 部分的均值，那么 的系数正如前文所释：在这个例子中，它将是在控制了模型中的任何其他变量后，处于婚姻满意度平均值的已婚人士所对应的 的期望值，减去未婚人士所对应的Y的期望值所产生的差值。

其实还可以再往回归模型中加入 和 的相互作用。不过，这样一来，如果不限制误差方差相等，就相当于分别估计了两个方程。

DVA 很好地估计了“不适用”情况下一个合理模型的系数。直观上，最大似然法或多重插补法难以做得比 DVA 更好。 Allison 在此处做了一个小型模拟，证明了这一点——相比较 DVA，最大似然法或多重插补法产生的估计都有很大的偏差。(详见附录 1)

第二种情况：随机对照试验中基线协变量数据缺失

在分析随机对照试验 (RCT) 的数据时，控制了基线测量的一个或多个协变量后，依据协方差分析 (ANCOVA) 框架估计处理效应的做法，是非常普遍的。即使更为简单的方差分析 (ANOVA) 可以得到处理效应的无偏估计，但纳入协变量以减少模型中的误差方差，继而减少效应估计的标准误，依然十分值得。此外，这一做法还增加了检验处理效应的统计效力。

然而，基线协变量常有缺失的数据。而删除缺失数据的样本背离了处理缺失值的预期目标，因为这往往会增加标准误，也违反了 意向性分析原则 (intention-to-treat，ITT)。在这种情况下，DVA 是一种有效的方法 (Puma et al. 2009, Groenwold et al. 2012)。正如 Jones 在 1996 年的论文中所指出的，只有当缺失数据的变量与模型中的其他变量相关时，误差才会成为 DVA 的问题。而在随机对照试验中，基线协变量与处理效应在假定上不相关。因此，DVA 不会对处理效应产生任何偏差。

而多重插补或最大似然，在随机对照试验中似乎也是合理的选择。DVA 与这些方法相比如何？Allison 也在此处做了一个小模拟，在 和 不相关的前提下，用这三种方法来估计 对 和 的线性回归模型。(详见附录 2)。

在模拟中， 是二分处理指标， 是缺失数据的基线协变量。结果表明，当 上的数据完全随机缺失时，三种方法对 和 的系数都是近似无偏的。

至于标准误差，根据 DVA 的估计， 的系数的标准误比最大似然和多重插补的标准误差大 17% 左右。然而这个系数在本文中并不重要。Allison 关心的是处理效应，即 的系数。而在这个系数上，DVA 产生的标准误小于最大似然或多重插补 (虽然这个差异并不大)。不过，当 非随机缺失时，差异会变得更大。

结论是，对于基线协变量数据缺失的 RCT，DVA 至少和最大似然法或多重插补法一样好，甚至可能更好。

4. 对原文附录内容的展示

附录1：针对“不适用”情况下的 DVA 效用评估的模拟

基于上述已婚和未婚人群的模型，生成回归方程，如下：

Married：

Unmarried：

其中 和 都服从标准正态分布，且二者相关性 。两个回归模型的 都服从正态分布且独立于 和 。

由于只关心偏差，而不关心抽样变异，可以生成一个大样本，其中有 10,000 个样本已婚，10,000 个样本未婚。 缺失 (假定为完全随机缺失) 的概率是 0.50。

下一步是使用三种不同的方法——DVA、多重插补法和完全信息下的最大似然法，估计 对 和 的回归方程。数据的生成和分析都是用 R 语言完成的。

Allison 使用 lm 函数进行 DVA 估计，使用 lavaan 包进行最大似然法估计，使用 jomo 包进行多重插补法估计 (25 个数据集通过多元正态 MCMC 方法进行插补)。

结果如下：

                DVA      Maximum Likelihood  Multiple Imputation
Variable  Estimate  S.E.  Estimate   S.E.      Estimate    S.E.
X          2.010    .020    1.846    .019        1.844     .019
Z          2.006    .015    1.780    .017        1.781     .017

由于真实系数都是 2.0，很明显，DVA 产生了近似无偏的结果。正如所预料到的，最大似然法和多重插补法产生的估计值彼此非常相似。然而，这两组估计值都比真实值低 8-10%。

Allison本 人还补充，这一模拟只能被视为尝试性的。准确地说，他只使用了一组参数值、一个协变量、一个线性模型，只探讨了一种数据缺失类型。但最起码，它直指多重插补法和最大似然法可能会对一个大致合理的“不适用”模型产生有偏的估计。

本节 R 代码：

#Appendix 1
library(lavaan)
library(jomo)
library(mitools)
library(mice)

#Generate data
set.seed(3495)
n <- 10000
z_mar <- rnorm(n)
x_mar <- .5*z_mar + sqrt(1-.5*.5)*rnorm(n)
y_mar <- 3 + 2*z_mar + 2*x_mar + 2*rnorm(n)
z_sing <- rnorm(n)
x_sing <- rep(NA, n)
y_sing <-  4 + 2*z_sing + 2*rnorm(n)
dat <- data.frame(y=c(y_mar, y_sing),
                 z=c(z_mar, z_sing),
                 x=c(x_mar, x_sing))
dat$d <- ifelse(is.na(dat$x), 0, 1)
dat$x.dva <- ifelse(is.na(dat$x), 0, dat$x)

#DVA method
mod.dva <- lm(y ~ x.dva + z + d, data=dat)
summary(mod.dva)

#Maximum likelihood
mod.ml <- sem('y ~ x + z', data=dat, missing='fiml.x')
summary(mod.ml)

#Multiple imputation
datsub <- subset(dat,select=-c(d,x.dva))
outjomo <-jomo(datsub, nimp=25)
mi_list <- imputationList(split(outjomo, outjomo$Imputation)[-1])
mi_results <- with(mi_list, lm(y ~ x + z))
summary(pool(as.mira(mi_results)),type="all")

附录2：针对 RCTs 情况下的 DVA 效用评估的模拟

假设一个回归方程，如下：

其中 为虚拟变量——当样本属于实验组时，取值为 1；属于对照组取值为 0，不同取值的概率均为 0.50。 和 都服从标准正态分布，彼此独立，也都与 无关。 缺失 (假定为完全随机缺失) 的概率是 0.50。

为了得到标准误的无模型估计，Allison 生成了 1000 个样本集，每个样本集大小为 200。对于每个样本集，采用 DVA、最大似然法和多重插补法对模型进行估计。

所有的数据生成和分析都是用 R 语言完成的。Allison 使用内置的 lm 函数来完成 DVA 估计，使用 lavaan 包来进行最大似然法估计，使用 jomo 包来进行多重插补法估计 (每个样本集取 10 个数据集合)。最大似然法和多重插补法均建立在多元正态假设的基础上，在设计上适用这些数据。

重复样本下的系数估计均值，以及这些估计的标准差 (估计的标准误) 如下：

                DVA       Maximum Likelihood  Multiple Imputation
Variable  Estimate  S.E.   Estimate   S.E.      Estimate    S.E.
X          1.987    .209     2.011    .177        1.975     .180
Z          2.004    .340     1.998    .344        1.979     .345

在 和 的系数上，这三种估计方法大致近似无偏。对于 的系数，最大似然法和多重插补法产生的标准误明显小于 DVA。但正如前面提到的，我们不关心这个系数。在处理效应上，DVA 的标准误差最小。可惜差异不大，并且很可能在蒙特卡罗误差范围内。

因此，Allison 进行了小调整：假定 不是随机缺失的。具体地说，假定 小于 0 的所有值都缺失，而大于 0 的所有值都能被观察到。结果如下：

                DVA       Maximum Likelihood  Multiple Imputation
Variable  Estimate  S.E.   Estimate   S.E.      Estimate    S.E.
X          1.997    .346     2.612    .333        2.519     .339
Z          1.987    .296     1.995    .377        1.981     .376

现在，DVA 估计下的 的系数是无偏的；但最大似然法和多重插补法下的系数估计明显比真实值高 (约 30%)。三种方法在处理效应上均近似无偏，但 DVA 的标准误比其他两种方法低 20% 以上。

同样，Allison认为，这一模拟只能被视为尝试性的。准确地说，他只使用了一组参数值、一个样本量、一个协变量、一个线性模型，且遗漏了两种数据缺失类型。但最起码，它表明了，DVA 是处理 RCTs 中基线协变量缺失值问题的有力竞争者。

本节代码

#Appendix 2
library(lavaan)
library(jomo)
library(mitools)
library(mice)
library(psych)
set.seed(23599)

#MCAR Scenario
#Create empty data frames to hold coefficients
dva <- data.frame(matrix(ncol = 4, nrow = 0))
   colnames(dva) <-c("(Intercept)", "x.dva","z","d")
ml <- data.frame(matrix(ncol = 4, nrow = 0))
   colnames(ml) <-c("y~x", "y~z","y~~y","y~1")
mult <- data.frame(matrix(ncol = 3, nrow = 0))
   colnames(mult) <-c("(Intercept)", "x","z")

for (i in 1:1000) {
  #Generate MCAR Data
  n <- 200
  z <- as.numeric(rnorm(n)>0)
  x <- rnorm(n)
  y <- 2*z + 2*x + 2*rnorm(n)
  xmiss <- ifelse(rnorm(n)<0,NA,x)
  dat <- data.frame(y=y,x=xmiss,z=z)

  #DVA method
  dat$d <- ifelse(is.na(dat$x), 0, 1)
  dat$x.dva <- ifelse(is.na(dat$x), 0, dat$x)
  mod.dva <- lm(y ~ x.dva + z + d, data=dat)
  dva[nrow(dva) + 1,] <- mod.dva$coefficients

  #Maximum likelihood
  mod.ml <- sem('y ~ x + z', data=dat, missing='fiml.x')
  ml[nrow(ml) + 1,] <- coef(mod.ml)

  #Multiple imputation
  datsub <- subset(dat,select=-c(d,x.dva))
  outjomo <-jomo(datsub, nimp=10, nburn=100, nbetween=20, output=2)
  mi_list <- imputationList(split(outjomo, outjomo$Imputation)[-1])
  mi_results <- with(mi_list, lm(y ~ x + z))
  mult[nrow(mult) + 1,] <- getqbar(pool(as.mira(mi_results)))
}

print(describe(dva, fast=T),digits=3)
print(describe(ml, fast=T),digits=3)
print(describe(mult, fast=T),digits=3)


#NMAR scenario
set.seed(4321)
#Create empty data frames to hold coefficients
dva <- data.frame(matrix(ncol = 4, nrow = 0))
  colnames(dva) <-c("(Intercept)", "x.dva","z","d")
ml <- data.frame(matrix(ncol = 4, nrow = 0))
  colnames(ml) <-c("y~x", "y~z","y~~y","y~1")
mult <- data.frame(matrix(ncol = 3, nrow = 0))
  colnames(mult) <-c("(Intercept)", "x","z")

for (i in 1:1000) {
  #Generate NMAR Data
  n <- 200
  z <- as.numeric(rnorm(n)>0)
  x <- rnorm(n)
  y <- 2*z + 2*x + 2*rnorm(n)
  xmiss <- ifelse(x<0,NA,x)
  dat <- data.frame(y=y,x=xmiss,z=z)

  #DVA method
  dat$d <- ifelse(is.na(dat$x), 0, 1)
  dat$x.dva <- ifelse(is.na(dat$x), 0, dat$x)
  mod.dva <- lm(y ~ x.dva + z + d, data=dat)
  dva[nrow(dva) + 1,] <- mod.dva$coefficients

  #Maximum likelihood
  mod.ml <- sem('y ~ x + z', data=dat, missing='fiml.x')
  ml[nrow(ml) + 1,] <- coef(mod.ml)

  #Multiple imputation
  datsub <- subset(dat,select=-c(d,x.dva))
  outjomo <- jomo(datsub, nimp=10, nburn=100, nbetween=20, output=2)
  mi_list <- imputationList(split(outjomo, outjomo$Imputation)[-1])
  mi_results <- with(mi_list, lm(y ~ x + z))
  mult[nrow(mult) + 1,] <- getqbar(pool(as.mira(mi_results)))
}

print(describe(dva, fast=T),digits=3)
print(describe(ml, fast=T),digits=3)
print(describe(mult, fast=T),digits=3)

5. 笔者补充佐证

蒙特卡洛模拟方法 (MC)，即从总体中抽取大量随机样本的计算方法。当根据总体的分布函数很难求出想要的数字特征时，可以使用蒙特卡洛模拟的方法，从总体中抽取大量样本，使用样本的数字特征估计总体的数字特征。

针对原文中 DVA 适用的第一种情况： 上的数据缺失是因为该变量对样本的这个子集“不适用”或没有意义，笔者沿用了 Allison 在附录 1 中的回归模型及相关假设，基于 Stata 生成数据并进行分析，模拟所得不完全变量的系数和标准误如下：

//DVA
    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
           b |     10,000    2.017304           0   2.017304   2.017304
          se |     10,000    .0214553           0   .0214553   .0214553

该结果与原文附录 1 相关分析相容，起到了进一步佐证的作用。

针对原文中 DVA 适用的第二种情况，笔者沿用了 Allison 在附录 2 中的回归模型及相关假设，基于 Stata 生成数据并进行分析，当假定 数据完全随机缺失时，模拟所得代表处理效应的系数估计及其标准误如下：

//DVA
    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
           b |      1,000    1.992375           0   1.992375   1.992375
          se |      1,000    .0283461           0   .0283461   .0283461

该结果与原文附录 2 相关分析相容，起到了进一步作证的作用。

本节 Stata 代码如下：

* 针对第一种情况

//生成数据集
clear
set obs 10000
set seed 123
gen z = rnormal()
gen x = 0.5 * z + sqrt(1-0.5*0.5) * rnormal()
gen y = 3 + 2 * z + 2 * x + 2 * rnormal() 
save "…\data00.dta", replace //生成未缺失了X的部分数据并保存
clear
set obs 10000
set seed 123
gen z = rnormal()
gen x = .
gen y = 4 + 2 * z + 2*rnormal() //生成缺失了X的部分数据
append using "…\data00.dta"     //合并两个数据集
save "…\data01.dta", replace    //获得符合 Alllison 在附录 1 中所提出的所有假定的数据集

//DVA-虚拟变量调整法
use data01.dta, clear
set seed 12345
gen dum = 0 if x == .
replace dum = 1 if dum == .
gen xdum = x
replace xdum = 0 if xdum == .
tempname sim //定义临时的内存区域的名字
tempfile results //定义临时的保存文件的名字
postfile `sim' b se using "`results'", replace
quietly{
  forvalues i = 1/10000 { //循环语句，做 10000 次蒙特卡洛模拟
    reg y xd z dum 
    post `sim' (_b[xd]) (_se[xd])
  }
}
postclose `sim'
use "`results'", clear
sum

* 针对第二种情况
//生成数据集
clear
set obs 100
set seed 123
gen z = 0
save "…\data02.dta", replace
clear
set obs 100
set seed 123
gen z = 1
append using "…\data02.dta"
save "…\data03.dta", replace
clear
set obs 100
set seed 123
gen x = .
save "…\data04.dta", replace
clear
set obs 100
set seed 123
gen x = rnormal()
append using "…\data04.dta"
save "…\data05.dta", replace
cross using "…\data03.dta" 
save "…\data06.dta", replace

//DVA-虚拟变量调整法
clear
use "…\data06.dta"
set seed 12345
gen dum = 0 if x == .
replace dum = 1 if dum == .
gen xdum = x
replace xdum = 0 if xdum == .
gen y = 2 * z + 2 * x + 2 * rnormal()
tempname sim     //定义临时的内存区域的名字
tempfile results //定义临时的保存文件的名字
postfile `sim' b se using "`results'", replace
quietly {
  forvalues i = 1/1000{
    reg y x z dum
    post `sim' (_b[z]) (_se[z])
  }
}
postclose `sim'
use "`results'", clear
sum

6. 相关推文

Note：产生如下推文列表的 Stata 命令为：
  lianxh 缺失值, m
安装最新版 lianxh 命令：
  ssc install lianxh, replace

• 专题：数据处理
• Stata数据处理：缺失值类型及应对方法
• Stata数据处理：缺失值填充-autofill-carryforward
• Stata：fillmissing-缺失值填充-数值和文字的前后填充！
• Stata数据处理：缺失值与多重补漏分析（一）
• Stata数据处理：缺失值与多重补漏分析（二）
• Stata数据处理：缺失值与多重补漏分析（三）
• Stata：缺失值与多重补漏-misstable-D204
• 缺失值能否用零代替？-L117
• Stata：让缺失值一览无余
• Stata缺失值专题：多重补漏分析
• Stata：缺失值的填充和补漏
• 专题：面板数据
• Stata：面板数据缺失值与多重补漏分析-twofold

New！ Stata 搜索神器：lianxh 和 songbl  GIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
  . ssc install lianxh
  . ssc install songbl
👉  使用：
  . lianxh DID 倍分法
  . songbl all

🍏 关于我们

• 连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。