加入交乘项后符号变了!?

2019-04-09 本文已影响4625人 stata连享会

作者：连玉君 (知乎 | 简书 | 码云)

Stata连享会精品专题 || 精彩推文

1. 简介

在前面的几篇推文中，我们对交乘项的基本设定、图示、边际效应分析等内容进行了较为细致的分析。最近适逢很多学生写毕业论文，有关交乘项的问题又涌上心头。其中，最突出的问题便是：为何加入交乘项后主变量变得不显著了，甚至符号都变掉了？

简单的解释是：此一时，彼一时！

因为，加入交乘项前后，主变量的系数含义发生了实质性的变化，二者不具可比性。本文的目的在于澄清这种差异，并介绍一种让主变量系数在加入交乘项前后不会发生大幅变化 (具有可比性) 的方法。

连享会-空间计量现场班 (2019年6月27-30日，西安)

2. 为何加入交乘项后主变量符号会变化？

对于模型

$y = \beta_0 + β_1x_1 + β_2x_2 + \epsilon_1 \qquad (1)$

系数 $\beta_1=\partial{y}/\partial{x_1} | x_2=\bar{x}_2$ ，也就是当 $x_2$ 取样本均值时， $x_1$ 变动一个单位对 $y$ 的影响。

当我们加入交乘项 $x_1 x_2$ 后， $x_1$ 的系数含义发生了很大的变化。

$y = \theta_0 + \theta_1x_1 + \theta_2x_2 + \theta_3 x_1 x_2 + \epsilon_2 \qquad (2)$

先看 $x_1$ 对 $y$ 的边际影响： $\partial{y}/\partial{x_1}=\theta_1+ \theta_3 x_2$ ，这是大家都了解的基本结论：包含交乘项时， $x_1$ 对 $y$ 的边际影响不再是常数，而是随着 $x_2$ 的取值不同而发生变化。

解释：此时，一阶项 $x_1$ 的系数为 $\theta_1 = \partial{y}/\partial{x_1} | x_2=0$ 。也就是说，在模型 (2) 中，一阶项 $x_1$ 的系数表示当 $x_2=0$ 时， $x_1$ 变动一个单位对 $y$ 的影响。显然，模型 (1) 中的 $\beta_1$ 与模型 (2) 中的 $\theta_1$ 估计值不同，甚至发生符号变化是很正常的事情。
举个简单的例子。假设 $y$ 表示收入； $x_1$ 表示丑陋程度； $x_2$ 表示教育年限，取值为 0, 1, 2, ……20，均值为 12。基本想法是想检验「教育能否扭转我在职场上的天生劣势？(你知道我为什么读 PhD 吗？)」假设估计模型 (2) 得到的参数为 $\theta_1=-1.6$ ， $\theta_3 = 0.2$ ，即 $\partial{y}/\partial{x_1}=-1.6 + 0.2x_2$ 。根据上面的数值，可以大致推断模型 (1) 中 $x_1$ 的系数约为 $\beta_1 \approx -1.6 + 0.2 \bar{x}_2 = -1.6 + 0.2\times 12 = 0.8$ 。大家可以自行分析一下 $\theta_1$ 和 $\beta_1$ 的经济含义。

3. 如何尽力保证模型间的系数可比性？

若想让加入交乘项前后的模型 (1) 和模型 (2) 中主变量 ( $x_1$ ) 的系数具有可比性，可以采用如下模型设定形式 (参见 Balli H O, Sørensen B E. Interaction effects in econometrics[J]. Empirical Economics, 2013, 45(1): 583-603. [PDF])：

$y = \gamma_0 + \gamma_1x_1 + \gamma_2x_2 + \gamma_3(x_1 − \bar{x}_1) (x_2 − \bar{x}_2) + \epsilon_3 \qquad (3)$

其中， $\bar{x}_1$ 和 $\bar{x}_2$ 分别表示 ${x}_1$ 和 ${x}_2$ 的样本均值。

此时，主变量 ( $x_1$ ) 的系数 $\gamma_1$ 会非常接近基于模型 (1) 得到的 $\beta_1$ ：

$\gamma_1=\frac{\partial{y}}{\partial{x_1}} \Big |_{x_2=\bar{x}_2} \qquad (3a)$

大家可能更加关心交乘项的系数是否会发生变化，答案是：不会！

因为，模型 (3) 相对于模型 (2) 无非是增加了一些一阶项和常数项，而交乘项并未发生变化。我们也可以用更为正式的方式来得到这一结论。对于模型 (2) 而言， $\partial[\partial{y}/\partial{x_1}]/\partial{x_2}=\theta_3$ ，而在模型 (3) 中 $\partial[\partial{y}/\partial{x_1}]/\partial{x_2}=\gamma_3$ 。

4. 部分离差还是全部离差？

在模型 (3) 中 $(x_1 − \bar{x}_1)$ 称为 $x_1$ 的离差形式，其实就是对 $x_1$ 的每个观察值都做去均值处理。

因此，文献中也会采用如下模型设定形式：

$y = \gamma_0 + \gamma_1(x_1 − \bar{x}_1) + \gamma_2 (x_2 − \bar{x}_2) + \gamma_3(x_1 − \bar{x}_1) (x_2 − \bar{x}_2) + \epsilon_3 \qquad (4)$

按照上面的分析逻辑不难看出，这个模型与模型 (3) 没有任何本质区别，因为展开后新增的项目 $\gamma_1\bar{x}_1$ 和 $\gamma_2\bar{x}_2$ 都是常数，即

$y = \delta_0 + \gamma_1x_1 + \gamma_2x_2 + \gamma_3(x_1 − \bar{x}_1) (x_2 − \bar{x}_2) + \epsilon_3 \qquad (4a)$

其中， $\delta_0 = (\gamma_0-\gamma_1\bar{x}_1-\gamma_2\bar{x}_2)$ 。简言之，相对于模型 (3)，由模型 (4) 中得到的 $x_1$ , $x_2$ 以及 $x_1x_2$ 的系数都是完全相同的，唯一差别在于常数项。

需要补充说明的是，无论是采用模型 (3) 还是模型 (4)，本意都是为了方便对系数的含义进行解释，并不是所谓的克服共线性之类的说辞。

5. 模拟分析

参考 Balli et al. (2013, [PDF]) 文中的做法进行模拟，发现在使用交乘项时，在模型中用 $(x_1 − \bar{x}_1) (x_2 − \bar{x}_2)$ 替换 $x_1x_2$ ，一次项的系数更容易解释一些。

Note： 这里的 $x$ 表示上文中的 $x_1$ ，这里的 $z$ 表示上文中的 $x_2$ 。

*-Source: 
/*
Balli, H. O., B. E. Sørensen, 2013, 
   Interaction effects in econometrics, 
   Empirical Economics, 45 (1): 583-603.
*/

/* Table 1
The true model is Y = 3X1 + 5X2 + 8X1X2 + e 
where X1 = 1 + e1 and X2 = 1 + e2, 
ei~N(0,1) for i = 1, 2 
(X1 and X2 are not correlated) and  e~N(0,100). 
A constant is included but not reported. 
The sample size is 500 and the number of simulations is 20,000. 
Averages of estimated t statistics are shown in parentheses
*/
clear 
set obs 500
set seed 135
local rhox = 0
gen x = 1 + rnormal()
gen z = 1 + rnormal() + `rhox'*x
gen e = rnormal(0,10)

gen y = 10 + 3*x + 5*z + 8*x*z + e

pwcorr y x z

center x z, prefix(c_)
*-模型 (0)
reg y x 
est store m0
*-模型 (1)
reg y x z     
est store m1
*-模型 (2)
reg y x z c.x#c.z
est store m2
*-模型 (3)
reg y x z c.c_x#c.c_z  // Balli2013, Eq.(3)
est store m3
*-模型 (4)
reg y c_x c_z c.c_x#c.c_z
est store m4

*-结果对比
  local m "m0 m1 m2 m3 m4"
  local m "m1 m2 m3 m4"
  esttab `m' `s', nogap replace order(x z c_x c_z)  ///
         b(%6.3f) s(N r2_a) drop(`drop')   ///
         star(* 0.1 ** 0.05 *** 0.01)      ///
         addnotes("*** 1% ** 5% * 10%")


----------------------------------------------------------------------------
Model                (1)             (2)             (3)             (4)   
----------------------------------------------------------------------------
x                   9.979***        2.904***       10.047***                
                  (17.66)          (4.48)         (21.68)                   
z                  12.898***        5.450***       13.101***                
                  (22.53)          (8.14)         (27.90)     
c_x                                                                10.047***
                                                                  (21.68)   
c_z                                                                13.101***
                                                                  (27.90)                 
x#z                                 7.479***                                
                                  (15.59)                                   
c_x#c_z                                             7.479***        7.479***
                                                  (15.59)         (15.59)   
_cons               3.024***        9.792***        2.485***       25.275***
                   (3.12)         (10.81)          (3.12)         (53.61)   
----------------------------------------------------------------------------
N                 500.000         500.000         500.000         500.000   
r2_a                0.629           0.751           0.751           0.751   
----------------------------------------------------------------------------
t statistics in parentheses
*** 1% ** 5% * 10%

.   sum y x z 
    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
           y |        500    25.55029    21.10401  -33.97629   103.5704
           x |        500       1.023    1.018471  -2.755543   3.808831
           z |        500    .9550672    1.005523  -1.766887   3.902355

结果分析：

在模型 (1)-(3) 中，交乘项的系数和 t 值都完全相同；
对比模型 (3) 和模型 (4)，除了常数项外，其他变量的系数和 t 值完全相同；
模型 (2) 和模型 (3) 中的一阶项系数之间的关系。在两个模型中 x 的系数分别为 2.904 和 10.047。z 的样本均值为 zm = 0.955。因此，对于模型 (2) 而言，当 z 取其样本时，dy/dx (z=zm) = 2.904 + 7.479*zm = 2.904 + 7.479*0.955 = 10.046。这与第三列中由模型 (3) 估计出的 x 的系数 (10.047) 非常接近。期间的微小差别主要源于四舍五入。
模型 (1) 和模型 (3) 中的一阶项系数之间的关系。这是我们最关心的事情。由上面的分析可知，虽然模型 (1) 和 (2) 中 x 的系数存在很大差异，但 (1) 和 (3) 中 x 的系数应该比较接近。我们看到的结果也的确如此。

事实上，在很多论文中，通常会先估计 y = a + b*x，而不是 y = a + b1*x + b2*z ，即本文的模型 (1)。如果 $corr(x, z) = 0$ ，这两个模型中得到的 x 的系数不会有明显差异，但如果 x 和 z 彼此相关，则简化模型 y = a + b*x 就会存在遗漏变量的问题，其系数是有偏估计。感兴趣的读者，可以把上述模拟分析代码中的 local rhox = 0 修改为 local rhox = 0.5 或 local rhox = -0.5 等数值，并在结果呈现部分也列示出 m0 的结果，看看系数估计值会发生哪些变化。

6. 结论和实操建议

若只关心交乘项的系数，则使用模型 (1)-(3) 中的任何一种形式，都不影响系数的统计推断。
若希望让加入交乘项前后一阶项的系数容易解释，也具有可比性，可以在论文中呈现模型 (1) 和模型 (3) 的结果。
如果 $x_2$ 是一个虚拟变量，则无需做中心化处理，在论文中呈现模型 (1) 和模型 (2) 的结果即可，此时基于模型 (2) 进行分析反而更为直观。

参考资料

Balli H O, Sørensen B E. Interaction effects in econometrics[J]. Empirical Economics, 2013, 45(1): 583-603. [PDF] 这篇讲的非常清楚，还做了 MC 分析。
Bun M J G, Harrison T D. OLS and IV estimation of regression models including endogenous interaction terms[J]. Econometric Reviews, 2018: 1-14. -PDF-
Ebbes, P., D. Papies, H. van Heerde. 2016, Dealing with endogeneity: A nontechnical guide for marketing researchers[C], Handbook of market research 1-37. 从 IV 讲起，很细致，进一步讨论了交乘项内生以及多个内生变量的情形。[PDF]
Papies, D., P. Ebbes, H. J. Van Heerde. 2017, Addressing endogeneity in marketing models[C], Advanced methods for modeling markets (Springer, 581-627. [PDF], [web]
交乘项 - 这一篇就够了，说得无比清楚！
交乘项 - 一系列文章，写的很清楚
Interpreting coefficients from interaction (Part 1)
Interpreting three-way interactions in R - 模拟分析

后续写作安排

离差形式的系数含义图示，边际效应图示
包含内生变量以及内生变量的交乘项的情形
包含两个以上内生变量
估计方法
- IV
- 2SLS
- 控制函数法 (Control Function Approach, see Ebbes, Papies, van Heerde, 2016, pp.28 讲的非常清楚)
- 潜工具变量法 (LIV, see Ebbes, Papies, van Heerde, 2016)

关于我们

【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。
公众号推文同步发布于 CSDN-Stata连享会、简书-Stata连享会和知乎-连玉君Stata专栏。可以在上述网站中搜索关键词Stata或Stata连享会后关注我们。
点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
Stata连享会精彩推文1 || 精彩推文2

联系我们

欢迎赐稿： 欢迎将您的文章或笔记投稿至Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。
招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
联系邮件： StataChina@163.com

往期精彩推文

Stata连享会推文列表
Stata连享会精品专题 || 精彩推文

连享会-Python爬虫与文本分析现场班 (2019年5月17-19日，太原)

欢迎加入Stata连享会(公众号: StataChina)