倒U型+RDD:利用断点回归检验 U 形关系

发布时间:2021-12-31 阅读 4526

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc, ihelp, rdbalance, gitee, installpkg

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

PDF下载 - 推文合集

作者:郑伊达 (中山大学)
邮箱zhengyida@mail2.sysu.edu.cn

备注:本推文翻译自以下论文,特此致谢!
Source: Simonsohn, 2018"Two lines: A valid alternative to the invalid testing of U-shaped relationships with quadratic regressions.". -Link-


目录


1. 检验方法介绍

1.1 应用背景

在实证分析中,有时我们能够观察到 x 与 y 呈现出 U 形关系,又或者是倒 U 形关系。例如,工资收入与年龄之间就呈现出倒 U 形关系,人们多数在年轻时收入水平较低,中年时期其收入水平达到最高值,随后又因退休而导致收入减少 (Deming, 2019)。这些结论都是我们通过可观测数据观察到的结果,那么这样的结果能反映事实吗?有多大程度能够反映事实?基于这样的问题,我们需要进一步检验两者之间的关系,由此出现了对 U 形关系的检验。

1.2 二次回归检验 U 形关系的不足之处

在实证过程中,当研究者认为 x 和 y 之间存在一个 U 形关系,那么使用二次回归会是一个方便又快捷的方法,来对数据进行拟合,从而判断它们是否存在 U 形关系。但是,这种方法存在一定的风险,如果真实的回归模型不是一个二次方程,最终结果可能会误导研究者。

举一个例子,我们先定义一个函数关系式:y=log(x)。很明显,x 和 y 之间不是 U 形关系,如果我们用二次回归来拟合它们就会出现错误的结果。接下来我们就用 Stata 做个实验。

preserve

clear 

set obs 500       // 设置观测数
set seed 552      // 设置种子数
gen x=runiform()  // 生成x值(均匀分布)
gen log_y=log(x)  // 生成y值;y=log(x)
gen er=1*invnormal(uniform())  // 生成残差项,N~(0, 1)
gen y=log_y+er    // 生成观测值
sort x

gen x_2=x^(2)  // x二次项
reg y x x_2    // 二次回归
predict y_hat  // 得到拟合值y_hat
twoway scatter y x || line log_y x || line y_hat x

restore

上面的代码构造了一系列的伪随机观测值,真实的函数关系为 y=log(x),其残差项服从正态分布,均值为 0,方差为 1。随后我们生成一个 x2,作二次回归,结果如下:

      Source |       SS           df       MS      Number of obs   =       500
-------------+----------------------------------   F(2, 497)       =    141.97
       Model |  317.282005         2  158.641002   Prob > F        =    0.0000
    Residual |   555.37353       497  1.11745177   R-squared       =    0.3636
-------------+----------------------------------   Adj R-squared   =    0.3610
       Total |  872.655534       499  1.74880869   Root MSE        =    1.0571

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |   7.181829   .6636716    10.82   0.000     5.877881    8.485777
         x_2 |  -4.659204   .6456834    -7.22   0.000    -5.927809   -3.390598
       _cons |  -3.082529   .1439533   -21.41   0.000    -3.365361   -2.799697
------------------------------------------------------------------------------

从二次回归的结果可以看到,二次项的回归系数为负数,而且显著。因此,结果带给我们的结论是,x 和 y 之间存在倒 U 形关系。从这个例子中,可以显示出用二次回归来检测 U 形关系存在的问题,其核心问题是我们错误的假设了函数形式。

如果想了解跟多关于二次回归所引发的问题,可以参见连享会往期推文 :平方项 = 倒U型 ?

1.3 检验方法的理念

在检验 U-shaped 关系之前,我们先对 U 形关系进行定义:x 存在一个中间值 xc,小于 xc 的 x 为低数值组,大于 xc 的 x 为高数值组;在低数值组的 f(x) 与高数值组的 f(x) 两者之间为异号。简单来说,在 U 形关系中,x 数值较低的部分其线段斜率是负的,x 数值较高的部分其斜率是正的。另外,在 U-shaped 的中还包含额外的特征,例如,对称 vs 不对称,连续 vs 不连续,有极值 vs 无极值,想要深入研究这些特征需要用额外的方法去检测,在本推文中并不涉及。

这里我们引入一个概念,线性回归所计算出来的斜率系数是平均斜率,无论 x 和 y 的真实函数形式是如何。因此将 (xxc) 和 (xxc) 作两组线性回归,其两个斜率系数如果是异号显著,可以判断 x 和 y 存在 U 形关系。利用两段线性回归来检测是否存在 U 形关系,其最大的好处是我们无需对回归模型进行假设。

为加深对平均斜率的理解,以下我们利用 Stata 展示一个小例子。

preserve

clear

range x 0 3 4  // value-x: 0, 1, 2, 3
gen y = x^(2)  // value-y: 0, 1, 4, 9
reg y x
predict y_hat

twoway scatter y x || line y_hat x

restore

在这个例子中,我们假设真实的函数形式为: y=x2 。通过简单的计算我们可以得知 x 从 0 到 3 其线段的平均斜率为 3 [(1+3+5)/3]。随后,我们利用 OLS 回归,得出的斜率系数也为 3。

      Source |       SS           df       MS      Number of obs   =         4
-------------+----------------------------------   F(1, 2)         =     22.50
       Model |          45         1          45   Prob > F        =    0.0417
    Residual |           4         2           2   R-squared       =    0.9184
-------------+----------------------------------   Adj R-squared   =    0.8776
       Total |          49         3  16.3333333   Root MSE        =    1.4142

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |          3   .6324555     4.74   0.042     .2787635    5.721237
       _cons |         -1   1.183216    -0.85   0.487    -6.090967    4.090967
------------------------------------------------------------------------------

2. 断点回归

在上面的描述中,提到要将数据分为两部分别进行线性回归,其断点为 xc。这里我们引入断点回归 (see e.g., Marsh & Cormier, 2001, p. 7),其通用的回归模型为:

Where:

  • xlow=xxc if xxc, and 0 otherwise
  • xhigh=xxc if xxc, and 0 otherwise
  • high=1 if xxc, and 0 otherwise
  • Z is the (optional) matrix with covariates, and
  • BZ its vextor of coefficients.

完成断点回归后,可以得到两条线性回归的拟合值,对于 (xxc) 组别的拟合值为:y^low=a^+b^xlow;对于 (xxc) 组别的拟合值为:y^high=a^+d^+c^xhigh

在进行断点回归时,需要设置断点。那么如何找到一个合适的断点 xc,使两段回归的斜率系数其显著水平最大化?这里需要用到 Robin Hood 算法 (see e.g., Simonsohn, 2018),来寻找最合适的断点。这种算法不用对函数形式 f(x) 进行假设,不用对 x 的分布进行假设,也不用对残差项的分布进行假设。其算法有三个核心概念:(i) 两段回归最终要求两个斜率系数其显著性最大化,为了能进一步提升整体的显著水平,我们的目标是提升显著性较弱的回归;(ii) 其方法是令它们变得陡峭一点;(iii) 实现方法是给予它们更多的观测值。总的来说,就是通过改变断点的位置,将观测值从显著水平较高的回归转移至显著性水平较低的回归中。

3. 实现过程

3.1 生成伪随机数

首先,我们需要生成一些伪随机数,来作为实践对象。在这里,我们定义一个真实的函数形式 f(x)

  • 若 x<0.5,其斜率为 1,即 f(x)=1
  • 若 0.5x0.7,则 f(x)=0.5,f(x)=0
  • 若 x>0.7,其斜率为 0.15,即 f(x)=0.15

随后,我们在 f(x) 中加入随机扰动项 eN(0,0.5σy),生成伪随机观测值。数据生成过程可以表示如下:

其中,xiU[0,1]eiN(0,0.5σy)

Stata 代码如下:

clear

set obs 500       // 设置观测数
set seed 552      // 设置种子数
gen x=runiform()  // 生成x值(均匀分布)

gen pre_y=x if x<0.5  // x小于0.5,斜率为1
replace pre_y=0.5*1-(x-0.7)*0.15 if x>0.7  // x大于0.7,斜率为0.15
replace pre_y=0.5 if 0.5<=x & x<=0.7       // x的中间部分,y=0.5

egen pre_y_sd=sd(pre_y) 
gen er=0.5*pre_y_sd*invnormal(uniform())  // 生成残差项,N~(0, 0.5*y.sd)
gen y=pre_y+er

//* 绘图
twoway scatter y x  // Robin_1_data

图像展示:

3.2 数据拟合

在得到观测值后,我们需要对数据进行拟合,拟合的目的是为了找到待定断点,Robin Hood 算法需要运用待选断点的位置来计算最终的断点。在拟合方面,可以运用不同的工具进行拟合,例如:多项式回归(polynomial regression)、局部建模(local regression)、核回归(kernel regression)、样条回归(spline regression)等等。这里我们使用限制三次样条回归(restricted cubic spline regression),这种方法的好处是不用对回归模型进行假设,同时在回归后可以得到预测值的标准误来构造预测值的置信区间。在 Stata 中。限制三次样条回归(RCSR)可以用 mkspline 命令实现。

ssc install mkspline, replace

语法介绍 (mkspline):

mkspline stubname = oldvar [if] [in] [weight] , 
         cubic [nknots(#) knots(numlist) displayknots]
  • stubname: 定义一个新变量
  • oldvarx 变量
  • cubic: Restricted cubic spline 模式
  • options:
    • nknots: 选择结点数量
    • knots: 选择结点的位置
    • displayknots: 显示结点位置

通过 mkspline 命令可以对数据作内部处理并生成新解释的变量 stubname,随后将旧的解释变量替换成成新的解释变量作回归,详细操作会在下面展示。如果想了解更多关于 RCSR,可以参见:(Wood, 2006.)。以下为 Stata 代码实现:

mkspline k=x, cubic displayknots  // Restricted cubic spline,新变量为k~
reg y k*, vce(r)    // 回归
predict y_hat, xb   // 生成y的拟合值
predict y_se, stdp  // 获得y的标准误

//* 绘图
sort x  // x排序,画图预处理
twoway scatter y x	|| line y_hat x  // Robin_2_mkspline

图像展示:

3.3 划定平滑区域范围和选出待选断点

在 U-shaped 检验中,我们假设 x 和 y 之间存在 U 形关系而不是 V 形关系,所以在拟合线中存在一段水平的平滑区域。在这个区域中存在待选断点以及最终的断点,所以我们先要划定这个平滑区域。首先,我们要找出拟合值的最大值 y_max(或最小值);随后计算拟合值 1 个标准误的上置信区间 y_ub(或下置信区间);接着将 y_ub 大于 y_max 对应的 x 称为 xflat,对应的 y 称为 yflat,这样我们就划出了平滑区域。其 Stata 代码如下:

gen y_ub=y_hat+y_se  // yhat置信区间上边界
qui sum y_hat
gen flat=1 if y_ub>abs(r(max))  // 若yhat置信区间上界大于yhat最大值,flat==1
gen xflat=x if flat==1      // 找出flat区域对应的x,为xflat
gen yflat=y_hat if flat==1  // 找出flat区域对应的y,为yflat

随后,我们将平滑区域的中点设置为待选中点,其 Stata 代码如下:

program define which_median
    version 15.0
	//* this programme was designed for finding x which is most closer to median
	
	syntax varlist 				   	    ///
		[,								///
		gen(string) 					///
		]
	
	if "`gen'"=="" {
		di as err `"Error: gen must to be define"'
	}
	
	cap drop merry_chrismas   // 后面用于排序用途,临时变量
	qui sum `varlist',detail  // 用sum命令找中位数
	qui gen merry_chrismas=abs(`varlist'-r(p50))  // 如果数据个数为偶数,那么中位数不存在于x中,所以我们选取离中位数最近的x作为中位数
	sort merry_chrismas
	qui gen `gen'=`varlist'[1]  // 成为中位数变量
	drop merry_chrismas 
	
	
end

which_median xflat,gen(mid_xflat)

这里以程序的形式运行,which_median 命令可以寻找离中位数最近的 x 作为最终的中位数,随后生成中位数变量(变量名须自行设定)。

语法介绍 (which_median):

which_median varlist, gen(string)
  • varlist: 输入需要寻找中位数的变量
  • gen: 创建新变量作为中位数

现在,我们用图像的形式将平滑区域与待选断点显示出来:

sort x
gen mid_yflat=y_hat if x==mid_xflat
twoway scatter y x	|| line y_hat x	|| area yflat xflat || ///
 dot mid_yflat mid_xflat  // Robin_3_flat

图像展示:

3.4 利用 Robin Hood 算法寻找最终断点

根据 Robin Hood 算法的思想理念,我们先要对数据进行断点回归,分别获得两段回归其斜率系数的 t-value,随后计算百分位数调整断点的位置,来提高整体的显著性水平。那么我们先进行断点回归,这里使用的断点是平滑区域的中点,Stata 代码如下:

//* 断点回归
program define reg2
    version 15.0
	//* this programme was designed for U-shaped test 
	//* base on Uri Simonsohn(2018) pp.15
	//* R-code website: https://osf.io/zdert/
	
	syntax varlist 				   	    ///
		[,								///
		xc(string) 						///
		family(string) 					///
		robin_hood(string)				///
		savedata(string)				///
		]
	//* syntax:
		// the regression is using glm
		// var1 var2 are y and x respectively
		// xc: where to set the breakpoint, a number
		// Robin_Hood: show the ratio: t2/(t1+t2)
		// savedata: whether save new var(Y or N)
		// link:
			// Gaussian for OLS 
			// binomial for probit
			
	tokenize `varlist'
	local var1 `1'
	local var2 `2'
	//* setting default value for option
	qui sum `var2'
	cap drop max
	cap drop min
	qui gen max = r(max)
	qui gen min = r(min)
	if "`xc'"=="" {
		di as err `"Error: xc need to be define"'
	}
	else if `xc'>max {
		di as err `"Error: xc is not be observed"'
	}
	else if `xc'<min {
		di as err `"Error: xc is not be observed"'
	}
	
	// family
	if "`family'"=="" {
		local infamily "gaussian"
	}
	else {
		local infamily "`family'"
	}
	
	// Robin_Hood
	if "`robin_hood'"=="" {
		local in_Robin_Hood "N"
	}
	else if "`robin_hood'"=="Y" {
		local in_Robin_Hood "Y"
	}
	else if "`robin_hood'"=="N" {
		local in_Robin_Hood "N"
	}
	else {
		di as err `"Error: Robin_Hood must be Y or N"'
	}
	
	// savedata
	if "`savedata'"=="" {
		local insavedata "Y"
	}
	else if "`savedata'"=="Y" {
		local insavedata "Y"
	}
	else if "`savedata'"=="N" {
		local insavedata "N"
	}
	else {
		di as err `"Error: savedata must be Y or N"'
	}

	//* breakpont include in the fist line
	cap drop xlow1
	cap drop xhigh1
	cap drop high1

	qui gen xlow1=`var2'-`xc' if `var2'<=`xc'
	qui replace xlow1=0 if xlow1==. 
	qui gen xhigh1=`var2'-`xc' if `var2'>`xc'
	qui replace xhigh1=0 if xhigh1==. 
	qui gen high1=1 if `var2'>`xc'
	qui replace high1=0 if high1==. 
	
	//* breakpont include in the second line
	cap drop xlow2
	cap drop xhigh2
	cap drop high2
	
	qui gen xlow2=`var2'-`xc' if `var2'<`xc'
	qui replace xlow2=0 if xlow2==. 
	qui gen xhigh2=`var2'-`xc' if `var2'>=`xc'
	qui replace xhigh2=0 if xhigh2==. 
	qui gen high2=1 if `var2'>=`xc'
	qui replace high2=0 if high2==. 
	
	//* glm1 (regression process for first line)
	qui glm `var1' xlow1 xhigh1 high1,family(`infamily') vce(r)
	scalar beta1=_b[xlow1]
	scalar t1=abs(_b[xlow1]/_se[xlow1])
	qui test xlow1
	scalar p1=r(p)
	
	//* glm 2 (regression process for second line)
	qui glm `var1' xlow2 xhigh2 high2,family(`infamily') vce(r)
	scalar t2=abs(_b[xhigh2]/_se[xhigh2])
	scalar beta2=_b[xhigh2]
	qui test xhigh2
	scalar p2=r(p)
	
	//* the ratio
	scalar ratiopp=t2/(t1+t2)
	
	//* dis 
	dis "*********** For first line ***********"
	dis ""
	dis "slope: "beta1
	dis "t-value: "t1
	dis "p-value: "p1
	dis ""
	dis "*********** For second line ***********"
	dis ""
	dis "slope: "beta2
	dis "t-value: "t2
	dis "p-value: "p2
	
	//* the result of Robin Hood algorithm
	if "`in_Robin_Hood'"=="Y" {
		dis ""
		dis "======================================"
		dis "Robin Hood result: "ratiopp
		dis "======================================"
	}
	
	//* data
	if "`insavedata'"=="N" {
		drop xlow* xhigh* high*
	}
	
end

reg2 y x, xc(mid_xflat) robin_hood(Y)

语法介绍 (reg2):

reg2 varlist, xc(string) family(string) robin_hood(string) savedata(string)
  • varlist: 依顺序输出变量,第一个变量应为 y,第二个变量为 x
  • options:
    • xc(string): 输入断点的位置
    • family(string): 回归采用广义线性模型,所以需要设置 “family”,默认值为 “Gaussian”
    • robin_hood(string): Robin Hood 算法的计算结果选项,“Y” 表示显示结果,“N” 表示隐藏结果
    • savedata(string): 变量保存选项,在计算过程中会生成多个变量,“Y” 选项代表保持,“N” 为不保存

以上的代码是以程序的形式运行,直接将代码复制在 DO 文档运行即可,其运行的结果如下:

*********** For first line ***********

slope: .86045729
t-value: 34.270446
p-value: 2.16e-257

*********** For second line ***********

slope: -.09809849
t-value: 1.8139696
p-value: .06968241

======================================
Robin Hood result: .05027017
======================================

结果分别显示两段回归线的斜率和斜率系数的 t 值以及 p 值。这里我们主要关注斜率系数和 p 值,两个斜率系数互为异号,说明 x 和 y 存在 U 形关系;再来看 p 值,第一段回归的斜率系数是显著的,第二段的斜率系数 p 值仅为 0.07 左右,这个结果似乎不能有力地说明 x 和 y 之间的 U 形关系。为此,我们用 Robin Hood 算法进一步寻找合适的断点。

根据 Robin Hood 算法的理念,我们应该给予第二段回归更多的观测值。算法的计算公式为 t2/(t1+t2)。计算的结果将会作为平滑区域的百分位数,对应的 x 值就是最终的断点。以下为 Stata 代码实现:

program define qq
    version 15.0
	// quantile
	syntax varlist 				   	    ///
		[,								///
		p(string) 						///
		gen(string)						///
		]
	
	if "`gen'"=="" {
		di as err `"Error: gen must to be define"'
	}
	
	qui sum `varlist'
	scalar quantilepoint=r(min)+(`p'*(r(max)-r(min)))
	gen `gen'=abs(quantilepoint)
	
end

qq xflat, p(.05027017) gen(breakpoint)

程序 qq 的功能是输出对应的百分位数可以寻找相应的 x

语法介绍 (qq):

qq varlist, p(string) gen(string)
  • varlist: 输入变量
  • p(string): 输入需要查找的百分位数
  • gen(string): 将结果输出为变量

我们用图像的形式观察结果,其代码如下:

sort x
gen y_breakpoint=y_hat if x==breakpoint
twoway scatter y x	|| line y_hat x	|| area yflat xflat || ///
 dot y_breakpoint breakpoint  // Robin_4_breakpoint

图像展示:

可以看到断点的位置在平滑区域的最左边,在下次的断点回归中,第二段回归会拥有更多的观测值,这能够提升回归的显著性水平。现在,我们利用计算得出的断点进行断点回归。

reg2 y x, xc(breakpoint)

其结果为:

*********** For first line ***********

slope: .91059031
t-value: 32.790951
p-value: 7.92e-236

*********** For second line ***********

slope: -.14647279
t-value: 3.2746319
p-value: .001058

可以看到,第二段回归的斜率系数其 p 值较之前有大幅度地下降,整体的显著性水平有所改善。现在,从结果来看我们能进一步确定,x 和 y 存在着 U 形关系。

4. 总结

最后我们来总结一下 U-shaped 检验的步骤:

  1. 拟合数据
  2. 寻找 y_hat, y_hat 的最大值 (yhat_max) 和 y_hat 的上置信区间 (yhat_ub)
  3. 划定平滑区域,将 yhat_ub 大于 yhat_max 的 y 划定为 yflat,对应的 x 为 xflat
  4. 在平滑区域寻找中点,作为待选断点 mid_xflat
  5. 利用待选断点进行断点回归,获得两段回归其斜率系数的 t 值
  6. 计算百分位数 (t2/(t1+t2))
  7. 在平滑区域中,对应百分数的 x 为最终断点
  8. 利用最终断点计算断点回归,观察两段回归其斜率系数的结果

在最后需要指出这种方法的两个缺点:

  • 其一,是如果真实的函数关系是 N 形、W 形或者是 X 形,那么采用这种检测方法会出现误导的结果。
  • 其二,断点回归是基于线性回归,任何影响线性回归的因素都会导致检测结果有所偏差,例如,有效性,可解释性,偏差,稳健性等。

5. 参考资料和相关推文

  • Deming, David J., and Kadeem Noray. "Earnings dynamics, changing job skills, and STEM careers." The Quarterly Journal of Economics 135.4 (2020): 1965-2005. -PDF-
  • Simonsohn, Uri. "Two lines: A valid alternative to the invalid testing of U-shaped relationships with quadratic regressions." Advances in Methods and Practices in Psychological Science 1.4 (2018): 538-555. -PDF-
  • Wood, Simon N. Generalized additive models: an introduction with R. chapman and hall/CRC, 2006. -PDF-

6. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh U型 断点, m
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxh 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh