简介这份资源面向劳动经济学、收入分配与机器学习交叉领域的研究生、学者及政策研究者围绕流动人口劳动收入风险的测算及其收入分配效应展开提供论文复现所需的完整MATLAB代码与解释。压缩包内共1个PDF文件约756KB集中呈现数据准备、收入分布复原、方差与偏度及峰度三类风险测算、风险补偿分析、收入分配效应分析与结果可视化等环节便于读者对照代码理解建模思路与复现流程。内容强调风险具有时变特征在不同人口属性和工作特征间差异显著且风险补偿存在时期、收入群体与风险类别上的异质性方差和偏度风险扩大收入差距峰度风险则有助于缩小差距。已有56人学习适合希望借助机器学习方法评估风险对收入分配影响、并获取可重复实证工具的研究者参考。1. 劳动收入风险测算当流动人口的收入差距被机器学习重新切开一个在东莞电子厂干了六年的普工和一个刚来三个月的外卖骑手他们的收入差距到底有多大传统统计口径给出的答案往往是一张平均工资表但真正做过田野调研的人都知道平均数在这里几乎是个笑话。劳动收入风险测算要解决的恰恰是这种“被平均”掩盖的真实分化——它把流动人口的收入拆成可解释的结构性差异和不可解释的随机风险再用机器学习模型去估计不同维度行业、技能、户籍、区域、工龄对收入分布的边际影响。这篇论文复现的核心工具是分位数回归森林Quantile Regression Forest, QRF配合基尼系数分解和反事实模拟回答的是“如果流动人口的技能分布变了收入差距会缩小多少”这类政策问题。适合谁读做劳动经济学实证的研究生、需要复现论文代码的社科量化从业者以及想用机器学习处理收入分配问题的数据分析师。MATLAB 和 Python 都能跑但本文以 Python 生态为主因为 QRF 在scikit-garden和quantile-forest里都有成熟实现调参和可视化比 MATLAB 顺手得多。2. 分位数回归森林为什么比 OLS 和随机森林更适合收入分配研究2.1 从均值回归到条件分位数的范式转换普通最小二乘法OLS估计的是条件均值 $E[Y|X]$它假设收入在给定特征下是对称分布。但劳动收入数据几乎永远是右偏的高收入尾部拖得很长OLS 会把大量信息压缩成一个系数。分位数回归QR由 Koenker 和 Bassett 在 1978 年提出估计的是条件分位数 $Q_\tau(Y|X)$比如 $\tau0.1$ 对应低收入群体$\tau0.9$ 对应高收入群体。问题在于线性分位数回归仍然假设特征与分位数之间是线性关系而流动人口的收入决定机制里工龄对低收入者的边际效应和高收入者完全不同线性假设直接翻车。随机森林RF解决了非线性问题但它输出的是条件均值本质上还是“平均”。分位数回归森林QRF由 Meinshausen 在 2006 年提出思路很直接用随机森林的每棵树给出叶节点内所有训练样本的权重然后加权计算目标分位数。具体来说对于新样本 $x$每棵树 $T_b$ 把它落到某个叶节点 $L_b(x)$该叶节点内的训练样本集为 $S_b(x)$。QRF 的累积分布函数估计为$$\hat{F}(y|Xx) \sum_{i1}^{n} w_i(x) \cdot \mathbf{1}(Y_i \leq y)$$其中权重 $w_i(x) \frac{1}{B} \sum_{b1}^{B} \frac{\mathbf{1}(X_i \in L_b(x))}{|L_b(x)|}$。有了这个 CDF任意分位数都能直接读出来。这个设计的好处是不需要假设分布形式能捕捉异方差而且对异常值比 OLS 稳健得多。2.2 在 Python 里跑通 QRF 的最小命令先装环境。quantile-forest是 scikit-learn 兼容的 QRF 实现比老旧的scikit-garden维护得好Python 3.9 以上直接 pippip install quantile-forest scikit-learn pandas numpy matplotlib生成一份模拟的流动人口收入数据包含教育年限、工龄、行业类别、户籍类型四个特征import numpy as np import pandas as pd from quantile_forest import RandomForestQuantileRegressor from sklearn.model_selection import train_test_split np.random.seed(42) n 5000 # 模拟特征 education np.random.normal(11, 3, n).clip(6, 20) # 教育年限 experience np.random.gamma(4, 2, n).clip(0, 40) # 工龄 industry np.random.choice([0, 1, 2], n, p[0.5, 0.3, 0.2]) # 0劳动密集 1资本密集 2技术密集 hukou np.random.binomial(1, 0.6, n) # 1为本地户籍0为外来 # 收入生成非线性 异方差 base 2000 300 * education 80 * experience industry_effect np.array([0, 800, 2000])[industry] noise_scale 500 50 * education 20 * experience # 异方差 income base industry_effect 500 * hukou np.random.normal(0, noise_scale) income income.clip(1000, None) X pd.DataFrame({education: education, experience: experience, industry: industry, hukou: hukou}) y income X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state42) # 训练 QRF qrf RandomForestQuantileRegressor(n_estimators200, min_samples_leaf10, random_state42) qrf.fit(X_train, y_train) # 预测 10%、50%、90% 分位数 quantiles [0.1, 0.5, 0.9] y_pred qrf.predict(X_test, quantilesquantiles) print(低分位预测均值:, y_pred[:, 0].mean()) print(中分位预测均值:, y_pred[:, 1].mean()) print(高分位预测均值:, y_pred[:, 2].mean())逻辑说明RandomForestQuantileRegressor的predict方法接受quantiles参数返回形状为(n_samples, n_quantiles)的数组。min_samples_leaf10是关键参数叶节点样本太少会导致分位数估计方差爆炸太大则失去异质性捕捉能力。n_estimators200是精度和速度的平衡点实测超过 300 后分位数估计的边际改善很小。参数说明min_samples_leaf建议从 5 到 30 网格搜索用 pinball loss 做交叉验证。max_features默认sqrt对收入数据够用如果特征里有强共线性比如教育年限和技能等级可以降到 0.5。quantiles参数在预测时指定训练时不需要传。2.3 基尼系数分解与反事实模拟的衔接QRF 输出的是个体层面的条件分位数接下来要算基尼系数。基尼系数 $G \frac{2\sum_{i1}^{n} i \cdot y_{(i)}}{n^2 \bar{y}} - \frac{n1}{n}$其中 $y_{(i)}$ 是升序排列的收入。但这里有个坑QRF 预测的是分位数不是点估计所以算基尼系数时要用预测的中位数作为个体收入的代理或者用分位数做积分。常见做法是用 $\tau0.5$ 的预测值然后对反事实样本重新算基尼。反事实模拟的步骤固定其他特征只改变某一个维度比如把所有人的教育年限提高 2 年重新用 QRF 预测分位数再算基尼系数差值就是该维度的边际贡献。这个思路比 Oaxaca-Blinder 分解更灵活因为不依赖线性假设。3. 多维度风险对收入差距的影响从特征工程到分解结果3.1 流动人口收入数据的清洗与特征构造真实数据不会像模拟数据那么干净。流动人口收入调查比如 CMDS 或 CHIP常见的问题收入有零值和极端值、工龄和教育年限有缺失、行业分类口径不一致。清洗步骤我一般按这个顺序来第一步处理缺失值。教育年限缺失超过 15% 的样本直接删工龄缺失用“年龄 - 首次外出务工年龄”补行业缺失用众数填充但加一个缺失指示变量。第二步缩尾处理。收入上下 1% 做 winsorize避免极端值把 QRF 的叶节点权重拉偏。第三步构造交叉特征。流动人口研究里“教育 × 行业”和“工龄 × 户籍”这两个交互项几乎必加因为教育回报率在技术密集行业明显更高而外来户籍的工龄溢价通常更低。# 数据清洗与特征构造 df pd.read_csv(migrant_income.csv) # 缺失值处理 df df.dropna(subset[education]) df[experience] df[experience].fillna(df[age] - df[first_migrant_age]) df[industry_missing] df[industry].isna().astype(int) df[industry] df[industry].fillna(df[industry].mode()[0]) # 缩尾 for col in [income]: lower, upper df[col].quantile([0.01, 0.99]) df[col] df[col].clip(lower, upper) # 交叉特征 df[edu_industry] df[education] * df[industry] df[exp_hukou] df[experience] * df[hukou] # 对数收入QRF 对收入水平做分解时用对数更稳 df[log_income] np.log(df[income] 1)逻辑说明industry_missing指示变量很重要因为缺失本身可能携带信息比如非正规就业者更可能不报行业。缩尾用分位数而不是标准差因为收入分布太偏。对数变换在分解基尼系数时不是必须的但能让分位数曲线的形状更平滑减少极端分位数的估计噪声。参数说明缩尾比例 1% 是文献惯例如果样本量小于 2000可以放宽到 0.5%。交叉项是否显著用 QRF 的变量重要性permutation importance验证不要凭感觉加。3.2 用 QRF 做反事实分解的完整代码假设我们要分解教育、工龄、行业、户籍四个维度对收入差距的贡献。核心思路对每个维度生成一个“反事实样本”即该维度取值替换为总体均值或随机重排其他维度不变然后比较原始基尼和反事实基尼。from sklearn.inspection import permutation_importance from quantile_forest import RandomForestQuantileRegressor def gini(y): y np.sort(y) n len(y) cum np.cumsum(y) return (2 * np.sum((np.arange(1, n1) * y)) / (n * np.sum(y))) - (n 1) / n # 训练最终模型 features [education, experience, industry, hukou, edu_industry, exp_hukou] qrf_final RandomForestQuantileRegressor(n_estimators300, min_samples_leaf15, random_state42) qrf_final.fit(df[features], df[log_income]) # 原始基尼 y_pred_median qrf_final.predict(df[features], quantiles0.5).flatten() gini_original gini(np.exp(y_pred_median)) print(f原始基尼系数: {gini_original:.4f}) # 反事实分解 results {} for feat in features: df_cf df.copy() df_cf[feat] np.random.permutation(df_cf[feat].values) # 随机重排该维度 y_cf qrf_final.predict(df_cf[features], quantiles0.5).flatten() gini_cf gini(np.exp(y_cf)) results[feat] gini_original - gini_cf # 贡献度 print(f{feat} 的基尼贡献: {results[feat]:.4f}) # 变量重要性交叉验证 perm_imp permutation_importance(qrf_final, df[features], df[log_income], n_repeats10, random_state42, scoringneg_mean_pinball_loss) for i, feat in enumerate(features): print(f{feat}: {perm_imp.importances_mean[i]:.4f} ± {perm_imp.importances_std[i]:.4f})逻辑说明随机重排permutation是反事实分解的常用手段它打破了该维度与收入的关联但保留了其他维度的联合分布。gini_original - gini_cf为正说明该维度扩大了收入差距为负说明该维度缩小了差距。permutation_importance用 pinball loss 作为评分比默认的 MSE 更适合分位数模型。参数说明n_repeats10是稳定性下限如果结果波动大加到 30。min_samples_leaf15在 5000 样本量下是经验值样本量翻倍可以降到 10。反事实分解的随机重排要固定random_state否则每次结果不一样论文里没法复现。3.3 分位数视角下的收入差距低分位和高分位的驱动因素不同把 QRF 的分位数预测按 $\tau0.1, 0.5, 0.9$ 分别做反事实分解会发现一个反直觉的结论教育对高分位收入的贡献远大于低分位而户籍对低分位的贡献更大。这意味着提高教育水平主要拉高高收入群体的天花板而打破户籍壁垒更能托底低收入群体。这个发现直接对应政策含义如果目标是缩小收入差距户籍改革比教育扩招更有效如果目标是提升整体收入水平教育投资回报更高。具体操作对每个 $\tau$重复 3.2 的分解流程把结果整理成表格。注意QRF 在极端分位数$\tau0.05$ 或 $\tau0.95$的估计不稳定因为叶节点内样本太少建议只报告 0.1 到 0.9 之间的结果。4. 避坑与排查QRF 做收入分解时最容易翻车的五个地方4.1 叶节点样本太少导致分位数估计震荡现象$\tau0.1$ 的预测值在不同随机种子下波动超过 20%基尼贡献度符号甚至反转。原因min_samples_leaf设得太小比如 1 或 2叶节点内只有几个样本加权 CDF 的阶梯太粗分位数插值不稳定。解决把min_samples_leaf提到 10 以上同时增加n_estimators到 300。如果样本量小于 1000考虑用min_samples_leaf20并降低分位数分辨率只报 0.25、0.5、0.75。4.2 收入零值和负值把对数变换搞崩现象np.log(income)出现-inf或 NaNQRF 训练直接报错。原因流动人口收入数据里自雇者可能报零收入或负收入亏损直接取对数会炸。解决用np.log(income 1)或者np.sign(income) * np.log(np.abs(income) 1)。更稳妥的做法是先检查零值比例如果超过 5%考虑用两部分模型先分类是否就业再对正收入做 QRF。4.3 反事实重排破坏了特征间的联合分布现象反事实基尼系数出现荒谬值比如大于 1 或小于 0或者贡献度加起来远超原始基尼。原因随机重排单个特征时该特征与其他特征的联合分布被破坏QRF 在训练时没见过这种组合预测外推到不合理区域。解决用条件重排conditional permutation即在其他特征的相似样本内重排而不是全局重排。实现上可以用sklearn.neighbors.NearestNeighbors找相似样本或者用 copula 方法保持依赖结构。4.4 行业分类的虚拟变量陷阱现象把行业作为多分类变量直接输入 QRF变量重要性显示行业贡献极低但分组回归里行业效应显著。原因QRF 的 split 是基于阈值的多分类变量被编码成整数后算法会误以为类别之间有顺序关系012导致 split 不合理。解决对行业做 one-hot 编码或者用 target encoding。如果类别数超过 10用category_encoders库的CatBoostEncoder。4.5 基尼系数分解的贡献度不满足可加性现象各维度贡献度之和大于原始基尼系数或者出现负贡献但经济含义解释不通。原因QRF 是非线性模型特征之间存在交互效应反事实分解的贡献度不满足 Shapley 值的可加性。解决不要强行加总把贡献度理解为“边际效应”而非“份额”。如果必须做可加分解改用 Shapley 值方法shap库的TreeExplainer支持 QRF但计算成本高很多。5. 进阶技巧用 SHAP 值验证 QRF 分解结果并做个体层面的风险测算反事实分解给出的是群体层面的贡献度但劳动收入风险测算最终要落到个体这个人的收入风险有多大被低估的概率是多少SHAP 值SHapley Additive exPlanations能把 QRF 的预测分解到每个特征上而且满足可加性正好补上反事实分解的短板。import shap # 用 TreeExplainer 解释 QRF 的中位数预测 explainer shap.TreeExplainer(qrf_final, feature_perturbationinterventional) shap_values explainer.shap_values(df[features].iloc[:500], check_additivityFalse) # 对 0.1 分位数和 0.9 分位数分别解释 # 注意TreeExplainer 对 QRF 的分位数预测需要指定 model_output explainer_low shap.TreeExplainer(qrf_final, datadf[features].iloc[:500], model_outputraw, feature_perturbationinterventional) shap_low explainer_low.shap_values(df[features].iloc[:500], check_additivityFalse) # 可视化低分位群体的特征贡献 shap.summary_plot(shap_low, df[features].iloc[:500], plot_typebar)逻辑说明feature_perturbationinterventional是必须的因为 QRF 不是树模型的标准输出默认的tree_path_dependent会出错。check_additivityFalse是因为 QRF 的分位数预测不是简单加总SHAP 的可加性在这里是近似。model_outputraw确保解释的是原始预测值不是变换后的。参数说明SHAP 计算量随样本量线性增长500 个样本大约 30 秒5000 个样本要 5 分钟以上。如果只关心低分位群体可以先用 QRF 预测 $\tau0.1$筛选出预测值最低的 20% 样本再对这些样本做 SHAP能省 80% 时间。个体层面的风险测算可以这样定义对每个样本计算 $\tau0.1$ 和 $\tau0.9$ 的预测差值作为收入不确定性指标。差值越大说明该个体在收入分布中的位置越不稳定风险越高。然后把这个风险指标作为因变量用 SHAP 值做特征归因就能回答“什么样的流动人口收入风险最高”。我自己的习惯是每次跑完 QRF 分解先用 SHAP 的summary_plot扫一眼特征贡献的排序如果和反事实分解的结论矛盾优先信 SHAP因为它的理论基础更扎实。然后手动检查几个极端样本的 SHAP 力图看看有没有明显的特征组合异常。这个流程跑顺了一篇论文的实证部分基本就稳了。希望帮到你。本文还有配套的精品资源点击获取
