★ FEATURED ARTICLE 补贴100万和补贴50万企业的减排效果会差多少我前阵子被一个做环境政策评估的朋友问住。他手里有几百家企业的补贴数据想用倾向得分匹配PSM跑因果效应结果卡住了——PSM要求处理变量是0/1而他的处理变量是连续的补贴金额把金额二分成“补贴/不补贴”又太浪费信息连补贴强度这个核心维度都丢了。这个问题很典型补贴强度、培训时长、研发投入、药物剂量政策世界里大量处理变量是连续的PSM这套标准工具根本接不住。这也正是广义倾向得分匹配GPSM要解决的问题。GPSM把经典PSM从“是否处理”推广到“处理多少”估计的是一条剂量响应曲线Dose-Response Function而不是一个单一的平均处理效应。在Stata里跑GPSM有不少实现路径也有不少坑。这篇就把我自己的使用经验完整写出来从原理到实操到避坑一次讲透。1. 补贴强度、培训时长这类“连续处理”PSM为什么接不住1.1 PSM的设计前提是二值处理这是它的底层边界经典PSM的定义起点是一个二值处理变量DD1表示参与政策D0表示未参与。倾向得分定义为e(X)P(D1|X)然后用这个得分在参与组和未参与组之间做匹配、分层或加权消除自选择偏差。这个框架在处理“是否参与”这类问题时非常成熟但它从骨子里依赖二值处理。你可以尝试把连续处理变量强行二分比如“补贴超过50万算1否则算0”但这么做的代价很大一是丢掉了强度差异里的增量信息50万和200万在减排效果上可能完全不同一刀切之后这种差异被并入误差项二是分组的临界值本身就是主观的换一个切点估计结果可能完全改变稳健性很难说服人。1.2 连续处理场景的三个核心难点当处理变量T是连续型时政策评估直接面临三个问题。第一个问题是反事实的维度暴增。二值处理下每个样本只有两个潜在结果Y(0)和Y(1)而连续处理下T可以取任意值理论上需要考虑Y(t)在整个取值区间上的变化。研究者想要的不是一个ATE而是整个剂量响应曲线当T从t变成t1时Y平均变化多少。第二个问题是自选择偏差更隐蔽。企业拿到的补贴金额不是随机分配污染强度大、清理成本高的企业可能倾向于申请更多补贴也可能符合更高额度的条件。直接画出“补贴金额—减排量”散点图拟合回归估计量是有偏的因为补贴金额和企业潜在减排能力相关。第三个问题是匹配逻辑不再成立。经典PSM可以在二值状态下找“得分相近”的对子但连续处理下每个样本的“处理组”和“对照组”边界模糊很难定义什么是精确匹配。GPSM的处理思路不是匹配而是利用广义倾向得分构造权重和回归调整这也就是它被称为“match”但又不完全等价于传统匹配的原因。1.3 GPSM具体能回答什么问题什么场景适合用GPSM直接估计的是平均剂量响应函数ADRF, Average Dose-Response Function。它回答的问题是如果所有人都在处理强度t下平均结果会是多大这个比ATE信息量大得多政策上尤其有用。补贴预算有限时决策者想知道的是“边际补贴增量还能带来多少减排回报”而不是“补贴与否的差距”。适合用GPSM的典型场景包括政府补贴金额对就业或产出的影响、培训时长对工资提升的作用、研发投入强度对创新产出的贡献、药物剂量对疗效的剂量反应关系。一句话处理变量必须是连续或至少有序多值且存在与协变量相关的自选择机制这才是GPSM的用武之地。2. GPSM估计框架拆解Hirano-Imbens三步法是怎么运转的2.1 广义倾向得分到底是什么定义GPSM的核心概念由Hirano和Imbens在2004年提出控制在统计学界的论文里算是非常可读的一篇。广义倾向得分被定义为在给定协变量X的条件下处理变量T取某个值t时的条件密度函数。用公式写就是G g(t, X) f_{T|X}(t | X)这个定义比传统倾向得分更灵活传统PSM的倾向得分是条件概率属于密度函数的一个特殊退化形式。GPS把一个概率推广成了密度自然就从二值处理扩展到了连续处理。GPSM依赖的识别假设叫弱无混淆假设weak unconfoundedness比传统PSM的强无混淆假设还要宽松一些。它的含义是在控制协变量X的前提下潜在结果Y(t)与处理变量T相互独立即Y(t) ⊥ T | X对任意t都成立。换成大白话所有会导致“处理强度不同”的混淆因素都被我们观测到了不存在遗漏变量。这是一个很强但无法直接检验的假设实证研究能做的就是尽量收集协变量并在论文中做敏感性讨论。GPS最漂亮的性质是它的均衡性给定GPS的取值后样本的处理强度与协变量无关。这和传统PSM的均衡性质完全平行。正是这个性质让我们可以在给定GPS的条件下对处理变量和结果做条件调整剥掉协变量带来的偏差。2.2 三步估计每一步的具体设定Hirano-Imbens框架的标准估计分三步。第一步估计GPS。实际操作中几乎都用参数模型。最常见的是假设处理变量的条件分布为正态T | X ~ N(Xβ, σ²)用普通最小二乘回归得到β和σ²的估计后对每个样本i在观测到的T_i处计算密度值G_i normalden(T_i, X_iβ̂, σ̂)注意这里的密度是在观测到的T_i处计算的每个样本只有一个GPS值。这和后面第三步在假设剂量t处重新计算密度值不是一回事初学者经常在这里绕晕。第二步建立结果变量的条件期望模型。用处理变量T和GPS记作G作为解释变量对结果Y做回归。Hirano-Imbens建议使用灵活的函数形式最常用的是二阶多项式加交叉项E[Y|Tt, Gg] α₀ α₁t α₂t² α₃g α₄g² α₅t·g为什么选多项式因为我们不知道真实的剂量响应函数长什么样多项式本质上是在做一个局部逼近用尽可能灵活的设定去接近未知函数。二次加交叉项已经可以捕捉U形、倒U形、边际效应随剂量变化这类非线性关系。项数不是越多越好后面实操部分我会讲怎么权衡。第三步计算ADRF。选择一个剂量网格t₁, t₂, ..., t_K覆盖处理变量的实际取值范围。对每个网格点t用第一步的模型重新计算每个样本在这个t处的GPS值ĝ(t, X_i)。再将这个值代入第二步的回归方程得到每个样本在t处的预测结果然后取平均μ(t) (1/N) Σ_i [α̂₀ α̂₁t α̂₂t² α̂₃ ĝ(t,X_i) α̂₄ ĝ(t,X_i)² α̂₅ t·ĝ(t,X_i)]这一步的意义是把每个样本都“设定”在相同的处理强度t下然后平均预测结果。所有样本统一放在t上直接用回归方程插值就是剂量响应的基本逻辑。2.3 为什么不能用普通回归直接估计必须走GPS这条路如果不用GPSM直接做Y对T的回归得到的系数是条件相关而不是因果效应。原因还是自选择偏差T取值高的样本可能本身X特征就不同这些X同时影响T和Y造成混杂。GPSM的逻辑是两步分离第一步把T对X的依赖关系显式建模得到GPS第二步在回归方程中把GPS也作为控制变量放进去相当于在给定T和GPS的条件下考察Y的变化。GPS作为协变量X的一个充分总结把混淆因素“压缩”成一个标量避免了高维X带来的稀疏性问题。这是GPSM在有限样本下比直接高维回归更实用的原因。用生活类比的话传统PSM像给两个条件相似的人配对比效果GPSM则更像把每个人在不同“剂量”下的预期结果全部拟合成一条曲线。这条曲线上的每个点都是全体样本在统一剂量下的平均水平所以才能谈“因果”。3. Stata实现的三条路径官方命令、用户命令与手写代码3.1 路径一Stata 18的官方doseresponse命令如果你的Stata版本在18及以上直接有官方命令doseresponse可以用。这是Stata官方实现的连续处理变量因果推断命令最大的优势是语法规范、文档齐全、结果输出干净适合正式论文使用。基本调用框架大致是doseresponse Y T X1 X2 X3, effect(response) dose_model(linear)具体选项在help doseresponse里有完整说明。这个命令的边界是必须有比较新版本的Stata版本不够的话只能往下看路径二和路径三。如果你的Stata版本还没有这个命令不要硬写直接走路径二。3.2 路径二Bia-Flores的gpscore与doseresponse命令GPSM在Stata社区里最经典的实现是Michela Bia和Carlos A. Flores开发的两个用户命令gpscore和doseresponse。它们通过SSC安装命令如下ssc install gpscore ssc install doseresponsegpscore负责第一步和平衡性检验doseresponse负责第二步和第三步的剂量响应估计。大致用法是gpscore T X1 X2 X3, gpscore(gscore) predict(hat) sigma(sig) cutpoints(3) doseresponse Y T, gpscore(gscore) dose(dtreat) sigma(sig) cutpoints(3) index(normal)cutpoints(3)选项会在GPS的分位数上分3层做平衡性检验。这两个命令的优点是针对GPSM专门开发检验功能也比较完整缺点是命令的选项有一些历史包袱返回结果的变量名需要看清楚再使用。3.3 路径三手写完整代码彻底搞明白每一步我个人最推荐初学者走的路径是第三步——手工实现。你不需要特别的命令包只用Stata最基础的regress、predict和normalden就能把整个流程跑通。手写一遍的好处是你对每一步的数学含义都会有肌肉记忆后面用任何现成命令都不会心虚。手写全流程的代码骨架如下。假设处理变量是T结果变量是Y协变量是X1到X4。* 第一步对处理变量做标准化提升数值稳定性 summarize T scalar mu_T r(mean) scalar sd_T r(sd) gen T_std (T - mu_T) / sd_T * 第二步估计GPS模型保存均方根误差和线性预测值 regress T_std X1 X2 X3 X4 scalar sigma_gps e(rmse) predict T_hat, xb gen G normalden(T_std, T_hat, sigma_gps) * 第三步剂量响应回归二阶多项式 交叉项 gen T2 T_std^2 gen G2 G^2 gen TG T_std * G regress Y T_std T2 G G2 TG * 保存回归系数 scalar b0 _b[_cons] scalar b1 _b[T_std] scalar b2 _b[T2] scalar b3 _b[G] scalar b4 _b[G2] scalar b5 _b[TG] * 第四步在剂量网格上计算ADRF summarize T_std local tmin r(min) local tmax r(max) local K 30 matrix results J(K, 2, .) forvalues k 1/K { local tval tmin (tmax - tmin) * (k - 1) / (K - 1) gen double G_t normalden(tval, T_hat, sigma_gps) gen double yhat_t b0 b1*tval b2*tval^2 /// b3*G_t b4*G_t^2 b5*tval*G_t quietly summarize yhat_t matrix results[k, 1] tval matrix results[k, 2] r(mean) drop G_t yhat_t } * 第五步把矩阵转为变量画剂量响应曲线 svmat results, names(dose_adrf) twoway line dose_adrf2 dose_adrf1, sort /// ytitle(预测结果) xtitle(标准化处理强度)这个代码里最需要注意的地方是循环内部重新计算G_t。很多初学者会误以为直接用第一步算出来的那一个GPS值就行其实不对。第三步的ADRF公式里对每一个候选剂量t都需要用第一步的模型重新算出“如果这个样本对应Tt它的GPS会是多少”。可以这样理解G是每个样本在自身实际T位置上的密度G_t则是全体样本统一移动到t位置上之后的密度。三种路径怎么选我做了个实际对比实现路径门槛灵活性适合场景官方doseresponseStata 18中正式研究、需要规范输出gpscore/doseresponse需联网安装中需要自带平衡性检验时手写代码无额外依赖高教学、理解原理、复杂设定扩展4. 从数据到结论一个补贴减排案例的完整Stata演示4.1 模拟数据的设计逻辑为了演示完整流程我构造一个模拟数据集。背景设定某地政府给企业发放减排补贴单位百万元我们关心的结果是企业当年单位产值的碳排放下降幅度单位%。数据生成过程中补贴金额与企业的污染强度、规模、利润率、是否位于重点控制区相关这是为了模拟真实的“选择性补贴”机制。同时真实的减排效果与补贴强度之间存在非线性关系补贴较低时边际减排效果递增补贴过高时边际递减。clear all set obs 600 set seed 20240826 * 协变量 gen size rnormal(50, 15) // 企业产值规模百万元 gen pollution rnormal(0.30, 0.08) // 行业污染强度 gen zone rbinomial(1, 0.35) // 是否位于重点控制区 gen roa rnormal(0.08, 0.03) // 利润率 * 连续处理变量补贴金额百万元受协变量影响存在自选择 gen subsidy 6 0.12*size - 1.5*pollution 0.8*zone 15*roa rnormal(0, 2.5) replace subsidy 1.5 if subsidy 1.5 * 结果变量碳排放下降幅度%对补贴的真实响应为非线性 gen co2_reduction 8 2.1*subsidy - 0.08*subsidy^2 /// - 0.3*pollution*subsidy 0.4*roa*subsidy rnormal(0, 2.5) * 标准化 summarize subsidy gen subsidy_std (subsidy - r(mean)) / r(sd)这个DGP里加入了污染强度与补贴的交互项意思是同样一笔补贴污染强度高的企业减排空间大、边际效果强。如果不用GPSM而直接回归这个交互项带来的异质效应会被协变量相关关系搞混。4.2 三步估计的Stata输出解读第一步估计GPS模型并生成GPSregress subsidy_std size pollution zone roa scalar sigma_gps e(rmse) predict subsidy_hat, xb gen G normalden(subsidy_std, subsidy_hat, sigma_gps)回归输出里size、pollution、zone、roa的系数基本都很显著说明补贴金额确实与企业特征相关自选择机制存在直接用OLS回归y会有偏。这一步的目的不是解释而是拿到准确的GPS。第二步剂量响应回归gen subsidy2 subsidy_std^2 gen G2 G^2 gen subG subsidy_std * G regress co2_reduction subsidy_std subsidy2 G G2 subG输出中subsidy_std和subsidy2的系数分别反映了剂量响应的线性项和曲率项。如果subsidy2的系数为负值且显著说明存在边际递减效应剂量响应曲线是倒U形。本例中应该能看到显著的负二次项。第三步计算ADRF并画图scalar b0 _b[_cons] scalar b1 _b[subsidy_std] scalar b2 _b[subsidy2] scalar b3 _b[G] scalar b4 _b[G2] scalar b5 _b[subG] summarize subsidy_std local tmin r(min) local tmax r(max) local K 30 matrix adrf_out J(K, 2, .) forvalues k 1/K { local tval tmin (tmax - tmin) * (k - 1) / (K - 1) gen double G_t normalden(tval, subsidy_hat, sigma_gps) gen double yhat_t b0 b1*tval b2*tval^2 /// b3*G_t b4*G_t^2 b5*tval*G_t quietly summarize yhat_t matrix adrf_out[k, 1] tval matrix adrf_out[k, 2] r(mean) drop G_t yhat_t } svmat adrf_out, names(dose) twoway line dose2 dose1, sort /// ytitle(预测碳减排幅度 (%)) xtitle(标准化补贴强度) /// title(补贴强度-减排效果的剂量响应曲线)4.3 结果解读的正确姿势画出来的剂量响应曲线就是GPSM的核心输出。横轴是标准化的补贴强度纵轴是调整协变量后预测的平均减排幅度。曲线形状可能和直接回归的拟合线有明显差异差异的来源正是GPS调整了自选择偏差。读ADRF曲线要看三点。第一整体水平曲线各点对应的Y均值就是“如果全样本都接受该补贴强度”时的平均减排。第二边际效应曲线的局部斜率随着t的变化斜率增长还是下降对应边际递增还是递减。第三最优区间如果曲线有峰值峰值对应的补贴强度就是政策上的最优参考点这在实际补贴政策制定中最实用。在论文或报告中ADRF曲线通常配合置信区间一起展示。手写代码里我没有做标准误如果想加可以用bootstrap包一层循环对样本有放回抽样后重复整个三步估计把2.5%和97.5%分位数作为置信区间端点。Stata中可以用bootstrap命令把整个流程包装成程序或者直接用gpscore命令自带的bootstrap选项。5. 那些文档里不写的实操细节与避坑经验5.1 先标准化再跑模型这是数值稳定性的关键GPS的计算依赖正态密度函数normalden。如果处理变量的原始尺度很大比如说补贴金额的单位是万元取值范围在几百万到几千万那么变量方差会非常大在均值附近的密度值会趋近于0计算机精度很容易出问题。标准化处理变量和连续型协变量后处理变量标准差为1密度值落在合理区间数值上稳定很多。这也是为什么我在案例里特意把subsidy标准化成subsidy_std。需要说明的是标准化只会影响回归系数的尺度不会影响模型的拟合优度和假设检验结论。ADRF结果的横轴坐标要用标准化数值标注但在解读政策含义时再换算回原始单位这样报告给政策部门看也不难理解。5.2 处理变量的条件分布假设不能无脑用正态GPS的估计依赖于对T | X分布形态的假设。线性回归加正态残差是最常用、最省事的选择但处理变量明显偏态时直接用正态会带来偏误。一个常见做法是先用直方图或分位数图看看处理变量的分布。如果明显右偏可以在第一步使用对数变换gen lnT ln(T) regress lnT X1 X2 X3然后再基于lnT计算GPS。这时GPS对应的是lnT的密度剂量网格也应该在对数尺度上均匀铺开解读结果时需要把横轴换回原始尺度。如果处理变量是取值有限的多分类有序变量比如0、1、2、3级严格的GPSM并不完全适用。更合适的做法是使用有序Probit/Logit模型作为第一步或者退一步用广义倾向得分在暴风算法下的扩展版本。我的建议是处理变量最好是真正连续且分布相对良好的变量这样GPSM的估计结果才可靠。5.3 共同支撑区间第一时间检查而不是最后才检查传统PSM里大家都习惯检查倾向得分的共同支撑区间GPSM同样需要。如果某些样本的GPS极端低说明在它们的协变量特征下几乎观察不到那种处理强度这些样本会对ADRF的预测值产生外推偏差。我常做的一件事是看GPS的分布直方图以及检查协变量取值范围的交集。如果发现GPS尾部有极长尾或者大量接近0的样本常规处理是修剪尾部。一种简单做法是删除GPS低于1%分位数的样本后再跑第二步回归。修剪多少、删除后结果是否稳定都应该在论文的敏感性分析里报告。这也验证了GPSM的结果不能只看一个设定。5.4 多项式阶数和交叉项不是固定答案第二步回归里用几阶多项式是否保留交叉项没有一个统一的正确值。Hirano-Imbens原文用的是二阶多项式加交互项但也有文献用三阶。我的经验是先用二阶加交互作为基准设定再跑一个三阶的扩展作为稳健性检查。怎么比较设定好坏可以直接比较回归的AIC或BIC。如果三阶项的系数不显著且AIC没有下降就说明二阶已经足够。另一个更朴素的标准是看ADRF曲线的形状是否对设定变化敏感。曲线在二阶和三阶设定下保持一致的形状结果的置信度才高。如果形状随设定剧烈变化说明数据对函数形式的依赖太强这时与其继续扩大多项式不如回到第一步检查GPS的模型设定。5.5 平衡性检验不能靠直觉要用分层检验说话GPSM的核心主张是在给定GPS后协变量不再与处理强度相关。这个主张需要验证。Bia-Flores的gpscore命令里cutpoints(3)会按GPS的三分位分层在层内检验处理变量与协变量的相关性并输出显著性检验结果。一个简化的手动检验思路是先把处理变量按分位数分成三层低强度、中强度、高强度再把GPS也分层在每一个“处理强度层 × GPS层”的子块里检验协变量X与处理变量T的秩相关性是否显著。如果大部分子块的相关性不显著就可以认为GPS调整有效。这个检验很反直觉很多人跑完GPS直接就出ADRF忽略了平衡性验证审稿人一问就心虚。把这个结果输出并贴进论文说服力会强很多。5.6 GPSM和传统PSM不是二选一联合使用更有说服力在实测项目里我很少只用GPSM单独出结果。通常的做法是先用二值PSM做一个基准分析确认“是否有政策效果”的方向再用GPSM刻画“政策强度如何调节效果”。前者回答“政策有没有用”后者回答“用多少最合适”两个问题在政策评估里缺一不可。另外GPSM的结果对协变量集合非常敏感。加一个关键协变量ADRF曲线可能整体平移甚至改变形状。我在做正式报告前会对协变量做敏感性分析依次删除每个协变量重新估计看ADRF是否保持稳定。这个做法不需要额外软件手写循环就能完成但能显著提升结果的公信力。5.7 样本量不够时不要碰GPSM最后说一个很少被提及但极其现实的问题样本量。GPSM的三步估计本质上是多个模型叠加对样本量的消耗比传统PSM大得多。系数估计、多项式交叉项、分层平衡性检验每个环节都会消耗自由度。我的底线经验是有效样本量低于500时GPSM的ADRF曲线会非常毛糙置信区间宽到无法得出政策结论。这种情况下不妨考虑将连续处理变量离散化为多分类处理配合传统PSM或IPW做稳健性分析虽然损失了强度信息但至少结果不至于散成一片。我在自己的项目里吃过一次亏只有200多个样本跑GPSM出来的剂量响应曲线呈锯齿状加bootstrap置信区间后几乎覆盖整个纵轴最后只能放弃GPSM降级成三分类处理分析。在那之后我给自己定了个规矩跑GPSM前先看样本量少于500坚决不硬上。
阅读完成 · 觉得有帮助？
