Stata如何做机器学习？——SHAP解释框架下的足球运动员价值驱动因素识别：基于H2O集成学习模型
SHAP解释框架下的足球运动员价值驱动因素识别——基于H2O集成学习模型⚽
欢迎关注 「阿水实证通」，前沿方法时刻看！🌟🌟🌟
文章目录
- SHAP解释框架下的足球运动员价值驱动因素识别——基于H2O集成学习模型⚽
- 聚焦：决策树的足球运动员集合⚽
- 数据探索热身♾️
- 准备H2O开始分析 👍
- 机器学习实战
- 利用变量重要性和SHAP解读洞见
- 终场哨声
- 对咱们实证分析有什么帮助？
- 一、数据预处理：适配经济变量特性与多源整合逻辑
- 二、集成模型与因果推断的协同建模策略
- 三、SHAP解释框架与经济理论的双向映射
- 核心价值：从“预测”到“解释”的经济学研究闭环
聚焦：决策树的足球运动员集合⚽
你是否拥有一个大型数据集？数据中是否存在复杂的、可能是非线性的关系？你是否不确定哪些预测变量最为重要？如果对上述任何一个问题的回答是“是”，那么机器学习可能是处理你的数据的合适选择。
幸运的是，你不必成为程序员也能有效使用这些方法。借助H2O的支持，现在你可以在Stata中使用新的h2oml命令套件进行机器学习。流线型的H2O集成和图形界面使执行用于回归、二元分类和多分类的梯度提升机（GBM）和随机森林（RF）变得轻松。
在这篇聚焦文章中，我们通过分析场上表现指标如何预测足球运动员的平均市场价值——这是球员转会、合同谈判和足球行业整体估值的关键因素——来展示基本的H2O工作流程。我们将训练两个模型进行预测比较，并使用选定的模型解释预测变量的影响。
数据探索热身♾️
我们的分析始于一个从两个不同来源（GitHub)和Kaggle）合并的数据集(https://www.stata.com/stata-news/news40-3/ensemble-decision-trees/fifa.dta)，其中包含2023年以来1000多名不同球员的观测数据。你可以使用以下命令获取这些数据：
. use https://www.stata.com/users/lil/fifa, clear
我们生成箱线图，按位置可视化以百万欧元为单位的市场价值分布：
. graph hbox mil_average_market_value, over(position) title("Market value distribution by position")ytitle("Market value (in millions)")
想知道最右边那个点是谁吗？那是法国前锋基利安·姆巴佩。在本文结尾，我们将探究他为何价值如此之高。
由于存在异常值和右偏分布，我们选择应用对数变换来稳定方差，这可能有助于我们的模型在后续更准确地泛化：
. generate ln_average_market_value = log(average_market_value)
准备H2O开始分析 👍
让我们演示如何从Stata开始使用H2O。首先，h2o init启动一个新的H2O集群，在Stata和H2O之间建立连接。接下来，我们将当前的Stata数据集传输到一个名为fifa的H2O框架中，并使其成为后续操作的当前活动框架：
. h2o init
(output omitted). _h2oframe put, into(fifa) currentProgress (%): 0 100
现在我们已经建立了H2O环境，开始为机器学习准备数据。我们使用_h2oframe toenum将字符串变量转换为分类（枚举）类型。为了验证这种编码是否正确，我们运行_h2oframe describe：
. _h2oframe toenum position nationality league_rank, replace. _h2oframe describeRows: 1108Cols: 21Column Type Missing Zeros +Inf -Inf Cardinalityname string 0 0 0 0
position enum 0 269 0 0 4
age int 0 0 0 0
height int 0 0 0 0
league_rank enum 0 281 0 0 5
average_marke~e real 0 0 0 0
highest_marke~e int 0 0 0 0
total_played_~s int 0 0 0 0
average_minut~d real 0 0 0 0
average_assis~e real 0 419 0 0 0
total_assists int 0 419 0 0 0
assist_per_mi~e real 0 419 0 0 0
average_goals~e real 0 414 0 0 0
total_goals int 0 414 0 0 0
goals_per_min~e real 0 414 0 0 0
total_yellow_~s int 0 163 0 0 0
team_win_ratio real 0 1 0 0 0
data_year int 0 0 0 0
nationality enum 0 6 0 0 73
mil_average_m~e real 0 0 0 0
ln_average_ma~e real 0 0 0 0最后一步是使用_h2oframe split将数据分为训练集和测试集。我们使用标准的训练-测试拆分，分配80%的数据用于训练，保留剩余20%用于测试，并设置随机种子以确保可重复性：
. _h2oframe split fifa, into(train test) split (0.8, 0.2) rseed(19)
机器学习实战
我们准备开始训练！首先，我们切换框架，使训练数据集成为工作框架：
. _h2oframe change train
我们定义全局宏predictors来存储模型将使用的变量，这些变量包括球员人口统计信息、联赛信息、表现指标以及其他行为和球队成功指标：
. global predictors position age height nationality league_rank average_minutes_playedaverage_goals_per_game average_assists_per_game total_yellow_cards team_win_ratio
我们首先实现随机森林回归，使用cv()选项应用三折交叉验证，并使用h2orseed()选项确保可重复性：
. h2oml rfregress ln_average_market_value $predictors, cv(3) h2orseed(19)Progress (%): 0 28.9 100Random forest regression using H2OResponse: ln_average_market_value
Frame: Number of observations:Training: train Training = 877Cross-validation = 877
Cross-validation: Random Number of folds = 3Model parametersNumber of trees = 50actual = 50
Tree depth: Pred. sampling value = -1Input max = 20 Sampling rate = .632min = 17 No. of bins cat. = 1,024avg = 18.9 No. of bins root = 1,024max = 20 No. of bins cont. = 20
Min. obs. leaf split = 1 Min. split thresh. = .00001Metric summaryCross-Metric Training validationDeviance .6068654 .6541131MSE .6068654 .6541131RMSE .7790156 .8087726RMSLE .0480274 .0497248MAE .6133148 .6394611R-squared .6423362 .6144901在这个例子中，我们主要关注均方误差（MSE）作为整体拟合优度的关键指标。使用这个模型，我们通过交叉验证得到的MSE为0.65，作为基准性能水平。让我们看看是否可以通过调优降低这个值。
调优是通过调整超参数来调整模型的过程。为了演示，我们将仅通过比较20到80棵树的模型来调整树的数量。有关RF的可调整选项的完整列表，请参阅 h2oml rf(https://www.stata.com/manuals/h2omlh2omlrf.pdf)。
. h2oml rfregress ln_average_market_value $predictors, cv(3) h2orseed(19)ntrees(20(10)80)Progress (%): 0 100Random forest regression using H2OResponse: ln_average_market_value
Frame: Number of observations:Training: train Training = 877Cross-validation = 877
Cross-validation: Random Number of folds = 3Tuning information for hyperparametersMethod: Cartesian
Metric: DevianceGrid valuesHyperparameters Minimum Maximum SelectedNumber of trees 20 80 70Model parametersNumber of trees = 70actual = 70
Tree depth: Pred. sampling value = -1Input max = 20 Sampling rate = .632min = 17 No. of bins cat. = 1,024avg = 18.9 No. of bins root = 1,024max = 20 No. of bins cont. = 20
Min. obs. leaf split = 1 Min. split thresh. = .00001Metric summaryCross-Metric Training validationDeviance .587266 .6434246MSE .587266 .6434246RMSE .7663328 .8021375RMSLE .0472597 .0493562MAE .6040548 .6320858R-squared .6538873 .6207895我们将MSE适度降低到0.64。在实践中，我们会继续调优以找到交叉验证MSE最低的模型。让我们继续存储这个模型：
. h2omlest store rf
接下来，我们使用相同的三折交叉验证、随机数种子和树网格训练梯度提升回归，并存储此模型。
. h2oml gbregress ln_average_market_value $predictors, cv(3) h2orseed(19)ntrees(20(10)80)Progress (%): 0 100Gradient boosting regression using H2OResponse: ln_average_market_value
Loss: Gaussian
Frame: Number of observations:Training: train Training = 877Cross-validation = 877
Cross-validation: Random Number of folds = 3Tuning information for hyperparametersMethod: Cartesian
Metric: DevianceGrid valuesHyperparameters Minimum Maximum SelectedNumber of trees 20 80 80Model parametersNumber of trees = 80 Learning rate = .1actual = 80 Learning rate decay = 1
Tree depth: Pred. sampling rate = 1Input max = 5 Sampling rate = 1min = 5 No. of bins cat. = 1,024avg = 5.0 No. of bins root = 1,024max = 5 No. of bins cont. = 20
Min. obs. leaf split = 10 Min. split thresh. = .00001Metric summaryCross-Metric Training validationDeviance .1025126 .6337925MSE .1025126 .6337925RMSE .3201759 .7961109RMSLE .0198991 .0492887MAE .23022 .6083482R-squared .9395829 .6264663. h2omlest store gbm
为了帮助我们选择模型，我们使用h2oml事后估计命令在测试样本上评估其预测性能。首先，我们使用h2omlest restore恢复每个模型。然后，我们使用h2omlpostestframe test激活测试数据集进行评估。准备好两个模型后，我们运行h2omlgof来比较它们在测试数据上的性能指标：
. h2omlest restore rf
(results rf are active now). h2omlpostestframe test
(testing frame test is now active for h2oml postestimation). h2omlest restore gbm
(results gbm are active now). h2omlpostestframe test
(testing frame test is now active for h2oml postestimation). h2omlgof rf gbmPerformance metrics for model comparison using H2O
Testing frame: testrf gbmTestingNo. of observations 231 231Deviance .5692276 .5008833MSE .5692276 .5008833RMSE .7544718 .7077311RMSLE .0458549 .0430169MAE .6169705 .567296R-squared .6051586 .6525652GBM在测试集中具有最低的MSE，证明了其优越的预测性能。如果我们有兴趣在测试数据中进行预测，可以在恢复我们选择的模型后使用h2omlpredict：
. h2omlest restore gbm
(results gbm are active now). h2omlpredict gbm_pred_valueProgress (%): 0 100
利用变量重要性和SHAP解读洞见
现在我们已经选择了性能最佳的模型，我们通过使用h2omlgraph varimp绘制变量重要性图来检查其预测变量的影响。在基于树的模型中，变量重要性通过每个预测变量的分裂导致的MSE总减少量来衡量：
. h2omlgraph varimp
变量重要性图确定市场价值的前三个预测变量为team_win_ratio、age和nationality。相比之下，球员position和height等变量的重要性较小。
为了进一步解释，我们现在转向SHAP（夏普利加性解释）值，以了解每个预测变量如何对每个球员的预测市场价值做出贡献。让我们使用h2omlgraph shapvalues查看我们最有价值的球员基利安·姆巴佩的SHAP贡献。他的预测平均市场价值（对数变换后）为18.65，比训练数据集的平均预测值15.59高出3.06个单位。我们使用SHAP值来解释这一差异：
. h2omlgraph shapvalues, obs(590) title("SHAP values for Kylian Mbappé") frame(fifa)
团队胜率（1.3）、平均每场进球数（0.47）和平均出场时间（0.42）共同贡献了从基线到姆巴佩预测市场价值增长的三分之二以上，表明它们是其高估值的关键驱动因素。
最后，为了在样本层面获得洞见，我们可以使用h2omlgraph shapsummary命令生成SHAP汇总图，也称为蜂群图：
. h2omlgraph shapsummary, frame(fifa)Progress (%): 0 100
在这个图中，每个观测值对应每个预测变量有一个点，其水平位置表示该观测值的预测变量的SHAP值。预测变量的高观测值用红色表示，低观测值用蓝色表示。预测变量在y轴上按其整体SHAP重要性排序。
团队胜率是最强的预测变量，较高的胜率始终提升球员的市场价值。年龄呈现相反的模式——年轻增加价值，而年长降低价值，反映了对年轻球员的溢价。
虽然每个球员的个体预测变量贡献不同，但SHAP汇总结果与变量重要性排名大致一致，表明球队成功、年龄和表现指标是球员估值的最有影响力的驱动因素。
终场哨声
通过H2O，Stata的机器学习工具通过对变量重要性的详细评估，既提供了高预测准确性，又提供了有价值的可解释性。
对咱们实证分析有什么帮助？
上文所提到的H2O集成学习与SHAP解释框架的应用，对于实证分析领域，有何指导意义呢？
一、数据预处理：适配经济变量特性与多源整合逻辑
- 分布转换与异常值控制：针对经济学中常见的右偏变量（如收入、资产价格），采用对数变换等方法稳定方差（如文章对球员市场价值的处理），并通过箱线图等可视化工具识别极端值（如球星高估值样本），避免模型被异常点主导。
- 多源数据的结构化整合：经济学研究常涉及宏观-微观跨层数据（如行业指标与企业财务数据），需参照文章合并GitHub与Kaggle数据的思路，重点校验时间频率一致性（如季度数据与年度数据对齐）、量纲统一性（如货币单位标准化），并利用H2O等工具的自动缺失值算法（如均值插补、树模型预测填充）提升数据完整性。
二、集成模型与因果推断的协同建模策略
- 非线性关系捕捉与模型优化：突破线性回归假设，采用梯度提升机（GBM）、随机森林（RF）等集成模型捕捉经济变量的高阶交互效应（如人力资本理论中“年龄-经验-工资”的非线性关系），通过交叉验证（如3折CV）与超参数网格搜索（如树深度、学习率调优）优化预测精度，以均方误差（MSE）、R²等指标作为模型筛选基准。
- 因果机制的双重验证：机器学习模型的预测相关性需与传统计量方法结合：先用集成模型筛选关键变量（如文章中团队胜率、年龄等），再通过工具变量法（IV）、双重差分法（DID）等验证因果关系（如解决“球队胜率与球员能力”的内生性问题），避免将统计关联误判为因果效应。
三、SHAP解释框架与经济理论的双向映射
- 变量重要性的经济机理解构：基于H2O输出的特征重要性排序（如球队胜率＞年龄＞国籍），结合经济学理论解读驱动因素：例如“团队绩效对球员估值的溢出效应”可关联产业组织中的“品牌外部性”理论，“年龄的倒U型影响”对应劳动经济学的“经验溢价衰减”规律。
- 个体与总体层面的边际效应分析：通过SHAP值实现双维度解释：
- 个体层面：如文章中姆巴佩的SHAP贡献图，可用于分析特定经济主体（如龙头企业、高收入群体）的特征影响路径（如胜率+1.3、进球效率+0.47的价值贡献）；
- 总体层面：借助蜂群图分析变量边际效应异质性（如年龄对价值的非线性影响），对应计量经济学中的“异质性处理效应”分析，为政策模拟（如最低工资调整对不同年龄段就业的影响）提供微观基础。
核心价值：从“预测”到“解释”的经济学研究闭环
上述框架通过集成学习提升复杂经济关系的捕捉能力，以SHAP解释连接数据模式与理论机制，最终实现“数据驱动建模—因果机制验证—经济理论创新”的研究闭环，尤其适用于劳动经济学（工资决定）、行为经济学（决策偏差）、金融经济学（资产定价）等领域的高维非线性问题。
