分位数回归和Stata这对组合我用了快十年了。最开始接触这个主题是因为做收入研究时发现一个让人抓狂的问题OLS回归跑了半天教育年限对收入的“平均效应”是正的、显著的但我想知道的明明是——教育对低收入人群、中等收入人群、高收入人群的影响是不是一样均值回归完全没法回答这个问题。后来摸到分位数回归像打开了另一扇门才发现过去做的很多实证模型其实只讲了故事的一半。今天把我这些年用Stata做分位数回归的经验整理出来从原理到实操再到图表呈现和写论文时的汇报规范一次说清楚。这篇文章适合正在做实证研究、写论文、做政策评估或者分析调查数据的朋友。不管你是刚接触Stata的初学者还是已经用过qreg命令但始终搞不懂结果怎么解读、图形怎么画的老手这篇都值得从头过一遍。我会尽量少堆公式多用大白话和真实案例来解释——分位数回归的本质、Stata里几种命令的区别、标准误到底选哪种、多分位数联合估计怎么做、结果图怎么画得让审稿人眼前一亮以及我踩过的那些坑全部交代清楚。1. 分位数回归到底在解决什么问题1.1 均值回归的“盲区”为什么OLS经常不够用OLS普通最小二乘回归是所有计量入门课的第一课它的目标函数是残差平方和最小化估计的是因变量的条件均值E(y|x)与自变量之间的关系。问题在于条件均值是一个“汇总统计量”它把整个条件分布的差异全部平均掉了。一旦数据存在偏态、异方差、截尾或者异常值均值回归的结果容易失真甚至会给出极具误导性的结论。我举个特别直观的例子研究“参加职业培训对收入的影响”。用OLS估计出来的平均处理效应可能是正的——平均来看培训确实提高了收入。但真实情况很可能是培训对原本收入就在低分位比如10%分位的人帮助极大对中等收入人群帮助一般对高收入人群几乎没作用甚至还有负作用。这三个群体的效应平均在一起得到一个“看起来还行”的正值然后政策制定者就拿着这个结果去推广项目了。可问题是培训资源是有限的如果只能服务一部分人你该优先给谁均值回归回答不了。分位数回归的思路完全不同——它估计的是因变量在特定分位点比如10%、50%、90%上的条件分位数与自变量之间的关系。还是同一个数据集分位数回归可以告诉你职业培训使低收入群体10%分位收入提高了30%但对高收入群体90%分位的影响统计上不显著。这才是真正能指导决策的信息。另一个非常典型的场景是金融领域的VaR在险价值研究。VaR本质上就是一个条件分位数——比如“5%分位数的日收益率”。你要是用OLS去估计均值收益率跟风险管理的需求完全脱节。分位数回归在金融风险度量、医疗费用分析、教育经济学、劳动经济学这些领域被广泛使用根本原因就是它们都关心分布的尾部、关心异质性效应而不是一个平均数字。1.2 分位数回归的直觉理解不看平均看分布分位数的概念大家中学就接触过中位数就是50%分位数。分位数回归可以理解为不把数据“压缩”成一个均值而是把整个条件分布切成一排截面每个截面都做一次回归。你可以选择任意分位点tautau在0到1之间估计的是y在x条件下的第tau个分位点。举一个生活化的类比。想象你要了解一条路上不同时段的车速与事故率的关系。均值回归告诉你的结果是“平均车速每提高10公里平均事故率增加X”——这是把早高峰、午间、深夜全部混在一起算的平均数。分位数回归则是把事故率从低到高排队分别看“事故率最低的10%路段”“事故率中等的50%路段”“事故率最高的10%路段”然后问车速对这三类路段的影响分别是多少很可能结论是车速对本来就安全的路段影响不大但对本来事故率就极高的路段影响巨大。知道这个区别你才知道治理重点应该放在哪。分位数回归还有一个巨大的优势它对异常值不敏感。因为它用的是绝对离差最小化而不是平方离差最小化即便因变量出现极大的极端值比如收入数据里的亿万富翁只要它们没有改变某个分位点的位置回归结果就不会像OLS那样被“拉偏”。所以当你面对厚尾分布、右偏严重的数据时分位数回归往往比OLS稳健得多。2. Stata实现分位数回归的核心原理2.1 从最小二乘到最小绝对离差目标函数的变化理解分位数回归的数学原理不需要啃艰深的优化理论抓住目标函数就行。OLS是让残差平方和最小即 min Σ(yi - xiβ)²它的核心特性是“对大的残差惩罚更重”这也是为什么离群值会严重影响估计结果。分位数回归的目标函数换成了对残差的加权绝对离差最小化min Σ ρτ(yi - xiβ)这里的ρτ是分位数损失函数也叫tick loss因为形状像钟表的指针当残差为正yi ≥ xiβ即实际值高于拟合值时权重是τ当残差为负yi xiβ即实际值低于拟合值时权重是1-τ。这个设计非常巧妙。当τ0.5时正负残差的权重相等此时目标函数等于最小化绝对离差之和得到的解就是条件中位数回归也叫LAD回归Least Absolute Deviations。当τ取0.9时模型给正残差的权重是0.9给负残差的权重是0.1——这意味着它更加在意“预测值低于实际值”的情况从而把拟合平面往上推估计的实际上是条件分布较高的位置。这个目标函数没有显式解需要借助线性规划方法求解。Stata的qreg命令内部用的就是基于单纯形法的优化算法。好消息是你完全不需要手动处理这些计算细节但理解目标函数的意义很重要——它能帮你想清楚为什么不同分位点的回归系数会不同以及为什么极端分位点如0.1、0.9的估计通常不如中位数的稳定。2.2 分位数回归的系数到底怎么解释分位数回归系数的解释比OLS多绕一个弯但一旦绕过去就一劳永逸。OLS的系数解释是“x每增加一个单位y的平均值变化多少”而分位数回归的系数解释是“在控制其他变量的条件下x每增加一个单位y的第τ个条件分位数变化多少”。用工资数据举例。假设你做了τ0.1和τ0.9两个分位点的分位数回归得到教育年限的系数分别是0.08和0.15。那么应该这样读教育年限每增加一年收入分布最底部10%人群的收入水平大约提高8%而收入分布最顶端10%人群的收入水平大约提高15%。这直接说明教育对高收入群体的边际回报更高——教育溢价在分布的不同位置存在显著差异。这里有一个初学者特别容易混淆的点分位数回归不是在把“人”分成高收入、低收入两组来分别回归。它是在利用全样本信息但给不同位置的观测赋予不同的权重从而刻画条件分布不同位置的边际效应。这一点要牢牢记在心里写论文时表述错了会闹笑话。另一个重要概念是“条件分位数”和“无条件分位数”的区别。qreg估计的是条件分位数即给定x的前提下y的分位数它回答的问题是“在可观测特征相同的人群中x对y分布位置的影响”。如果你的研究想要回答“对整个人群收入分布的某个分位点的影响”那是无条件分位数回归UQR要用RIF回归recentered influence function来实现。Stata里有rifreg命令需要ssc install rifreg但那是另一个专题了这篇文章先聚焦qreg这条主线。2.3 为什么要同时做多个分位点系数差异检验的价值很多人做分位数回归只挑一个分位点跑比如只跑中位数回归这其实浪费了方法的核心优势——异质性分析。分位数回归最有价值的用法是同时估计多个分位点比如0.1、0.25、0.5、0.75、0.9然后检验不同分位点的系数是否存在显著差异。假设教育年限在0.1分位点的系数是0.08在0.9分位点的系数是0.15这两个系数表面上看差挺多但它们各自的置信区间可能很宽重叠严重——统计上可能根本无法拒绝“两者相等”的原假设。如果直接下结论说“教育回报率在高低收入群体间存在显著差异”审稿人一定会质疑。Stata里进行这个检验有两种常规做法。第一种是使用sqreg命令这个命令通过bootstrap同时估计多个分位数方程并存储分位数回归系数的联合方差-协方差矩阵然后用test命令检验跨分位点的系数相等性。第二种是用iqreg命令估计分位数区间直接看置信区间是否重叠。我实际使用中更推荐sqregtest的组合因为操作直接、输出规范受到的质疑最少。3. Stata中分位数回归的完整实操流程3.1 基础命令qreg从安装到第一个回归先说明一下qreg、sqreg、bsqreg、iqreg这几个命令都是Stata官方自带的不需要额外安装装好Stata直接能用。如果你的Stata版本比较老比如12以下建议升级到新版因为旧版的分位数回归在速度和标准误计算上有不少限制。做任何回归之前先看数据概况这一步不能省。我习惯用describe、summarize、histogram来摸清数据的分布形态。尤其是因变量的分布——如果它明显偏态比如收入的右偏长尾那就是用分位数回归的强信号。下面用Stata自带数据集auto演示最基础的操作sysuse auto, clear summarize price weight mpg * 标准分位数回归默认估计中位数 qreg price weight mpg * 指定分位点估计10%分位 qreg price weight mpg, quantile(0.1)qreg输出结果里你可以看到和OLS基本一样的结构系数、标准误、t值、P值、置信区间。区别在于系数含义不同——这是第τ个条件分位数的偏效应。还有一个细节值得注意qreg默认使用解析标准误基于稀疏性矩阵估计它假设残差密度函数在分位点附近是平滑的。如果样本量不够大、或者分布形态怪异这个标准误可能偏得厉害。稳妥的做法是结合bootstrap标准误对比看。从实操手感来说qreg最让人满意的一点是速度——即使面对几万条数据单纯形法也稳得住。但跑大型调查数据几十万观测时建议用qreg的替代命令sqreg配合bootstrap或者优化数据量后再跑不然单纯形法的迭代次数会拉满等待时间感人。3.2 核心命令sqreg与bsqreg标准误的正解qreg一次只能估计一个分位点而且标准误计算方式单一。实际研究中更常用的是sqreg——它一次可以跑多个分位点并且使用bootstrap获得标准误和分位点之间的协方差矩阵。协方差矩阵的意义在于你可以检验不同分位点的系数是否存在显著差异这是分位数回归异质性分析的基石。sqreg的标准语法* 同时估计5个分位点的回归bootstrap重复200次 sqreg price weight mpg, quantile(0.1 0.25 0.5 0.75 0.9) reps(200)这个命令跑出来后你会看到5个分位点各自的回归结果依次排列。和qreg不同sqreg的置信区间基于bootstrap因此对数据分布形态的假设更少结果更稳健。reps(200)是bootstrap的重复次数官方推荐至少100但我个人习惯设置200到500。在最终论文里我会用reps(500)保证结果稳定可复现因为bootstrap的随机性会导致每次跑出来的标准误略有不同重复次数越多波动越小。bsqreg是sqreg的单分位点版本它只做一个分位点但使用bootstrap标准误。如果你只需要报告一个分位点的结果bsqreg比sqreg更轻量、更快。它的语法是bsqreg price weight mpg, quantile(0.5) reps(200)还有一个容易被忽略的命令是iqreg它估计一个分位数区间。用法是* 估计25%分位与50%分位之间的系数区间 iqreg price weight mpg, quantile(0.25 0.5)输出的结果是两个分位点之间系数差异的估计和置信区间。如果你想直接回答“0.25分位和0.5分位的效应是否相同”iqreg给你的就是最直接的回答连test命令都省了。3.3 多分位数联合估计与系数差异检验的具体操作假设我要研究“汽车重量对价格的影响在不同价位段是否有差异”完整操作流程如下sysuse auto, clear * 第一步画价格分布直方图确认是否需要分位数回归 histogram price, frequency normal * 第二步sqreg同时估计5个分位点 sqreg price weight mpg, quantile(0.1 0.25 0.5 0.75 0.9) reps(200) * 第三步查看完整结果 estat bootstrap * 第四步检验weight在0.1分位点和0.9分位点的系数是否相等 test [q10]weight [q90]weight这里test命令的写法要注意sqreg估计后每个分位点的结果都存放在对应的方程里方程名默认是q10、q25、q50、q75、q90对应你指定的分位点。如果我把quantile选项里的分位点改成0.2、0.4、0.6那么方程名就变成q20、q40、q60这种命名规律一定要记住不然test命令会找不到对象。检验结果如果P值小于0.05就说明weight在0.1分位和0.9分位的系数存在显著差异异质性效应统计上成立。反之如果P值大于0.1说明两个分位点的系数差异在统计上与零无异——即使点估计值看起来天差地别你也不能脑补出一个“显著差异”的故事。从我个人经验来看这个检验是分位数回归论文中最容易被审稿人追问的地方。很多作者只展示不同分位点的系数表却从来不做跨分位点的差异检验这等于讲了一个没有统计推断支撑的故事。把这个检验加进去论文的严谨度立刻提升一个档次。4. 结果的可视化让分位数回归的结论一目了然4.1 系数随分位数变化的图形grqreg与coefplot表格可以精确地报告数字但要让读者快速抓住“系数随分位点变化”的趋势图形比表格高效得多。Stata里最常用的画图工具是grqreg和coefplot两个都是外部命令需要安装ssc install grqreg ssc install coefplotgrqreg专门为分位数回归设计用法非常简单。以之前的auto数据为例如果你已经用sqreg估完了5个分位点直接运行grqreg weight, ci就能得到一张以分位点为横轴、weight系数为纵轴的折线图图形中还会画出置信区间的阴影带。如果系数随分位点上升而上升图形就是一条向上倾斜的线这种视觉冲击力是表格给不了的。coefplot的功能更通用一些它可以把多次估计的系数画在一起对比。如果你想展示不同分位点上所有变量的系数变化可以循环估计后用coefplot把结果叠起来quietly sqreg price weight mpg, quantile(0.1 0.25 0.5 0.75 0.9) reps(200) coefplot, keep(weight mpg) /// vertical yline(0) /// xlabel(1 { tau}0.1 2 { tau}0.25 3 { tau}0.5 4 { tau}0.75 5 { tau}0.9) /// title(Coefficients across quantiles)coefplot输出结果中每个变量在不同分位点的系数估计值和置信区间会以“系数点置信区间线”的形式排在同一张图中一眼就能看出哪个变量的效应在哪个分位段显著、哪个不显著。唯一要注意的是coefplot对中文标题的支持在旧版本中不太好建议直接用英文标题或者设置好中文字体后再绘图。4.2 从图形中读出的三个核心信息看分位数回归图的时候有三个信息点值得重点留意。第一中位数τ0.5的系数和OLS系数的差异。如果两者差别很大说明数据分布偏态严重均值回归的表征能力存疑这时候用分位数回归就不只是“补充分析”而是“必要分析”。如果两者非常接近说明数据分布相对对称分位数回归的价值更多体现在对尾部异质性的挖掘上。第二系数曲线的形状和置信带宽窄。曲线单调上升意味着x对y的边际效应随分位点提高而增强曲线呈U型意味着两端的效应高于中间。同时要注意置信带的宽度。通常来说中位数附近的置信带最窄两端的置信带会变宽这是正常的统计现象——因为极高分位和极低分位的估计本身就有更大的不确定性。第三看系数曲线是否穿越零线。如果weight在0.1分位的系数置信区间包含0而在0.9分位的系数置信区间不包含0那么这个变量的效应就存在“从无到有”的转变这是非常有故事性的研究发现写作时值得浓墨重彩。我常跟人强调一个习惯画好图之后不要急着截图放进论文先盯着图看30秒问自己三个问题——变量的效应格局是否如预期有没有哪个分位点的结果非常反直觉置信区间是否过宽导致结论不具备说服力这套自问流程帮我避免过好多次“先入为主”的解读错误。5. 真实案例复盘用分位数回归做一份完整的实证分析5.1 案例背景、变量选择与模型设定为了把上面的方法串起来我用一个经典的虚构案例来演示研究“教育年限educ和工作经验exper对个人收入wage的影响”这是劳动经济学里最常见的模型。数据设定为1000个观测的截面数据wage呈现明显的右偏分布少数高收入者拉高了均值这正是分位数回归的用武之地。模型的基本设定wage β0(τ) β1(τ)×educ β2(τ)×exper ετ这里β1(τ)表示教育年限在第τ个条件分位点上的边际效应。我的研究问题是教育回报率在不同收入层次的人群中是否存在差异如果有差异有多大数据处理阶段我需要先构造变量。如果原始数据里没有现成的工资变量而只是报告了收入等级区间那需要用一些方法转换。Stata里需要注意分位数回归对因变量的要求是连续变量或至少是有序数值变量类别变量如“高收入/中收入/低收入”的三个分类标签不适合直接作为因变量。如果遇到收入取对数后分布仍然偏态严重的情况可以考虑同时跑“对数收入”和“收入水平”两个模型看结论是否一致这叫稳健性检验。5.2 模型估计与结果解读一个教科书式的输出跑分位数回归的命令如下* 先跑ln(wage)的OLS作为基准 reg ln_wage educ exper * 再跑5个分位点的sqreg sqreg ln_wage educ exper, quantile(0.1 0.25 0.5 0.75 0.9) reps(500) * 检验教育年限系数的跨分位差异 test [q10]educ [q90]educ典型的结果大概是OLS给出的educ系数大约在0.09也就是教育年限每增加一年收入平均上升9%。但分位数回归的结果是educ在0.1分位的系数只有0.05在0.9分位的系数高达0.13。这意味着什么教育回报率在整个收入分布中增长强劲——最底层10%的人多读一年书带来收入增长5%最顶层10%的人多读一年书带来收入增长13%。用通俗的话讲书读得越多收入越高而且对本来就收入高的人教育的“溢价”更大。这个结果如果通过了test检验P值小于0.01就是一篇实证论文里非常有价值的“核心发现”。它告诉我们平均效应9%实际上掩盖了两个极端群体的巨大差异政策上若要缩小收入差距光靠普及教育可能不够还需要配套措施帮助低收入群体把教育转化为实际收入增长。5.3 论文汇报规范系数表、附图与稳健性检验的呈现把分位数回归结果写进论文有几个行文规范我建议严格遵守。第一主回归表建议用“OLS结果多分位点分位数回归结果”合并的形式。常见的做法是第一列放OLS系数后面几列分别放0.1、0.25、0.5、0.75、0.9分位点的系数。每行对应一个变量每列顶部标注分位点。标准误一律放在括号里星号标注显著水平。第二图形要配合表格使用。表格报告精确数字图形展示趋势。如果杂志社对图片数量有限制优先保留“系数随分位点变化”的那张核心图比如grqreg画出来的因为它最直观地证明了异质性效应的存在。第三稳健性检验不能少。我通常会做三个方向的稳健性处理一是改变分位点的选取比如改用0.05/0.15/0.5/0.85/0.95看结论方向是否一致二是换标准误的计算方式比如从bootstrap换成解析标准误或者增加reps次数到1000看显著性是否稳定三是增加控制变量或换核心解释变量的度量方式比如把educ换成最高学历虚拟变量看结论是否改变。这三个方向都经得起考验论文的说服力就上去了。6. 常见问题与排查技巧实录6.1 报错与异常结果的常见原因实战中我遇到过不少报错和怪象这里盘点几个最高频的附上排查思路错误消息“no observations”。分位数回归对每个分位点的观测数量有要求如果样本量太小或者某些变量的缺失严重Stata会自动剔除缺失观测出现“no observations”提示。解决办法先跑describe和misstable summarize检查缺失值用drop if missing或显式地处理缺失值确保有效样本量充足。估计结果出现“coefficient not estimable”。这个常见于自变量之间存在完全共线性或者某个变量在某个分位点附近取值过于稀疏。排查方式用estat vif检查多重共线性删除共线性强的变量如果变量本身是虚拟变量且某个类别占比极低如少于1%考虑合并类别。标准误异常大bootstrap结果每次跑都不一样。这多半是reps(200)这种小重复次数加上数据本身方差大导致的。解决方法是增加reps我一般至少500或者改用“bca”置信区间。Stata的sqreg支持bca选项但要注意bca区间计算在极端分位点可能不稳定需要结合数据情况判断。极端分位点如0.05、0.95的系数置信区间宽得离谱。这不一定是错误而是数据稀疏导致的正常现象。但如果你发现0.95分位的拟合值出现“阶梯状跳跃”即相邻观测值之间的系数突然跳变这通常意味着在那个分位点附近可用的数据点太少优化算法难以稳定求解。这时候与其报一个不可靠的极端分位数结果不如在论文里只报告0.1到0.9范围的结果并在脚注里说明理由。我自己的经验是遇到任何异常输出第一步永远是“回到数据”——用summarize、histogram、centile查看分布形态用scatter查看变量关系绝大多数问题都是数据问题而不是方法问题。方法本身不会骗人数据会。6.2 分位数回归的一些独门避坑经验除了报错我还想分享几条浅一点但非常实用的经验这些都是靠时间和踩坑换出来的。bootstrap随机数种子一定要设置。sqreg和bsqreg的bootstrap过程依赖随机数如果不事先用set seed设置种子那么你每次跑出来的标准误都会略有不同读者以及你自己将无法精确保存和复现结果。我通常固定为set seed 12345并在论文的方法部分注明种子值。这一条看似微不足道但在结果复核时能省掉无数麻烦。分位点不是越多越好。有人为了“全面”一下子跑0.01到0.99每隔0.01一个分位点共99个回归结果图形乱成一团估计也不稳定。我自己的经验是常规研究报告5到9个分位点就够了比如0.1、0.25、0.5、0.75、0.9或者再加密到0.05到0.95既能反映分布全貌又能保证每个分位点的估计相对可靠。面板数据能不能用分位数回归能但不要直接用qreg。面板分位数回归需要专门的命令和模型设定比如xtqreg需要ssc install xtqreg或者广义分位数回归genqreg命令。如果直接把面板数据当混合截面跑qreg固定效应的控制就会缺失估计结果有偏。这是我的一个血泪教训早期做面板数据时偷懒用了qreg审稿人一眼看穿直接打回。分位数回归样本量要求比OLS更高。因为分位数估计是“局部”的每个分位点实际利用的有效信息量小于全样本。如果样本量只有一两百跑出来的结果可能极不稳定。我建议少于500观测的数据尽量少报极端分位点如0.05和0.95如果非要报加上bootstrap置信区间并明确提醒读者谨慎解读。末尾再分享一个小技法。grqreg画出来的图默认是连续折线式的但很多人喜欢用“点置信区间”的方式呈现——这时候可以用coefplot实现。我自己从踩坑中总结出一个个人偏好先用sqreg跑完5个分位点再跑2到3个额外的分位点作为稳健性参考然后用grqreg画一张主图用coefplot画一张全变量对比图最后在附录里放回归系数完整表格。这套组合拳打下来审稿人一般不会再对分位数回归部分提出“分析不够深入”的意见。研究这东西方法不怕基础怕的是用了方法却讲不清楚结论——分位数回归恰好是一个让你“讲清楚”的绝佳工具。