线性分类 方匡南 教授 博士生导师 耶鲁大学博士后 厦门大学数据挖掘研究中心 副主任 厦门大学经济学院统计系 中华数据挖掘协会(台湾） 顾问 厦门思创优合信息技术有限公司 顾问 厦门数析信息技术有限公司 顾问 Email: xmufkn@xmu.edu.cn 线性分类 Linear Classification 方匡南 教授 博士生导师

线性分类 第5章我们探讨的线性回归模型假设因变量Y是定量（quantitative）的，但在很多实 际问题中，因变量却是定性的（qualitative）。所谓定性变量，是指这些量的取值 并非有数量上的变化，而只有性质上的差异。定性变量也称为分类（categorical） 变量，预测一个观测的定性响应值的过程也称为分类（classification）。大部分的 分类问题都是先从预测定性变量取不同类别的概率开始，进而将分类问题作为概率 估计的一个结果，所以从这个角度看，分类问题与回归问题有许多类似之处。 根据定性因变量取值的特点，我们又可将其分为二元变量（binary variable）和多分 类变量（multinomial variable）。二元变量的取值一般为1和0，当取值为1时表示某 件事情的发生，取值为0则表示不发生。对于二元因变量，我们可考虑用Logistic模 型或判别分析法来处理，而对于多分类因变量，判别分析法的应用会更为广泛。本 章将分别对这两种方法进行介绍。

目录 Contents 一、问题的提出 二、Probit与Logistic模型 三、判别分析 四、分类问题评价准则 五、R语言实现

问题的提出 The Raise of Problems 1 问题的提出 The Raise of Problems

例1 新教学方法对学生成绩的影响 obs GRADE GPA TUCE PSI   1 2.66 20 17 2.75 25 2 2.89 22 18 2.83 19 3 3.28 24 3.12 23 .. … 为了考察一种新的经济学教学方法对学生成绩的影响，进行了调查，共得到了32 个样本数据。数据见上表。GRADE取1表示新近学习成绩提高，0表示其他； GPA是平均积分点；TUCE是以往经济学成绩；PSI取1表示受到新的经济学教学 方法的指导，0表示其他。假如想要了解GPA，TUCE和PSI因素对学生成绩是否 有影响，以及根据学生的GPA，TUCE和PSI能否预测学生成绩是否会提高？ 注：本例及例中的数据引自Greene（2000）第19章例19.1

Probit与Logistic模型Probit and Logistic Model 2 Probit与Logistic模型Probit and Logistic Model

（一）线性概率模型（Linear Probability Model, LPM） 考虑二元选择模型： (6.1) 其中， 是包含常数项的𝑘元设计矩阵， 是二元取值的因变量： 若假定 ，则总体回归方程为: (6.2)

进一步地，假设在给定 的时候，某一事件发生的概率 , 不发 生的概率为 ,即 ， ， 则因为 只取1和0两个值，所以其条件期望为： 进一步地，假设在给定 的时候，某一事件发生的概率 , 不发 生的概率为 ,即 ， ， 则因为 只取1和0两个值，所以其条件期望为： (6.3) 综合（6.2）式和（6.3）式可以得到： (6.4)

因此(6. 1)式拟合的是当给定自变量 的时候，某一事件发生（即 取值为1）的平均概率。在(6 因此(6.1)式拟合的是当给定自变量 的时候，某一事件发生（即 取值为1）的平均概率。在(6.4)式中，这一概率体现为线性的形式 ，因此(6.1)式称为线性概率模型（Linear Probability Model, LPM）。这 实际上就是用普通的线性回归方法对二元取值的因变量直接建模。 对于线性概率模型，我们也可以采用普通的最小二乘法进行估计， 但是会存在如下三个问题： (1)我们对（6.1）式进行的拟合，实际上是对某一事件发生的平均 概率的预测，即 。但是，这里的 值并不能保证 在0和1之间，完全有可能出现大于1和小于0的情形。

比如开头例中，GRADE关于GPA做一元线性回归，拟合的结果如图： (2)由于Y是二元变量，因此扰动项: 也应该是二元变量，它应该服从二 项分布，而不是我们通常假定的正 态分布。但是当样本足够多时，二 项分布收敛于正态分布。 线性概率模型拟合图

(3)在LPM中，扰动项的方差为： 因此，扰动项是异方差的。 由于存在着上述的诸多问题，因此对于二元定性因变量，一般 不推荐使用LPM，而是需要其他更为科学的方法。

（二）Probit与Logistic模型 在LPM中，通过适当的假设可以使得 的概率 与 是线性关系，即： 同时为了保证估计的概率的取值范围能在[0, 1]区间上，一个 直接的想法就是在外面套上分布函数 。 接下来介绍两种常用的分布函数：

(1)如果分布函数在 用标准正态分布函数 ，即 其中 是正态分布的分布函数，取值范围是[0, 1]。这时的概 率模型我们就称为Probit模型。 (2)如果分布函数在 取Logistic分布函数 ，则产生的概率 模型为Logit模型：

特别地对于Logit模型来说还可以写成： 其中 称为称为链接函数（link function） 称之为赔率 （odds ratio），即 取1的概率与取0的概率比值。

（三）基于潜变量模型的理解 二元选择模型也可以从潜变量回归模型的角度去解释，首先考察 以下模型： 其中， 是潜变量或隐变量（Latent Variable），它无法获得实 际观测值，但是却可以观测到 还是 。因此，我们实 际上观测到的变量是 而不是 。上式称为潜变量响应函数 （Latent Response Function）或指示函数（Index Function）。

如果我们假设： A1 ； A2 是i.i.d.的正态分布； A3 。 根据潜变量响应函数可得 的概率特征： 则当 为标准正态分布的概 率密度函数： 时， 可以写成： 这正是Probit模型。

如果我们假设： A1 ； A2 是i.i.d.的Logistic分布； A3 。 根据潜变量响应函数可以表示 的概率特征： 这正是Logit模型。

（四）最大似然估计（MLE） Probit和Logit模型的参数估计常用最大似然法。对于Probit或Logit 模型： 所以似然函数为： 对数似然函数为：

最大化这个对数似然函数𝑙𝑜𝑔𝐿的一阶条件为： 上式不存在显示解或封闭解，所以要用非线性方程的迭代的 方法进行求解。常用的有Newton-Raphson法或二次爬坡法 （Quadratic hill climbing）。

对于Probit模型来说，其边际效应为： （五）边际效应分析 对于Probit模型来说，其边际效应为： 对于Logit模型，其边际效应为： ，

从上面两式可以看出，在Probit和Logit模型中，自变量对 取值为 1的概率的边际影响不是常数，会随着自变量取值的变化而变化。所 以对于Probit和Logit模型来说，它们的边际影响不能像线性回归模型 直接等于系数。 这两个模型的边际效应分析常用的方法，是计算其平均边际效应 ，即对于非虚拟的自变量，一般是用其样本均值代入到两式，估计平 均边际影响。但是，对于虚拟自变量而言，需要先分别计算其取值为 1和0时的值，二者的差即为虚拟自变量的边际影响。 如果 则 随 的增加而增加，若 则 随 的增加而减少。

（六）似然比检验 似然比检验类似于检验模型整体显著性的F检验，原假设为全部自 变量的系数都为0，检验的统计量LR为： 其中𝑙𝑛𝐿为对概率模型进行MLE估计的对数似然函数值， 𝑙𝑛𝐿 0 为 只有截距项的模型的对数似然函数值，往往也称为空模型，即模型 中不包含任何自变量。当原假设成立时，LR的渐近分布是自由度 为𝑘−1（即除截距项外的自变量的个数）的 分布。

（七）预测 如果我们得到了系数的估计值 ，那么我们就可以预测出在给定 下， 的概率预测值。即， Probit: Logit: 也可以直接求链接函数的值：

Discriminant Analysis 3 判别分析 Discriminant Analysis

对于二元因变量，Logistic模型直接对 进行建模，即在给定自变量 下建立因变量 的条件分布模型。相反，判别分析采取的方法是先对每一个给定 的 建立自变量 的分布，然后使用贝叶斯定理反过来再去估计 。 何时使用判别分析？一方面，当类别的区分度较高的时候，或者当样本量 𝑛较小 且自变量 近似服从正态分布时，Logistic模型的参数估计会相对不够稳定，而判 别分析就不存在这样的问题；另一方面，在现实生活中，有很多因变量取值超过 两类的情形，虽然我们可以把二元Logistic模型推广到多元的情况，但这在实际应 用中并不常用。实际中对于因变量取多类别的问题，我们更常使用的是判别分析 法。 本节介绍判别分析的三种常用方法，包括朴素贝叶斯判别分析、线性判别分析和 二次判别分析。

（一） Naive Bayes判别分析 (1) 贝叶斯分类器 对于分类模型，我们的目的是构建从输入空间（自变量空 间） 到输出空间（因变量空间）的映射（函数）： ， 它将输入空间划分成几个区域，每个区域对应一个类别。区域的 边界可以是各种函数形式，其中最重要且最常用的一类就是线性 的。对于第 类，记 则第 类和第 类的判 别边界为 ，也就是所有使得 成立 的点。需要说明，实际上我们只需要 个边界函数。

为了确定边界函数，在构造分类器时，我们最关注的便是一 组测试观测值上的测试错误率，在一组测试观测值上的误差计算 具有以下形式： 其中 是用模型预测的分类变量， 是示性变量，当 时值等于1，说明测试值被误分，当 时值等于0，说明测试 值被正确分类。一个好的分类器应使上式表示的测试误差最小。

一个非常简单的分类器是将每个观测值分到它最大可能所在 的类别中，即给定 的情况下，将它分到条件概率最大的 类 中是比较合理的： 这类方法称为贝叶斯分类器，这种分类器将产生最低的测试 错误率，称为贝叶斯错误率，在 这一点的错误率为 ，整个分类器的贝叶斯错误率为：

(2) 贝叶斯定理 假设观测分成𝐾类，即因变量𝑌的取值为 ，取 值顺序对结果并无影响。设 为一个随机选择的观测属于因变量 𝑌的第𝑘类的概率，即先验概率； 表示第𝑘类观测的 𝑋的密度函数。则根据贝叶斯定理，我们可得观测 属于第𝑘类 的后验概率为:

我们记 ，因此只要估计出 和 就可以计 算 。通常 的估计比较容易，分别计算属于第𝑘类的样本占总样 本的比例，即 。 的估计则要复杂一些，除非假设它们 的密度函数形式很简单。 我们知道，贝叶斯分类器是将一个观测分到 最大的类 中，它在所有分类器中测试错误率是最小的（注意：只有当后验 概率中的各项假设正确，该结论才是对的）。如果我们找到合适 的方法估计 ，便可构造一个与贝叶斯分类器类似的分类方法。

（二） 线性判别分析（linear discriminant analysis, LDA） (1) 一元线性判别分析 我们在进行分类时，首先要获取 的估计，然后代入 第𝑘类的后验分布估计 ，并根据 的值，将观测分到值 最大的一类中。为获取 的估计，要先对其做一些假设。 通常假设 的分布是正态的，当 ，密度函数为一 维正态密度函数：

其中 和 分别为第𝑘类的均值和方差。再假设 ， 即所有𝐾个类别方差相同，均为 。将 带入 ： 贝叶斯分类器是将观测值 分到上式中 最大的一类。我 们对上式取对数，并将对 大小无影响的项去掉，整理式子， 不难看出贝叶斯分类器也将观测分到 最大的 一类。

例如 ，且 ，当 时，贝叶斯分 类器把观测分到第一类，反之分到第二类。此时贝叶斯决策边界 对应的点为 在实际中，即使确定𝑋服从正态分布，但总体参数 ， ， 我们仍然不知道，需要进行估计。线性判别分析方法 与贝叶斯分类器类似，计算 的值，常用的参数估计如下： ，

其中𝑛为随机抽取的样本数， 为属于第𝑘类的样本数， 的估计 值为第𝑘类观测的均值； 为第𝑘类观测的样本方差， 的估计值可以看作𝑘类样本方差的加权平均。 实际中，有时我们可以掌握每一类的先验概率， ，但当 信息不全时需要用样本进行估计，LDA是用属于第𝑘类的观测的比 例作为 的估计，即： 将 , 和 的估计值代入 ，即得线性判别分析的判别函数

LDA分类器将观测 分入 值最大的一类中。由于判别函 数 是关于𝑥的线性函数，所以称该方法为线性判别分析。 与贝叶斯分类器比较，LDA分类器是建立在观测都来自于均值 不同、方差相同的正态分布假设上的，将均值、方差和先验概率 的参数估计代入贝叶斯分类器便可得到LDA分类器。

(2) 多元线性判别分析 若自变量维度 ，假设 服从一个均值不同、 协方差矩阵相同的多元正态分布，即假设第𝑘类观测服从一个 多元正态分布 ，其中 是一个均值向量， 为所有𝐾 类共 同的协方差矩阵，其密度函数形式为： 通过类似于一维自变量的方法，我们可以知道贝叶斯分类 器将 分入 最大的一类。

同样，需要估计未知参数 ， 和 ，估计方法 与一维情况类似。同样的，LDA分类器将各个参数估计值带 入判别函数 中，并将观测值 分入 最大的一类。我 们发现多元情况下判别函数关于𝑥也是线性的，即可以写作 可以看到，对于类别𝑘和𝑙 ，决策边界 是关于 𝑥的线性函数，如果我们将 空间分成𝐾个区域，这些分割将 是超平面。

(3) 二次判别分析(quadratic discriminant analysis, QDA) 正如前面讨论的，LDA假设每一类观测服从协方差矩阵相 同的多元正态分布，但现实中可能很难满足这样的假设。二 次判别分析放松了这一假设，虽然QDA分类器也假设每一类 观测服从一个正态分布，并把参数估计代入贝叶斯定理进行 预测，但QDA假设每一类观测有自己的协方差矩阵，即假设 第𝑘类观测服从的分布为 ，其中 是第𝑘类观测的协 方差矩阵。此时，二次判别函数为：

QDA分类器把 、 、 的估计值代入上式然后将观测分入使 值最大的一类。我们发现判别函数 是关于𝑥 的二次函数，类别𝑘和𝑙 的决策边界也是一条曲线边界，这也是二次判别分析名字的由来。

面对一个分类问题时，我们该如何在LDA和QDA中做出选择呢？这其实是一个偏差-方差权衡的问题。这里不妨假设我们有𝑝 个自变量，并且因变量包含𝐾个不同类别，于是上述问题我们可以这样分析：

其次，若我们假设𝐾类的协方差矩阵相同，那么LDA模型对𝑥来说是线性的，这时候就只需要估计𝐾𝑝个线性系数。从这个角度看，LDA没有QDA分类器光滑，因此拥有更低的方差，所以说LDA模型有改善预测效果的潜力。 但是如果我们假设𝐾类的协方差矩阵相同不成立，那么LDA就会产生很大的偏差，需要在方差与偏差之间权衡。一般而言，训练数据相对较少时倾向于降低模型的方差，LDA比QDA更好而如果训练集非常大，我们更倾向于使用QDA，因为这时候𝐾类的协方差矩阵相同的假设是站不住脚的。

Criteria for Evaluation of Classification Problems 4 分类问题评判准则 Criteria for Evaluation of Classification Problems

在分类问题中需要基于一定的准则对分类器进行评价，以二分类问题为例，设 ，例如在信用卡用户的违约问题中，我们可以将“+”理解为“违约”，将“-”理解为“未违约”。如果与经典的假设检验进行结合，那么就可以将“-”看作零假设，而将“+”看作备择（非零）假设。

对于二分类问题，预测结果可能出现四种情况：如果一个点属于阴性（-）并被预测到阴性（-）中，即为真阴性值（True Negative，TN）；如果一个点属于阳性（+）但被预测到阴性（-）中，称为假阴性值（False Negative，FN）；如果一个点属于阳性（+）并且被预测到阳性中，即为真阳性值（True Positive，TP）；如果一个点属于阴性（-）但被预测到阳性（+）中，称为假阳性值（False Positive，FP）。可用下表的混淆矩阵来表示这四类结果。

预测分类 真实分类 于是，模型整体的正确率可表示为 ，相应的，整体错误率即为 。 混淆矩阵 真阴性值（TN） 假阳性值（FP） N - 或零 + 或非零 总计 真实分类 真阴性值（TN） 假阳性值（FP） N 假阴性值（FN） 真阳性值（TP） P N* P* 于是，模型整体的正确率可表示为 ，相应的，整体错误率即为 。

不过，很多时候我们更关心的其实是模型在每个类别上的预测能力，尤其是在不平衡分类（imbalance classification）问题下，模型对不同类别点的预测能力可能差异很大，如果只关注整体预测的准确性，模型很有可能将所有数据都预测为最多类别的那一类，而这样的模型是没有意义的。所以需要如下表所示评价指标来综合判断模型的准确性。

名称 定义 相同含义名称 分类和诊断测试中重要的评价指标 假阳性率 FP/N 第I类错误，1-特异度（1-specificity） 真阳性率 TP/P 1-第II类错误，灵敏度（sensitivity），召回率（recall） 预测阳性率 TP/P* 精确度，1-假阳性率 预测阴性率 TN/N*  

对于二分类模型，很多时候不是直接给出每个样本的类别预测，而是给出属于其中一类的概率，因此需要选取一个阈值，当预测概率大于阈值时，将观测预测为此，否则为另一类。不同的阈值对应不同的分类预测结果。对于不平衡数据可以通过ROC曲线来比较模型优劣。ROC曲线可以同时展示出所有可能阈值对应的两类错误，通过将阈值从0到1移动，获得多对FPR（1-specificity）和TPR （sensitivity）以FPR为横轴TPR为纵轴，将各点连接起来而得到的一条曲线。下面我们详细地介绍一下ROC曲线的原理。

第二个点(1,0) ，即FPR=1，TPR=0，类似地分析可以发现这是一个最差的分类器，因为它对所有样本的分类都是错误的。 首先考虑图中的四个点和一条线。 第一个点(0,1)，即FPR=0，TPR=1，这意味着FP（false positive）=FN（false negative）=0，表示我们得到了一个完美的分类器，因为它将所有的样本都正确分类。 第二个点(1,0) ，即FPR=1，TPR=0，类似地分析可以发现这是一个最差的分类器，因为它对所有样本的分类都是错误的。 ROC曲线

类似的，第四个点(1,1) ，表示此时的分类器将所有的样本都预测为正样本 （positive）。 第三个点 (0,0)，即FPR=TPR=0，这时FP（false positive）=TP（true positive）=0，可以发现此时的分类器将所有的样本都预测为负样本（negative）。 类似的，第四个点(1,1) ，表示此时的分类器将所有的样本都预测为正样本 （positive）。 ROC曲线

经过上述分析我们就可以明确，ROC曲线越接近左上角，说明该分类器的性能是越好的。另外，图中用虚线表示的对角线上的点其实表示的是一个采用随机猜测策略的分类器的结果，例如 (0.5,0.5)，它表示该分类器随机的将一半的样本猜测为正样本，另外一半的样本猜测为负样本。

接下来，我们讨论如何绘制出ROC曲线。我们知道，对于一个特定的分类器和数据集，我们显然只能得到一个分类结果，即一对FPR和TPR的结果，而要得到一条ROC曲线，实际上是需要多对FPR和TPR的值的。但是我们前面提到了，很多分类器的输出其实是一个概率值，即表示该分类器认为某个样本具有多大的概率属于正样本。对于后面章节中将介绍的其他分类器，例如第9章的基于树的分类器，第10章的支持向量机等，也都是可以通过某种机理得到一个概率输出的。

当我们得到了所有属于正样本的概率,将它们从小到大排列，从最小的概率值开始，依次将概率值作为阈值，当样本的输出概率大于等于这个值时，判定为正样本，反之为负样本。每次我们都可以得到一对FPR和TPR，即ROC曲线上的一点，最终我们就可以得到与样本个数相同的FPR和TPR的对数。注意，当我们将阈值设置为1和0时，就可以分别得到ROC曲线上的(0,0)和(1,1)这两个点，将这两个点与之前得到的这些（FPR, TPR）对连接起来，就得到了ROC曲线。选取的阈值个数越多，ROC曲线就越平滑。

不过，很多时候由ROC曲线并不能清晰地看出哪个分类器的效果更好，所以可以采用ROC曲线下面的面积（area under the ROC curve, AUC）来代表分类器的性能，显然这个数值不会大于1。前面提到，一个理想的ROC曲线会紧贴左上角，所以AUC值越大，分类器的效果就越好。另外不难理解，任何一个有效的分类器的AUC值都应该大于0.5。

5 R语言实现 The Realization of R

现在对本章开头例1的数据进行分析。在建模之前，求出数据的基本描述统计量，R程序如下（整理后的结果见下页）： （一）描述统计 现在对本章开头例1的数据进行分析。在建模之前，求出数据的基本描述统计量，R程序如下（整理后的结果见下页）： > grade <- read.table ( file = "grade.txt" , header = T ) > summarys <- function ( x ) { list ( mean = mean ( x ) , max = max ( x ) , min = min ( x ) , sd = sd ( x ) ) } # 自编一个求基本描述统计量简单函数 > summarys ( subset ( grade , PSI == 0 ) $ GRADE ) # subset ( )筛选PSI = 0的数据 > summarys ( subset ( grade , PSI == 1 ) $ GRADE ) > summarys ( grade $ GRADE ) > summarys ( subset ( grade , PSI == 0 ) $ GPA ) > summarys ( subset ( grade , PSI == 1 ) $ GPA )

数据的基本描述 变量 均值 最大值 最小值 标准差 GRADE PSI=0 0.166667 1.000000 0.000000   PSI=0 0.166667 1.000000 0.000000 0.383482 PSI=1 0.571429 0.513553 全部 0.343750 0.482559 GPA 3.101111 4.000000 2.630000 0.422985 3.137857 2.060000 0.533511 3.117188 0.466713 TUCE 21.55556 29.00000 12.00000 4.003267 22.42857 28.00000 14.00000 3.857346 21.93750 3.901509 PSI 0.437500 0.504016

（二）Logistic模型 R中可以用glm()函数拟合广义线性模型，包含Logistic回归中的probit模型和logit模型。glm()的形式与lm()类似，只是多了一些参数。函数的基本形式为： glm ( formula , family = family ( link = function ) , data = ) 其中，formula是模型表达式，与lm()的表达式一致。family参数用于设置模型的连接函数对应的分布族，比如gaussian分布，Poisson分布等，详见下表：

glm()的family参数代表的分布族 连接函数 binomial （link = "logit" 或 "probit" 或 "cauchit"） gaussian （link = "indentity"） gamma （link = "inverse" 或 "identity" 或 "log"） inverse.gaussian （link = "1/mu^2"） poisson （link = "log" 或 "identity" 或 "sqrt"） quasi （link = "identity" , variance = "constant"） quasibinomial （link = "logit"） quasipoisson （link = "log"）

与分析线性回归模型时lm()连用的许多函数在glm()中也有对应的形式，其中常用的函数见下表： 用途 summary() 给出拟合模型的信息摘要 coefficients() / coef() 列出拟合模型的参数 confint() 给出模型参数的置信区间 residuals() 列出拟合模型的残差值 anova() 生成两个拟合模型的方差分析表 plot() 生成评价拟合模型的诊断图 Predict() 用拟合的模型对原有数据进行拟合或者对新数据进行预测 aic() 计算拟合模型的AIC值

接下来，就开始建立模型。首先设定以下线性概率模型： 其中GRADE取1表示新近学习成绩提高，0表示其他；GPA是平均积分点；TUCE是以往经济学成绩；PSI取1表示受到新的经济学教学方法的指导，0表示其他。

用OLS估计这一线性概率模型的R程序和结果如下： > lpm <- lm ( GRADE ~ GPA + TUCE + PSI , data = grade ) > summary ( lpm ) Call: lm(formula = GRADE ~ GPA + TUCE + PSI, data = grade) Residuals: Min 1Q Median 3Q Max -0.781530 -0.277314 0.005306 0.210891 0.811449 Coefficients: Estimate Std. Error t -value Pr(>|t|) (Intercept) -1.49802 0.52389 -2.859 0.00793 ** GPA 0.46385 0.16196 2.864 0.00784 ** TUCE 0.01050 0.01948 0.539 0.59436 PSI 0.37855 0.13917 2.720 0.01109 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.3881 on 28 degrees of freedom Multiple R-squared: 0.4159, Adjusted R-squared: 0.3533 F-statistic: 6.646 on 3 and 28 DF, p-value: 0.001571

从分析的结果看，在5%的显著性水平上，PSI对GRADE的影响是显著的。也就是说，当GPA和TUCE都一样的情况下，接受过新的教学方法的学生与没有接受过新的教学方法的学生相比，学习成绩提高的概率要多0.3786。此外，GPA对成绩提高的边际影响是0.46，也就是说，在其他条件相同的情况下，GPA每增加1,学习成绩提高的概率是46%。

用Probit模型估计的R程序和结果如下： > grade.probit <- glm ( GRADE ~ GPA + TUCE + PSI , family = binomial ( link = "probit" ) , data = grade ) # Probit模型 > summary ( grade.probit ) Call: glm ( formula = GRADE ~ GPA + TUCE + PSI , family = binomial ( link = "probit" ) , data = grade ) Deviance Residuals: Min 1Q Median 3Q Max -1.9392 -0.6508 -0.2229 0.5934 2.0451 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -7.45231 2.57152 -2.898 0.00376 ** GPA 1.62581 0.68973 2.357 0.01841 * TUCE 0.05173 0.08119 0.637 0.52406 PSI 1.42633 0.58695 2.430 0.01510 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 41.183 on 31 degrees of freedom Residual deviance: 25.638 on 28 degrees of freedom AIC: 33.638 Number of Fisher Scoring iterations: 6

从分析结果看Probit模型的设定形式是： 其中 是标准正态分布的累积分布函数。将系数的估计结果代入得到估计的模型为：

下面我们进行模型的似然比检验和拟合优度: > install.packages ( "lmtest" ) > library ( lmtest ) # 需首先安装lmtest包，并载入包 > lrtest ( grade.probit ) # LR检验 Likelihood ratio test Model 1: GRADE ~ GPA + TUCE + PSI Model 2: GRADE ~ 1 #Df LogLik Df Chisq Pr(>Chisq) 1 4 -12.819 2 1 -20.592 -3 15.546 0.001405 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’

上面的检验结果中还给出了有关模型的似然比检验和拟合优度的信息。𝐿为Model1 LogLik的值，为-12 上面的检验结果中还给出了有关模型的似然比检验和拟合优度的信息。𝐿为Model1 LogLik的值，为-12.81880， 𝐿 0 为Model2 LogLik的值，为-20.59173，LR值即为Chisq，为15.546，它对应的𝑝值只有0.001405，表明模型整体是显著的。

(3)Logit模型估计 Logit 模型的R程序和估计结果： > grade.logit <- glm ( GRADE ~ GPA + TUCE + PSI , family = binomial (link = "logit" ) , data = grade) # 注意link设为logit > summary ( grade.logit ) Call: glm ( formula = GRADE ~ GPA + TUCE + PSI , family = binomial ( link = "logit" ) , data = grade ) Deviance Residuals: Min 1Q Median 3Q Max -1.9551 -0.6453 -0.2570 0.5888 2.0966 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -13.02135 4.93127 -2.641 0.00828 ** GPA 2.82611 1.26293 2.238 0.02524 * TUCE 0.09516 0.14155 0.672 0.50143 PSI 2.37869 1.06456 2.234 0.02545 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 41.183 on 31 degrees of freedom Residual deviance: 25.779 on 28 degrees of freedom AIC: 33.779 Number of Fisher Scoring iterations: 5

Logit模型的估计结果也给出了与Probit模型估计结果相似的似然比和拟合优度指标。但是对应的Logit模型为：

下面我们进行模型的似然比检验和拟合优度: > library ( lmtest ) > lrtest ( grade.logit ) # LR检验 Likelihood ratio test Model 1: GRADE ~ GPA + TUCE + PSI Model 2: GRADE ~ 1 #Df LogLik Df Chisq Pr(>Chisq) 1 4 -12.890 2 1 -20.592 -3 15.404 0.001502 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

𝐿为Model 1 LogLik 的值，为-12. 890， 𝐿 0 为Model 2 LogLik的值，为-20 𝐿为Model 1 LogLik 的值，为-12.890， 𝐿 0 为Model 2 LogLik的值，为-20.592，LR值即为Chisq，为15.404，它对应的𝑝值只有0.001502，因此，它是显著的，表明模型整体是显著的。 Probit和Logit模型中的回归系数与线性概率模型不同，并没有实际的经济意义。但可以计算自变量GPA 和TUCE对GRADE的平均边际影响。接下来我们进行计算：

有关的计算结果整理后见表: > coe <- coef ( grade.probit ) # 提取probit模型系数 > probit <- dnorm ( coe [ 1 ] + coe [ 2 ] * mean ( grade $ GPA ) + coe [ 3 ] * mean ( grade $ TUCE ) + coe [ 4 ] * mean ( grade $ PSI ) ) # 求probit模型平均边际影响 > ( m.gpa = coe [ 2 ] * probit ) # 求GPA平均边际影响 > ( m.tuce = coe [ 3 ] * probit ) > ( m.PSI = coe [ 4 ] * probit ) > coe.l <- coef ( grade.logit ) # 提取logit模型系数 > logit <- dlogis ( coe.l [ 1 ] + coe.l [ 2 ] * mean ( grade $ GPA ) + coe.l [ 3 ] * mean ( grade $ TUCE ) + coe.l [ 4 ] * mean ( grade$PSI ) ) # 求logit模型平均边际影响 > ( m.gpa.l = coe.l [ 2 ] * logit ) > ( m.tuce.l = coe.l [ 3 ] * logit ) > ( m.PSI.l = coe.l [ 4 ] * logit ) 有关的计算结果整理后见表:

Probit 和Logit 模型边际影响分析对比 变量 回归系数 平均边际影响 GPA 1.625810 0.5333 2.826113 0.5339 TUCE 0.051729 0.0170 0.095158 0.0180 PSI 1.426332 0.4644 2.378688 0.4493

上表中自变量GPA和TUC对因变量GRADE的边际影响是通过将相应的回归系数乘以 的值得到的。例如，对于Probit模型，GPA对GRADE的边际影响等于1.625810 × 0.3281 = 0.5333。但是，这一算法不适用于像PSI这一离散的自变量。对于Logit模型，PSI对GRADE的平均边际影响是PSI分别取值为1和0时，GRADE取值为1的概率差，即：

上表的边际影响分析取的是自变量的均值。但实际上自变量对因变量的影响是非线性的。例如，在Probit模型当中，PSI对GRADE的影响是随着GAP和TUCE取值的不同而不同的。假设TUCE取均值，则这一边际影响的函数为：

不同GPA水平下PSI对GRADE的影响 用各样本的GPA值代入这一边际影响函数，可以得到在不同的GPA水平下，PSI对GRADE的边际影响，见图。 不同GPA水平下PSI对GRADE的影响 （假设TUCE取均值）

上图横轴为GPA，纵轴为概率。图中折断线表示的是接受过新方法训练的学生（PSI = 1）成绩提高的概率曲线；点线是没有接受过新方法训练的学生（PSI = 0）成绩提高的概率曲线。可以看到，二曲线之间的差距并不是不变的，它开始随着GPA的提高而增大，但是，当GPA高于一定水平后，这一差距又开始缩小。图中的实线即为这两曲线的差值，它先是上升的，后又开始下降。

（三）判别分析 在这部分，我们同样用开头例1的数据进行判别分析。R中的MASS包提供了lda()和qda()函数分别做线性判别分析和二次判别分析，e1071包提供了naiveBayes()函数做朴素贝叶斯分类。 (1) 线性判别分析 我们首先使用MASS包中的lda()函数来拟合一个LDA模型。

> library ( MASS ) > lda.fit <- lda ( GRADE ~ GPA + TUCE + PSI , data = grade ) > lda.fit Call: lda(GRADE ~ GPA + TUCE + PSI, data = grade) Prior probabilities of groups: 0 1 0.65625 0.34375 Group means: GPA TUCE PSI 0 2.951905 21.09524 0.2857143 1 3.432727 23.54545 0.7272727 Coefficients of linear discriminants: LD1 GPA 1.91853948 TUCE 0.04340893 PSI 1.56574254

LDA输出表明 ， ，也就是说，约34.4%的观测对应着新近学习成绩提高。同时，LDA也输出了类平均值，即每类中每个自变量的平均，用来估计 。另外，线性判别系数输出给出了线性判别函数中GPA、TUCE和PSI的组合系数，用来形成LDA的决策准则。换句话说，该决策函数是由 各变量相乘表示的，若1.91853948 × GPA + 0.04340893 × TUCE + 1.56574254 × PSI很大，则LDA分类器预测新进学习成绩提高；若很小，则预测为其他。

> plot ( lda.fit ) 还可以通过plot()函数生成线性判别图像，可以通过对每个观测计算1.91853948 × GPA + 0.04340893 × TUCE + 1.56574254 × PSI获得: 线性判别图像

另外，可以使用predict()函数对LDA模型进行预测。predict()函数返回一个三元列表，第一个元素class存储了LDA关于新进学习成绩是否提高的预测；第二个元素posterior是一个矩阵，其中第𝑘列是观测属于第𝑘类的后验概率；最后，𝑥包含着线性判别。

例如这里我们对LDA模型的训练误差进行估计： > lda.pred <- predict ( lda.fit, grade ) > names ( lda.pred ) [1] "class" "posterior" "x" > lda.class <- lda.pred $ class > table ( lda.class , grade $ GRADE , dnn = c ( "Prediction" , "Actual" ) ) Actual Prediction 0 1 0 18 3 1 3 8 > mean ( lda.class != grade $ GRADE ) [1] 0.1875

我们对上面的数据用QDA模型进行拟合，可以通过MASS库中的qda()函数实现，它的语法与lda()函数一样。 (2)二次判别分析 我们对上面的数据用QDA模型进行拟合，可以通过MASS库中的qda()函数实现，它的语法与lda()函数一样。 > qda.fit <- qda ( GRADE ~ GPA + TUCE + PSI , data = grade ) > qda.fit Call: qda(GRADE ~ GPA + TUCE + PSI, data = grade) Prior probabilities of groups: 0 1 0.65625 0.34375 Group means: GPA TUCE PSI 0 2.951905 21.09524 0.2857143 1 3.432727 23.54545 0.7272727

从结果可以看出，QDA模型的输出同样包含类平均值，但不再包含线性判别系数，这是有因为QDA分类器是一个二次函数，不是自变量的线性函数了。 同样可以用predict()函数进行预测，使用语法与LDA一样。

估计得到的QDA模型训练误差为0.15625，略好于LDA模型。 > qda.pred <- predict ( qda.fit , grade ) > names ( qda.pred ) [1] "class" "posterior" > qda.class <- qda.pred $ class > table ( qda.class , grade $ GRADE , dnn = c ( "Prediction" , "Actual" ) ) Actual Prediction 0 1 0 18 2 1 3 9 > mean ( qda.class != grade $ GRADE ) [1] 0.15625 估计得到的QDA模型训练误差为0.15625，略好于LDA模型。

(3) Naive Bayes判别分析 最后，使用e1071包中的naiveBayes()函数来做朴素贝叶斯分类，它的语法与lda()和qda()均类似，不过这里要注意的是，使用naiveBayes()函数时，因变量必须是因子型的格式。

> library ( e1071 ) > bayes.fit <- naiveBayes( as.factor ( GRADE ) ~ GPA + TUCE + PSI , data = grade ) > bayes.fit Naive Bayes Classifier for Discrete Predictors Call: naiveBayes.default(x = X, y = Y, laplace = laplace) A-priori probabilities: Y 0 1 0.65625 0.34375 Conditional probabilities: GPA Y [,1] [,2] 0 2.951905 0.3572201 1 3.432727 0.5031320 TUCE Y [,1] [,2] 0 21.09524 3.780275 1 23.54545 3.777926 PSI Y [,1] [,2] 0 0.2857143 0.4629100 1 0.7272727 0.4670994

估计得到的Naïve Bayes模型的训练误差为0.15625，与QDA模型相同。 该函数有两项输出，其中A-priori probabilities给出该样本中每个类别出现的频率，Conditional probabilities给出了每个自变量在各个类别上的服从正态分布下的均值和标准差。同样可以用predict()函数进行预测。 > bayes.pred <- predict ( bayes.fit , grade ) > table ( bayes.pred , grade $ GRADE , dnn = c ( "Prediction" , "Actual" ) ) Actual Prediction 0 1 0 19 3 1 2 8 > mean ( bayes.pred != grade $ GRADE ) [1] 0.15625 估计得到的Naïve Bayes模型的训练误差为0.15625，与QDA模型相同。

（四）模型比较 我们将从两个方面，即模型的错分率和ROC曲线来分别对前面介绍的四种模型（Logit、LDA、QDA、Naïve Bayes）的分类效果进行比较。 (1)错分率 前面我们已经计算了LDA、QDA和Naïve Bayes三种模型在训练集上的错分率，这里我们再算一下Logit模型的错分率，并将结果整理见下表。

> logit. pred <- predict ( grade > logit.pred <- predict ( grade.logit , grade , type = "response" ) > logit.class <- rep ( 0 , nrow ( grade ) ) > logit.class [ logit.pred > 0.5 ] <- 1 #阈值设为0.5 > table ( logit.class , grade $ GRADE , dnn = c ( "Prediction" , "Actual" ) ) Actual Prediction 0 1 0 18 3 1 3 8 > mean ( logit.class != grade $ GRADE ) > [1] 0.1875

四种模型的错分率比较 模型 Logit LDA QDA Naïve Bayes 错分率 0.1875 0.15625 根据上表的结果可以发现，若单纯从错分率来看，模型的优劣情况为，QDA与Naïve Bayes相同，Logit与LDA相同，且前面一组一致优于后面一组。

接下来，加载ROCR包，并且编写一个rocplot()函数，使得每次只需要输入拟合概率值和相对应的真实类别就能自动画出ROC曲线。 R中的ROCR包可以用于生成ROC曲线。首先，为了能画出ROC曲线，需要将所有输出变为概率值，不同的模型概率输出的方式不一样，详见下面的代码： 接下来，加载ROCR包，并且编写一个rocplot()函数，使得每次只需要输入拟合概率值和相对应的真实类别就能自动画出ROC曲线。 > logit.pred2 <- predict ( grade.logit , grade , type = "response" ) > lda.pred2 <- predict ( lda.fit , grade ) $ posterior [ , 2 ] > qda.pred2 <- predict ( qda.fit , grade ) $ posterior [ , 2 ] > bayes.pred2 <- predict ( bayes.fit , grade , type = "raw") [ , 2 ]

> library ( ROCR ) > rocplot <- function ( pred , truth , ... ) { predob <- prediction ( pred , truth ) perf <- performance ( predob , "tpr", "fpr" ) plot ( perf , ... ) auc <- performance ( predob , "auc" ) auc <- unlist ( slot ( auc , "y.values" ) ) auc <- round ( auc , 4 ) #保留4位小数 text ( x = 0.8 , y = 0.1 , labels = paste ( "AUC =" , auc ) ) }

在上面所编写的函数中，prediction()和performance()是ROCR包中自带的函数，其中prediction()用于创建一个预测的对象，即将输入值转换为某种特定的格式；而对于performance()，它的第一个参数是前面prediction()生成的对象，后面若还设置了两个参数，则它们共同决定了得到的是什么图形，若只设置了一个参数，且是“auc”时，表示我们要输出的是AUC值。

最后，就可以通过调用rocplot()函数画出每种模型的ROC曲线了，四种模型的ROC曲线如图所示。

> par ( mfrow = c ( 2 , 2 ) ) > y <- grade $ GRADE > rocplot ( logit.pred2 , y , main = "Logit" ) > rocplot ( lda.pred2 , y , main = "LDA" ) > rocplot ( qda.pred2 , y , main = "QDA" ) > rocplot ( bayes.pred2 , y , main = "Naive Bayes" ) 由图形中给出的AUC值就可以得到各个模型的优劣排序为QDA大于Naïve Bayes大于LDA大于Logit，这与单纯从错分率的角度得到的结论有所不同。在实际应用中，我们更多的是使用ROC曲线来对分类器的性能进行评价。

谢 谢！ Thank You 方匡南 xmufkn@xmu.edu.cn 请扫描加微信