温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
周洋 (清华大学, yangzhou0327@outlook.com)
李森林 (中南财经政法大学, 2289286835@qq.com)
连玉君 (中山大学, arlionn@163.com)
目录
在实证分析中我们常常假设解释变量和被解释变量存在线性关系。然而,在很多情况下,解释变量和被解释变量可能存在非线性关系。在模型中加入平方项,甚至是高阶项 (断点回归分析 RDD 中常用的处理方法),是文献中常用的模型设定方法。典型的模型设定方式为:
其中,
本文的目的在于分析加入平方项后我们该如何解释系数的含义?有哪些需要特别注意的地方?
浏览现有的期刊论文,我们经常发现不少作者只要发现平方项的系数
在社会科学研究中,我们必须考虑变量的经济含义和取值范围,这就使得我们不能像做初中代数习题那样来分析包含平方项的模型。比如,研究「收入 (Income) 与年龄 (Age) 」的非线性关系时,Age 就不可能取负值,甚至不可能小于 18。同时,Age 的最大值也有限制,取决于研究对象及样本特征。
一旦考虑到
图 1 呈现了一些典型的状况,我们可以可以通过在线性回归模型中添加二次项来拟合每一种情形。在 Stata 中,使用 regress
命令即可为所有以上曲线提供了线性和二次系数的最佳估计。
下面,举一些例子来具体说明。
图 1-A 呈现了经济学中反复强调的边际效用递减规律。这种例子很多。比如,随着收入的增加,食物消费支出也会增加,但增幅会不断降低。又如,喝咖啡可以提神,但第五杯的效果远远低于第一杯。「倒 U 形关系」的典型特征是模型 (1) 中的平方项系数
其中,
「U 形关系」或「倒 U 形关系」也很普遍,此时,典型的特征是,
例如,图 1-B 呈现了夫妻冲突程度与夫妻快乐程度之间可能存在的「倒 U 形关系」关系。夫妻经常有冲突吗? 如果他们一直有冲突,他们会不开心。 但是一种理论认为,从来没有或很少有冲突的夫妻也是不快乐的,那些有适度冲突的人是最快乐的。注意正的线性关系系数正在将
图 1-C 中呈现了一种「衰减递减型曲线」关系。假设被解释变量是各国的婴儿死亡率 (
此时,平方项的系数
类似的例子还有很多:止疼药可以减缓疼痛,但随着用药剂量的增加,药效逐渐下降 (产生了耐药性);经济援助可以降低一国的失业率,但政策效果可能会不断降低,等等。
图 1-D 和 图 1-E 分别呈现了另外两种情形,有兴趣的读者可以自行举例分析。
图 1: 加入平方项可能刻画的非线性关系
注意: 在加入二次项系数后,我们需要关注转折点:
,同时,还要作如下几个分析才能确保转折点是有意义的:
在最简单的情形中,
使用 OLS 估计出参数后,可以将
我们关注的是
这说明,
当
在(2)式中,如若估计出的
加入二次项的方程也可以具有
如果水平项和二次项的系数具有相同的符号(同为正或者同为负),而解释变量必须非负,那么,便不存在
例如,若
对于任意二次方程的转折点的一般公式为
因此,在计算转折点时,考虑
加入二次项之后,一般很难避免多重共线性。由于多重共线性具有方差膨胀(variance inflation)的作用,故加入二次项后一般会使得估计量的方差增大,导致回归系数的显著性下降。这当然不是我们想看到的效果。
那么,究竟是否应该在回归模型加入二次项呢?这就涉及到如何在遗漏变量偏差与多重共线性之间进行权衡。当然,如果线性模型已经是对于现实世界的足够好近似,那么就可以忽略遗漏变量偏差,而不必加入二次项或高次项了。为此,可以进行「回归方程设定误差检验」(Ramsey's RESET 检验,即 Regression Equation Specification Error Test)。比如,如果完整的方程为:
则可对原假设
进行 F 检验。如果拒绝
如果在模型中加入二次项,则一般应在论文中同时汇报仅包含一次项的简洁模型,以及包含二次项的完整模型之估计结果,这是所谓「稳健性检验」(robustness checks) 的一种形式。如果两种模型的定性结果类似 (qualitatively similar),或者不影响你感兴趣变量的显著性与符号,则也很容易处理。
困难之处在于,有时简洁模型与完整模型的结果并不一致,甚至影响了统计显著性或回归系数的符号。而产生这种现象的原因依然是遗漏变量偏差或多重共线性。如果存在遗漏变量偏差,则简单的线性模型并不一致,而包含二次项的完整模型才是一致估计,故二者的估计结果大相径庭,也在情理之中。此时,应该选择后者。
另一方面,即使遗漏变量偏差不存在或较微弱,加入二次项所导致的多重共线性,也可能通过「方差膨胀因子」(variance inflation factor),增大估计量的标准误,使得原来显著的项变得不再显著。(译者注:对于多数微观面板数据而言,如上市公司数据、劳动力调查数据等,共线性通常都不会构成太严重的威胁)。
下面我们使用 Stata 自带的 nlsw88.dta 数据为例,介绍回归分析中二次项在 Stata 中的具体用法。
其中, wage
表示个体的工资水平, ttl
表示个体的工作经验,而 ttl
的二次项。
在正式分析之前,我们可以使用如下命令绘制 ttl 与 wage 之间的散点图,以便初步确认二者之间是否存在非线性关系。
. sysuse "nlsw88.dta", clear
. twoway (scatter wage ttl) (qfit wage ttl, lcolor(red))
其中,twoway () ()
是 Stata 绘制二维图的标注架构,两对括号中分别绘制两个图形,分处于两个图层中 (每个图层可以视为一张透明的胶片)。qfit
用于绘制两个变量之间的非线性拟合曲线。 输出的图形如下:
图 2: 散点图及非线性拟合图
可以看出,工资 wage 与工作经验 ttl_exp 可能存在倒 U 型关系。然而,我们可能更为关注转折点的具体数值,下面进行计算:
. sysuse "nlsw88.dta", clear
. gen ttl_exp2 = ttl_exp*ttl_exp
. reg wage ttl_exp ttl_exp2 // 也可以添加一些控制变量
Source | SS df MS Number of obs = 2,246
-------------+---------------------------------- F(2, 2243) = 87.46
Model | 5379.78611 2 2689.89305 Prob > F = 0.0000
Residual | 68988.1813 2,243 30.7571027 R-squared = 0.0723
-------------+---------------------------------- Adj R-squared = 0.0715
Total | 74367.9674 2,245 33.1260434 Root MSE = 5.5459
------------------------------------------------------------------------------
wage | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
ttl_exp | 0.557 0.109 5.10 0.000 0.343 0.771
ttl_exp2 | -0.009 0.004 -2.12 0.034 -0.018 -0.001
_cons | 2.462 0.639 3.85 0.000 1.209 3.716
------------------------------------------------------------------------------
. global tp: dis %3.1f -_b[ttl_exp]/(2*_b[ttl_exp2]) //转折点
. dis "Turn point = " $tp //能够得到 Turn point = 29.6
*---------------------------图示-----------------------------
. local b0 = _b[_cons]
. local b1 = _b[ttl_exp]
. local b2 = _b[ttl_exp2]
. sum ttl_exp
. local min: dis %3.1f r(min)
. local max: dis %3.1f r(max)
. #delimit
. twoway ( function y = `b2'*x^2 + `b1'*x + `b0',
range(`min' `max') )
( function y = `b2'*x^2 + `b1'*x + `b0',
range(`max' 40) lp(dash) )
,
ytitle("wage") xtitle("ttl_exp")
xline(`min',lc(red) lw(*1.5))
xline($tp, lp(dash) lc(green))
xline(`max',lc(red) lw(*1.5))
xlabel(`min' "Min (`min')" 5(5)25
`max' "Max (`max')", angle(50))
legend(off)
;
. #delimit cr
. graph export "fig_Ushape.png", replace //保存图形
下图绿色虚线所在位置即为转折点(turn point,TP),等于 29.6:
图 3: 图示转折点
为了便于查看,我们还在图中用红色实线标注了解释变量 ttl_exp 的最小值和最大值。显然,本例中计算所得到的转折点在 ttl_exp 的最大值的右侧。也就是说,虽然表面上存在倒 U 型关系 (ttl_exp2) 的系数显著为负,但在本例的样本区间并不能覆盖转折点及其右侧的观察值。这意味,虽然 ttl_exp 与 wage 之间存在非线性关系,但这种关系主要表现为 ttl_exp 对 wage 的边际影响逐渐递减,但二者的关系始终是单调增的。简言之,此时只能讨论「倒 U 型关系」的左半支,而转折点本身并不具有实际的经济意义。
我们也可以使用 count
命令统计一下转折点两侧的观察值的分布情况:
. count if ttl_exp>$tp & e(sample)==1
0
在上述模型表示,工作经验 (ttl_exp) 及其二次项 (ttl_exp2) 会影响工资水平,我们要计算的是工作经验对工资的边际效应。
我们先用 regress
命令 (需要使用因子变量形式来表示二次项,help fvvarlist
,详情参阅「Stata:因子变量的全攻略」 估计模型,进而使用 margins
命令计算边际效应,并以此为基础使用 marginsplot
图示之。
sysuse "nlsw88.dta", clear
reg wage c.ttl_exp##c.ttl_exp $xx
sum ttl_exp
global min: dis %4.1f r(min)
global max: dis %4.1f r(max)
margins, dydx(ttl_exp) at(ttl_exp=($min 1(3)28 $max))
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
ttl_exp |
_at |
1 | .5551562 .1084481 5.12 0.000 .3424871 .7678252
2 | .5382228 .1007071 5.34 0.000 .3407339 .7357117
3 | .4817783 .0752654 6.40 0.000 .3341811 .6293755
4 | .4253338 .0510196 8.34 0.000 .3252833 .5253843
5 | .3688893 .0309232 11.93 0.000 .3082482 .4295304
6 | .3124448 .0269191 11.61 0.000 .2596558 .3652338
7 | .2560003 .0436827 5.86 0.000 .1703376 .3416629
8 | .1995557 .0671331 2.97 0.003 .0679062 .3312053
9 | .1431112 .0923075 1.55 0.121 -.0379059 .3241283
10 | .0866667 .1181086 0.73 0.463 -.1449469 .3182803
11 | .0302222 .1442004 0.21 0.834 -.252558 .3130024
12 | .0132889 .1520623 0.09 0.930 -.2849086 .3114864
------------------------------------------------------------------------------
下面进行绘图,附带 95% 置信区间。
marginsplot, xlabel(, format(%3.1f) angle(60)) //绘图,附带 95% 置信区间
graph export "$Out\fig_Wage_marginplot_A.png", replace
marginsplot, addplot(function y=0, range($min $max) ///
lcolor(red) lpattern(dash) legend(off)) ///
xlabel(, format(%3.1f) angle(60)) ///
ylabel(, format(%2.1f) angle(0))
graph export "$Out\fig_Wage_marginplot_B.png", replace
图 4: 边际效应图示
图 4 采用另一种方式展示了
在 图 3 中,我们展现的是
转而对比观察 图 4,可以清晰地看出上述变化过程。由于本例中,转折点 (TP = 29.6)
与样本中
在模型中加入变量的二次项是寻找 U 形关系的典型处理方法之一,然而,文献中也经常见到一些其他的处理方法,比如,分组回归、面板门槛模型等。他们之间有何关系?能否相互替代呢?
在模型进行异质性分析时,通常的做法是按照核心解释变量 (
文献中常用的分组标准包括:样本均值,中位数,第 33 和 66 百分位数 (将样本分成三组)。当然,在某些情况下,可以根据理论分析或数据的具体情况,选择更有意义的分组界点 (cut-point)。比如,研究退休对健康的影响时,往往会以 55 岁 或 60 岁 作为分组依据。
对于 图 5 中的散点图而言,我们可以分别以
采用二次项的方式有一些独特的好处:
margins
和 marginsplot
命令分析边际效应的变化。相比之下,分组回归只能看到两个组内的平均边际效应 (分组回归假设每个组内的边际效应均为常数)。当然,使用二次项时也有一些局限,主要在于
图 5a:分组回归的拟合图
图 5b:分组回归与加入二次项
为了让读者对上述分析有一点直观的感受,我们可以延续 nlsw88.dta 的例子,来对比一下分组回归和加入平方项两种设定下的估计结果有何差异。
在分组回归中,出于稳健性考虑,我们分别以 ttl_exp 的样本均值 (分组 1) 和第 75 百分位数 (分组 2) 为分组界点 ,将样本分成两组。从 图 3 的结果来看,「分组 2」的结果应该更接近于加入平方项的结果。 下表中的会结果证实按上述猜测,具体分析留给各位读者吧。
*-分组回归
. sysuse "nlsw88.dta", clear
* 根据 ttl_exp 的样本均值分组
qui sum ttl_exp, detail
gen G1 = (ttl_exp>r(mean)) // r(mean) 表示 ttl_exp 的样本均值
*gen G2 = (ttl_exp>r(p75)) // r(mean) 表示 ttl_exp 的样本均值
global cx "age hours i.race i.industry"
reg wage ttl_exp $cx if G1==0 // Group 1
est store Lmean
reg wage ttl_exp $cx if G1==1 // Group 2
est store Rmean
* 根据 ttl_exp 的第 75 百分位数分组 (p75)
qui sum ttl_exp, detail
gen G2 = (ttl_exp>r(p75)) // r(75) 表示 ttl_exp 的 p75
reg wage ttl_exp $cx if G2==0 // Group 1
est store Lp75
reg wage ttl_exp $cx if G2==1 // Group 2
est store Rp75
* 对比:加入二次项
reg wage c.ttl_exp##c.ttl_exp $cx
est store x2
*-----表:回归结果-------
local m "Lmean Rmean Lp75 Rp75 x2"
esttab `m', mtitle(`m') nogap compress replace ///
b(%6.3f) s(N r2_a) drop(`drop') ///
star(* 0.1 ** 0.05 *** 0.01) ///
addnotes("*** 1% ** 5% * 10%") ///
indicate("种族效应 =*.race" "行业效应 =*.industry" )
------------------------------------------------------------------------------
分组 1 分组 2 平方项
----------------------- ------------------------- ------------
(1) (2) (3) (4) (5)
Lmean Rmean Lp75 Rp75 x2
------------------------------------- ------------------------- ------------
ttl_exp 0.295*** 0.152** 0.348*** 0.118 0.423***
(5.17) (2.08) (8.89) (1.00) (3.85)
age -0.063 -0.133** -0.107** -0.042 -0.107***
(-1.21) (-2.34) (-2.36) (-0.57) (-2.77)
hours 0.048*** 0.055*** 0.052*** 0.044* 0.051***
(3.37) (3.00) (3.98) (1.77) (4.47)
ttl_exp^2 -0.006
(-1.45)
_cons 3.260 8.520*** 4.427* 6.449 4.600**
(1.20) (2.67) (1.80) (1.60) (2.11)
种族效应 Yes Yes Yes Yes Yes
行业效应 Yes Yes Yes Yes Yes
------------------------------------------------------------------------------
N 1019.000 1209.000 1671.000 557.000 2228.000
r2_a 0.098 0.065 0.112 0.096 0.123
------------------------------------------------------------------------------
t statistics in parentheses
*** 1% ** 5% * 10%
* p<0.1, ** p<0.05, *** p<0.01
上文介绍的「分组回归」方法存在一个明显的局限:我们需要事先人为设定分组界点。这显然是最容易被诟病的地方。Hansen (1999, JoE) 提出的「面板门槛模型」可以克服这一问题。其基本思想是:假设分组界点可以发生在
由此,可以获得当
上面介绍的是「单一门槛模型」的估计思路,若模型中有多个门槛,我们可以先按上述方法搜索第一个门槛,进而固定之 (相当于在模型中加入一个已知界点的交乘项),继而搜索第二个、第三个门槛值。图 6 展示了一个典型的「双重门槛模型」的估计结果。
由此看来,所谓的「门槛模型」本质上就是一个「分组界点未知时的分组回归模型」。它本着「让数据说话」的基本原则,预先通过模型的拟合情况 (RSS 最小) 来寻找门槛值,进而以门槛值为基础构造寻变量及其与
当然,具体应用过程中,还需要基于 Bootstrap 执行「门槛效应检验」,以及「门槛值是否等于真实值」等检验。详情可以查阅 Hansen (1999) 原文,或用以实现该模型的 Stata 命令的帮助文件:help xthreg
。
图 6:面板门槛模型图示
margins
和 marginsplot
命令进行边际效应及其图示分析;
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟 Stata 33 讲 - 连玉君, 每讲 15 分钟. 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看,所有课程可以随时购买观看。
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 ⭐ | DSGE, 因果推断, 空间计量等 | |
⭕ Stata数据清洗 | 游万海 | 直播, 2 小时,已上线 |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD