IV和Matching老矣, “弹性联合似然法”成新趋势

Original 因果推断研究小组计量经济圈 2021-10-23

凡是搞计量经济的，都关注这个号了

邮箱：econometrics666@sina.cn

所有计量经济圈方法论丛的code程序, 宏微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.

圈子在上一日，引荐了"Queens"女王喜欢对别人开战? 已婚vs未婚。这篇文章发表在Top 5上的JPE，文章标题为“Queens”，不仅霸气而且非常有趣。里面对工具变量的应用和对机制分析，值得各位学者认真研习。

”因果推断研究小组“引荐了很多因果推断文章，受到倾向于做实证研究的社科学者青睐。毕竟内生性问题是每一位学者都要面对和处理的难题，它关乎着因果关系推断的命运。下面，咱们列举10篇代表性文章，有进一步研习需求的学者，可以到公众号搜索关键词。这40个微观数据库够你博士毕业了, 反正凭着这些库成了教授

1.广义PSM,连续政策变量因果识别的不二利器

2.PSM倾向匹配Stata操作详细步骤和代码

3.看完顶级期刊文章后, 整理了内生性处理小册子

4.处理效应模型选择标准，NNM和PSM

5.因果推断中的匹配方法:最全回顾和前景展望

6.内生性和倾向得分匹配, 献给准自然试验的厚礼

7.倾向值匹配与因果推论,史上最全面精妙的锦囊

8.匹配还是不匹配？这真是个值得考虑的问题

9.匹配方法(matching)操作指南, 值得收藏的16篇文章

10.倾向匹配分析深度(Propsensity matching)

今天，我们将为学者们引荐一个处理内生性问题的无冕之王，尤其是在处理”不可观测混淆变量“所带来的内生性问题方面。

内生性问题通常由三种因素导致：遗漏变量、测量误差和双向因果。如果不理解这些问题，可以读一读这篇文章”计量分析中的内生性问题综述,一篇不得不读的经典作品“。关于怎么处理测量误差问题所导致的内生性问题，从而让估计的结果得到修正，可以看看”测量误差消除, 直接和间接效应计算, 多数据和指标使用的方法“和”测量误差和新回归法则，引子“。对于内生性问题，在当前最普遍的处理方法是”工具变量法“和”匹配法“(实际上遗漏变量难以彻底根除)。

如果想要进一步对IV方法处理内生性问题有所理解，可以看看”工具变量IV与内生性处理的精细解读“、我的"工具变量"走丢了，寻找工具变量思路手册、”内生性处理的秘密武器－工具变量估计“、”如何寻找工具变量？和”得工具者得实证计量“、工具变量在社会科学因果推断中的应用”、”AEA期刊的IV靠不靠谱？“、”IV和GMM相关估计步骤，内生性、异方差性等检验方法、面板数据“、”面板数据、工具变量选择和HAUSMAN检验的若干问题“。看了这么多文章后发现，工具变量法在使用的过程中会出现很多问题，比如工具变量与内生变量的相关性微弱所导致的”弱工具变量问题“(小组会尽快推出关于这个的解决方法，社群里也有推荐的重要文献)，又比如工具变量不满足”排除约束“，即工具变量只通过内生变量去影响结果变量。除了这两大问题，工具变量方法在使用中经常被忽视的另一大问题，是工具变量进行因果推断时的效力(power)往往是与方程形式(function form)紧密关联的，文献表明即使微小的方程形式差异也会带来巨大的估计偏差。

正是因为严重依赖于方程形式的IV估计方法所带来的问题，Matching方法就成为了最近非常火的因果推断新方法。Matching依靠observable covariates在处理组与控制组进行匹配，然后获得一个不依赖于模型而是强调数据本身特征的政策处理效应。如果想要进一步对Matching方法处理内生性问题有所理解，可以看看“广义PSM,连续政策变量因果识别的不二利器”、“PSM倾向匹配Stata操作详细步骤和代码”、“PSM, RDD, Heckman模型的操作程序”、“处理效应模型选择标准，NNM和PSM”、“因果推断中的匹配方法:最全回顾和前景展望”、“内生性和倾向得分匹配, 献给准自然试验的厚礼”、“倾向值匹配与因果推论,史上最全面精妙的锦囊”、“匹配还是不匹配？这真是个值得考虑的问题”、“匹配比OLS究竟好在哪里？这是一个问题”、”倾向匹配分析深度(Propsensity matching)“和”PSM和马氏匹配已淘汰, '遗传匹配'成因果推断匹配之王“。

但是，匹配方法所依赖的”selection on observables“实际上是一个非常强的假设，在很多情况下表现得并不是很好。在社会科学研究中，有很多变量不能直接衡量，因此我们倾向于使用一些proxy变量作为近似代理变量。而且，对于有些变量而言，他们的数据收集成本过高，因此我们很可能把他们在问卷中就直接省略掉了。在这种情况下，如果只是使用观测到的变量(而忽视未观测到的混杂变量)进行匹配，那么得到的政策处理效应必定存在着很大的问题。

这就出现了经常所说的“遗漏变量问题”。在实证研究中，我们能够做的就是加入足够多的变量去逐一消除遗漏变量所带来的估计偏误问题，尽管不可能得到一个纯粹的无偏误估计量。但是，遗漏变量问题中最难以解决的是“不可观测的混淆变量”部分，毕竟，变量的不可观测性使我们没有办法去很好地衡量并有效地控制他们的影响。举一个例子，劳动经济学当中经常实证检验“教育回报率”，其中教育的内生性问题在经济学家中是得到公认的。内生性的来源是，有一个不可观测到的混淆变量——能力，会同时影响个体的受教育程度和将来的收入。不过，能力这个指标是很难去客观衡量的，如果此时贸然用诸如“工作经验”、“年龄”和“IQ”这些指标进行量化，那有很大一部分未能控制的“不可观测的能力”将会进入到残差项，从而导致"教育"严重的内生性问题(残差项与教育相关)。

Heckman表明Heckman-style多方程模型能够解决“非观测到的混淆因素”所带来的问题。不过，后面的研究表明，他所依赖的方程形式和相关假设会让其估计值出现很大偏差。比如，在学者经常使用的bivariate probit模型中，尤其是遇到所谓的“样本选择偏差时”，通常假定选择方程和结果方程的变量都服从边际正态概率分布，且两个方程中的残差项e1和e2服从二元正态分布。此时由二个方程所组成的联立方程组看起来就比较僵化呆板，但即将引荐的“弹性联合似然模型”(flexible joint likelihood model)就能通过一个Copula函数放松对方程形式的各种假设。

至于具体怎么放松各种假设形式，“弹性联合似然模型”主要表现在三点，①允许联立方程组中的各个方程里的变量服从不一样的边际概率分布，②允许“dependence参数”(联立方程中的残差项之间相关关系)服从包括但不限于Gaussian正态分布，③允许各个方程里出现半参数回归形式，即y是x1和f(x2)的函数。我们知道，之前的文献中一般都假定两个方程中的残差项服从二元正态分布，此时，若联立方程中残差项的相关系数θ≠0，那就表明需要通过联立方程去估计各个方程中的系数。然而，Copula函数(包括Clayton, Frank, Gumbel and Joe和 Gumbel and Joe等)能够把两个残差项的关系扩展到很多种不同形式。

下面可以看看这个关于Copula函数各种形式的解释"Copula函数"：Contour plots of some classic copula functions with standard normal margins for data simulated using association parameters 2, 5.74, 0.71, 2, 2.86, and 0.71, respectively (these values are consistent with a medium positive correlation). The Gaussian, Student-t (here with three degrees of freedom) and Frank copulas allow for equal degrees of positive and negative dependence. Gaussian and Frank show a weaker tail dependence as compared to Student-t, and Frank exhibits a slightly stronger dependence in the middle of the distribution. Clayton is asymmetric with a strong lower tail dependence but a weaker upper tail dependence. Vice versa for the Gumbel and Joe copulas.

下面，我们就从一个实际操作中理解一下，“弹性联合似然模型”是如何处理“不可观测的混淆变量”所带来的内生性问题(注: 下方操作是由R软件完成，Stata软件中没有对应的软件包)。这下面只展示了一部分的运行程序，完整的运行程序和软件包放在社群和因果推断研究小组(白天上传)。

以后，若遇到可能存在的不可观测的混淆变量时，那么推荐使用以下方法去解决其带来的内生性问题。

内生性问题

例一：0-1内生性变量和0-1结果变量的情况

set.seed(0)
n <- 400
Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u <- rMVN(n, rep(0,2), Sigma)
cov <- rMVN(n, rep(0,2), Sigma)
cov <- pnorm(cov)
x1 <- round(cov[,1]); x2 <- cov[,2]
f1 <- function(x) cos(pi2x) + sin(pix)
f2 <- function(x) x+exp(-30(x-0.5)^2)
y1 <- ifelse(-1.55 + 2x1 + f1(x2) + u[,1] > 0, 1, 0)
y2 <- ifelse(-0.25 - 1.25y1 + f2(x2) + u[,2] > 0, 1, 0)
dataSim <- data.frame(y1, y2, x1, x2)
检验y1是否是内生变量, 如果值小于0.05则具有内生性
LM.bpm(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), dataSim, Model = "B")

1.标准的RECURSIVE BIVARIATE PROBIT, y1是0-1内生变量, y2是0-1结果变量(结果变量和处理变量都是0-1虚拟变量)

out <- gjrm(list(y1 ~ x1 + x2,y2 ~ y1 + x2),data = dataSim, margins = c("probit", "probit"), Model = "B")  ##B=bivariate mode，连接函数是probit
conv.check(out)   ##收敛性检测，无关紧要
summary(out)  ##显示回归结果, theta是dependence参数,Kendall's tau系数AIC(out); BIC(out)   ##信息准则指标
AT(out, nm.end = "y1", hd.plot = TRUE) ##平均处理效应ATE

关于Kendall's tau系数可以看看下面的解释

这里，我们关注的是内生政策变量y1对结果变量y2的影响，因此，需要把焦点放在第二个方程中的y1系数和显著性。从回归结果来看，y1给y2带去了负面的政策效应(三颗星显著)。得到的AIC和BIC分别为870.0191和897.9593。

得到的平均处理效应为：Treatment effect (%) with 95% interval: -44.1 (-54.4,-29.5)

2.弹性的RECURSIVE BIVARIATE PROBIT, 半参数回归, 其中x2和x3替换成他们各自的函数形式

out <- gjrm(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), data = dataSim, margins = c("probit", "probit"), Model = "B")
conv.check(out)
summary(out)
AIC(out); BIC(out)
估计后再做y1是否是内生性变量的检验, 小于0.05则具有内生性
gt.bpm(out)
处理效应ATE, 风险比例RR和odds比例
mb(y1, y2, Model = "B")
AT(out, nm.end = "y1", hd.plot = TRUE)
RR(out, nm.end = "y1")
OR(out, nm.end = "y1")
AT(out, nm.end = "y1", type = "univariate")

这里，我们关注的是内生政策变量y1对结果变量y2的影响，因此，需要把焦点放在第二个方程中的y1系数和显著性。从回归结果来看，y1给y2带去了负面的政策效应(三颗星显著)。得到的AIC和BIC分别为788.7169和858.1936。根据信息准则，此种形式的回归是最优的。

得到的平均处理效应为：Treatment effect (%) with 95% interval: -49.9 (-58.8,-41.7)。

3.弹性的RECURSIVE BIVARIATE PROBIT，dependence参数服从Clayton copula分布

outC <- gjrm(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), data = dataSim, BivD = "C0", margins = c("probit", "probit"), Model = "B") ##BivD=bivariate error服从clayton copula分布
conv.check(outC)
summary(outC)
AT(outC, nm.end = "y1")

这里的结果部分就省略了，只展示平均处理效应ATE。得到的AIC和BIC分别为870.0191和897.9593。

得到的平均处理效应为：Treatment effect (%) with 95% interval:-50.4 (-61.5,-38.0)

4.弹性的RECURSIVE BIVARIATE PROBIT，dependence参数服从Joe copula分布

outJ <- gjrm(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), data = dataSim, BivD = "J0", margins = c("probit", "probit"),Model = "B") ##BivD=bivariate error服从Joe copula分布
conv.check(outJ)
summary(outJ)
AT(outJ, "y1")
VuongClarke(out, outJ) ##Vuong's test和Clarke's test

这里的结果部分就省略了，只展示平均处理效应ATE。得到的AIC和BIC分别为870.0191和897.9593。

得到的平均处理效应：Treatment effect (%) with 95% interval:-47.1 (-59.7,-35.5)

例二：0-1内生变量y1和连续结果变量y2的情形

set.seed(0)
n <- 1000
Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u <- rMVN(n, rep(0,2), Sigma)
cov <- rMVN(n, rep(0,2), Sigma)
cov <- pnorm(cov)
x1 <- round(cov[,1]); x2 <- cov[,2]
f1 <- function(x) cos(pi2x) + sin(pix)
f2 <- function(x) x+exp(-30(x-0.5)^2)
y1 <- ifelse(-1.55 + 2x1 + f1(x2) + u[,1] > 0, 1, 0)
y2 <- -0.25 - 1.25y1 + f2(x2) + u[,2]
dataSim <- data.frame(y1, y2, x1, x2)

1.标准的RECURSIVE Bivariate Model(这里不是bivariate probit)

rc <- resp.check(y2, margin = "N", print.par = TRUE, loglik = TRUE) ##产生因变量y2的柱状图和QQ图
out <- gjrm(list(y1 ~ x1 + x2, y2 ~ y1 + x2), data = dataSim, margins = c("probit","N"), Model = "B") ##B=bivariate mode，y1为0-1变量因此连接函数为probit，y2为连续变量因此连接函数为identity
conv.check(out)
summary(out)
post.check(out)Treatment effect with 95% interval: -1.211 (-1.457,-0.998)

这里，我们关注的是内生政策变量y1对结果变量y2的影响，因此，需要把焦点放在第二个方程中的y1系数和显著性。从回归结果来看，y1给y2带去了负面的政策效应(三颗星显著)。得到的AIC和BIC分别为3914.293和3953.555。

得到的平均处理效应：Treatment effect with 95% interval: -1.211 (-1.457,-0.998)

2.弹性的RECURSIVE Bivariate Model, 半参数回归，其中x2和x3替换成他们各自的函数形式

eq.mu.1 <- y1 ~ x1 + s(x2)
eq.mu.2 <- y2 ~ y1 + s(x2)
eq.sigma2 <- ~ 1 ##sigma2是dispersion参数(方差)
eq.theta <- ~ 1 ##theta是dependence参数
fl <- list(eq.mu.1, eq.mu.2, eq.sigma2, eq.theta)
out <- gjrm(fl, data = dataSim, margins = c("probit","N"), gamlssfit = TRUE,
Model = "B")
conv.check(out)
summary(out)
post.check(out)
jc.probs(out, 1, 1.5, intervals = TRUE)[1:4,]
AT(out, nm.end = "y1")
AT(out, nm.end = "y1", type = "univariate")

这里的结果部分就省略了，只展示平均处理效应ATE。得到的AIC和BIC分别为3692.566和3766.468。根据信息准则，此种形式的回归是最优的。

得到的平均处理效应：Treatment effect with 95% interval: -1.23 (-1.41,-1.04)。

例三：内生连续变量y1, 0-1虚拟结果变量的情形

set.seed(0)
n <- 1000
Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u <- rMVN(n, rep(0,2), Sigma)
cov <- rMVN(n, rep(0,2), Sigma)
cov <- pnorm(cov)
x1 <- round(cov[,1]); x2 <- cov[,2]
f1 <- function(x) cos(pi2x) + sin(pix)
f2 <- function(x) x+exp(-30(x-0.5)^2)
y1 <- -0.25 - 2x1 + f2(x2) + u[,2]
y2 <- ifelse(-0.25 - 0.25y1 + f1(x2) + u[,1] > 0, 1, 0)
dataSim <- data.frame(y1, y2, x1, x2)
eq.mu.1 <- y2 ~ y1 + s(x2)
eq.mu.2 <- y1 ~ x1 + s(x2)
eq.sigma2 <- ~ 1
eq.theta <- ~ 1
fl <- list(eq.mu.1, eq.mu.2, eq.sigma2, eq.theta)
out <- gjrm(fl, data = dataSim, margins = c("probit","N"), Model = "B") ##注意y1和y2的方程顺序，这里sigma2和theta参数都为常数(其实你可以去掉eq.sigma2和eq.theta).
conv.check(out)
summary(out)
post.check(out)
AT(out, nm.end = "y1") ##average treatment effect
AT(out, nm.end = "y1", type = "univariate")
RR(out, nm.end = "y1", rr.plot = TRUE) ##risk ratios
RR(out, nm.end = "y1", type = "univariate")
OR(out, nm.end = "y1", or.plot = TRUE) ##odds ratio
OR(out, nm.end = "y1", type = "univariate")

这里，我们关注的是内生政策变量y1对结果变量y2的影响，因此，需要把焦点放在第一个方程中的y1系数和显著性。从回归结果来看，平均而言，y1给y2带去了负面的政策效应(三颗星显著)。

因为内生政策变量y1是连续变量，因此得到的平均处理效应如下表和图所示，在y1=2之前，y1每增加一个单位则政策效应变得更强了(绝对值)，而在y1=2之后，则y1每增加一个单位则政策效应变得更弱了。但在这个区间里，有一个共同点是，政策效应始终是负的。

下面这部分属于衍生，不属于这篇文章的主题“不可观测混淆变量”的内生性处理主题。

鉴于“弹性的联合似然模型”的灵活性，该框架在样本选择问题方面也具有很大的优势。正如上面所说，他允许dependce参数是协变量的函数，还允许半参数回归来拟合选择和结果变量等。

样本选择问题

例一：样本选择偏误

set.seed(0)
n <- 2000
rh <- 0.5
sigmau <- matrix(c(1, rh, rh, 1), 2, 2)
u <- rMVN(n, rep(0,2), sigmau)
sigmac <- matrix(rh, 3, 3); diag(sigmac) <- 1
cov <- rMVN(n, rep(0,3), sigmac)
cov <- pnorm(cov)
bi <- round(cov[,1]); x1 <- cov[,2]; x2 <- cov[,3]
f11 <- function(x) -0.7(4x + 2.5x^2 + 0.7sin(5x) + cos(7.5x))
f12 <- function(x) -0.4( -0.3 - 1.6x + sin(5x))
f21 <- function(x) 0.6(exp(x) + sin(2.9x))
ys <- 0.58 + 2.5bi + f11(x1) + f12(x2) + u[, 1] > 0
y <- -0.68 - 1.5bi + f21(x1) + u[, 2]
yo <- y(ys > 0)
dataSim <- data.frame(ys, yo, bi, x1, x2)

看一下这个里面的数据结构，选择变量ys就是0-1结构，对于ys=0的数据，结果变量yo=0，因此，我们倾向于使用两个方程来估计bi对yo的影响。

1.标准的样本选择偏差方程，第一个方程ys必须是选择方程

resp.check(yo[ys > 0], "N")
out <- gjrm(list(ys ~ bi + x1 + x2, yo ~ bi + x1), data = dataSim, Model = "BSS",margins = c("probit", "N")) ##BSS=bivariate model with sample selection
conv.check(out)
post.check(out)
summary(out)
AIC(out); BIC(out)

这里，我们关注的是bi对yo的影响，因此，需要把焦点放在第二个方程中的bi系数和显著性。从回归结果来看，bi对yo有负面的效应(三颗星显著)。得到的AIC和BIC分别为4562.817和4613.225。

2.弹性的样本选择偏差方程，半参数样本选择偏误

out <- gjrm(list(ys ~ bi + s(x1) + s(x2), yo ~ bi + s(x1)), data = dataSim, Model = "BSS", margins = c("probit", "N"))
conv.check(out)
post.check(out)
AIC(out); BIC(out)

summary(out) ##考虑了选择偏差
summary(out$gam2) ##没有考虑选择偏差

考虑了样本选择偏差后所得到的估计结果：

下面是没有考虑样本选择偏差时所得到的估计结果。通过与上面考虑了样本选择偏差的情形相比，没有考虑样本选择偏差会导致overestimation。

3.弹性的样本选择偏差，半参数样本选择偏差：相关和离差参数是协变量的函数

eq.mu.1 <- ys ~ bi + s(x1) + s(x2)
eq.mu.2 <- yo ~ bi + s(x1)
eq.sigma2 <- ~ bi
eq.theta <- ~ bi + x1
fl <- list(eq.mu.1, eq.mu.2, eq.sigma2, eq.theta)
out <- gjrm(fl, data = dataSim, Model = "BSS",margins = c("probit", "N"))
conv.check(out)
post.check(out)
summary(out)
out$sigma2
out$theta
jc.probs(out, 0, 0.3, intervals = TRUE)[1:4,]
outC0 <- gjrm(fl, data = dataSim, BivD = "C0", Model = "BSS", margins = c("probit", "N")) ##bivariate error服从Clayton copula分布
conv.check(outC0)
post.check(outC0)
AIC(out, outC0); BIC(out, outC0)

通过比较信息准则指标AIC和BIC，我们发现outc0这个联立方程系统不如out表现得好。

4. 二次样本选择(double sample selection)，总共有三个方程，其中二个是0-1选择方程

set.seed(0)
n <- 5000
Sigma <- matrix(c(1, 0.5, 0.4,
0.5, 1, 0.6,
0.4, 0.6, 1 ), 3, 3)
u <- rMVN(n, rep(0,3), Sigma)
f1 <- function(x) cos(pi2x) + sin(pix)
f2 <- function(x) x+exp(-30(x-0.5)^2)
x1 <- runif(n)
x2 <- runif(n)
x3 <- runif(n)
x4 <- runif(n)
y1 <- 1 + 1.5x1 - x2 + 0.8x3 - f1(x4) + u[, 1] > 0
y2 <- 1 - 2.5x1 + 1.2x2 + x3 + u[, 2] > 0
y3 <- 1.58 + 1.5*x1 - f2(x2) + u[, 3] > 0
dataSim <- data.frame(y1, y2, y3, x1, x2, x3, x4)
f.l <- list(y1 ~ x1 + x2 + x3 + s(x4),y2 ~ x1 + x2 + x3,y3 ~ x1 + s(x2))
out <- gjrm(f.l, data = dataSim, Model = "TSS", margins = c("probit", "probit", "probit"))
conv.check(out)
summary(out)
plot(out, eq = 1)
plot(out, eq = 3)
prev(out)
prev(out, type = "univariate")
prev(out, type = "naive")