Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:郑伊达 (中山大学)
邮箱: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-
目录
在实证分析中,有时我们能够观察到
在实证过程中,当研究者认为
举一个例子,我们先定义一个函数关系式:
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
上面的代码构造了一系列的伪随机观测值,真实的函数关系为
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-shaped 关系之前,我们先对 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
在这个例子中,我们假设真实的函数形式为:
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
------------------------------------------------------------------------------
在上面的描述中,提到要将数据分为两部分别进行线性回归,其断点为 xc。这里我们引入断点回归 (see e.g., Marsh & Cormier, 2001, p. 7),其通用的回归模型为:
Where:
完成断点回归后,可以得到两条线性回归的拟合值,对于
在进行断点回归时,需要设置断点。那么如何找到一个合适的断点 xc,使两段回归的斜率系数其显著水平最大化?这里需要用到 Robin Hood 算法 (see e.g., Simonsohn, 2018),来寻找最合适的断点。这种算法不用对函数形式
首先,我们需要生成一些伪随机数,来作为实践对象。在这里,我们定义一个真实的函数形式
随后,我们在
其中,
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
图像展示:
在得到观测值后,我们需要对数据进行拟合,拟合的目的是为了找到待定断点,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
图像展示:
在 U-shaped 检验中,我们假设
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
图像展示:
根据 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 值,两个斜率系数互为异号,说明
根据 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 形关系。
最后我们来总结一下 U-shaped 检验的步骤:
在最后需要指出这种方法的两个缺点:
Note:产生如下推文列表的 Stata 命令为:
lianxh U型 断点, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh