Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:黄海嫦 (中山大学)
邮箱:huanghch23@mail2.sysu.edu.cn
编者按:本文部分整理自 PAUL ALLISON 的博客 Is Dummy Variable Adjustment Ever Good for Missing Data?,特此致谢!
目录
社会调查数据难免会有缺失值,常见情形之一是问卷题目跳转。比如,未婚群体受访时会跳过「婚姻幸福度」等类似题目。一般而言,缺失值有实际值,但由于受访者拒答、隐瞒或谎报造成答案不合理、追踪调查数据后期未能成功追访等原因无法观测。目前,数据缺失大致分为以下三种类型:
在数据缺失类型中,完全随机缺失和随机缺失属于可忽略的缺失,而非随机缺失属于不可忽略的缺失。总的来说,由于数据缺失原因复杂、影响难测,对缺失值进行合理处理就成为一项重要工作。目前,缺失数据的处理方法常被分为三类:加权法、删除法和插补法。
在处理数据缺失问题上,虚拟变量调整法 (DVA, Dummy variable adjustment) 曾在许多年间大行其道。它主要针对回归分析中不完全的自变量,适用于任何类型的回归——线性、logistic、Cox 等,易于理解与实现。
然而, DVA 并不正确。Michael Jones (1996) 证明了:即使数据的缺失是完全随机的,DVA 也常常会产生有偏的系数估计。由此,目前在回顾缺失值处理方法时,DVA 常被弃置一旁。
本文将重新介绍 DVA。尽管 Jones 一再强调 DVA 的缺陷,但 Allison 本人认为,仍然存在两种情况,可以使 DVA 成为处理缺失值的有效方法。这种方法,甚至会优于 多重插补法和 (完全信息下的) 最大似然法。
假设现有因变量
从本质来看,常数 C 的选择不影响回归模型本身。准确地说,无论 C 是多少,
需要注意的是,Allison 本人在发表于 2010 年的一篇文章 中曾较正式地介绍了虚拟变量调整法中将缺失值替换为 0 这一特殊情形,即把上文中的常数 C 设为 0。
实际上,当有足够理由说明缺失数据的实际值为 0 时,可以将缺失值替换为 0。当理由不充分时,这一做法不可行。详见 When should missing data, in numerical variables, be replaced by zeros?。
在本文中,Allison 进一步指出,
模型给出了方程中系数的估计值。其中,
如果将 C 设为
如果变量
虚拟变量调整法优点众多:
但正如 Jones 所指出的,DVA 所得系数往往是有偏的。这不难理解。如果用任意常数值替换掉
不过,Allison 在 Blog 文章中提出了两种例外情况。在这两种情况下,Jones 的结论不足以完全否定 DVA。
第一种情况:
Jones 的证明有一个隐含假设:存在未被观察到的
假设一个回归模型,其中,
Married:
Unmarried:
这两个方程之间的明显区别是
分别估计已婚和未婚人士这两个方程是可行的。但
为了解决这个问题,本文直接假设
这个方程包含了
还要注意,对于
其实还可以再往回归模型中加入
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 也在此处做了一个小模拟,在
在模拟中,
至于标准误差,根据 DVA 的估计,
结论是,对于基线协变量数据缺失的 RCT,DVA 至少和最大似然法或多重插补法一样好,甚至可能更好。
附录1:针对“不适用”情况下的 DVA 效用评估的模拟
基于上述已婚和未婚人群的模型,生成回归方程,如下:
Married:
Unmarried:
其中
由于只关心偏差,而不关心抽样变异,可以生成一个大样本,其中有 10,000 个样本已婚,10,000 个样本未婚。
下一步是使用三种不同的方法——DVA、多重插补法和完全信息下的最大似然法,估计
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 效用评估的模拟
假设一个回归方程,如下:
其中
为了得到标准误的无模型估计,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
在
因此,Allison 进行了小调整:假定
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 估计下的
同样,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)
蒙特卡洛模拟方法 (MC),即从总体中抽取大量随机样本的计算方法。当根据总体的分布函数很难求出想要的数字特征时,可以使用蒙特卡洛模拟的方法,从总体中抽取大量样本,使用样本的数字特征估计总体的数字特征。
针对原文中 DVA 适用的第一种情况:
//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
Note:产生如下推文列表的 Stata 命令为:
lianxh 缺失值, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh