github-actions[bot]
commented
4 years ago 广义线性模型应用举例之logistic回归及R计算 by 优雅R
logistic回归概述
前文曾提到的一系列线性或非线性回归类型，如线性回归、指数回归、泊松回归、负二项回归等，响应变量均是数值型变量（连续型或计数型）。事实上很多回归方法并不局限于对数值型响应变量的分析，还可应用于类别型响应变量的建模。
当回归中的响应变量是类别型时，与其说是回归，更普遍的名称是分类，归因于此时的模型实际上执行了带监督的分类器功能。在前文也已经介绍过类似的例子，回归和分类的联系与区别。如用于数值型响应变量的随机森林回归和用于类别型响应变量的随机森林分类，用于数值型响应变量的偏最小二乘回归和用于类别型响应变量的偏最小二乘判别分析等。
在广义线性模型中，logistic回归就是一类对类别型响应变量进行建模的方法，其假设响应变量服从二项分布或多项分布，进而通过极大似然估计获得模型参数。
二项分布和二分响应变量的logistic回归
二分响应变量可以用取值为0或1的变量来表示，用伯努利分布或二项分布来模拟。例如，为了调查某物种的分布，在每个地点只调查一次，可记录为1（该物种存在）或0（该物种不存在）。如果对每个地点调查n次，所得到的数据就记录了观察到该物种的次数。从统计学角度看，n次调查可视为试验次数，每次试验中如果观测到该物种就视为成功，反之视为失败，则可以通过二项分布描述成功次数k的概率模型：
式中参数p代表了试验成功的概率，y的期望值是np，方差是np(1-p)。
伯努利分布可认为是二项分布在n=1时的特例，因此下文统称为二项分布。
二项分布是指数分布族的一员，推广到广义线性模型中就是logistic回归，其似然度函数通过上述二项概率分布函数来定义。若响应变量Y服从二项分布，则对于给定的自变量X，广义线性模型（此时为logistic回归）的拟合式为：
多项分布和多类别响应变量的logistic回归
上述是当响应变量只有两种取值（0或1）时的情况，最初的logistic回归基于响应变量的二项分布应用于对二分响应变量的建模。
而在很多实际情况中，响应变量并非简单的0-1关系，具有更多的类别属性。在此基础上，拓展了应用于多类别响应变量的logistic回归。例如：
多项logistic回归，响应变量包含两个以上的无序类别时，可使用基于多项分布的多项logistic回归。例如，寻找临床上具有预后价值的标志物RNA或蛋白分子时，响应变量包含多种组织类型，胃、肝、肠等。
序数logistic回归，响应变量是一组有序的类别时，可使用序数logistic回归。例如，寻找临床上具有预后价值的标志物RNA或蛋白分子时，响应变量是患者的肿瘤分期，早期、中期、晚期。
R语言执行logistic回归
以二分响应的logistic回归为例
临床研究中经常使用到logistic回归，因为经常涉及到有关生存、死亡或患病与健康等非常受关注的问题，这些临床事件就是类别型的变量。下文主要以一个简单示例，展示logistic回归在R语言中的计算过程，及对结果的解读。
R语言中执行logistic回归的方法很多，如基础包函数glm()提供了对二分响应变量的方法，mlogit包mlogit()函数可实现多项logistic回归，ms包lrm()函数可实现序数logistic回归等。下文主要以二分响应变量为例，展示最常见的二分响应的logistic回归。对于包含多类别响应变量时的logistic回归就不再举例了，需要时结合相关文献的方法描述学习对应的R包即可。
测试数据，R代码等的百度盘链接（提取码，ow0p）：
https://pan.baidu.com/s/1wY0QUmDGOB_sZS5LO3w9zA
若百度盘失效，也可在GitHub的备份中获取：
https://github.com/lyao222lll/sheng-xin-xiao-bai-yu
示例数据集
威斯康星州乳腺癌数据集，包含699个细胞组织样本，由空心细针在皮下肿瘤组织中抽得。其中458个（65.5%）为良性样本，241个（34.5%）为恶性样本。
#威斯康星州乳腺癌数据集
breast <- read.csv('http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data',na.strings = '?', header = FALSE)
names(breast) <- c('ID', 'clumpThickness', 'sizeUniformity', 'shapeUniformity', 'maginalAdhesion',
'singleEpithelialCellSize', 'bareNuclei', 'blandChromatin', 'normalNucleoli', 'mitosis', 'class')
str(breast)
#网盘中也备份了一份数据
#breast <- read.csv('breast.csv')
数据集中共有11个变量，包括：
ID：样本ID；
clumpThickness：肿瘤厚度；
sizeUniformity：细胞大小的均匀性；
shapeUniformity：细胞形状的均匀性；
maginalAdhesion：边际附着力；
singleEpithelialCellSize：单个上皮细胞大小；
bareNuclei：裸核，其中有16个样本中有缺失数据并用“?”表示；
blandChromatin：乏味染色质
normalNucleoli：正常核；
mitosis：有丝分裂；
class：肿瘤状态，2代表良性样本（共458个），4代表为恶性样本（共241个）；读取数据后对该列作个转换，可以修改为0-1类型变量，或者直接定义为类别型或因子型变量。
#预先将 class 列转换为类别或因子变量，或者重新定义为 0-1 类型
#为了方便下文部分作图过程，以重新定义的新列（0，良性；1，恶性）数据进行后续分析
breast[which(breast$class == 2),'class2'] <- 0
breast[which(breast$class == 4),'class2'] <- 1
除样本ID和肿瘤状态外，其余9个变量均可视为与判别恶性肿瘤相关的细胞特征，特征以得分为1（最接近良性）至10（最接近病变）之间的整数来表示。现在期望对细胞特征得分与肿瘤状态的关系进行建模，寻找最能代表肿瘤状态的细胞特征，logistic回归就是一种备选方案。
尽管用于分类型响应变量的建模方法有很多，某些机器学习算法如随机森林、支持向量机、判别分析等或许比logistic回归更优越，但本篇主讲logistic回归，故不讨论或比较其他方法的性能。
细胞特征得分与肿瘤状态的单变量logistic回归
首先展示最简单的单变量回归类型，帮助初步了解函数使用和结果解读。
单变量logistic回归及回归系数的生物学意义解读
假如只关注某个自变量，如clumpThickness（肿瘤厚度得分）对响应变量（肿瘤状态）的独立关系，就可以使用单变量logistic回归拟合自变量与响应变量的关系。
R语言中拟合广义线性模型的函数有很多，这里直接使用基础包函数glm()执行logistic回归。
#首先展示肿瘤厚度得分与肿瘤状态的一元 logistic 回归示例
#通过广义线性回归函数 glm() 执行 logistic 回归，详情 ?glm
#通过 family 参数指定了二项回归，即二分响应的 logistic 回归
fit_logistic1 <- glm(class2~clumpThickness, data = breast, family = binomial())
summary.glm(fit_logistic1) #展示拟合回归的简单统计
输出结果列出了回归系数、回归系数的标准误和参数为0的检验。
统计显著性是用检验统计量（z值）和相应的p值来表达的。首先根据检验的p值，可知clumpThickness（肿瘤厚度得分）是显著的，表明肿瘤厚度水平与肿瘤状态密切相关。
logistic回归中响应变量是Y=1的对数优势比（log），在单变量logistic回归中，回归系数表示了该自变量的变化可引起的响应变量对数优势比的变化。clumpThickness（肿瘤厚度得分）的回归系数是0.935，代表了肿瘤厚度的特征得分每升高1时，肿瘤状态恶性的对数优势比将相应增加0.935，即更高的肿瘤厚度得分对应了肿瘤恶性更大的可能性。而对于截距项，由于自变量无法为0（最低取值1），所以截距项尽管在统计上显著，但无实际的生物学意义。
由于对数优势比解释性差，不妨对回归系数进行指数化。如下所示，clumpThickness（肿瘤厚度得分）的回归系数的指数转化值为2.548，代表了当肿瘤厚度的特征得分每升高1时，肿瘤状态恶性的优势比将乘以2.548。
logistic_coef <- coef(fit_logistic1) #提取回归系数
exp(logistic_coef) #指数化回归系数
如果存在困惑，不妨结合公式理解logistic回归中自变量和响应变量的关系。根据上述单变量logistic回归获得的模型参数，代入公式可得clumpThickness（肿瘤厚度得分）与肿瘤状态的关系为：
尽管截距项无生物学意义，但数学上显著，因此需要在公式中保留，毕竟clumpThickness取值不会为0，故模型不会产生无意义的拟合值。p为二项分布试验成功的概率，在此处的logistic回归中，概率p的生物学意义代表了肿瘤为恶性的概率，根据公式可继续推算出其值：
可知，概率p（肿瘤为恶性的概率）的值将随clumpThickness（肿瘤厚度得分）值的增加而呈现非线性增长，随着肿瘤厚度的提升，代表了肿瘤是恶性可能性更高。
R中可使用predict()函数获得各肿瘤厚度得分下对应的肿瘤恶性概率。
#predict() 默认输出肿瘤为恶性的对数概率
#指定参数 type='response' 即可得到预测肿瘤为恶性的概率
prob <- predict(fit_logistic1, breast, type = 'response')
head(prob)
plot(breast$clumpThickness, prob, pch = 20, col = 'red',
xlab = 'clumpThickness（肿瘤厚度得分）',
ylab = 'probability（肿瘤为恶性的可能性）')
#一元 logistic 回归的散点图
library(arm)
plot(jitter(class2, factor = 0)~clumpThickness, data = breast,
type = 'n', xlab = 'clumpThickness（肿瘤厚度得分）',
ylab = 'probability（肿瘤为恶性的可能性）')
points(jitter(class2, factor = 0.3)~clumpThickness, data = breast,
col = '#00000022', pch = 20)
points(breast$clumpThickness, prob, pch = 20, col = 'red')
curve(invlogit(logistic_coef[1] + logistic_coef[2] * x), add = TRUE, col = 'red')
logistic回归将二项分布试验成功的概率p（肿瘤为恶性的概率）作为clumpThickness（肿瘤厚度得分）的函数来进行估计。
图中灰色的点呈现了原数据集中699个细胞组织样本的肿瘤厚度得分与肿瘤良性或恶性状态的对应关系（肿瘤良性时值为0，恶性恶性时值为1，作图时对点添加了少量上下抖动效果仅为减少点的重叠）。通过数据集真实的观测值可以看出，随着肿瘤厚度的提升，更多肿瘤表现为恶性。
红色曲线为logistic回归获得的概率拟合曲线，形象呈现了肿瘤为恶性的概率随肿瘤厚度得分增加的非线性响应。
模型预测
相比解释变量的生物学意义，对于类别型响应变量的建模方法而言，更多是关注模型的预测能力。上文通过单变量logistic回归既知clumpThickness（肿瘤厚度得分）更高时肿瘤为恶性的概率更大，不妨使用该回归对训练集样本作个预测，评估肿瘤厚度得分对于划分肿瘤状态的功效。
#predict() 默认输出肿瘤为恶性的对数概率
#指定参数 type='response' 即可得到预测肿瘤为恶性的概率
prob <- predict(fit_logistic1, breast, type = 'response')
#将概率大于 0.5 判定为恶性肿瘤，不大于 0.5 为良性
pred <- ifelse(prob > 0.5, 1, 0)
perf <- table(breast$class2, pred, dnn = c('Actual', 'Predicted'))
perf
#通过单变量 clumpThickness（肿瘤厚度得分）的预测准确性
(perf[1,1]+perf[2,2]) / sum(perf)
预测与实际情况对比的交叉表（即混淆矩阵，confusion matrix）显示，通过clumpThickness（肿瘤厚度得分）可以将86%的样本成功划分归类（肿瘤良性或恶性）。
由于没有提供测试集或验证集样本，未再进行更进一步准确度评估。
细胞特征得分与肿瘤状态的多变量logistic回归
上文首先以一个简单的单变量logistic回归，初步帮助理解了对R函数的计算过程及logistic回归参数意义的解读。现在继续延伸至多变量logistic回归的情形。
在上文结果中也不难看到，尽管较低的clumpThickness（肿瘤厚度得分）对应了肿瘤大概率是良性的可能性，但不可忽略少数肿瘤样本是恶性的事实。通过肿瘤厚度得分的单变量logistic回归可以将86%的样本正确分类，但仍存在14%的样本错误划分。
实际中，任一单变量其实都很难单独作为判别良性或恶性肿瘤的标准。如果以模型的预测性能作为主要的出发点，考虑更多自变量（细胞特征得分）的效应将有利于提升精度。
多变量logistic回归及回归系数的生物学意义解读
因此接下来综合考虑多种细胞特征并通过多变量logistic回归进行建模，目的是找到一组显著的变量组合区分良性或恶性肿瘤状态。与单变量的回归区别在于，回归中存在多个自变量时，各自变量之间互为协变量对待。
下文多变量logistic回归的示例中，暂且忽略多重共线性的影响。
#首先不妨使用全部细胞特征变量拟合与肿瘤状态的多元 logistic 回归
#通过广义线性回归函数 glm() 执行 logistic 回归，详情 ?glm
#通过 family 参数指定了二项回归，即二分响应的 logistic 回归，其余参数暂且使用默认值
fit_logistic_multi <- glm(class2~clumpThickness+sizeUniformity+shapeUniformity+maginalAdhesion+singleEpithelialCellSize+bareNuclei+blandChromatin+normalNucleoli+mitosis,
data = breast, family = binomial())
summary.glm(fit_logistic_multi) #展示拟合回归的简单统计
#上式可以简化为下式，其中“.”代表了数据集中除 class2 列外的其它所有列
#使用该式时，需保证除 class2 列外，其与所有列都作为自变量看待
#因此需要将数据集中的无效列，如代表样本编号的 ID 列和原始肿瘤状态的 class 列剔除
breast <- breast[ ,-which(colnames(breast) %in% c('ID', 'class'))]
fit_logistic_multi <- glm(class2~., data = breast, family = binomial())
summary.glm(fit_logistic_multi) #展示拟合回归的简单统计
输出结果列出了回归系数、标准误和参数为0的检验。
首先根据检验的p值，可知clumpThickness（肿瘤厚度）、maginalAdhesion（边际附着力）、bareNuclei（裸核）以及blandChromatin（乏味染色质）更为显著，表明它们的水平与肿瘤状态密切相关。同时normalNucleoli（正常核）接近显著性水平。相对于其它不显著的细胞特征，这些细胞特征或许可以作为判别良性或恶性肿瘤的候选标准。
对于截距项，由于所有自变量无法为0（所有细胞特征得分最小取值1），所以截距项没有生物学意义，尽管它数学上显著。
正回归系数代表了当该细胞特征得分增加时，肿瘤状态恶性的可能性增高；负回归系数则表示该细胞特征得分增加时，肿瘤状态良性的可能性增高。这里所有显著的细胞特征得分都对应了正回归系数，可知它们的发展与肿瘤恶性状态密不可分，代表了一组恶性肿瘤的特征。类似地，多变量logistic回归中响应变量是Y=1的对数优势比（log），由于各自变量之间互为协变量，因此回归系数的含义是当其它自变量不变时，一单位自变量的变化可引起的响应变量对数优势比的变化。例如，clumpThickness（肿瘤厚度）的回归系数是0.535，代表了在当其它自变量不变的情况下，肿瘤厚度的特征得分每升高1时，肿瘤状态恶性的对数优势比将相应增加0.535。
类似地，对回归系数指数转化后进行理解，正值的回归系数将转化为＞1的值，负值的回归系数将转化为<1的值。clumpThickness（肿瘤厚度）的回归系数的指数转化值为1.707，代表了在当其它自变量不变的情况下，肿瘤厚度的特征得分每升高1时，肿瘤状态恶性的优势比将乘以1.707。
logistic_coef <- coef(fit_logistic_multi) #提取回归系数
exp(logistic_coef) #指数化回归系数
大于3个自变量的多元回归无法在三维空间中绘制散点图，故不再展示肿瘤恶性概率与多变量间的响应关系，参考上文的单变量logistic回归理解即可。
模型简化和变量选择
由于sizeUniformity（细胞大小的均匀性）、shapeUniformity（细胞形状的均匀性）、singleEpithelialCellSize（单个上皮细胞大小）和mitosis（有丝分裂）的贡献不显著，不妨将它们从原回归模型中去除，只使用显著的clumpThickness（肿瘤厚度）、maginalAdhesion（边际附着力）、bareNuclei（裸核）、blandChromatin（乏味染色质）以及接近显著性的normalNucleoli（正常核）重新拟合logistic回归以简化模型。
此时，前后的两个logistic回归（使用所有9种细胞特征的logistic回归，和只使用5种细胞特征的logistic回归）可视为嵌套模型。所谓嵌套模型，即其中之一的一些项完全包含在另一模型中，这种关系有利于进行评估，简化后模型的综合性能是否比先前的要好。
#只使用观察到的显著和接近显著的 5 种细胞特征拟合与肿瘤状态的 logistic 回归
#这里除了更改协变量外，其余参数使用默认值，和先前一致
fit_logistic_multi2 <- glm(class2~clumpThickness+maginalAdhesion+bareNuclei+blandChromatin+normalNucleoli,
data = breast, family = binomial())
summary.glm(fit_logistic_multi2) #展示拟合回归的简单统计
#anova() 检验两个嵌套模型的拟合优度
#卡方检验可用于比较广义线性模型
anova(fit_logistic_multi, fit_logistic_multi2, test = 'Chisq')
#AIC（Akaike Information Criterion，赤池信息准则）评估两个回归复杂性与拟合优度的关系
#AIC 值较小的回归优先选择，表明较少的预测变量已经获得了足够的拟合度
AIC(fit_logistic_multi, fit_logistic_multi2)
检验结果和AIC值均显示，前后的两个logistic回归在性能上没有很明显的区别。也就是说，sizeUniformity（细胞大小的均匀性）、shapeUniformity（细胞形状的均匀性）、singleEpithelialCellSize（单个上皮细胞大小）和mitosis（有丝分裂）可有可无，不怎么影响模型精度。相比之下更倾向选择后者，即去除了这4个不显著自变量后的多变量logistic回归，因为它更简洁。（两个嵌套模型的比较，一般而言：2优于1时，选2；1优于2时，选1；1和2无区别时，选更简洁的那个）
对于新的logistic回归，由于排除了先前部分协变量的影响，所以自变量的回归系数发生了微小的变动。若要尝试解释各细胞特征的效应时，需要重新结合回归系数解释。
此外，除了像上述这种手动比较选择变量外，也可以通过逐步回归进行自动选择，如下是一个简单的后向选择例子。关于逐步回归的概念及对结果的解读，可见前文“多元回归中常见的变量选择方法”，不再多说。
#后向逐步回归自动选择变量
fit_logistic_reduced <- step(fit_logistic_multi, direction = 'backward')
fit_logistic_reduced
模型预测
现在不妨使用上述模型对对训练集样本作个预测，评估选择的细胞特征得分组合是否能够准确区分良性或恶性肿瘤状态，以实现对恶性肿瘤的准确预测。
#predict() 默认输出肿瘤为恶性的对数概率
#指定参数 type='response' 即可得到预测肿瘤为恶性的概率
prob_multi <- predict(fit_logistic_multi, breast, type = 'response')
prob_multi2 <- predict(fit_logistic_multi2, breast, type = 'response')
#将概率大于 0.5 判定为恶性肿瘤，不大于 0.5 为良性
pred_mult <- ifelse(prob_multi > 0.5, 1, 0)
perf_mult <- table(breast$class2, pred_mult, dnn = c('Actual', 'Predicted'))
perf_mult #全部的 9 个细胞特征得分对肿瘤恶性状态的预测性能
(perf_mult[1,1]+perf_mult[2,2]) / sum(perf_mult)
pred_mult2 <- ifelse(prob_multi2 > 0.5, 1, 0)
perf_mult2 <- table(breast$class2, pred_mult2, dnn = c('Actual', 'Predicted'))
perf_mult2 #去除 4 个不显著细胞特征得分后对肿瘤恶性状态的预测性能
(perf_mult2[1,1]+perf_mult2[2,2]) / sum(perf_mult2)
与上文仅使用clumpThickness（肿瘤厚度得分）的单变量logistic回归（预测准确度86%）相比，考虑多种细胞特征得分的多变量logistic回归的准确度明显提升，96%以上的样本被准确划分了肿瘤良性或恶性状态。
使用全部9个细胞特征得分与使用5个显著或接近显著的细胞特征得分的多变量logistic回归在预测精度上没有观察到区别（后者仅1个预测值与前者不同），选择使用后者更合适，至于那4个不显著的细胞特征得分可视为噪声来看待了。
偏大离差问题及解决方案
在前文“广义线性模型之泊松回归”中，提到当响应变量观测值的方差比回归模型预测值的方差大时，回归可能发生偏大离差（overdispersion，或称过度离势）。类似地，logistic回归中，当响应变量观测值的方差大于期望的二项分布的方差时也可能会带来偏大离差问题。
在二分响应变量的logistic回归中，假定响应变量服从二项分布。如果成功的概率p已知，响应变量（m次试验中成功的次数）的方差就是已知的（mp(1-p)）。如果二项试验不是独立的，或者各二项响应的p是不同的，或者重要的自变量没有包含在关于p的模型中，那么响应变量的方差通常会比二项模型条件下的期望方差要大，产生偏大离差。
偏大离差的存在将导致奇异的标准误检验和不精确的显著性检验，对结果的可解释性造成负面影响，因此在建模时是不可忽视的。
以上述获得的第二个多变量logistic回归（只包含5个显著或接近显著细胞特征得分的多变量logistic回归）为例，在R中评估偏大离差的过程如下。
#检测偏大离差的一种方法是比较二项分布模型的残差偏差与残差自由度
#比值偏离 1 程度较大提示存在偏大离差
deviance(fit_logistic_multi2)/df.residual(fit_logistic_multi2)
#检验，零假设不存在偏大离差
#p 值显著则拒绝零假提示存在偏大离差
fit <- glm(class2~clumpThickness+maginalAdhesion+bareNuclei+blandChromatin+normalNucleoli,
data = breast, family = binomial())
fit.od <- glm(class2~clumpThickness+maginalAdhesion+bareNuclei+blandChromatin+normalNucleoli,
data = breast, family = quasibinomial())
pchisq(summary(fit.od)$dispersion * fit$df.residual, fit$df.residual, lower = FALSE)
提示上述讨论的回归模型实际上是存在偏大离差问题的，可能会产生不可信的生物学意义解读。
如果存在偏大离差，可考虑重新调整对拟合模型的统计推断，通常是将估计出的模型系数的标准误乘上偏大离差的平方根。此举不会改变原模型的回归系数，但扩大了回归系数标准误和置信区间，增加显著性检验的严格度，以缓解由偏大离差导致的检验不精确的问题。
存在偏大离差的二分响应数据可以用考虑了偏大离差问题的logistic回归来拟合。在R 函数glm()中，可以通过指定参数family='quasibinomial'代替先前的family='binomial'。
#考虑了偏大离差问题的logistic回归
#通过广义线性回归函数 glm() 执行时， family 参数指定 quasibinomial，详情 ?glm
fit_quasibinomial <- glm(class2~clumpThickness+maginalAdhesion+bareNuclei+blandChromatin+normalNucleoli,
data = breast, family = quasibinomial())
summary.glm(fit_quasibinomial) #展示拟合回归的简单统计
#先前的多变量 logistic 回归
summary.glm(fit_logistic_multi2)
输出结果列出了回归系数、标准误和参数为0的检验，前后回归的唯一区别在回归系数标准误的估计值上。
能够看到，各自变量的回归系数和先前相比没有改变。但回归系数的标准误变大了，此举扩大了标准误和置信区间，增加了显著性检验的严格度，也很容易注意到这里的p值也普遍比先前回归中的大。
尽管如此，在考虑了偏大离差问题的logistic回归中，5个细胞特征得分的自变量仍然全部是显著的。并且由于回归系数并未发生改变，也就是说这5个自变量的生物学意义解释和先前相比不会改变，预测精度也和先前一样。因此，尽管上文讨论了大篇幅的回归模型存在偏大离差问题，但它并没有产生不可信的统计检验，大可以放心地用于生物学解释或预测中。
但若和先前的统计检验区别较大，以考虑了偏大离差问题的logistic回归结果为准，并重新解释变量的效应，或者构建可靠的预测模型。
参考资料
关于多元回归中需考虑的问题：
一般线性模型（LM）：
常见的非线性模型：
回归树：
生存分析：
基于相似或相异矩阵的回归或降维思想的回归：
基于距离的冗余分析（db-RDA），或称典范主坐标分析（CAP）
https://mp.weixin.qq.com/s/lVt_1UZLs0VMlo0keGfHLg