逻辑回归

逻辑回归是我们第一个用于进行分类的模型，是学习更先进的处理回归问题的技术的先导。其输出是定性的，由一组有限的值构成，称这些值为类别（class）。以二元分类为起点，然后拓展到多个类别的区分。

3.1 利用线性回归进行分类

假设利用线性模型，使预测值在[0, 1]，套用一个阈值，使得当模型输出值小于0.5，就预测类别为0，否则为1。
这样的模型不完美。

• 线性回归的原始输出会超出[0, 1]区间
• 线性回归解决的是最小化MSE的问题，对找到两个划分类作用不大。
  因此，不再通过最小化MSE来拟合函数。

3.2 逻辑回归入门

逻辑函数（logistic function）,对输入进行了非线性的变换，确保输出的范围（可以解释为输入属于类别1的概率）处于[0, 1]内。形式：

图形如下：

image.png

x =	y =
0	0.5
	1
	0

广义线性模型

逻辑回归 广义线性模型（generalized linear model, GLM)
GLM三个特点：

• 包含了输入特征的一个线性组合
• 输出被认为具有属于指数分布一类的相关概率分布, 有正态分布，泊松分布和二项分布。
• 输出分布的均值和输入特征的线性组合是通过一种函数的方式关联起来的，这种函数称为联系函数（link function)

在逻辑回归中，这些特点的结合展示如下。
输入特征的线性组合，创建一个x项：

注意，在逻辑回归里，建模对象是输出属于类别1的概率，而不是像线性回归里那样直接对输出建模。于是，我们不需要为误差项建模，因为输出就是一个概率， 它直接包含了模型中内在的随机性。

下面对这个项运用逻辑函数，以便产生模型的输出：

左边项直接告诉我们，根据输入特征的值，我们要计算输出属于类别1的概率。对于逻辑回归，它的输出的概率分布是伯努利分布。对于单次实验，它和二项分布是相同的，也是在一个具有两个常熟概率的输出的实验（如抛硬币）中会得到的分布。伯努利分布的均值就是成功的结果（本例中是类别1）的概率。因此前面的方程的左边也是相关输出的概率分布的均值。对输入特征的线性组合进行变换的函数，有时候称为均值函数，我们刚看到的这个函数就是逻辑回归的逻辑函数。
现在为了确定逻辑回归的联系函数，我们可以进行某些简单的代数运算，以便把输入特征的线性组合分离出来。






左边项被称为对数几率（log-odds)或逻辑函数（logit function），也就是逻辑回归的联系函数。在对数里面的部分，它的分母是给定数据的输出类别为0的概率。因此，对数里面这部分分别代表了类别1和类别0出现的概率比，也称几率比（odds ratio）。

3.2.2 解释逻辑回归中的系数

在逻辑回归里，特征的单位增量会导致几率比以的数量倍增。当为正，就是给几率比乘以一个大于1的数，因而我们就知道增大特征将有效地增大输出被标记为类别1的概率。类似的为负值，会让平衡向类别0
的预测转移。最后，当我们改变某个输入特征的值时，产生的效果是对几率比而非模型输出本身的倍乘，而输出被看作是预测类别为1的概率。从绝对数量上来说，由于输入的变化而导致模型输出的变化并非恒定，而依赖于输入特征的当前值。这一点和线性回归不同。在线性回归中，回归系数总是代表输入特征的单位增量在输出中所对应的固定增量。

3.2.3 逻辑回归的假设

相比于线性回归，我们不需要：

• 残差的正态性假设
• 同方差性假设
• 误差项互相独立， 也就是特征本身不再需要互相独立

但是如果特征表现出较高程度的多元共线性，我们的模型还是会遇到问题。
特征的比例缩放不会影响逻辑回归模型。逻辑回归具有最大似然的不变性（invariance property of maximum likelihood）的性质。

3.2.4 最大似然估计

线性回归的目标是把误差项的平方和最小化来得到系数。

逻辑回归是通过把数据的似然（likelihood）最大化来进行的。某条观测数据的似然就是在特定模型下看到该观测数据的概率。

在示例中，看到某条观测数据X属于类别1的似然就是由概率P(Y = 1|X)给出的，这个公式的形式是在本章前面的部分将结果的。因为只有两个类别，看到某条观测数据属于类别0的似然可由1- P(Y = 1|X) 给出。看到观测数据的完整数据集的整体似然则是每个数据点的所有单个似然的乘积，因为观测数据可以看作是独立获得的。由于每条观测数据的似然是以回归系数为参数，因此整体数据集的似然函数可以以这些系数为参数。可以把似然函数表示为一个方程，如下面的方程所示：

对数似然（log likelihood）

可以构建一个单个求和形式，进一步简化：

下一节会用R语言训练一个逻辑回归模型时，对此做一些练习。把似然最大化等同于把对数似然最大化，这两种方法会得到相同参数。

最大似然估计是参数拟合的基本技术，尽管在很多模型中流行，但是最大似然估计并非万能。

3.3 预测心脏病

数据来自UCI机器学习知识库的真实数据集Statlog(Heart)，其包含270个可能由心脏问题的病人的观测数据。其中120个病人心脏有问题，因此两个类别划分较为平均。我们的任务是根据病人的基本情况和一系列的医学检查来预测病人是否患有心脏病。首先是加载数据和对列名重命名：

heart <- read.table('YuCeFenXiRawData/heart.data', quote = '\"')
names(heart) <- c('AGE', 'SEX', 'CHESTPAIN', 'RESTBP', 'CHOL', 'SUGAR', 'ECG', 'MAXHR', 'ANGINA',
                  'DEP', 'EXERCISE', 'FLOUR', 'THAL', 'OUTPUT')

数据中各列名的含义：

列名	类型	定义
AGE	数值	年龄（以岁为单位）
SEX	二元类别	性别
CHESTPAIN	类别	胸痛类型，分4类
RESTBP	数值	静息血压(mmHg)
CHOL	数值	血清胆固醇（以mg/dl为单位）
SUGAR	二元类别	空腹血糖是否大于120mg/dl？
ECG	类别	静息心电图结果，分3个类别
MAXHR	数值	达到的最大心率（每分钟心跳次数）
ANGINA	二元类别	心绞痛是否由运动诱发？
DEP	数值	ST压低症状在运动中诱发和静息时诱发的比例
EXERCISE	有序类别	运动ST曲线的峰值斜率
FLUOR	数值	透视着色的主要血管数量
THAL	类别	地中海贫血的情况，分3类别
OUTPUT	二元类别	是否患有心脏病（输出）

要注意到把分类变量转化为因子类型, 把输出变成1表示由心脏病， 0表示没有心脏病：

heart$CHESTPAIN <- factor(heart$CHESTPAIN)
heart$ECG <- factor(heart$ECG)
heart$THAL <- factor(heart$THAL)
heart$EXERCISE  <- factor(heart$EXERCISE)
heart$OUTPUT <- heart$OUTPUT - 1

把数据分成训练集和测试集：

library(caret)
set.seed(987954)
heart_sampling_vector <- createDataPartition(heart$OUTPUT, p = 0.85, list = F)
heart_train <- heart[heart_sampling_vector,]
heart_train_labels <- heart$OUTPUT[heart_sampling_vector]
heart_test <- heart[-heart_sampling_vector,]
heart_test_labels <- heart$OUTPUT[-heart_sampling_vector]

现在测试集230条观测数据，测试机40条观测数据。在R中，用glm()函数训练逻辑回归模型，它代表了广义线性模型（generalized linear model）。该函数可以用来训练多种广义线性模型，但我们专注于逻辑回归模型，调用如下：

heart_model <- glm(OUTPUT ~ ., data = heart_train, family = binomial('logit'))
summary(heart_model)
[out]Call:
glm(formula = OUTPUT ~ ., family = binomial("logit"), data = heart_train)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.7137  -0.4421  -0.1382   0.3588   2.8118  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -7.946051   3.477686  -2.285 0.022321 *  
AGE         -0.020538   0.029580  -0.694 0.487482    
SEX          1.641327   0.656291   2.501 0.012387 *  
CHESTPAIN2   1.308530   1.000913   1.307 0.191098    
CHESTPAIN3   0.560233   0.865114   0.648 0.517255    
CHESTPAIN4   2.356442   0.820521   2.872 0.004080 ** 
RESTBP       0.026588   0.013357   1.991 0.046529 *  
CHOL         0.008105   0.004790   1.692 0.090593 .  
SUGAR       -1.263606   0.732414  -1.725 0.084480 .  
ECG1         1.352751   3.287293   0.412 0.680699    
ECG2         0.563430   0.461872   1.220 0.222509    
MAXHR       -0.013585   0.012873  -1.055 0.291283    
ANGINA       0.999906   0.525996   1.901 0.057305 .  
DEP          0.196349   0.282891   0.694 0.487632    
EXERCISE2    0.743530   0.560700   1.326 0.184815    
EXERCISE3    0.946718   1.165567   0.812 0.416655    
FLOUR        1.310240   0.308348   4.249 2.15e-05 ***
THAL6        0.304117   0.995464   0.306 0.759983    
THAL7        1.717886   0.510986   3.362 0.000774 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 315.90  on 229  degrees of freedom
Residual deviance: 140.36  on 211  degrees of freedom
AIC: 178.36

Number of Fisher Scoring iterations: 6

评估回归模型

glm()函数产生的逻辑回归模型摘要部分的格式类似于lm()函数为线性回归模型产生的摘要。它显示，对于二分类变量，二元特征比原始变量的水平数少一个，例如具有3个特征值的输入特征THAL产生了2个二元变量，分别标记为THAL6和THAL7。我们首先检查模型预测出的回归系数。它们和对应的z统计量一起显示。z统计量类似于t统计量。z的绝对值越大，对应的特征和输出变量相关的可能性越高。z统计量旁边的p值用概率形式表达了这个概念，并用星号和点标记，如同在线性回归里的形式一样，表明了包含对应p值的最小置信区间。
根据逻辑回归模型按照最大似然准则进行训练的事实，我们利用标准正态分布对系数进行显著性检验。例如要重现THAL7特征列出的z值3.362对应的p值，我们用以下代码：

pnorm(3.362, lower.tail = F) * 2
[1] 0.0007738012

安利了Larry Wasserman的一本参考书《统计学大全》，Springer出版社出版

根据模型摘要，我们看到FLOUR、CHESTOAIN4和THAL7对于心脏病是最强的特征预测因子。很多输入特征是有相对较高的p值。这表明了在有其他特征存在的情况下，它们很可能不是号的心脏病指标。我们再次强调正确解释这张表的重要性。比如，该表格并不是说心脏年龄就不是心脏病的良好指标，相反，模型摘要表示在有其他输入特征存在的情况下，年龄其实没有给模型带来什么贡献。此外，要注意，因为年龄的回归系数是负数，而我们期望发生心脏病的可能性应该随着年龄的增长而增长，所以几乎可以肯定在特征里存在某种程度的共线性。当然这个假设只有在其他所有输入特征都不存在的情况下才成立。实际上，如果重新训练一个只有AGE变量的逻辑回归模型我们会得到一个正的回归系数和一个较低的p值，着两者都支持我们对于特征存在共线性的判断：

heart_model2 <- glm(OUTPUT ~ AGE, data = heart_train, family = binomial('logit'))
summary(heart_model2)
[out]Call:
glm(formula = OUTPUT ~ AGE, family = binomial("logit"), data = heart_train)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5027  -1.0691  -0.8435   1.2061   1.6759  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -2.71136    0.86348  -3.140  0.00169 **
AGE          0.04539    0.01552   2.925  0.00344 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 315.90  on 229  degrees of freedom
Residual deviance: 306.89  on 228  degrees of freedom
AIC: 310.89

Number of Fisher Scoring iterations: 4

注意这个更简单的模型的AIC(赤池信息准则)值比我们在完整模型里得到的更大，可以预期这个简单模型的表现会更差。

模型的偏差

为了理解模型摘要的其它部分，我们就需要介绍一个重要的概念，它叫做偏差（deviance）。在线性回归里，残差的定义就是我们尝试预测的输出的预测值和实际值之差。逻辑回归是利用最大似然进行训练的。因此我们会很自然的认为，逻辑回归里和残差相似的概念会与似然有关。关于偏差的概念有一些紧密关联的定义。这里，我们要用glm()函数使用的一些定义来解释模型的输出。一条观测数据的偏差可以用-2乘以该观测数据的对数似然进行计算。一个数据集的偏差就是所有观测数据的偏差之和。一条观测数据的偏差残差（deviance residual）是从偏差本身派生出来的，和线性回归的残差是类似的。它们可以用如下公式进行计算：

对于观测数据，代表了偏差残差，而代表偏差。注意，对于一个偏差残差取平方有效的消除了正负号函数，并恰好产生了该观测数据的偏差。因此，偏差残差的平方和就是该数据集的偏差。它正好是数据集对数似然用常数-2进行比例缩放的结果。因此，把数据的对数似然最大化也就等同于把偏差残差的平方和最小化，这样我们和线性回归的类比就完成了。

为了重现我们在模型摘要里显示的结果，并且理解偏差是如何运算的。我们要在R里编写我们自己的一些函数，我们要先运用本章之前看到的对数似然方程对我们的数据集计算对数似然。根据方程创建两个函数。log_likelihoods()函数会根据模型预测的概率和实际的目标标签，为数据集的所有观测数据计算一个对数似然的向量，dataset_log_likelihood()会给这些值求总和，产生整个数据集的对数似然：

log_likelihoods <- function(y_labels, y_probs){
  y_a <- as.numeric(y_labels)
  y_p <- as.numeric(y_probs)
  y_a * log(y_p) + (1 - y_a) * log(1 - y_p)
}
dataset_log_likelihood <- function(y_labels, y_probs){
  sum(log_likelihoods(y_labels, y_probs))
}

下一步，利用偏差的定义计算两个类似的函数：deviances()和dataset_deviance()。第一个函数计算观测数据偏差的向量，第二个对整个数据集进行求和：

deviances <- function(y_labels, y_probs){
  -2 * log_likelihoods(y_labels, y_probs)
}
dataset_deviance <- function(y_labels, y_probs){
  sum(deviances(y_labels, y_probs))
}

有了这些函数之后，我们现在可以创建一个能计算模型偏差的函数。为此，我们需要利用predict()函数对训练集里的观测数据计算模型的概率预测。它的工作原理和在线性回归中是一样的，只是在默认的情况下，返回的概率是洛基量尺上的。为了确保获得实际的概率，我们需要给type参数指定response值：

model_deviance <- function(model, data, output_column){
  y_labels = data[[output_column]]
  y_probs <- predict(model, newdata = data, type = 'response')
  dataset_deviance(y_labels, y_probs)
}

要检查函数是否管用，让我们计算心脏病模型的偏差，它也被称为残余偏差（residual deviance）:

model_deviance(heart_model, data = heart_train, output_column = 'OUTPUT')
[1] 140.3561

image.png

如上图所示，这个值确实和模型摘要里面的值是相同的。评价逻辑回归模型的一种方法是计算模型偏差和空模型（即不采用任何特征训练出来的模型）偏差的差值。空模型偏差被称为空偏差（null deviance）。空偏差因为没有特征，会以一个常数概率预测类别1.这个概率的估算是根据训练集中类别为1的观测数据的比例而得到的，我们可以通过直接取OUTPUT列的平均值而得到：

model_deviance(heart_model, data = heart_train, output_column = 'OUTPUT')

null_deviance <- function(data, output_column){
  y_labels <- data[[output_column]]
  y_probs <- mean(data[[output_column]])
  dataset_deviance(y_labels, y_probs)
}

null_deviance(data = heart_train, output_column = 'OUTPUT')
[1] 315.9023

重现了模型摘要的Null deviance的值。残余偏差和空偏差的两个概念和我们在线性回归里看到的残差平方和与总平方和是类似的。如果两者之间的差值很小，类似于线性回归里残差平方和（RSS）主导解释了在输出变量中观测到的方差。要延续这个类比，可以给我们的模型定义一个伪（pseudo ）值，采用用来计算线性回归的统计量的相同公式，但是替换成用偏差来计算。在R语言里面可以用如下代码实现：

model_pseudo_r_squared <- function(model, data, output_column){
  1 - (model_deviance(model, data, output_column) / null_deviance(data, output_column))
}

model_pseudo_r_squared(heart_model, data = heart_train, output_column = 'OUTPUT')
[1] 0.5556977

这个结果说明，逻辑回归解释了大约56%的空偏差。这个数字不是特别大，最有可能的原因是我们的特征集还不够丰富，不足以用逻辑回归得出精确的预测。和线性回归的不同，伪是有可能大于1的，意味着残余偏差大于空偏差的疑难情况。如果出现了这种模型，就不应该信任模型，而是应该采用特征选择的方法，或尝试其他模型。

除了伪之外，我们还需要一种统计学检验来检查空偏差和残余偏差之间的差值是否显著。在模型摘要里面，残余偏差旁边缺少了p值说明R语言没有创建任何检验。实际上，残余偏差和空偏差之间的差值是近似并渐近的服从分布（即卡方分布）。我们要定义一个函数来计算这个差值的p值，不过要说明的是，这只是一个近似值。

首先，我们需要空偏差和残差之间的差值。还需要关于这个差值的自由度，他是通过空偏差的自由度数减去我们模型的自由度数得到的。空模型只有一个截距，因此它的自由度就是我们数据集里的观测数据总数减1。对于残余偏差，我们会计算一组回归系数，再加上截距，所以我们要把这些数量从观测数据的总数中减去。最后，利用pchisq()函数获得一个p值，注意，我们要创建的是上尾的计算值，因此需要设置lower.tail参数为FALSE。代码如下所示：

model_chi_squared_p_value <- function(model, data, output_column){
  null_df <- nrow(data) - 1
  model_df <- nrow(data) - length(model$coefficients)
  difference_df <- null_df - model_df
  null_deviance <- null_deviance(data, output_column)
  m_deviance <- model_deviance(model, data, output_column)
  difference_deviance <- null_deviance - m_deviance
  pchisq(difference_deviance, difference_df, lower.tail = F)
}
model_chi_squared_p_value(heart_model, data = heart_train, output_column = 'OUTPUT')
[1] 7.294219e-28

得到p值很小，所以我们感觉能确定模型产生的预期比平均预测的结果更好。在原来的模型摘要里面，我们看到了偏差残差的摘要。利用之前给出的偏差残差的定义，可以定义一个函数来计算偏差残差的向量：

model_deviance_residuals <- function(model, data, output_column){
  y_labes <- data[[output_column]]
  y_probs <- predict(model, newdata = data, type = 'response')
  residual_sign <- sign(y_labels - y_probs)
  residuals <- sqrt(deviances(y_labels, y_probs))
  residual_sign * residuals
}

最后对调用model_deviance_residuals()函数得到的偏差残差再调用summary()函数，得到一个数据表：

model_deviance_residuals(heart_model, data = heart_train, output_column = 'OUTPUT')
summary(model_deviance_residuals(heart_model, data = heart_train, output_column = 'OUTPUT'))
[out]    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
    -2.71369 -0.44214 -0.13823 -0.02765  0.35875  2.81178 

同样验证结果是正确的。模型摘要再最后还给我们提供了一个诊断，即费雪评分算法的迭代次数（Number of Fisher Scoring iterations），我们目前还没有将结果。这个数字通常再4到8之间，是一种收敛诊断。如果R语言用来训练逻辑模型的优化过程没有收敛，我们会预期得到一个较大的数字。如果发生了这种情况，模型是由疑问的，我们可能不能采用它进行预测。在示例中，这个数据在正常范围内。

测试集的性能

我们已经看到如何利用predict()函数来计算模型的输出。这个输出就是属于类别1的概率。我们可以运用一个阈值来进行二元分类。我们可以针对训练集和测试集进行预测，并用我们的预期输出与之比较，来衡量其分类的精确度：

train_predictions <- predict(heart_model, newdata = heart_train, type = 'response')
train_class_predictions <- as.numeric(train_predictions > 0.5)
mean(train_class_predictions == heart_train$OUTPUT)

test_predictions <- predict(heart_model, newdata = heart_test, type = 'response')
test_class_predictions <- as.numeric(test_predictions > 0.5)
mean(test_class_predictions == heart_test$OUTPUT)

与原著不一样

原著里面测试集精确度是0.9，与程序运行结果不一样

训练集里精确度接近88.7%，测试集里是77.5%。模型的系数表显示，由几个特征看起来不显著，我们也看到一定程度的共线性，这意味着现在可以接着进行变量选择，并通过计算或获得病人的附加数据来找到更多特征。伪的计算让我们看到，在模型中并没有解释足够多的偏差，也对前面的结论提供了证据。