凡是搞计量经济的，都关注这个号了
邮箱：econometrics666@126.com
所有计量经济圈方法论丛的code程序, 宏微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.
前面，推荐了一些Stata高效使用Tips，例如：23年的时候Stata 18出来了，参看：1.Stata18版本对外发布了, 里面提供了30个让学者激动的新程序板块！2.群友在计量社群里分享了一份Stata18, 确实可以使用
25年Stata 19终于出来了，并且详细介绍了 Stata 19 的新功能（参看：Stata19新功能独家揭秘, ML,CS等全面升级, 附详细教程和案例助轻松上手, 绝对独一份指南.）。
下面就详细看看“Machine learning via H2O: Ensemble decision trees”。
Stata 19 推出的全新 h2oml 套件，使得咱们能够利用强大的 H2O 机器学习平台，即使在传统统计模型表现不佳时，也能深入挖掘数据中的复杂模式。该套件专注于集成决策树模型，如梯度提升机 (Gradient boosting machine, GBM) 和随机森林 (Random forest)，支持灵活的分类与回归分析。学者不仅可以方便地调整模型超参数、使用验证集或交叉验证 (CV) 评估性能，还能有效解释模型预测结果，执行更多高级机器学习任务。许多学者习惯使用线性回归、逻辑回归等传统统计模型来分析数据。然而，现实世界的数据关系往往更为复杂，呈现显著的非线性特征。在这些情况下，经典方法可能难以捕捉数据中的精细模式，导致预测能力受限。是否遇到过以下挑战？预测变量间的交互作用复杂，超出了传统模型的表达能力？模型的预测准确率遭遇瓶颈，难以进一步提升？需要模型能有效处理缺失值、多重共线性，并具备良好的泛化能力以预测新数据？希望在追求高预测精度的同时，不牺牲模型的可解释性（即避免预测能力与可解释性之间的矛盾）？这些正是梯度提升机 (GBM) 和随机森林等先进机器学习方法展现优势之处。Stata 19 通过 h2oml 套件，将 H2O 的高性能算法无缝集成到 Stata 环境中，让学者能够直接运用这些强大的预测模型。现在，只需使用熟悉的 Stata 命令语法，即可训练出性能可能超越传统技术的复杂集成学习模型，即使是机器学习新手也能快速上手。更重要的是，h2oml 打破了“机器学习是难以理解的‘黑箱’”这一固有印象。借助套件内置的 Shapley 加性解释 (Shapley additive explanations, SHAP) 值、部分依赖图 (Partial dependence plots, PDP) 和变量重要性排序等强大解释工具，GBM 和随机森林模型在提供卓越预测能力的同时，依然保持了良好的可解释性——无需在预测能力与模型理解之间进行妥协。下面概览下使用 h2oml 的标准机器学习工作流程。如何在 Stata 中使用 H2O 进行机器学习分析：初始化 H2O 集群，并将当前 Stata 内存中的数据集导入 H2O，命名为 dat_h2oframe put, into(data)
将 H2O 中的 data数据框分割为训练集（80%）和验证集（20%），并设置随机种子以确保可重复性：_h2oframe split data, into(train valid) split(0.8 0.2) rseed(19)
将 train数据框设为 H2O 当前活动的数据框，用于后续模型训练：
使用默认超参数，在训练集 (train) 上拟合 GBM 二分类模型，并在验证集 (valid) 上监控性能：h2oml gbbinclass response predictors, h2orseed(19) validframe(valid)
或者，不使用单独验证集，而是通过 3 折交叉验证 (3-fold CV) 来评估模型：h2oml gbbinclass response predictors, h2orseed(19) cv(3)
存储拟合好的基准模型结果，命名为 gbm_ref：
可以直接指定超参数值。例如，设置树的数量为 200，学习率为 0.2：h2oml gbbinclass response predictors, h2orseed(19) cv(3) ntrees(200) lrate(0.2)
存储这个使用自定义超参数的模型，命名为 gbm_user：执行超参数自动调优，搜索树的数量 (20到200，步长10) 和学习率 (0.1到1，步长0.1) 的最佳组合：h2oml gbbinclass response predictors, h2orseed(19) cv(3) ntrees(20(10)200) lrate(0.1(0.1)1) tune()
(注意：tune() 默认激活网格搜索，并使用默认评估指标)
存储调优后找到的最佳模型，命名为 gbm_tuned：自定义调优过程：例如，使用 accuracy (准确率) 作为评估指标，并采用随机网格搜索策略：h2oml gbbinclass response predictors, h2orseed(19) cv(3) ntrees(20(10)200) lrate(0.1(0.1)1) tune(metric(accuracy) grid(random))
绘制得分历史图，检查模型在训练过程中的收敛情况，诊断过拟合或欠拟合：查看交叉验证性能摘要统计（如果使用了 cv() 选项）：查看超参数网格搜索（或随机搜索）的结果摘要，按性能排序：从搜索结果中手动选择某个特定 ID 的模型作为当前活动模型：h2omlselect id = # // 将 # 替换为所需的模型 ID
h2omlgof gbm_tuned gbm_user gbm_ref
假设已完成随机森林模型的训练和调优，并将最佳模型存储为 rf_tuned。现在比较最佳 GBM 模型 (gbm_tuned) 和最佳随机森林模型 (rf_tuned)：h2omlgof gbm_tuned rf_tuned
在同一张图上比较它们的精确率-召回率 (PR) 曲线：h2omlgraph prcurve, models(gbm_tuned rf_tuned)
在同一张图上比较它们的受试者工作特征 (ROC) 曲线：h2omlgraph roc, models(gbm_tuned rf_tuned)
-
首先，恢复之前选定的最佳模型，例如 gbm_tuned：h2omlest restore gbm_tuned
将需要进行预测的数据框（这里以原始的 data数据框为例）设为 H2O 当前活动的数据框：执行预测，并将结果存储在 Stata 当前数据集的一个新变量中（变量名由 h2omlpredict 自动生成或通过选项指定）：h2omlpredict newvarname // newvarname 是你希望存储预测结果的变量名
(可以通过 class 或 probabilities 等选项控制预测类型)
生成部分依赖图 (PDP)，展示一个或多个变量对平均预测结果的边际效应：h2omlgraph pdp predictors // predictors 是一个或多个变量名
生成个体条件期望 (ICE) 图，展示单个变量变化对每个观测预测值的影响，揭示异质性：h2omlgraph ice predictor // predictor 是单个变量名
使用 SHAP 值分析预测变量对个体或整体预测结果的贡献：为单个观测（例如第 # 号观测）生成 SHAP 值贡献图：h2omlgraph shapvalues, obs(#)
生成 SHAP 摘要图（蜂群图），展示全局变量重要性、变量值与贡献方向的关系等：h2omlest: 管理（存储 store、列出 dir、删除 drop、恢复 restore/use`）H2O 模型估计结果h2omlpostestframe: 指定用于后续后估计命令的默认 H2O 数据框架
h2omlestat metrics: 显示当前模型的详细性能指标h2omlgof: 表格化比较多个已存储模型的拟合优度指标h2omlestat cvsummary: 显示交叉验证结果摘要 (若使用 CV)h2omlestat gridsummary: 显示超参数搜索（网格/随机）结果摘要 (若使用 tune())h2omlexplore: 浏览超参数搜索中特定模型的详细结果h2omlselect: 从超参数搜索结果中选择一个模型作为当前活动模型h2omlgraph scorehistory: 绘制模型训练过程中的得分历史图
h2omlestat threshmetric: 显示基于不同分类阈值的性能指标 (Precision, Recall, F1 等)h2omlestat confmatrix: 显示混淆矩阵h2omlgraph prcurve: 绘制精确率-召回率 (PR) 曲线图h2omlgraph roc: 绘制受试者工作特征 (ROC) 曲线图
h2omlpredict: 使用当前活动模型生成预测（连续值、类别概率或类别标签）
h2omlgraph varimp: 生成变量重要性图h2omlgraph pdp: 生成部分依赖图 (PDP)h2omlgraph ice: 生成个体条件期望 (ICE) 图h2omlgraph shapvalues: 为单个观测生成 SHAP 值贡献图h2omlgraph shapsummary: 生成 SHAP 值摘要图 (蜂群图)
h2omltree: 保存集成模型中单棵决策树的结构为 DOT 文件格式，并可显示规则。
为了具体展示 h2oml 套件的功能，本节以一个电信客户流失数据集为例，重点演示如何使用梯度提升机 (GBM) 进行二元分类预测。请注意，使用随机森林模型或进行回归分析的工作流程与此类似。分析一家名为 Telco 的虚构电信公司的数据。该公司在加利福尼亚州提供家庭电话和互联网服务，数据集（由 IBM 提供）包含了 7,043 位客户的 26 个变量信息。核心目标是构建一个预测模型，用于识别哪些客户具有较高的流失风险——即他们可能终止与 Telco 的服务。模型的响应变量是一个二元变量 churn（流失），它标识了客户在上个月是否已停止服务（流失）或仍是活跃用户。人口统计信息: 如性别、是否为老年人 (seniorcitizen)、是否有伴侣 (partner)、是否有家属 (dependents) 等。账户详细信息: 如合同类型 (contract)、在网时长 (tenuremonths)、账单类型 (paperlessbill)、支付方式 (paymethod) 等。服务订阅情况: 如是否开通电话服务 (phoneservice)、多条线路 (multiplelines)、互联网服务 (internetserv)、在线安全 (onlinesecurity)、在线备份 (onlinebackup)、技术支持 (techsupport)、设备保护 (deviceprotect) 等。希望通过构建这个预测模型，揭示可能导致客户流失的关键因素和行为模式，从而帮助 Telco 公司采取更具针对性的主动措施，以有效挽留有价值的客户。
首先，加载数据集到 Stata 内存，并初始化本地 H2O 集群：接着，将 Stata 当前的数据集传输到 H2O 集群中，创建一个名为 churn 的 H2O Frame（H2O 使用的数据结构），并将其设定为 H2O 当前的活动数据框架：_h2oframe put, into(churn)
现在，将 churn 数据框架分割为训练集和测试集，按 80%/20% 的比例划分。后续模型训练将在训练集上进行，并使用交叉验证 (CV) 来评估和控制过拟合。测试集将用于最终的模型比较。设置 rseed(19) 以确保数据分割过程可重复：_h2oframe split churn, into(train test) split(0.8 0.2) rseed(19)
将 train 设置为当前活动的数据框架，用于接下来的模型训练：为了方便后续命令的书写，创建一个名为 predictors 的 Stata 全局宏，用于存储所有预测变量的名称：global predictors latitude longitude tenuremonths monthlycharges totalcharges gender seniorcitizen partner dependents phoneservice multiplelines internetserv onlinesecurity onlinebackup streamtv techsupport streammovie contract paperlessbill paymethod deviceprotect
首先拟合一个基准的梯度提升机 (GBM) 模型。这里使用 3 折分层交叉验证 (stratified 3-fold cross-validation, SCV)，并通过 h2orseed(19) 确保结果可重复。对于分类任务，特别是当类别不平衡（如此例中流失与未流失客户数量差异较大）时，分层交叉验证非常有用，它能确保每个交叉验证折叠中类别比例与原始数据大致相同。
h2oml gbbinclass churn $predictors, h2orseed(19) cv(3, stratify)
执行命令后，Stata 的输出会提供模型的详细信息，包括：模型类型和使用的损失函数（对于 h2oml gbbinclass，是伯努利损失函数）。训练集样本量以及交叉验证的设置。模型参数（超参数）: 显示用户指定的超参数值以及算法实际使用的值。性能指标摘要: 表格展示模型在训练集和交叉验证集上的各项二元分类性能指标。由于本例中响应变量 churn 存在类别不平衡，应重点关注精确率-召回率曲线下面积 (Area under the precision-recall curve, AUCPR) 作为主要的性能评估指标。AUCPR 对少数类（即更关心的“流失”客户）的预测性能更为敏感。AUCPR 的取值范围是 0 到 1，值越接近 1 表示模型性能越好。虽然交叉验证的指标（如交叉验证 AUCPR）是评估模型泛化能力的核心依据，但同时也要关注训练集指标。比较训练集 AUCPR 和交叉验证 AUCPR 的差距有助于判断模型是否过拟合。通常训练集指标会略好于交叉验证指标（例如，训练集 AUCPR=0.8024，交叉验证 AUCPR=0.6585，差距 0.1439），但如果差距过大，则可能意味着模型过度拟合了训练数据，在新数据上的表现可能不佳。将这个基准模型的估计结果保存起来，命名为 gbm_default，以备后续比较：h2omlest store gbm_default
此外，可以使用 h2omlestat cvsummary 命令查看交叉验证过程中各折性能指标的均值和标准差。较小的标准差通常意味着模型性能对数据划分不敏感，模型更稳定。
基准模型的交叉验证 AUCPR 为 0.6585。为了尝试提升模型性能，进行超参数调优。GBM 有多个可调参数，为简化演示，这里仅在一个较小的搜索空间内调整 ntrees()（树的数量）和 predsamprate()（每次分裂时使用的预测变量抽样比例）这两个参数。请注意，超参数调优通常是一个需要多次尝试和迭代的过程。h2oml gbbinclass churn $predictors, h2orseed(19) cv(3, stratify) ntrees(50 100 200) predsamprate(0.15 0.3 0.6) tune(metric(aucpr))
(这里添加了 tune(metric(aucpr)) 选项来激活调优，并指定 AUCPR 为优化目标。默认使用网格搜索穷尽所有组合。)
Stata 的输出会显示调优过程的设置（如调优方法、优化指标、搜索网格）以及最终选定的最佳超参数组合（例如，ntrees()=100, predsamprate()=0.15），该组合对应了在交叉验证中获得最高 AUCPR 的模型。输出的其余部分则展示了这个最优模型的详细参数和性能指标。通过本次调优，交叉验证 AUCPR 从基准的 0.6585 提升至 0.6739。由于本示例仅探索了超参数空间的一小部分，提升幅度可能有限；实际应用中可能需要探索更广泛的参数组合。保存这个经过调优找到的最佳模型，命名为 gbm_tuned：可以使用 h2omlestat gridsummary 命令查看调优过程的详细摘要。该命令会列出所有尝试过的超参数配置及其对应的性能（按 AUCPR 排序）。h2omlestat gridsummary // 查看网格搜索摘
如果需要基于其他性能指标比较排名前列的模型，可以使用 h2omlexplore 命令：h2omlexplore id = 1/2 // 比较网格搜索中性能最好的前两个模型
若基于某些考虑（例如，希望模型更简单，使用更少的树）想选择一个并非最优指标但性能尚可的模型，可以使用 h2omlselect 命令指定其 ID 来激活。接下来，使用 h2omlgof 命令，更全面地比较调优后的最佳模型
gbm_tuned 和基准模型 gbm_default 在各项指标上的表现：h2omlgof gbm_tuned gbm_default
输出结果会分别展示训练集和交叉验证集的性能对比。关注交叉验证结果，可以看到调优是否提升了模型的整体性能（例如，是否有更低的对数损失、错误率、MSE，以及更高的 AUC、AUCPR、Gini 系数）。此外，还可以通过检查变量重要性图来考虑是否可以简化模型（移除重要性低的预测变量）：根据变量重要性图，研究者可能会决定移除某些贡献较小的预测变量（如 onlinebackup），然后重新训练模型。
假设已经按照类似步骤（拟合、调优、评估）训练了一个随机森林分类模型（使用 h2oml rfbinclass），并将其最佳调优结果存储为 rf_tuned。现在，需要在独立的测试集（即之前分割数据时创建的 test 数据框架）上比较最终选定的 GBM 模型 (gbm_tuned) 和随机森林模型 (rf_tuned) 的性能，以做出最终的方法选择。首先，使用 h2omlpostestframe 命令指定 test 数据框架作为后续所有后估计命令（如性能评估、预测）使用的默认数据源：(Stata 会提示：测试框架 test 现在已激活用于 h2oml 后估计)如果当前活动的模型不是要比较的模型，需要先恢复它们：h2omlest restore gbm_tuned // 确保 gbm_tuned 是活动模型之一h2omlest restore rf_tuned // 确保 rf_tuned 也是活动模型之一（如果需要分别评估）
使用 h2omlgof 命令可以方便地在测试集上并排比较两种方法的性能指标：h2omlgof gbm_tuned rf_tuned
比较结果（例如 AUCPR 或 AUC），选择在测试集上表现更优的方法。本例中，假设 GBM 的 AUCPR 更高，因此选择 GBM 作为最终模型。还可以使用 ROC 曲线进行可视化比较（曲线越靠近左上角越好）：h2omlgraph roc, models(gbm_tuned rf_tuned)
图形结果应与指标结果一致，进一步确认哪种方法性能更优。比较分类预测的另一种常用方法是查看各自的混淆矩阵 (h2omlestat confmatrix)，它可以详细展示正确和错误预测的数量。
假设 Telco 公司收集了一批新客户的数据，存储在 newchurn.dta 文件中。公司希望使用最终选定的模型 (gbm_tuned) 来预测这些新客户的流失概率。首先，加载新数据到 Stata，并将其传输到 H2O，命名为 newchurn 框架：webuse newchurn.dta, clear // 加载新数据_h2oframe put, into(newchurn)
恢复（如果需要）最终选定的模型 gbm_tuned：h2omlest restore gbm_tuned
使用 h2omlpredict 命令进行预测。重要：由于之前用 h2omlpostestframe 将 test 设置为了默认的后估计框架，现在要对 newchurn 数据进行预测，必须使用 frame(newchurn) 选项明确指定目标数据框架。同时预测类别标签和每个类别的概率：h2omlpredict churnhat, frame(newchurn) classh2omlpredict churnprob*, frame(newchurn) probabilities
预测结果（类别 churnhat 和概率 churnprob1, churnprob2）会作为新变量存储在 H2O 的 newchurn 框架中。如果需要将这些预测结果传回 Stata，可以使用 _h2oframe get 命令。可以查看 H2O 框架中的预测结果（例如，使用 _h2oframe list 或相关命令）。churnprob2 列显示了每个新客户被预测为“流失”（假设 "Yes" 是类别 2）的概率。churnhat 列则显示了基于某个阈值预测的最终类别（"Yes" 或 "No"）。这个阈值通常是基于优化某个指标（如 F1 分数）得到的，可以通过 h2omlestat threshmetric 命令找到模型对应的最优阈值（例如，之前提到的 0.2378）。
机器学习模型的一个关键优势在于其预测能力，但理解模型为何做出特定预测同样重要（即可解释性）。h2oml 提供了多种工具来帮助理解模型行为。解释方法可分为：局部可解释性: 解释模型对单个观测样本做出特定预测的原因。已经看过了变量重要性图 (h2omlgraph varimp)。这里介绍另外两种全局方法：全局代理模型 (Global surrogate model): 用一个更简单、本身可解释的模型（如单棵决策树）来近似复杂机器学习模型（如 GBM）的预测行为。首先，确保最佳 GBM 模型 (gbm_tuned) 已恢复，并在整个数据集 (churn框架) 上生成预测类别 churnhat：h2omlest restore gbm_tunedh2omlpredict churnhat, class frame(churn)
现在，训练一个简单的随机森林模型（设置 ntrees(1)即为单棵决策树，maxdepth(3)限制树深度）来拟合churnhat（GBM 的预测结果），以此作为 GBM 的代理模型：h2oml rfbinclass churnhat $predictors, h2orseed(19) ntrees(1) maxdepth(3)
这个简单的决策树可以被可视化，帮助理解复杂 GBM 模型做出预测的大致规则。使用 h2omltree命令可以保存树结构为 DOT 文件，再用 Graphviz 等工具绘制：h2omltree, dotsaving(churntree.dot, title(Surrogate tree for class "No"))
(通过分析代理树的路径和终端节点概率，可以了解哪些特征组合大致对应低/高流失风险）部分依赖图 (Partial Dependence Plot, PDP): 展示一个或多个预测变量如何平均地影响模型的预测结果（在控制其他变量后）。确保 GBM 模型已恢复，并将 churn 数据框架设为后估计分析的活动框架（但不是作为测试集）：h2omlest restore gbm_tunedh2omlpostestframe churn, notest
为几个重要的预测变量生成 PDP 图，并将它们合并显示：h2omlgraph pdp contract tenuremonths onlinesecurity techsupport, combine
(观察图中曲线的趋势，可以了解变量值变化与平均流失概率预测的关系，例如，长期合同、长在网时间的客户平均流失概率更低)SHAP (SHapley Additive exPlanations) 值: 基于博弈论，SHAP 值将单个预测与平均预测之间的差异，公平地分配给每个预测变量，量化每个变量对该特定预测的贡献（正向或负向）。使用 h2omlgraph shapvalues 为某个特定观测（如下例中的第 19 号观测）生成 SHAP 值贡献图，通常显示最重要的前几个变量：h2omlgraph shapvalues, obs(19) top(10) xlabel(-2.5(0.5)2)
(图中条形表示各变量对该观测预测结果的推力：蓝色增加流失可能，红色降低。结果有助于理解为何该客户被预测为会流失，例如按月合同、短在网时长是主要推高风险的因素)
注意：对于二元分类，SHAP 值通常在 logit 尺度上报告。它们解释的是预测 logit 值与基准 logit 值之间的差异。若要理解其对概率的影响，需结合逆 logit 变换。
使用 h2omlgraph shapsummary 生成 SHAP 摘要图（蜂群图），它将数据集中许多观测的 SHAP 值点绘制在一起，提供全局视角：h2omlgraph shapsummary, top(4) rseed(19)
(此图能同时展示：变量的全局重要性（点分布的宽度）、变量值高低与预测贡献方向的关系（点的颜色模式）、效应的异质性（点的垂直散布情况）。
例如，图中可能显示合同类型变量整体重要性高，其中代表短合同的值（蓝色点）普遍对应正 SHAP 值，即增加流失风险)通过综合运用这些模型评估和解释工具，研究者不仅能够构建出高性能的预测模型，还能深入理解模型的决策机制，从而更有信心地应用模型结果并采取相应行动。
source: https://www.stata.com/new-in-stata/machine-learning-h2o-ensemble-decision-trees/
群友可在社群中下载该文PDF版本和对应的code和代码。
关于Stata, 1.Stata19新增功能有哪些? 满满干货拿走不谢，2.Stata资料全分享，快点收藏学习，3.Stata统计功能、数据作图、学习资源，4.Stata学习的书籍和材料大放送, 以火力全开的势头
，5.史上最全Stata绘图技巧, 女生的最爱，6.把Stata结果输出到word, excel的干货方案，7.编程语言中的函数什么鬼？Stata所有函数在此集结，8.世界范围内使用最多的500个Stata程序，9.6张图掌握Stata软件的方方面面, 还有谁, 还有谁? 10.LR检验、Wald检验、LM检验什么鬼？怎么在Stata实现，11.Stata15版新功能，你竟然没有想到，一睹为快，12."高级计量经济学及Stata应用"和"Stata十八讲"配套数据，13.数据管理的Stata程序功夫秘籍，14.非线性面板模型中内生性解决方案以及Stata命令，15.把动态面板命令讲清楚了，对Stata的ado详尽解，16.半参数估计思想和Stata操作示例，17.Stata最有用的points都在这里,无可替代的材料，18.PSM倾向匹配Stata操作详细步骤和代码，干货十足，19.随机前沿分析和包络数据分析 SFA,DEA 及Stata操作，20.福利大放送, Stata编程技巧和使用Tips大集成，21.使用Stata进行随机前沿分析的经典操作指南，22.Stata, 不可能后悔的10篇文章, 编程code和注解，23.用Stata学习Econometrics的小tips, 第二发礼炮，24.用Stata学习Econometrics的小tips, 第一发礼炮，25.广义合成控制法gsynth, Stata运行程序release，26.多重中介效应的估计与检验, Stata MP15可下载，27.输出变量的描述性统计的方案，28.2SLS第一阶段输出, 截面或面板数据及统计值都行，29.盈余管理指标的构建及其Stata实现程序, 对应解读和经典文献
，30.Python, Stata, R软件史上最全快捷键合辑！，31.用Stata做面板数据分析, 操作代码应有尽有，32.用Stata做面板数据分析, 操作代码应有尽有，33.没有这5个Stata命令, 我真的会活不下去！，34.第一(二)卷.Stata最新且有趣的程序系列汇编，35.第三卷.Stata最新且急需的程序系列汇编，36.第四卷.Stata最新且急需的程序系列汇编
下面这些短链接文章属于合集，可以收藏起来阅读，不然以后都找不到了。
