温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
谢雁翔 (南开大学),xyxmask1995@163.com
钟舜斌 (北京工商大学),1315829300@qq.com
编者按: 本文部分参考了游万海老师的分位数回归讲义和陈强老师的《高级计量经济学及Stata应用》,特此致谢!
目录
在此前的推文中,我们对分位数回归和面板分位数回归都做过简单介绍,参见:
本文通过几篇论文的实操对分位数回归进行更为全面的介绍,内容涉及:分位数回归的基本思想、面板分位数回归、边际效应估计及图示等。
一般回归模型中着重考察的是解释变量
从实际来看,分位数信息也十分重要,比如,在评估某项干预对受众群体的影响时,我们不但希望了解干预的「平均」影响,更希望掌握干预对位于特征分布不同位置 (分布末端或顶端) 人群的「异质性」影响。
针对样本数据「异质性」特征,常用做法是根据数据特征进行「分组回归」,但这样的做法会导致「样本数据的损失」。为此,Koenker 和 Bassett (1978) 提出「分位数回归 (Quantile Regression, QR)」,并使用残差绝对值的加权平均 (如,
由此可知,通过设置不同的分位点,分位数回归模型可以全面的刻画解释变量与被解释变量之间的关系。此外,相比于普通的线性回归 (均值回归),分位数回归的估计结果对「偏态、多峰和异常值数据」更为稳健。
对一个随机变量
有一组数据
即:
对
可得:
将标量
令
对应参数
(1) 如果损失函数定义为二次函数
(2) 如果损失函数定义为
(3) 如果损失函数定义为:
则:
即为对应下式的解:
损失函数就是上述例子中模型所表现的误差,最小化损失函数的过程其实是通过损失函数反过来优化模型参数。设分位点
当
在实际运用中,选择不对称的分位数损失函数往往可以反映因变量的分位数分布,不同分位数的分布可能由于各种原因存在差异。目前,分位回归模型的估计算法主要有单纯形法、内点算法和平滑算法。
clear
set seed 1 //设置种子值,保证结果可重复
qui set obs 100 //设置100观测值
egen x = seq(), from(1) to(100) //生成1到100序列
gen sig = 0.1 + 0.05*x //sig代表方差,方差随着x在变化,不是同方差,异方差不对称
scalar b_0 = 6
scalar b_1 = 0.1
gen e = rnormal(0,sig) //生成误差项
gen y = b_0 + b_1*x + e
twoway scatter y x, plotregion(fcolor(white)) graphregion(fcolor(white))
*散点图呈现喇叭状,非均匀分布考虑分位数回归,若均匀分布不考虑分位数回归
从上图可以看出,随着
分位数回归模型的基本步骤:
绘制散点图,判断数据分布 (对称/非对称) 是否适合于用分位数回归模型;
建模估计 qreg
,sqreg
,bsqreg
,qreg2
(聚类稳健标准误);
系数差异性检验 wald 检验;
系数可视化 grqreg
。
在 Stata 的进行分位数回归估计的命令中,sqreg
较为常用,详见 help sqreg
。
*计算单个分位点
sysuse auto, clear
qreg price weight length foreign, quantile(.25)
est store m0
*多个分位点,从0.3-0.9,可用qreg的循环语句也可使用sqreg命令
forvalues i = 0.3(0.3)0.9{
qreg price weight length foreign, quantile(`i')
local j = `i' * 10
local j = int(`j')
est store m`j'
}
local m "0.25RQ 0.3RQ 0.6RQ 0.9RQ"
esttab m0 m3 m6 m9, mtitle(`m') nogaps star(* 0.1 ** 0.05 *** 0.01)
----------------------------------------------------------------------------
(1) (2) (3) (4)
0.25RQ 0.3RQ 0.6RQ 0.9RQ
----------------------------------------------------------------------------
weight 1.832*** 1.705** 5.975*** 9.666***
(2.89) (2.59) (4.00) (4.53)
length 2.846 5.090 -100.1* -248.8***
(0.13) (0.23) (-1.96) (-3.40)
foreign 2209.9*** 2148.6*** 3731.0*** 4164.4***
(5.24) (4.90) (3.75) (2.93)
_cons -1879.8 -1843.9 5901.2 25191.0***
(-0.76) (-0.72) (1.01) (3.02)
----------------------------------------------------------------------------
N 74 74 74 74
----------------------------------------------------------------------------
t statistics in parentheses
* p<0.1, ** p<0.05, *** p<0.01
不同分位数下系数差异检验用 Wald 检验,具体如下:
假设
相应的 wald
检验统计量为
这个统计量服从卡方分布
其中:q
为待检验的线性假设个数。
set seed 1 //设置种子值
sysuse auto, clear
*sqreg不用循环语句,默认情况下reps为20,bootstrap100次
sqreg price weight length foreign, quantile(.25 .5 .75) reps(100)
*wald检验
test [q25]weight=[q75]weight //检验0.25分位点处与0.75分位点处是否存在显著差异,拒绝原假设则存在差异
(1) [q25]weight - [q75]weight = 0
F( 1, 70) = 11.62
Prob > F = 0.0011
*ssc install grqreg,all replace //安装外部命令
sysuse auto, clear
qreg price mpg headroom
grqreg, cons ci ols olsci title(Fig.1a Fig.1b Fig.1c) //cons表示常数项绘图,ci表示置信区间,ols表示绘出普通ols,olsci表示ols回归置信区间
Median regression Number of obs = 74
Raw sum of deviations 71102.5 (about 4934)
Min sum of deviations 64549.83 Pseudo R2 = 0.0922
------------------------------------------------------------------------------
price | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg | -127.3333 60.78234 -2.09 0.040 -248.5299 -6.136791
headroom | -130.6667 415.6721 -0.31 0.754 -959.4933 698.16
_cons | 8272.333 2158.205 3.83 0.000 3968.995 12575.67
------------------------------------------------------------------------------
上图分别为 cons (常数项) 、mpg 与 headroom 的分位数估计结果,横轴的虚线为 ols 估计结果及其置信区间。但也存在问题,比如纵轴坐标重叠和背景非白色。
面板分位数回归,顾名思义即将截面分位数回归模型扩展至面板数据框架下,考虑了研究单元之间的个体异质性。特别注意的是,面板分位数回归为非线性关系,对于个体效应的估计不能通过 i.id 虚拟变量的形式控制。
惩罚分位 (Penalized quantile regression with fixed effects),应用了LASSO的思想:
缺点: (1) 当
实现:R 中的 rqpd
包进行估计,rqpd
不支持 R 3.6.3 及以后版本,建议直接在 R 中运行,更多内容见「rqpd-package: Regression quantiles for panel data」,也可以使用 Stata 外部命令rcall
运行,详细请参考 「Rcall:Stata 与 R 的无缝对接」。
温馨提示: 文中链接在微信中无法打开,请点击「阅读原文」。
## rcall 可以从Stata中运行R命令,为简便运算可以直接复制粘贴进R运算
install.packages("rqpd", repos="http://R-Forge.R-project.org")
install.packages("plm", repos="http://R-Forge.R-project.org")
library(rqpd) ##加载面板分位数回归包
library(plm) ##加载面板数据回归包
data(Produc,package="plm") ##调用plm包自带数据集进行面板估计
mydata = Produc
m <- 17 ##17年
n <- 48 ##美国48州
s <- as.factor(rep(1:n,rep(m,n))) ##声明面板数据结构变量
x <-cbind(mydata$pcap,mydata$hwy) ##数据中取出解释变量x
y <- mydata$unemp ##数据中取出被解释变量y
##rqpd直接运算公式,个体效应惩罚系数lambda为1,系统默认分位点是0.25,0.5,0.75
fit <- rqpd(y ~x | s, panel(lambda = 1))
sfit <- summary(fit)
sfit
sfit$coefficients ##能一并取出估计的个体效应系数
##依次定义九个分位点同时定义tau和tauw
fit1 <- rqpd(y ~x | s, panel(taus=seq(0.1,0.9,by=0.1),tauw=rep(1/9,9),lambda = 1))
sfit1 <- summary(fit1)
sfit1
两阶段估计法 (Two-step estimator),分位数回归为非线性运算,面板数据不能通过组内去心的方法解决。
思想:先消除个体效应,再进行估计。与 Koenker (2004) 的相同点:将
温馨提示: 文中链接在微信中无法打开,请点击「阅读原文」。
优点:(1) 个体效应允许随着分位点在变动 (location-scale);(2) 可以纳入其他内生解释变量;(3) 计算速度较快。
缺点: 假定变量系数与分位点之间关系是单调的(要不单调上升,要不单调下降),即不存在分位交叉 (No Quantile-Crossing),画图无意义。
实现: Stata
中 xtqreg
进行估计,详细参见 help xtqreg
。
关于 xtqreg
系数可视化,请参考以下代码:
clear //清空数据
set obs 100 //设置100观测值
egen country_id = seq(),b(10) //每10个重复构造整数序列,生成变量country_id
bysort country_id: gen time = _n //根据分组country_id生成时间变量
set seed 12345 //设置种子值确保可重复
gen y = uniform() //生成被解释变量y服从正态分布
forv i = 1/8{
gen z`i' = uniform()
} //循环语句生成8个解释变量
local x1 "z1 z2 z3 z4 z5 z6 z7 z8" //局部宏定义所有解释变量
*list `x1' //列出所有被解释变量
local v_num : word count `x1' //利用扩展函数 word count string 得到x1中变量的个数
di `v_num' //v_num表示最大的循环个数,8个
est clear
mat result1 = J(`v_num',1,.) //生成空矩阵result1,8行1列
matrix rownames result1 = `x1' //矩阵行命名为对应的被解释变量
mat list result1 //列出空矩阵
tempname a b b1 //设定临时矩阵变量名
forvalues quantile = 0.1(0.1)0.9 { //从0.1-0.9没隔0.1开始循环
local wanted : di %2.1f `quantile' //列出循环分位点长度设置为2,保留一位小数,定义为`wanted'
local wanted = `wanted'*10 //`wanted'为每一分位点重复十次
xtqreg y `x1', i(country_id) q(`quantile') //面板分位数回归
matrix `a' = r(table) //对应各个分位点估计的`x1'系数
mat list `a' //列出面板分位数回归`a'矩阵
matrix `b' = `a'[1,1..`v_num'] \ `a'[2,1..`v_num'] //仅取出面板分位数回归`a'矩阵的前两行构成`b'矩阵,第一行是系数。第二行是标准误,第四行是p值
mat list `b' //列出`b'矩阵
mat `b1' = `b''
matrix colnames `b1' = b_q`wanted' p_q`wanted' //`b1'为`b'矩阵的转置
mat list `b1' //列出`b1'矩阵
matrix result1 = (result1, `b1') //`b1'数据导入列出空矩阵result1
}
mat list result1
matselrc result1 cresult, c(2/19) //19列 = 2*分位个数+1列无用项(第一列),保留有用数据列
mat list cresult
local newnames : colfullnames cresult
di "`newnames'"
local cn : word count `newnames'
di `cn'
svmat cresult, names(reg) //矩阵数据转换为变量
keep reg*
format reg* %4.3f //设定格式为保留三位小数
drop if reg1==. //删除缺失值,行数为8,8个解释变量:变量为18,对应不同分位点的系数和标准误;
forv i = 1/9 { /*9:分位个数*/
local k1 = 2*`i' - 1 //奇数列为系数
local k2 = 2*`i' //偶数列为标准误
gen up_`i' = reg`k1' + 1.96*reg`k2' //上界
gen low_`i' = reg`k1' - 1.96*reg`k2' //下界
}
preserve
keep reg1 reg3 reg5 reg7 reg9 reg11 reg13 reg15 reg17 //保留系数列
rename (reg1 reg3 reg5 reg7 reg9 reg11 reg13 reg15 reg17) (reg1 reg2 reg3 reg4 reg5 reg6 reg7 reg8 reg9) //重命名
xpose,clear //转置数据
renvars v*, prefix(coef_) //重命名加前缀coef_
save coef.dta,replace //储存系数数据
restore //preserve+restore运行不保存
keep up* low* //保留上下界数据
xpose,clear //转置数据
preserve
keep if mod(_n,2)==1 //保留奇数行作为上界数据
renvars v*, prefix(up_) //重命名加前缀up_
save up.dta,replace //储存上界数据
restore
preserve
keep if mod(_n,2)==0 //保留偶数行作为下界数据
renvars v*, prefix(low_) //重命名加前缀low_
save low.dta,replace //储存下界数据
restore
use up.dta,clear //打开上界数据
merge 1:1 _n using low.dta //一对一合并数据
keep if _merge==3 //合并匹配成功为_merge==3,保留合并成功行
drop _merge //删除合并信息
merge 1:1 _n using coef.dta //一对一合并数据
keep if _merge==3
drop _merge //删除合并信息
matrix input myrvec = (0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90)
mat list myrvec
mat qq = myrvec' //qq矩阵为myrvec的转置
mat list qq
svmat qq, names(qqp) //矩阵数据转换为变量
forv i = 1/8{
twoway (line up_v`i' qqp,lc(black*1.4) lpattern(dash) xlabel(0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90)) ///
(line low_v`i' qqp,lc(black*1.4) lpattern(dash) xlabel(0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90)) ///
(line coef_v`i' qqp, lc(black*2) lw(thick)) ,xtitle("Quantile") ytitle("Coefficient(z`i')") plotregion(margin(zero)) graphregion(color(white)) legend(off)
graph save g_`i'
}
graph combine g_1.gph g_2.gph g_3.gph g_4.gph g_5.gph g_6.gph g_7.gph g_8.gph
qregpd
估计是实现广义分位数估计 genqreg
的特例,其解决了由其他固定效应分位数估计引起的一个基本问题:包含个体固定效应会改变对处理变量的估计系数的解释。由于标准差的计算有时可能收极值的影响,可以使用 MCMC 方法或网格搜索方法来估计广义分位数回归。关于该命令更多详细介绍,请参考 help qregpd
。相关代码如下:
*ssc install qregpd, replace
*ssc install moremata, replace
sysuse nlswork, clear
qregpd ln_wage tenure union, id(idcode) fix(south)
estimates store QREG1
matrix list e(gamma) //显示south固定效应估计值
qregpd ln_wage tenure union, id(idcode) fix(year)
estimates store QREG2
matrix list e(gamma) //显示year时点固定效应估计值
*带工具变量的情形--采用MCMC优化
*tenure为内生变量,union为外生变量做自己的工具变量,ttl_exp wks_work为工具变量
qregpd ln_wage tenure union, id(idcode) fix(year) optimize(mcmc) draws(1000) burn(100) arate(.5) noisy instruments(ttl_exp wks_work union)
mat list e(gamma)
*估计1/4分位数的结果--默认是0.5
qregpd ln_wage tenure union, q(0.25) id(idcode) fix(year) optimize(mcmc) draws(1000) burn(100) arate(.5) noisy instruments(ttl_exp wks_work union)
mat list e(gamma)
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟,课程主页 Stata 33 讲 - 连玉君, 每讲 15 分钟. Stata 小白的取经之路 - 龙志能,时长:2 小时,课程主页 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 因果推断, 空间计量,寒暑假班等 | |
⭕ 数据清洗系列 | 游万海 | 直播, 88 元,已上线 |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD