逻辑回归及其他分类变量分析方法

理论学习 / 线性模型
字数统计: 6.3k   |   阅读时长≈ 23 分钟

ISL 真是好书！

资料

An Introduction to Statistical Learning，下文简称 ISL

概述

线性回归的因变量都是连续变量（例如人的身高），但是生活中我们也有很多因变量是分类变量的情况。比如：

  1. 一位来到急诊室的病人，具有很多的临床症状，这些症状可能是由三种疾病导致的，那么到底这个病人是患上了什么疾病呢？
  2. 邮箱区分哪些邮件是垃圾邮件，哪些邮件不是垃圾邮件
  3. 全基因组关联分析：确定哪些基因突变与疾病有关联

在 ISL 书中，分析的例子是基于一个人的年收入和每月信用卡余额来预测这个人是否会拖欠信用卡付款。数据如下图所示

1

为什么不用线性回归？

假设我们想要基于一个急诊室的病人的症状来预测他/她患了什么病。在这个简化的例子中，假设我们有三个可能的疾病：stroke, drug overdose, epileptic seizure 。 我们可以将这些重编码为数字，如下：

1

如果这里我们用最小二乘的方式来预测 Y，并且我们认为这三者之间的差异是相同的。实际上我们这么重编码数字是没有道理的。比如，有些人还可能用另外一种编码方式：

1

如果换了一种重编码方式，那么这三者之间的关系就完全变了。

如果反应变量的值是有序的，例如 mild, moderate, severe ，这样我们把这几个变量编码为1,2,3 还有一点道理。但是这么做还不是一种自然的方式。

对于一个二变量的分类性状，情况就好一点。比如如果病人的病症只有两种情况：stroke 和 drug overdose 。 我们可以将这两个变量转换为哑变量。

1

我们可以对这个二变量进行线性回归，如果预测值大于0.5， 我们认为是 drug overdose；反之我们认为 stroke 。在这种情况，即使颠倒了重编码方法，预测结果还是一样的。这种情况下，线性回归是有意义的，我们得到的预测值可以理解为 drug overdose 的可能性，但是线性回归的预测值可能会超出 [0,1] 区间，使得这个结果难以解释。

衡量分类变量模型的准确性

training error rate : 公式如下，其中 I() 为哑变量，如果预测值和真实值不同则为1，相同则为0。这个式子其实就是计算 traing data 有多少比例的预测值预测错了。

但是，实际上我们更关心的是 test data 的结果，同样的我们计算 test error rate 。一个好的模型，我们得到的 test error rate 应该是很低的。

逻辑回归

逻辑回归（logistic 回归）采用的因变量是某一个水平的概率。

线性回归与逻辑回归比对的结果如下：

1

对于 ISL 的 Default 数据集（信用卡是否违约）， 逻辑回归会预测违约的概率。例如，基于 balance 因子的预测可以写成

1

简写为 p(balance)，如果这个值大于0.5，我们可能预测会违约（这个阈值也可以调整）。

逻辑回归模型

常规的线性模型如下：

1

线性模型的问题在于预测值会超出[0,1]区间，这是不合理的。为了避免预测值超出[0,1]区间，我们设定逻辑回归的预测函数为：

1

个人理解：β0+β1X 就是原来的线性模型的预测函数，范围是负无穷到正无穷。因此，eβ0+β1X 的范围就是 (0,+∞] 。整个 p(X) 在这个范围内是递增函数（由求导得到，易知导数永远大于0），整个p(x) 的两个极值在自变量的两个端点取得，范围为 (0,1)

为了求解这个预测函数，我们需要用到最大似然法。

对上面的函数进行变换，我们得到下面这个式子。这个式子 p(x)/[1-p(x)] 称为 odds ，这个式子的值的范围为 (0, ∞)

1

如果我们对左右两边取对数，我们可以得到下面这个式子：下式中左侧称之为 log-odds 或者 logit 。我们可以看到 logit 是符合线性模型的。

1

我们这里估计得到的 β1 的值虽然和 p(x) 不是线性相关，但是如果β1 是正数，说明增加 X 的值，会提高 p(x)。

估计回归系数

逻辑回归中的 β0 和 β1 是未知的，必须通过已有的 training data 来进行估计。在线性模型中，我们采用最小二乘法来估计线性回归的系数。虽然这里我们也可以用最小二乘来拟合模型，当我们偏向采用更一般的方式 - 最大似然法。最大似然法用于拟合逻辑回归的基本思路是：我们想寻找一组 β0 β1 的值，使得每个个体p(x) 的预测值和观测值尽可能的接近。在这个例子中，也就是对于存在违约的个体的预测值要接近1，对于不存在违约的个体的预测值要接近于0。

我们可以把上面这个思路用一个 似然函数 来表示，如下图所示：

1

使得这个似然函数最大的一组 β0 和 β1 的值就是我们选中的估计值。

在ISL书中，最大似然法是一种估计非线性模型参数的一种普遍的方法。在线性模型中，最小二乘法实际是最大似然法的一个特例（理由呢？）。

回到信用卡的例子，如果我们采用逻辑回归，结果如下图。我们可以看到 β1 的估计值为正数，说明 随着 balance 的增加，违约的概率在增加。

1

预测

一旦我们得到了系数的估计值，我们就可以带入逻辑回归的公式进行预测。例如，加入我们要预测 balance 为 1000 时的违约概率：

1

如果自变量是分类变量也是可以的，不过需要先转变为哑变量。比如下面是“是否是学生” 的分析结果，我们可以看到学生违约的风险高于非学生群体。

预测值如下

1

多重逻辑回归

如果我们多个自变量，我们可以将逻辑函数扩展如下：

1

1

之后我们还是采用最大似然法来估计这些系数。

比如，我们用三个自变量去预测信用卡违约，结果如下，我们可以看到 balance 和 student[Yes] 的 P 值都很小，说明这这两个自变量都是显著的。但是stdent[Yes] 的系数是负数，说明学生违约的概率更小，这和上面的简单逻辑回归的结论相反。

1

下图就说明了这个矛盾的结论：左图说明了，如果固定了 Balance 因子，既在相同的 Balance 水平下，学生的违约风险更低；右图中是不管Balance 和 income 水平，计算所有的学生和非学生的结果，可以看到学生的违约风险更高。这对于信用卡发放的决策者而言这是一个重要的区别，如果不知道balance信息，学生的违约风险比非学生更高； 但是如果在相同的balance水平下，学生比非学术的违约风险要低。

1

发生这种矛盾的原因还是自变量之间的相关，这种现象称之为 confounding 。是否是学生 与 balance 绝对是存在相关关系的，或者说学生的balance 的 分布 和 非学生的 balance 分布是不一样的。

我这里可以通过画图来看一下，代码和图片结果如下。我们可以看到学生和非学生的 balance 的分布基本还是类似于正态分布，但是学生的均值更高，而且右侧的尾巴更肥。我们可以通过这个图简单地总结一下，就是学生群体balance的值倾向于比非学生群体更高。

1
2
library(ggplot2)
ggplot(Default, aes(x=balance, fill=student)) + geom_histogram(position = 'identity', alpha = 0.4, bins=60)

1

大于2个水平的因变量

如果有多个水平的话，一般就不用逻辑回归了。我们会用 判别分析 方法。

线性判别分析

对于大于2个水平（K>2）的因变量，我们假定我们获取因变量的先验分布（即一个随机挑选的观察值来自于其中某个水平的概率），记为下面的符号。

另外一个式子就是基于某个因变量水平的自变量分布，见下式。如果这个值很高，说明当Y=k时，X=x 的概率很大。

我们通过贝叶斯定理，就可以得到基于某个X的Y的条件概率，即后验概率，我们用 pk(x) 来表示

在这个公式里，因变量的先验分布比较容易获得，我们先对总体进行随机抽样，之后算一下每个水平的比例即可。但是 fk(x) 就不是很容易获得。如果我们可以估计 fk(x) ，我们就可以计算不同因变量水平的后验概率，我们认为后验概率最大的水平就是预测值。

p=1的线性判别分析

如果只有一个自变量，我们需要首先估计 fk(x) ，我们需要先做一些假设。

我们假设 fk(x) 是符合正态分布的 （假设1），在一个自变量的情况下，即符合下面这个式子：

其中 uk 和 σ2k 分别为第K个水平的均值和方差。我们进一步假设所有的水平的方差均相同（假设2），即如下式，可以简化为 σ2 。

将这个式子带入到上面计算 pk(x) 的式子中，得到

计算所有水平的后验概率，选择后验概率最大的水平。我们可以看到对所有的水平，分母部分都是一样的，因此我们只需要关注分子部分，只要找对分子最大的水平。对分子部分采用 log 函数，剔除常数部分，我们得到下面这个式子。这个式子值最大的水平，就是后验概率最大的水平

假设 K=2，并且 π1 = π2 。那么如果满足下面这个式子，则最终分配为水平1；否则则为水平2

在这种情况下，贝叶斯决策线 为

但是在实际中， πk , μk ， σ2 这些参数都是未知的，我们需要先估计这些常数，再估计 pk(x)

这里 n 是所有的观测值数目，nk 是观测值为水平k 时的数目。至于 πk ，有先验就用，没有就从 training data 计算一下

将上面的参数估计值带入下式，计算得到值最大的水平

为什么判别分析称为线性判别分析呢？因为 discriminant functions δk(x) 是 x 的线性函数。

p>1的线性判别分析

我们现在看多个自变量的情况，并且这些自变量服从多元正态分布（每个变量均服从一维正态分布）， 下图为两个变量时服从多元正态分布的概率密度图。

1

为了表示 p维的随机变量X 服从多元正态分布，我们写作

其中，μ 为均值，是一个p维的向量。

方差分布为 Σ , 为 p × p 的矩阵。

最后，多元正态分布定义为

类似于p=1的判别分析，这里假定因变量所有水平的观测值均服从方差相同的多元正态分布。

p>1的判别函数如下

以下图为例，作图是三种水平所有观测值结果，圆圈为95%的概率的范围，虚线为画出的贝叶斯决策线；右图为各随机抽取了20个样本的结果，实线为抽样结果的贝叶斯决策线。

个人理解：这张图给了我一个提示，分类变量的每个水平必须是能几乎截然分开的，才能进行分析。不然对于判别分析，无法画出决策线。所以我感觉进行分析前有必要先进行一下可视化，看一下分类变量的各个水平的设置是否合理，有没有两个水平实际是差不多同一个水平的情况。

1

这里我们同样需要估计这些参数，类似于p=1的情况

从判别函数的形式看，判别函数依然是 X 的线性转换，这就是线性判别分析中线性的由来。

如果对信用卡违约的数据进行判别分析，因变量是是否违约，自变量是 credit card balance 和 student status ，通过对 10,000 个 traing samples 进行建模，我们得到的 training error rate 的结果是 2.75% 。这个结果看上去很好，但是有两个注意事项：

  1. 首先，training error rate 通常会低于 test error rate，也就是说 test data 的结果通常会更差一点，这是因为我们是根据 training data 来拟合的模型参数。如果参数 p 除以 样本 n 的比例越高，我们可能越可能遇到过拟合问题。
  2. 其次，由于 training data 中只有3.33%的个体违约，因此一个最简单的 classifier 可以将所有个体都预测为不会违约，无论这个人的 credit card balance 和 student status 是什么情况，这个 classifier 的 trainging error rate 就是 3.33%。这种情况下仅仅比上面的 2.75% 略高。

在实际情况中，一个二分类的 classifier 会产生两种错误：它可能将一个不违约的个体判定为违约，或者将违约的个体判定为不违约。这两种错误发生的原因，通常是需要关注的地方。我们一般会通过 confusion matrix 来查看一下，我们可以看到在真实为NO的个体中（9667），总共有23个个体是错误预测的，这看上去错误概率很低。但是，在总共333个违约个体中，总共有252（75.7%）的个体被错误预测为No。因此，虽然总的 training error rate 很低，但是针对违约的个体的error rate 很高。由于信用卡公司做预测的目的是为了确定高风险可能违约的用户，这种预测情况肯定是无法接受的。

1

sensitivity 和 specificity

有两个概念可以说明一个 classifier 的表现， sensitivity 和 specificity 。在这个例子中，sensitivity 就是违约者被正确预测的概率，这里是24.3%；specificity 就是非违约者被正确预测的概率，这里是 99.8%。

为什么判别分析在违约个体的表现如此糟糕呢？或者说，为什么 sensitivity 如此低？如我们所见，判别分析会采用贝叶斯方法，会含有最低的 total error rate ，也就是说会产生最少的总的错误分类的数目，无论这些错误分类原本来自于哪个分类。但是基于信用卡公司的角度，肯定是希望尽量避免将会违约的个体判定为不会违约的错误，相比之下，将不会违约的个体判定为违约问题则没有那么严重。因此我们需要需改这个方法，来新建一个 classifier 来更好地符合信用卡公司的需求。

判别分析会将个体判定为后验概率最大的分类。在只有两种分类的情况下，一个观测值会判定为违约的公式如下（后验概率大于0.5）

这里默认阈值就是0.5，但是，如果我们特别关心将违约个体预测为非违约的情况，那么我们就可以降低这个阈值。例如如果后验概率超过0.2，我们就判定为违约，公式如下

修改阈值后，预测结果如下，我们可以看到相比于上面的结果，违约个体预测错误下降，但是非违约个体的预测错误上升。总的预测错误轻微上升至3.73%。

1

下图说明了修改违约后验概率阈值的 trade-off 结果。如果阈值设定为0.5，那么总的错误概率是最低的，如下图的黑线，但是违约个体的错误概率却是最高的，如蓝色虚线所示。

1

ROC 曲线是用来展现不同的阈值下两种错误的错误概率的变化。下图就是 ROC 图 ，纵坐标是 sensitivity ，就是违约者被正确预测的概率；横坐标为 1-specificity ，为非违约者被错误预测的概率。这幅图不显示阈值。

1

一型错误和二型错误

我们把上面的表格进行一般化，用 - 或 Null 表示零假设（本例为不违约），+ 或 Non-null 表示备择假设（本例为违约）。

1

根据上面的表格，我们得到一些重要的统计量。

1

二次判别分析

线性判别分析假设所有水平的观测值均假设服从一个方差矩阵相同的多元正态分布。二次判别分析则假设所有水平的观测值均服从一个正态分布，但是不像线性判别分析，二次判别分析假设所有水平的方差矩阵各不相同。即对第k个水平的观测值的分布为：

在这个假设下，判别函数如下。这个函数是 x 的二次函数，这就是二次判别分析名字的来源。

我们应该在什么情况下采用线性判别分析，什么情况下采用二次判别分析呢？这个答案在 bias-variance trade-off 中。二次判别分析需要估计的参数比线性判别分析多得多，如果我们有 p 个预测因子，那么估计线性分析的方差矩阵需要 p(p+1)/2 个参数；而二次判别分析需要对每个水平均要估计一个单独的方差矩阵，则总共需要 Kp(p+1)/2 个参数。因此二次判别分析更灵活，bias 比较低，但 variance 比较高。一般来说，当样本数量比较少时，推荐使用线性判别分析来降低 variance；如果样本数量非常大或者严重偏理所有水平同方差的假设，更推荐使用二次判别分析。

下图说明了线性判别分析和二次判别分析在两种情况下的表现。

1

不同方法的比较

逻辑回归和线性判别分析虽然思路不同，但却是紧密相关的。在线性判别分析的框架里，log odds 计算公式可简化如下（就是两个判别函数相减）：

根据逻辑回归公式，log odd 可写为：

这两种方法得到的 log odd 都是 x 的线性函数，因此这两种方法的决策边界均为直线。这两种方法唯一的区别在于 β0 和 β1 是通过最大似然得到的，而 c0 和 c1 是通过估计正态分布的均值和方差得到的。

由于这两种方法仅仅是拟合参数的过程有区别，因此这两种方法的结果往往是非常相似的，但偶尔可能情况不一样。线性判别分析由于存在同方差正态分布的假设，因此如果数据符合这个假设，那么结果可能比较好。相反，如果数据不符合这个假设，那么逻辑回归的结果会更好。

而 KNN 方法但预测 X = x 的值时，会提取x附近的多个观测值，然后分配X为大多数观测值的结果。因此 KNN 方法是一种完全的非参数方法，没有对决策界限的任何假设。因此，如果只有当决策线是高度非线性的情况下，KNN 方法才可能由于线性判别方法和逻辑回归方法。另一方面，KNN 方面没有告诉我们哪些预测是重要的，哪些是不重要的。

最后，二次判别分析是 KNN 方法 和 LDA/逻辑回归方法 之间的折中方法，它的决策线是二次函数。

为了说明这四种方法的表型，ISL 书中模拟了六种情况，在前三种情况下，真实的决策线就是线性的，另外三种情况则是非线性的。对于每种情况下，均随机生成 100 个训练集，我们在每个训练集中拟合模型，然后计算一个很大的 test data 数据集的 test error rate 。

真实决策线是直线的三种情况画图如下：

1

真实决策线不是直线的三种情况画图如下：

1

其中，KNN 方法做了两次，一次 K = 1，一次 K 值通过 cross validation 获得。在所有的情况中，p 均等于2。

我们可以看到，如果真实的预测线是线性的，那么 LDA和逻辑回归方法较好；如果是非线性的，那么 QDA 方法更好。最后，对于非常复杂的预测线，那么非参方法，例如KNN方法更好，但是需要慎重选择平滑度。

最后，在线性回归中，我们可以通过对预测因子进行转换，来实现对非线性关系的模拟。我们同样可以在分类变量的分析中采用相同的操作，例如采用 X2 , X3 等加入到逻辑回归 和 LDA分析中，这样可以实现类似与 QDA 的效果。

R代码

ISl 使用的数据是股票数据，Smarket data ，来自于 ISLR 包。

1
2
3
4
library(ISLR)
names(Smarket)
dim(Smarket)
summary(Smarket)

使用 cor() 函数可以计算所有数据集的相关矩阵，所有变量必须为连续变量。

我们可以看到这里只有 Year 和 Volume 这两个变量之间存在实质相关。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
> cor(Smarket)
Error in cor(Smarket) : 'x'必需为数值
> cor(Smarket[,-9])
             Year         Lag1         Lag2         Lag3         Lag4         Lag5      Volume
Year   1.00000000  0.029699649  0.030596422  0.033194581  0.035688718  0.029787995  0.53900647
Lag1   0.02969965  1.000000000 -0.026294328 -0.010803402 -0.002985911 -0.005674606  0.04090991
Lag2   0.03059642 -0.026294328  1.000000000 -0.025896670 -0.010853533 -0.003557949 -0.04338321
Lag3   0.03319458 -0.010803402 -0.025896670  1.000000000 -0.024051036 -0.018808338 -0.04182369
Lag4   0.03568872 -0.002985911 -0.010853533 -0.024051036  1.000000000 -0.027083641 -0.04841425
Lag5   0.02978799 -0.005674606 -0.003557949 -0.018808338 -0.027083641  1.000000000 -0.02200231
Volume 0.53900647  0.040909908 -0.043383215 -0.041823686 -0.048414246 -0.022002315  1.00000000
Today  0.03009523 -0.026155045 -0.010250033 -0.002447647 -0.006899527 -0.034860083  0.01459182
              Today
Year    0.030095229
Lag1   -0.026155045
Lag2   -0.010250033
Lag3   -0.002447647
Lag4   -0.006899527
Lag5   -0.034860083
Volume  0.014591823
Today   1.000000000

逻辑回归

glm() 函数使用广义线性模型，包括一系列模型，其中就有逻辑回归。这里我们通过指定 family = binomial 来指定使用逻辑回归。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
> glm.fit = glm(Direction ~ Lag1+Lag2+Lag3+Lag4+Lag5+Volume, data=Smarket, family=binomial)
> summary(glm.fit)

Call:
glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
    Volume, family = binomial, data = Smarket)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.446  -1.203   1.065   1.145   1.326  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.126000   0.240736  -0.523    0.601
Lag1        -0.073074   0.050167  -1.457    0.145
Lag2        -0.042301   0.050086  -0.845    0.398
Lag3         0.011085   0.049939   0.222    0.824
Lag4         0.009359   0.049974   0.187    0.851
Lag5         0.010313   0.049511   0.208    0.835
Volume       0.135441   0.158360   0.855    0.392

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1731.2  on 1249  degrees of freedom
Residual deviance: 1727.6  on 1243  degrees of freedom
AIC: 1741.6

Number of Fisher Scoring iterations: 3

查看系数

1
2
3
4
5
6
7
8
9
10
11
12
> coef(glm.fit)
 (Intercept)         Lag1         Lag2         Lag3         Lag4         Lag5       Volume 
-0.126000257 -0.073073746 -0.042301344  0.011085108  0.009358938  0.010313068  0.135440659
> summary(glm.fit)$coef
                Estimate Std. Error    z value  Pr(>|z|)
(Intercept) -0.126000257 0.24073574 -0.5233966 0.6006983
Lag1        -0.073073746 0.05016739 -1.4565986 0.1452272
Lag2        -0.042301344 0.05008605 -0.8445733 0.3983491
Lag3         0.011085108 0.04993854  0.2219750 0.8243333
Lag4         0.009358938 0.04997413  0.1872757 0.8514445
Lag5         0.010313068 0.04951146  0.2082966 0.8349974
Volume       0.135440659 0.15835970  0.8552723 0.3924004

线性判别分析

我们可以用 MASS 包的 lda() 函数进行LDA分析。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
> library(MASS)
> lda.fit = lda(Direction ~ Lag1 + Lag2, data=Smarket)
> lda.fit
Call:
lda(Direction ~ Lag1 + Lag2, data = Smarket)

Prior probabilities of groups:
  Down     Up 
0.4816 0.5184 

Group means:
            Lag1        Lag2
Down  0.05068605  0.03229734
Up   -0.03969136 -0.02244444

Coefficients of linear discriminants:
            LD1
Lag1 -0.7567605
Lag2 -0.4707872

这里提供了先验概率的估计值，每个水平均值的估计值。最后的 Coefficients of linear discriminants 提供了两个因子的线性组合，用于形成决策线。在本例中，就是下式

如果这个式子很大，那么就会预测为 Up，否则就预测为 down 。

二次判别分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
> library(MASS)
> qda.fit = qda(Direction ~ Lag1 + Lag2, data=Smarket)
> qda.fit
Call:
qda(Direction ~ Lag1 + Lag2, data = Smarket)

Prior probabilities of groups:
  Down     Up 
0.4816 0.5184 

Group means:
            Lag1        Lag2
Down  0.05068605  0.03229734
Up   -0.03969136 -0.02244444

这里没有提供系数。

打赏
  • 版权声明： 本博客所有文章除特别声明外，著作权归作者所有。转载请注明出处！
  • Copyrights © 2019-2024 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信