Stata：机器学习分类器大全

发布时间：2021-01-07 阅读 4065

Stata连享会主页 || 视频 || 推文 || 知乎

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
. ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
. help lianxh

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

作者: 樊嘉诚 (中山大学)

E-Mail: fanjch676@163.com

1. 引言

“如何根据西瓜的色泽、根蒂、敲声等特征分辨出好瓜和坏瓜”，这是我们日常生活中经常面临的分类问题 (classification problem) 。而在学术研究中，诸多研究都离不开分类的影子：识别经济周期，判断未来经济形势；研究上市公司财务信息，对其财务困境或危机进行预警……此外，计算机视觉、垃圾邮件分类、医学诊断等也与分类问题密切相关。

回想经典的计量模型 Logit 回归，为解决离散选择问题提供思路，本质上也可视作一种分类算法。随着大数据时代的到来，许多分类任务面临着数据维度过高、数据质量较低、样本不平衡等诸多问题。机器学习 (Machine Learning, ML) 算法作为近年来炙手可热的方法，为解决这些问题开辟了新的思路与途径。

本推文将要介绍的命令 c_ml_stata 利用 Python 语言在 Stata 中实现了机学习分类算法，不仅囊括了众多分类算法，如支持向量机、决策树、神经网络等；也支持交叉验证 (Cross Validation, CV) ，是利用 Stata 处理分类问题有力工具。本推文的余下部分安排如下：在第二部分，对该命令的部分机器学习分类算法进行简单的理论介绍；在第三部分，对该命令的简要介绍和安装方法进行说明；在第四部分，利用该命令及其提供的数据集使用 Stata 进行分类问题的处理；在第五部分，对本推文主要内容进行总结。

2. 理论介绍

机器学习分类算法众多，由于篇幅有限，现结合 c_ml_stata 命令中提供的部分分类算法进行简要的理论介绍，以便对机器学习分类问题、算法及后续命令使用有更清楚的认识。熟悉这些算法的读者可以快速跳过。该部分主要介绍的机器学习算法包括：支持向量机 (Support Vector Machine, SVM) 、决策树 (Decesion Tree) 和神经网络 (Neural Network, NN) 。

2.1 支持向量机

2.1.1 支持向量与间隔

支持向量机是一种二分类器，它的基本思想是基于训练集 $D$ 在样本空间中寻找一个划分超平面，将不同类别的样本划分开。

支持向量机学习方法包括由简至繁的一系列模型：当训练数据线性可分时，通过硬间隔最大化 (hard margin maximization) ，学习一个线性分类器，即线性可分支持向量机；当训练数据近似线性可分时，通过软间隔最大化 (soft margin maximization) ，也学习一个线性分类器，即线性支持向量机；当训练数据线性不可分时，通过使用核方法 (kernel method) 及软间隔最大化，学习非线性支持向量机。

我们从最简单的线性可分支持向量机作为引入，假定给定一个特征空间训练数据集 $D = {(x_{i}, y_{i})}_{i = 1}^{N}$ ，其中， $x_{i} \in R^{n}$ 表示具有 $n$ 个特征的特征向量 (feature vector) ， $y_{i} \in {+ 1, - 1}$ 表示因变量的两个不同的分类标签 (class label) 。并且假定训练数据集是线性可分的。我们的目的是找到一个划分超平面 $w^{'} x + b = 0$ 将实例划分到不同的类，其中， $w$ 为法向量， $b$ 为截距。故可定义相应的分类函数为

y_{i} = f (x_{i}; w, b) = s i g n （ w^{'} x_{i} + b ）= {\begin{matrix} + 1, i f w^{'} x_{i} + b \geq + 1 \\ - 1, i f w^{'} x_{i} + b \leq - 1 \end{matrix}

特别地，若 $w x_{i} + b = 0$ 表示样本点位于划分超平面上，被称之为“支持向量” (support vector) 。

一般地，对于线性可分的数据，存在无穷多个划分超平面可以将两类数据正确地分开，那么如何获得一个唯一的最优划分超平面呢？在下图中，有 $A, B, C$ 三个点，表示 $3$ 个实例，且预测分类时均在划分超平面的一侧。因为 $A$ 点距分类超平面较远，就比较确信 $A$ 点被正确分类的可信度较高； $C$ 点离划分超平面较近，因此可能怀疑划分是否准确 (可信度较低) ；而 $B$ 点介于点 $A$ 和 $C$ 之间，预测其分类结果的可信度介于 $A$ 与 $C$ 之间。上述直觉引导我们如何确定最优的划分超平面——一般来说，一个点离划分超平面的远近可以表示分类结果的可信度，而最优的划分超平面应使得所有点的分类结果可信度尽可能的高。由基本的空间几何知识可以知道， $| y_{i} (w^{'} x_{i} + b) |$ 恰好可以相对地表示第 $i$ 个实例到划分超平面的距离。此外，若 $y_{i} (w^{'} x_{i} + b) > 0$ 则表示分类正确。在此，引入一个非常关键的概念—— 间隔 (margin)，定义为

γ = min_{i = 1, \dots, N} y_{i} (w^{'} x_{i} + b)

基于最初的直觉，我们的目标自然是希望即使是离划分超平面最近的点，其划分结果的可信度也较高。因此，最大化间隔是一个不错的想法。然而，注意到如果我们成比例地改变 $w$ 和 $b$ ，虽然我们的间隔变化了，但是超平面本身并未改变。因此，为求得唯一解，对法向量施加某种约束是必要的 (如规范化，约束 $∥ w ∥ = 1$ ) 。最终，我们可以写出如下最优化问题

\begin{array}{cc} max_{w, b} γ \\ s.t. ∥ w ∥ = 1 \end{array}

略微遗憾的是，以上最优化问题是“非凸” (Non-convex) 的，求解过程较为复杂。幸运的是，上述最优化问题可以通过等价变换，转换为以下凸 (Convex) 问题：

\begin{array}{l} min_{w, b} \frac{1}{2} ∥ w ∥^{2} \\ s.t. y_{i} (w^{'} x_{i} + b) \geq 1, \forall i = 1, \dots, N \end{array}

求得最优解 $w^{*}, b^{*}$ 即可得到最优划分超平面。可以证明，该最优划分超平面是存在且唯一的。

2.1.2 软间隔

以上的讨论中，我们假定了训练数据集可以被划分超平面准确、完全地分开。然而，在实际问题中，我们往往难以确定数据是否线性可分；或者，即使数据线性可分，也很难断定这个貌似线性可分的结果不是由于过拟合导致的。因此，为了减少数据中”噪声“的干扰，我们允许支持向量机在一些样本上的分类结果出错。为此，我们需要引入软间隔 (soft margin) 的概念。之前的约束要求所有样本均划分正确，即满足 $y_{i} (w^{'} x_{i} + b) \geq 1, \forall i = 1, \dots, N$ ，这可以理解为“硬间隔”。而软间隔则允许某些样本不满足该约束，当然这些不满足样本的约束应该尽可能的少，因此优化目标可以写为

min_{w, b} \frac{1}{2} {∥ w ∥}^{2} + C \sum_{i = 1}^{N} I (y_{i} (w^{'} x_{i} + b) - 1)

其中， $C$ 是正则化系数， $C$ 越大，表示我们对错误分类的容忍程度越低，间隔越小；反之， $C$ 越小，表示我们对错误分类的容忍程度越大，间隔越大。 $I (z)$ 为示性函数，即

I (z) = {\begin{matrix} 1, & i f z > 0 \\ 0, & o t h e r w i s e \end{matrix}

引入松弛变量 (slack variables) $ξ_{i} \geq 0$ ，可以将上式重写为

\begin{array}{ll} min_{w, b} \frac{1}{2} ∥ w ∥^{2} + C \sum_{i = 1}^{N} ξ_{i} \\ s.t. y_{i} (w^{'} x_{i} + b) \geq 1 - ξ_{i} \\ ξ_{i} \geq 0, \forall i = 1, \dots, N \end{array}

2.1.3 核方法

由于现实中的许多问题并非是线性可分的，对于非线性可分的数据，常采用核方法将样本从原始空间映射到更高维的特征空间 (如下图所示) ，使得样本在这个特征空间内线性可分。

令 $ϕ (x)$ 表示 $x$ 映射后的特征向量，于是，在特征空间中划分超平面所对应的模型可以表示为

w^{'} ϕ (x) + b = 0

因此，最优化问题可以写为

\begin{array}{ll} min_{w, b} \frac{1}{2} ∥ w ∥^{2} \\ s.t. y_{i} (w^{'} ϕ (x_{i}) + b) \geq 1, \forall i = 1, \dots, N \end{array}

在求解过程中，由于需要计算样本 $x_{i}$ 与 $x_{j}$ 映射到特征空间之间的内积 $ϕ (x_{i})^{'} ϕ (x_{j})$ 。因为特征空间的维度可能很高，直接计算内积通常十分困难，因而设想如下函数

κ (x_{i}, x_{j}) = ⟨ ϕ (x_{i}), ϕ (x_{j}) ⟩ = ϕ (x_{i})^{'} ϕ (x_{j})

可以大大简化运算过程。而函数 $κ (x_{i}, x_{j})$ 称为核函数 (kernel function) 。核函数的选取并非是任意的，需要满足一些条件，受限于篇幅我们不作详细讨论。以下列出一些常用的核函数：

线性核： $κ (x_{i}, x_{j}) = x_{i}^{'} x_{j}$ ；
多项式核： $κ (x_{i}, x_{j}) = (x_{i}^{'} x_{j})^{d}$ ，参数 $d \geq 1$ 为多项式的次数；
高斯核： $κ (x_{i}, x_{j}) = e x p (- \frac{| | x_{i} - x_{j} | |^{2}}{2 σ^{2}})$ ，参数 $σ > 0$ 为高斯核的带宽 (width) ；
拉普拉斯核： $κ (x_{i}, x_{j}) = e x p (- \frac{| | x_{i} - x_{j} | |^{2}}{σ})$ ，参数 $σ > 0$ ；
Sigmoid 核： $κ (x_{i}, x_{j}) = t a n h (β x_{i}^{'} x_{j} + θ)$ ，其中 $t a n h (\cdot)$ 为双曲正切函数， $β > 0, θ < 0$ 。

2.1.4 补充

支持向量机可以构造对偶问题，利用拉格朗日乘子法求解。支持向量机最终可转化为一个二次规划问题，使用诸如 SMO ( Sequential Minimal Optimization ) 等高效算法求解 (由于我们本着浅显地了解支持向量机的基本理论，便不详细介绍其求解优化过程)。

此外，支持向量机有以下优劣势，我们在使用该分类方法时需额外注意：

优势：
- 支持向量机的最优化函数仅由少数的支持向量确定，计算的复杂性取决于支持向量的个数，而非样本空间的维数，因此在某种程度上可以避免“维数灾难”。
- 支持向量机拥相对于其他黑箱 (Black box) 的机器学习方法有更多的更有理论，并且其结果具有较好的稳健性 (Robust) 。
劣势：
- 支持向量机算法难以用于大规模的训练样本的计算，对缺失数据和核函的选取较为敏感。
- 传统的支持向量机算法仅适用于二分类问题，而实际应用中常常面临多分类任务。

2.2 决策树

2.2.1 基本模型

决策树是一种基于的分类和回归方法，顾名思义，决策树呈现树形结构 (见下图) 。一颗决策树由结点 (node) 和有向边 (directed edge) 或分枝组成，结点一般包括根结点、内部结点和叶结点，可以形象类比为“树根”和“树叶”，有向边可以理解为“树枝”。用决策树分类的基本思想是，从根结点开始，对样本的某一特征 (划分依据) 进行测试，根据测试结果将样本分配到子结点；如此递归地对样本进行测试并分配，直至到达叶结点。

2.2.2 特征选择

决策树学习的关键问题之一是特征选择，即在每次分类时选择什么特征进行测试和划分。因此，我们需要确定选择特征的准则。直观上，如果某一个特征具有更好的分类能力，那么决策树按这一特征分类后的各个子类应尽可能地属于同一类别，结点的“纯度”越高。接下来，我们会依次引入一些概念：信息熵、条件熵、信息增益以及信息增益比，来理解如何进行特征选择。

信息熵 (information entropy) 是度量样本集合纯度的一种指标，对于一个有 $K$ 个离散取值的随机变量 $X$ ，其概率分布为

P (X = x_{k}) = p_{k}, k = 1, \dots, K

其信息熵定义为

H (p) = - \sum_{k = 1}^{K} p_{k} l o g p_{k}

特别地，若 $p_{k} = 0$ ，定义 $0 l o g 0 = 0$ ；上式的对数常以以 $2$ 为底或者以 $e$ 为底。熵越大，随机变量的不确定性越大。可以证明 $0 \leq H (p) \leq l o g K$ 。

对于随机变量 $(X, Y)$ ，其联合概率分布为

P (X = x_{i}, Y = y_{j}) = p_{i j}, i = 1, \dots, n, j = 1, \dots, m

条件熵 (conditional entropy) $H (Y | X)$ 表示在已知随机变量 $X$ 的条件下随机变量 $Y$ 的不确定性，定义为 $X$ 给定条件下 $Y$ 的条件概率分布的熵对 $X$ 的数学期望

H (Y | X) = \sum_{i = 1}^{n} p_{i} H (Y | X = x_{i})

这里， $p_{i} = P (X = x_{i}), i = 1, \dots, n$ 。当熵和条件熵中的概率由数据估计获得时，所对应的熵和条件熵分别成为经验熵和经验条件熵。

信息增益 (information gain) 表示得知特征 $X$ 的信息而使得类 $D$ 的信息不确定性减少的程度，因此我们定义特征 $X$ 对训练数据集 $D$ 的信息增益 $g (D, X)$ 为数据集 $T$ 的经验熵 $H (D)$ 与特征 $X$ 给定下其条件熵 $H (D | X)$ 的差，即

g (D, X) = H (D) - H (D | X)

显然，对于具有较强分类能力的特征，其信息增益更高。因此，我们利用信息增益选择特征的方法是，对于训练数据集 $D$ ，计算其每个特征的信息增益，选择信息增益最大的特征。

但是，使用信息增益作为划分标准存在偏向于选择特征取值较多的特征的问题，这样是不公平的。因此，引入了信息增益比 (information gain ratio) 的概念。定义特征 $X$ 对训练数据集 $D$ 的信息增益比 $g_{R} (D, X)$ 为其信息增益 $g (D, X)$ 与训练数据集 $D$ 关于特征 $X$ 的值的熵 $H_{X} (D)$ 之比，即

g_{R} (D, X) = \frac{g (D, X)}{H_{X} (D)}

2.2.3 树的生成

决策树的生成有多种算法，如 ID3 、C4.5 等经典的生成算法。为了理解决策树的生成过程，我们还是选择介绍其中的一种生成算法：ID3 ，其核心思想是在树的各个结点用信息增益作为特征选择准则，递归地构建决策树。具体方法是：

从根结点开始，对结点计算所有可能特征的信息增益，选择信息增益最大的特征作为结点的特征，并由该特征的不同取值建立不同的子结点；
对子结点递归地使用以上方法，构建决策树；
直到所有特征的信息增益很小或者所有特征均选择完毕为止。

C4.5 算法与 ID3 算法相似，不同之处在于， C4.5 使用信息增益比作为特征选择的依据。此外还有诸如 CART 算法等等多种多样的生成树的方法。

2.2.4 补充

生成决策树后，往往还需要对其进行剪枝 (pruning) ，顾名思义，就是从已生成的树上裁剪一些子树或者叶结点，对树的结构进行简化以防止其过拟合。总结一下决策树的优缺点：

优点：
- 易于理解和实现，并且能够理解决策树所表达特征的含义；
- 对缺失值不敏感，能处理大规模数据。
缺点：
- 对有时间顺序的数据，需要很多数据预处理。
- 忽略属性之间的相关性。

因此，在其基础上也有诸多拓展模型：与袋装法 (Bagging) 思想结合的随机森林 (Random Forest) ，与提升法 (Boosting) 结合的梯度提升树 (Gradient Boosting Decesion Tree) 、极端梯度提升树 (Extreme Gradient Boosting Decesion Tree) 等等。

2.3 神经网络

2.3.1 神经元

神经网络是现在比较流行的机器学习算法，可以处理回归、分类等多种问题。神经网络中最基本的结构是神经元 (neuron) ，其结构见下图。

一个最基本的神经元由输入 (input) 、权重 (weight) 、偏置 (bias) 或阈值 (threshold) 、激活函数 (active function) 和输出 (output) 组成。以一个有多个输入只有一个输出的神经元为例，其接受了 $n$ 个输入信号 $x_{1}, x_{2}, \dots, x_{n}$ ，即 $x = (x_{1}, x_{2}, \dots, x_{n})^{'}$ ，这些输入信号通过带权重的连接进行传递，其加权后的总输入与神经元的阈值比较，通过激活函数处理产生输出。若将第 $i$ 个输入信号的权重写为 $w_{i}$ ，则加权总输入为 $\sum_{i = 1}^{n} w_{i} x_{i} = w^{'} x$ 。设阈值为 $θ$ ，权重函数为 $f (\cdot)$ ，则神经元产生的输出为

y = (w^{'} x - θ) = f (\sum_{i = 1}^{n} w_{i} x_{i} - θ)

激活函数 $f (\cdot)$ 往往是非线性的，有以下常用多种激活函数可以使用：

Sigmoid: $σ (x) = \frac{1}{1 + e x p (- x)}$ ；
tanh: $t a n h (x) = \frac{e^{x} - e^{- x}}{e^{x} + e^{- x}}$
ReLU: $m a x (0, x)$

一般而言，常常选取 ReLu 激活函数，原因是该函数形式较为简单，计算快、收敛快且可以避免诸如梯度消失等问题。

2.3.2 前馈神经网络

神经网络结构一般由**输入层 (input layer) ** 、**隐层 (hidden layer) ** 和 输出层 (output layer) 组成。

输入层：由多个无激活函数结构的神经元构成；
隐层和输出层：由多个有激活函数结构的神经元组成，且同层间神经元的激活函数往往相同。

最经典也是最常见的前馈神经网络正是由输入层、隐层和输出层构成，其特点是含有多个隐层，每层神经元与下一层神经元全互连，神经元之间不存在同层连接，也不存在跨层连接。如下图所示。

2.3.3 神经网络分类算法

一般地，对于一个多分类问题，我们可以定义神经网络的输出 $K \times 1$ 维列向量 $y = (y_{1}, y_{2}, \dots, y_{K})^{'}$ ，其中 $y_{j}, \forall j = 1, \dots, K$ 既可以是连续取值，也可以是离散取值如 $y = (0, 1, \dots, 0)^{'}$ 。该网络共有 $L$ 个隐层，简单起见我们各个隐层及输出层神经元的激活函数均为 $f (\cdot)$ ，输入层接受的输入即分类的特征 $n \times 1$ 维列向量 $x = (x_{1}, x_{2}, \dots, x_{n})^{'}$ 。该前馈神经网络可以写为如下递归形式：

\begin{array}{l} x^{(0)} = x \\ x^{(t)} = f (W^{(l - 1)} x^{(l - 1)} θ^{(l - 1)}), l = 1, \dots, L \\ \hat{y} = f (W^{(L)} x^{(L)} θ^{(L)}) \end{array}

其中， $W^{(l)}$ 和 $θ^{(l)}$ 分别为各个隐层神经元权重与偏置组成的矩阵或向量， $\hat{y}$ 即为神经网络的输出。特别地，若输出 $\hat{y} = ({\hat{y}}_{1}, {\hat{y}}_{2}, \dots, {\hat{y}}_{K})^{'}$ 中 ${\hat{y}}_{k}, \forall k \in [1, K]$ 为连续取值，想要得到分类结果或离散的取值，我们可以使用 softmax 方法计算其“概率”，即

{\tilde{y}}_{k} = \frac{e x p (- {\hat{y}}_{k})}{\sum_{k = 1}^{K} e x p (- {\hat{y}}_{k})}

设定某一阈值，或将其最大的 ${\tilde{y}}_{k}$ 对应为 $1$ ，其余为 $0$ ，即可得到分类结果。

2.3.4 补充

前馈神经网络中的权重 $W$ 和偏置 $θ$ 均是未知参数，常用误差逆传播算法 (Error BackPropagation, BP) 进行确定。神经网络算法有以下需要注意的地方：

神经网络的结构 (即设定多少个隐层、包含多少个神经元) 需要事先设定，其结构与算法的收敛性、是否过拟合 (over fitting) 等问题密切相关；
为了便于计算、消除量纲的影响，一般对输入数据进行标准化 (normalization) 处理。
由于神经网络具有强大的表示能力，经常需使用一些策略避免过拟合，如早停 (early stop) 、正则化 (regularization) 、丢包法 (dropout) 等。

3. 命令介绍和安装

熬过以上略微繁琐的算法理论介绍，正是来到本次推文的主角命令——c_ml_stata 。虽然机器学习算法十分复杂 (远比以上的理论介绍要复杂许多) ，但是该命令较为简洁而直接地集成了多种算法，同时十分容易调用。该部分主要对命令进行基本介绍，说明其安装方法以及语法选项的使用及一些注意事项。

3.1 基本介绍

c_ml_stata 由 Giovanni Cerulli 编写，是在 Stata 16 中实现机器学习分类算法的命令，该命令使用 Python 中的 Scikit-learn 接口实现模型训练、预测等功能，主要有以下特点：

支持多种分类算法： 该命令提供了分类树 (Classification tree) 、袋装树和随机森林 (Bagging and random forests) 、提升算法 (Boosting) 、正则化多项式 (Regularized multiomial) 、K 近邻算法 (K-Neareast Neighbor) 、神经网络 (Neural network) 、朴素贝叶斯 (Naive Bayes) 和支持向量机 (Support vector machine) 多种分类算法，便于模型的尝试和比较。
支持交叉验证： 该命令提供了 cross_validation 选项，利用“贪婪搜索” (greed search) 实现 K 折交叉验证 (K-fold cross validation) 选择最优超参数 (hyper parameters) ，调优分类模型。

3.2 安装方法

c_ml_stata 需要在 Stata 16.0 及以上版本使用。在 Stata 的命令行中输入 ssc install c_ml_stata 即可下载，或者使用如下命令打开下载页面：

. veiw net describe c_ml_stata  // 命令包简介

查看与之相关的完整程序文件和相关附件，主要包括以下文件：

example_c_ml_stata.do ：作者提供示例的 do 文件。
c_ml_stata_data_example.dta ：作者提供示例的样本内训练数据。
c_ml_stata_data_new_example.dta ：作者提供示例的样本外预测数据。

如果想完全使用命令来安装，可以执行如下两条命令：

· net install  c_ml_stata  // 安装命令包
. net get      c_ml_stata  // 下载附件：dofile, .dta 等
                           // 存储在当前工作路径下

3.3 语法及选项

c_ml_stata 主要命令的语法格式为：

c_ml_stata outcome [varlist], mlmodel(modeltype) out_sample(filename)
		 in_prediction(name) out_prediction(name) cross_validation(name)
		 seed(integer) [save_graph_cv(name)]

输入变量的含义如下：

outcome：是一个数值型、离散的因变量 (或标签) ，表示不同的类别。若因变量有 $M$ 种类别，建议对其类别进行编码 (recode) ，取值范围为 $[1, 2, \dots, M]$ 。例如，对于一个二元变量 (取值为 $0$ 或 $1$ ) ，则应编码成 $[1, 2]$ 。注意：outcome 不接受缺失值。
varlist：是代表自变量 (或特征) 的数值型变量列表，属于可选项。若某一特征也是类别变量，则需先生成相应的数值型、离散的虚拟变量。注意：varlist 不接受缺失值。

各个选项的含义如下：

mlmodel(modeltype)：指定使用的机器学习分类算法 (模型) ，有以下几种选择：
- tree : Classification tree (分类树)
- randomforest : Bagging and random forests (袋装树和随机森林)
- boost : Boosting (提升算法，提升树)
- regularizedmultionmial : Regularized multinomial (正则化多项式)
- nearestneighbor : Nearest Neighbor (K 近邻算法)
- neuralnet : Neural network (神经网络)
- naivebayes : Naive Bayes (朴素贝叶斯)
- svm : Support vector machine (支持向量机)
out_sample(filename)：要求指定一个样本外的新数据集 (测试集) ，该数据集仅包含各个特征 (无因变量) ，用于样本外测试。filename 表示存放该数据集的文件名。
in_prediction(name)：保存样本内训练数据 (训练集和验证集) 的拟合结果，name 为文件名。
out_prediction(name)：保存样本外数据 (测试集) 的预测结果，name 为文件名。样本外数据从 out_sample 中获得。
cross_validation(name)：将 name 设定为 "CV" 可以执行交叉验证，默认为 10 折交叉验证。
seed(integer)：随机种子 (整数) 。
[save_graph_cv(name)]：可选项，保存交叉验证中模型在训练集和验证集上分类结果的准确性 ( Accuracy ) ，用于确定最优的超参数和模型。

c_ml_stata 的返回值：

输出测试集的预测结果和预测为各个类别的概率。
储存通过交叉验证确定的最优超参数 (不同模型储存的超参数不同) 、最优训练集准确率、最优验证集准确率，可以通过 ereturn list 命令查看 (数值型返回值储存在 scalars 中，字符型结果储存在 macros 中) 。

3.4 注意事项

运行 c_ml_stata 程序需要拥有 Stata 16 及 Python (2.7 及以上版本) ，同时需安装 Python 的 Scikit-learn 和 Stata Function Interface (SFI) 两个依赖库。
outcome 和 varlist 均不允许出现缺失值，因此在使用该命令前需检查数据集是否出现缺失值 (并删除缺失值)。
建议安装该命令的最新版本并及时更新，命令为：ssc install c_ml_stata, replace。
更多详细介绍可使用 help c_ml_stata 命令获取。

4. Stata 实操

我们使用 c_ml_stata 提供的数据集 (储存在 c_ml_stata_data_example.dta 文件中) 进行实操。由于作者并未提供过多有关该数据集的其他解释信息，所以该例子仅作为操作和演示提供，实际意义不大。

4.1 数据结构描述

该数据集共有 74 个样本，包含 4 个解释变量 (分别命名为 x1, x2, x3, x4 ) 和 1 个因变量 (命名为 y ) 。其中，4 个解释变量均为连续变量，而因变量为分类变量 (离散取值) ，因此我们分别使用 summary 和 tab 命令对变量进行描述性统计。可以看出，根据 c_ml_stata 命令的要求，输入数据集并无缺失值。待分类变量 (作为 outcome 的输入) y 进行数值化编码， $1, 2, 3$ 分别代表三种类别。

. use "c_ml_stata_data_example.dta", clear
. tab y

     y |  Freq.    Percent        Cum.
-------+------------------------------
     1 |     42      56.76       56.76
     2 |     22      29.73       86.49
     3 |     10      13.51      100.00
-------+------------------------------
 Total |     74     100.00

. sum x1-x4

 Variable |  Obs       Mean    Std. Dev.    Min     Max
----------+--------------------------------------------
       x1 |   74   6165.257    2949.496    3291   15906
       x2 |   74    21.2973    5.785503      12      41
       x3 |   74   3019.459    777.1936    1760    4840
       x4 |   74   187.9324    22.26634     142     233

4.2 模型训练和结果

由于 c_ml_stata 提供了多种有监督学习分类算法，我们使用变量 x1, x2, x3, x4 作为解释变量，y 作为标签进行模型训练。以下部分我们以支持向量机为例，详细介绍该命令的调用方法和输出结果。

4.2.1 支持向量机

将选项 mlmodel 设定为 svm 即可使用支持向量机算法进行分类。样本内预测结果保存在文件 in_pred_svm.dta 中，使用 c_ml_stata_data_new_example.dta 文件读取样本外数据 (仅包含特征 x1, x2, x3, x4 ) ，样本外预测结果保存在文件 out_pre_svm.dta 中。在 cross_validation 选项中使用 CV 即可进行交叉验证，交叉验证结果自动保存在 CV.dta 文件中；若有需要可使用 save_graph_cv 选项可视化交叉验证结果并保存。具体代码及部分输出结果如下：

. c_ml_stata y x1-x4, mlmodel(svm)   ///
    out_sample("c_ml_stata_data_new_example")  	 ///
	in_prediction("in_pred_svm")     ///
    out_prediction("out_pred_svm")   ///
	cross_validation("CV")           ///
    seed(10) save_graph_cv("graph_cv_svm")

-------------------------------------------
CROSS-VALIDATION RESULTS TABLE
-------------------------------------------
The best score is:
0.5678571428571428
-------------------------------------------
The best parameters are:
{'C': 1, 'gamma': 0.1}
1
0.1
-------------------------------------------
The best estimator is:
SVC(C=1, break_ties=False, cache_size=200, class_weight=None, coef0=0.0,
    decision_function_shape='ovr', degree=3, gamma=0.1, kernel='rbf',
    max_iter=-1, probability=True, random_state=None, shrinking=True, tol=0.001,
    verbose=False)
-------------------------------------------
The best index is:
0
-------------------------------------------

in_pred_svm.dta 部分数据如下图所示，其中 index 表示观测值样本标号，与原数据样本标号相对应；label_in_pred 表示样本内标签的预测结果；Prob_1, Prob_2, Prob_3 可能表示预测结果不是第 $1 / 2 / 3$ 类的概率。

out_pred_svm.dta 数据结果与 in_pred_svm.dta 类似。其中，label_out_pre 表示样本外标签的预测结果，可以看出 SVM 将样本外结果分类为第 $3$ 类。

对于支持向量机算法，根据第二部分理论部分的介绍，我们的主要超参数为正则化系数 $C$ 和核函数参数 $σ$ (常记作 GAMMA ) 。使用 ereturn list 可返回最优的超参数选取结果，可以看出最优的正则化系数应设定为 $1$ ，和函数参数选取为 $0.1$ ；交叉验证中训练集的正确率为 $100 %$ ，测试集 (即验证集) 正确率为 $56.79 %$ 。

. ereturn list

scalars:
              e(OPT_C) =  1
          e(OPT_GAMMA) =  .1
      e(TEST_ACCURACY) =  .5678571428571428
     e(TRAIN_ACCURACY) =  1
         e(BEST_INDEX) =  94.5

详细的交叉验证结果可在 CV.dta 中查看。此外，通过 save_graph_cv("graph_cv_svm") 可视化交叉验证结果并保存为 graph_cv_svm.gph 文件，结果如下所示。可以发现，随着 index 的改变 (不同超参数组合的设定) ，模型在训练集、验证集的表现均不改变。值得注意的是，这种情况在实际应用中极为少见。

4.2.2 决策树

将 mlmodel 选项设定为 tree 即可使用决策树进行分类。决策树的超参数主要为叶子节点个数 (leaves) ，通过交叉验证的可视化结果可以发现，随着 index 的增加 (叶子节点个数的增加) ，模型在训练集的表现不断提高，而在验证集的表现先上升后呈下降趋势 (出现过拟合问题) ，因此可以得到最优的叶子节点个数 (对应于验证集最高的分类正确率) 。使用 ereturn list 查看最有超参数的选取结果。可以发现，最优的叶子节点个数为 3 ；交叉验证中训练集的正确率为 $81.08 %$ ，测试集 (即验证集) 正确率为 $63.75 %$ 。

. c_ml_stata y x1-x4, mlmodel(tree)    ///
    out_sample("c_ml_stata_data_new_example")  	///
    in_prediction("in_pred_ctree")     ///
    out_prediction("out_pred_ctree")   ///
	cross_validation("CV")             ///
    seed(10) save_graph_cv("graph_cv_ctree")

. ereturn list

scalars:
         e(OPT_LEAVES) =  3
      e(TEST_ACCURACY) =  .6375
     e(TRAIN_ACCURACY) =  .8108095884215288
         e(BEST_INDEX) =  2

4.2.3 神经网络

将 mlmodel 选项设定为 neuralnet 即可使用神经网络进行分类。神经网络的主要超参数为神经网络层数 (layers) 和神经元个数 (neurons) 。通过查看神经网络在样本内外的预测结果发现，该分类算法在该数据集上的表现较差，“暴力”地将所有标签均分类为 $1$ ，分类正确率仅为 $56.75 %$ 。通过交叉验证可视乎结果也可以看出，该神经网络模型对超参数选取较为敏感，而训练集准确率均为 $56.75 %$ ，没有任何改善。

4.3 结果汇总

c_ml_stata 命令共提供了 8 种分类模型，逐一使用各模型对所提供的样本内数据集进行训练后，样本内最优分类结果的准确率及验证集集准确率如下所示：

模型	最优训练集准确率	最优验证集准确率
决策树	$81.08 %$	$63.93 %$
随机森林	$76.42 %$	$71.25 %$
提升数	$75.53 %$	$68.75 %$
K 近邻	$72.50 %$	$75.21 %$
神经网络	$56.76 %$	$56.76 %$
朴素贝叶斯	$68.33 %$	$68.93 %$
支持向量机	$100 %$	$56.79 %$

Note： 正则化多项式结果无法收敛，暂不参与比较。

5. 总结

本推文的内容即将进入尾声，简要回顾我们上述的主要内容：我们简单了解了什么是分类问题以及机器学习的分类算法，介绍了 c_ml_stata 这一 Python 和 Stata 结合的命令及其使用方法。我们可以看到随着时代的发展，机器学习算法的多样性以及其广阔的应用范围；但是，我们始终不能将方法做为最终的目的，其背后的经济学含义依然值得我们思考。

6. 参考资料

周志华，机器学习，清华大学出版社，2016
李航，统计学习方法，清华大学出版社，2012
连享会推文专题：
- Stata: 交叉验证简介
- 人工神经网络与Stata应用

专题	嘉宾	直播/回看视频
⭐ 最新专题		因果推断, 空间计量，寒暑假班等
⭕ 数据清洗系列	游万海	直播, 88 元，已上线
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]