(倒)U型+RDD：利用断点回归检验 U 形关系
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 2022 面板数据因果推断专题
作者：郑伊达 (中山大学)
邮箱：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 应用背景
1.2 二次回归检验 U 形关系的不足之处
1.3 检验方法的理念
2. 断点回归
3. 实现过程
3.1 生成伪随机数
3.2 数据拟合
3.3 划定平滑区域范围和选出待选断点
3.4 利用 Robin Hood 算法寻找最终断点
4. 总结
5. 参考资料和相关推文
6. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
1. 检验方法介绍
1.1 应用背景
在实证分析中，有时我们能够观察到 与 呈现出 U 形关系，又或者是倒 U 形关系。例如，工资收入与年龄之间就呈现出倒 U 形关系，人们多数在年轻时收入水平较低，中年时期其收入水平达到最高值，随后又因退休而导致收入减少 (Deming, 2019)。这些结论都是我们通过可观测数据观察到的结果，那么这样的结果能反映事实吗？有多大程度能够反映事实？基于这样的问题，我们需要进一步检验两者之间的关系，由此出现了对 U 形关系的检验。
1.2 二次回归检验 U 形关系的不足之处
在实证过程中，当研究者认为 和 之间存在一个 U 形关系，那么使用二次回归会是一个方便又快捷的方法，来对数据进行拟合，从而判断它们是否存在 U 形关系。但是，这种方法存在一定的风险，如果真实的回归模型不是一个二次方程，最终结果可能会误导研究者。
举一个例子，我们先定义一个函数关系式：。很明显， 和 之间不是 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
上面的代码构造了一系列的伪随机观测值，真实的函数关系为 ，其残差项服从正态分布，均值为 0，方差为 1。随后我们生成一个 ，作二次回归，结果如下：
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
------------------------------------------------------------------------------
从二次回归的结果可以看到，二次项的回归系数为负数，而且显著。因此，结果带给我们的结论是， 和 之间存在倒 U 形关系。从这个例子中，可以显示出用二次回归来检测 U 形关系存在的问题，其核心问题是我们错误的假设了函数形式。
如果想了解跟多关于二次回归所引发的问题，可以参见连享会往期推文 ：平方项 = 倒U型 ？。
1.3 检验方法的理念
在检验 U-shaped 关系之前，我们先对 U 形关系进行定义： 存在一个中间值 xc，小于 xc 的 为低数值组，大于 xc 的 为高数值组；在低数值组的 与高数值组的 两者之间为异号。简单来说，在 U 形关系中， 数值较低的部分其线段斜率是负的， 数值较高的部分其斜率是正的。另外，在 U-shaped 的中还包含额外的特征，例如，对称 vs 不对称，连续 vs 不连续，有极值 vs 无极值，想要深入研究这些特征需要用额外的方法去检测，在本推文中并不涉及。
这里我们引入一个概念，线性回归所计算出来的斜率系数是平均斜率，无论 和 的真实函数形式是如何。因此将 和 作两组线性回归，其两个斜率系数如果是异号且显著，可以判断 和 存在 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
在这个例子中，我们假设真实的函数形式为: 。通过简单的计算我们可以得知 从 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:
if , and 0 otherwise if , and 0 otherwise if , and 0 otherwise is the (optional) matrix with covariates, and its vextor of coefficients.
完成断点回归后，可以得到两条线性回归的拟合值，对于 组别的拟合值为：；对于 组别的拟合值为：。
在进行断点回归时，需要设置断点。那么如何找到一个合适的断点 xc，使两段回归的斜率系数其显著水平最大化？这里需要用到 Robin Hood 算法 (see e.g., Simonsohn, 2018)，来寻找最合适的断点。这种算法不用对函数形式 进行假设，不用对 的分布进行假设，也不用对残差项的分布进行假设。其算法有三个核心概念：(i) 两段回归最终要求两个斜率系数其显著性最大化，为了能进一步提升整体的显著水平，我们的目标是提升显著性较弱的回归；(ii) 其方法是令它们变得陡峭一点；(iii) 实现方法是给予它们更多的观测值。总的来说，就是通过改变断点的位置，将观测值从显著水平较高的回归转移至显著性水平较低的回归中。
🍎 连享会 · 2022 面板数据因果推断专题课程
👇 点击下方海报可了解课程详情！
👇 也可扫码获取课程介绍
3. 实现过程
3.1 生成伪随机数
首先，我们需要生成一些伪随机数，来作为实践对象。在这里，我们定义一个真实的函数形式 ：
若 ，其斜率为 1，即 ； 若 ，则 ； 若 ，其斜率为 0.15，即 。
随后，我们在 中加入随机扰动项 ，生成伪随机观测值。数据生成过程可以表示如下：
其中，，。
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: 定义一个新变量oldvar: 变量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 检验中，我们假设 和 之间存在 U 形关系而不是 V 形关系，所以在拟合线中存在一段水平的平滑区域。在这个区域中存在待选断点以及最终的断点，所以我们先要划定这个平滑区域。首先，我们要找出拟合值的最大值 y_max（或最小值）；随后计算拟合值 1 个标准误的上置信区间 y_ub（或下置信区间）；接着将 y_ub 大于 y_max 对应的 称为 xflat，对应的 称为 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: 依顺序输出变量，第一个变量应为 ，第二个变量为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 值，两个斜率系数互为异号，说明 和 存在 U 形关系；再来看 p 值，第一段回归的斜率系数是显著的，第二段的斜率系数 p 值仅为 0.07 左右，这个结果似乎不能有力地说明 和 之间的 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 检验的步骤：
拟合数据 寻找 y_hat, y_hat 的最大值 (yhat_max) 和 y_hat 的上置信区间 (yhat_ub) 划定平滑区域，将 yhat_ub 大于 yhat_max 的 y 划定为 yflat，对应的 x 为 xflat 在平滑区域寻找中点，作为待选断点 mid_xflat 利用待选断点进行断点回归，获得两段回归其斜率系数的 t 值 计算百分位数 (t2/(t1+t2)) 在平滑区域中，对应百分数的 x 为最终断点 利用最终断点计算断点回归，观察两段回归其斜率系数的结果
在最后需要指出这种方法的两个缺点：
其一，是如果真实的函数关系是 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
专题：交乘项-调节 追本溯源，U型关系你用对了么？ utest：检验U型和倒U形关系 平方项 = 倒U型 ？ 专题：Stata绘图 Stata绘图：世行可视化案例-条形图-密度函数图-地图-断点回归图-散点图 专题：断点回归RDD Stata：基于大带宽的断点分位数回归稳健推断-rdqte Stata：RDD-DID-断点回归与倍分法完美结合 RDD断点回归：多个断点多个分配变量如何处理 Stata+R：一文读懂精确断点回归-RDD RDD：离散变量可以作为断点回归的分配变量吗？ RDD：断点回归可以加入控制变量吗？ 断点回归RDD：样本少时如何做？ Stata：时间断点回归RDD的几个要点 Stata：断点回归RDD简明教程 RDD：断点回归的非参数估计及Stata实现 Stata: 两本断点回归分析 (RDD) 易懂教程 Stata: 断点回归 (RDD) 中的平滑性检验
课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。