[R 语言]以白血病结局预测为例说明逻辑回归建模（下）
一例白血病临床数据的逻辑回归建模和优化（下篇）
前言
临床场景中，很多诊断和预后的结局是二分类的（是或非），比如阳性结果的“发生”和“未发生”。医生需要根据病人身上得到的各类有价值的线索来综合预判最有可能的结果，这就是逻辑回归最常应用的场景。本课时用一组临床的白血病数据集为范例，目标是建立逻辑回归预测模型，预测疾病复发 (remission)。
复发与否就是一个二分类的结局，我们要从数据中找规律，即各个因素（“自变量”）与结局（“因变量”）之间是否存在联系，和联系的强弱。比如，某一个自变量在大部分未复发的病例中都等于1，在大部分复发的病例中都等于5，那这就呈现了一定的关联。这就是逻辑回归通过处理样本数据而帮我们归纳的信息。如果存在联系，逻辑回归就认为这个自变量应该成为预测因子，并根据整个数据集，赋予这个自变量一个系数来表示联系的强弱，那么预测因子的取值乘以相应系数就是对结局产生的作用。对于数据中提供的每个病例，其所有预测因子的作用总和会产生一个最终的预测值（即结局发生的概率，在 0 到 1 之间），逻辑回归尝试最大程度地让预测值符合实际情况，而实际情况也就是现有数据告诉我们的复发和未复发，我们需要将预测值与实际值对比，看模型表现。如果把这个模型拿到当前数据之外的其他数据集，也可以查看其预测表现，这叫做验证。
可以用我们的视频 来理解上述内容，会更加浅显易懂。
回到这节课的任务。当拿到一个现实世界来源的数据，第一步要先做一些清理，比如查找缺失和重复，发现不合理的数值，规范格式等基本操作，才能使数据提供分析价值。然后，根据数据内容做一些浅显的分析（如单因素和多因素分析），进一步查看数据合理性与数据特征，这也叫做探索性数据分析。在这之后，分析人员心中才会对接下来的模型建立有了初步的想法，于是就开始了建模的过程。建好的逻辑回归模型会对，会呈现一系列参数来说明它的表现。接下来讲解这些内容。
共分三篇：
学习目标：
a. 根据原始数据，以 “REMISS” （复发）为结局，用其他变量预测结局，优化模型
b. 批判地理解并尝试逐步回归选择法、LASSO 回归
c. 评价模型的主流方法，熟悉模型诊断
d. 理解并比较几种拟合优度检验
本课时相关资源：
高校资源
经典书籍
人工智能 R 语言工具
其他见文中附超链接
以上资源并非生硬搬运，而是以一种注重读者理解的方式融合、发散性地呈现。
下面开始正文：
1，2，3（数据导入、数据清理；探索性数据分析）见上篇 ；
4，5，6（逻辑回归建模，变量筛选）见中篇 ；直接阅读本篇不会影响理解
7. 模型诊断之拟合优度
得到模型后，如何从统计学的角度查看它的“好坏”，叫做模型诊断 (model diagnostics)，无论何种诊断，都是算出一个统计量进行假设检验来给出结论，这也是统计学的精髓。比如，衡量模型与数据的拟合程度就是模型诊断的首要方面，下面会讲述拟合优度 (goodness-of-fit, GoF) 的衡量指标，有似然比检验，伪 R 方，残差平方和，Hosmer-Lemeshow 检验。
7.1 似然比检验 (likelihood ratio test, LRT)
似然比 LR 其实是 AIC（和 BIC）的运算基础。与 AIC 一样，也是似然比越小、模型越好。
LR 的公式体现了它的作用和优势。它衡量的每个数据点距离拟合曲线的距离总和，也衡量这些点到 null model（零模型，通常是一条没有斜率的直线）的距离总和。然后将两个 “总和” 做 log 变换再相减（变换是为了让结果符合卡方分布）。其中，每个数据点连起来就是饱和模型 (saturated model)，因为 "saturated" 代表模型有着与数据点相同数量的参数，这种模型其实就是数据本身而已。
另外，LR 也可以将两个非零模型做对比，但其中一个模型的变量应该是另一个模型变量的子集 (nested model)，才能实现有效的比较。为加深理解，举例两个应用场景：我们想要知道年龄是否应该纳入模型，就分别计算有年龄和没有年龄的 LR 然后做比较；如果想知道几个人口相关因素 (demographics) 是否重要（比如，用来决定不同地区的数据建模时是否纳入这些因素），就可以分别计算有和没有这些因素时的 LR。
似然比公式：
**LR = -2 * [logLik(proposed_deviance) - logLik(null_deviance)] = -2 * {logLik[(proposed_deviance)/(null_deviance)]};
"proposed" 是指要衡量的模型，"null" 是零模型，两者的偏差 (deviance) 相除，结果取 logLik() 计算。
无论是非零与零模型对比，还是两个非零模型对比，原理是相同的，LRT 的计算结果（称为统计量 G2）符合卡方分布，所以只要通过查看其相应 p 值，就可以给出结论，**提示更复杂的模型是否有更好的拟合优度。在 R 语言中，只要将两个相比较的模型传给 anova() 函数即可。
只写一个 "simple_model"，则默认是与 null_model 相比（否则需要写两个 model 名称，不论顺序）：
anova(simple_model, test=“Chisq”)
返回结果：
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 26 34.372
## LI 1 8.2988 25 26.073 0.003967 **
以上就是用 LRT 对模型 REMISS ~ LI 做检验得到的结果，说明 LI 所代表的单个变量的模型比零模型拟合度更好。要记住 LR 对比的是更复杂的模型和更简单的模型，p 值小于设定值则说明更复杂的模型拟合度更好（这时 simple_model 虽然取名 simple，但还是比 null mode 更复杂）。
注意，之前（“中篇” ）用 LI 建立的单变量逻辑回归模型，当中也有对 LI 的检验，这个检验叫 Wald Test，检验的是相应变量对结果是否有显著影响，它的数值与前面的 LRT 不同，因为这两种检验内容本质上不同。
似然比检验用于检验 β 的任何子集是否等于 0 的零假设。Wald 检验是逻辑回归中单个回归系数的显著性检验。
7.2 伪 R 方 (pseudo R-squared)
伪 R 方公式里的要素与 LR 相同，但计算方式不同：
R-squared = 1 - logLik(reduced_deviance)/logLik(null_deviance)
LR 与（伪）R 方还有以下不同：
本案例的 “R-squared” = 0.2414424 即 24.14%，伪 R 方数值越大，模型拟合度越好。LR 则是越小拟合度越好；
我们可以使用 LR 来估计模型所解释的变异性（R 方）。在理想情况下，R 方可以直接比较不同的预测因子 (predictors)，且无论预测因子是分类的还是连续的；
从公式来看R 方也与样本量无关，LR 的绝对值取决于 n（数据量），类似于线性回归中的平方和；
此外，我们还需要定义一个相对于零模型的 R 方数值才能计算下去，LR 则不依赖于零模型。
7.3 残差平方和与 Pearson 卡方检验
Pearson 卡方拟合优度检验 和残差平方和拟合优度检验 （类似于多元线性回归 F 检验）需要数据集包含重复观测值。这两种检验的统计量都近似于具有 c - k - 1 自由度的卡方分布，其中 c 是预测变量的不同组合数。如果拒绝接受H0检验，则在统计上存在显著的不拟合。否则，只能说还没有证据表明缺乏拟合度。
这两者的算式背后的原理可以参考上述超链接，我们只需要知道他们的 R 代码很简单：
# simple_model 表示只有自变量 “LI” 的模型
sum(residuals(simple_model, type=“deviance”)^2) # 结果为 26.07296
sum(residuals(simple_model, type=“pearson”)^2) # 结果为 23.93298
由于 simple_model 只有 LI 一个自变量，所以 LI 的似然比检验和该模型的残差平方和均为 26.07296， 用代码 simple_model$deviance 也可以得到这一结果。
然后用卡方统计量来检验：
# 公式 pchisq(q, df, lower.tail = TRUE) lower.tail 默认为 TRUE
pchisq(26.07269, 25, lower.tail = FALSE)
# 或者
1-pchisq(26.07269, 25)
如果 p 值较大，接受假设，没有证据证明缺乏拟合度。
用残差图 (Residual Plots) 评估上述两种残差：
使用 Pearson 残差和偏差残差绘制残差图，看看这些图中是否呈现不好的趋势。理想的残差图称为空残差图（null residual plots），显示的是随机的点散布在特征线周围，形成一个宽度近似恒定的带状区域，否则，表明模型存在问题。
残差图用下列方式作图：
# 先用 length() 函数先确定 residual 的个数为 27
length(residuals(simple_model, type=“pearson”))
length(residuals(simple_model, type=“deviance”))
作图的要求是横坐标为n（表示第n个residual），1:27 表示 “从 1 到 27” 个横坐标点，type=“b”是作图格式：
plot(1:27, residuals(simple_model, type=“pearson”), type=“b”, main = “Pearson”)
plot(1:27, residuals(simple_model, type=“deviance”), type=“b”,main = “Deviance”)
plot(residuals(simple_model,type=“pearson”),residuals(model.2, type=“deviance”))
下面两张图均显示数据点在特征线上下随机分布：
7.4 Hosmer-Lemeshow test (HL 检验)
另一个常用的模型拟合度检验是 HL 检验。它的原理是，预测频率 (predicted frequency) 和观测频率 (observed frequency) 匹配越密切，拟合度越高。HL 拟合度统计量是这两组数值的列联表 (contingency table) 计算出的 Pearson 卡方数值。这种列联表的关联检验与双向表 (two-way table，即四格表加上总和统计) 的关联检验类似，较大的 p 值代表拟合度良好。
当模型中有连续变量 (continuous variable) 时，会导致一个很长的列联表，这往往容易错误地产生显著的结果（较小的 p 值），从而认为模型拟合度差。因此，通常的做法是将预测变量分为 g 个组。组数的确定没有理论依据可循，一般范围在5-10。例如当 g=10，形成 2×10 的列联表，第一组数据就是预测值最小的那 10% 的样本。HL 检验得到的检验统计量近似于具有 g-2 自由度的卡方分布。
需要注意，对于样本量小于 500 的数据, HL 检验的效力较低。且自定义不同的组数会计算出不同的 p 值大小，从而可能得出不同的结论。
下面来运行 HL 检验用于本例的 R 代码。从加载 R 包 ResourceSelection 开始。HL 函数默认 g=10，即 hoslem.test(x, y, g = 10)，本例设置 g=9：
library(ResourceSelection)
hoslem.test(simple_model$y, fitted(simple_model), g=9)
HL 检验还可以用来校准 (calibration)，因为可以看到每组样本范围内预测频率和观测频率的一致程度，g 个组的一致程度连在一起可以作为一种校准曲线 (calibration curve)。但是也有反对这种做法的。首先，样本越大，伪显著性越高（HL 检验越容易认为校准度差）；也有认为模型校准应该直接避免用 HL 检验 。
由于 Pearson 卡方和残差平方和这两种拟合优度评价都要求存在重复观测值 ("multiple observations with the same values for all the predictors" ，有了重复观测值就会有相同的预测值，但实际结果值不同) 相比之下，HL 检验适用于重复观测值没有或很少的数据。
所以本例应采取 HL 检测的结果，而非另两种 GoF 的检验结果：
## Goodness-of-Fit Tests
## Test DF Chi-Square P-Value
## Deviance 25 26.07 0.404
## Pearson 25 23.93 0.523
## Hosmer-Lemeshow 7 6.87 0.442
与大多数拟合优度检验一样，这里的 H0 假设本质是将模型与 saturated model 相比，SSR (sum of squared residuals) 数值越大，p 值越小。小的 p 值（通常低于 5%）意味着模型拟合得不好。但是，大的 p 值并不一定意味着模型拟合得很好，只是没有足够的证据表明它拟合得不好。此外，很多其他情况也都会导致 p 值偏大，包括样本量小检验效能差。
注意：GoF 衡量的是fit的程度，所以只用这一指标衡量的话，可能导致模型过拟合 (overfitting)。
至此全部讲解都结束了。恭喜你已经学了非常多的知识，现在你已经可以尝试做一个逻辑回归建模！！
我们详细的项目解析，除了融合最实用的 R 语言代码和最前沿的统计学解释，更强调相应的医学背景！！我们希望让不同的学科背景的读者都能够“知其然，知其所以然”，真正明白每一步“为什么那么做”、”下次我自己会如何实现“。
同时为了大家能够更加习惯学术英文的使用，确保关键词汇遵循学术使用习惯，并且能直接用在学术写作中，我们在一些关键部分使用了中英双语注解。
学习建议
关于 R 语言，大家不必有焦虑，更不要啃下一本教程 - Preparation isn't readiness!
无论学习还是使用 R，应该在遇到实际问题时用高效的人工智能辅助，可以广泛利用开发好的资源，比如找代码模板，甚至 “开口要” —— 用大语言模型获取基本代码框架。
又比如公益在线工具 RTutor：将数据上传至 RTutor，用自然语言（中英文皆可）告诉它想要做什么数据处理，实现什么样的统计、作图，它会返回代码和结果，并且可以根据结果追问。你可以拿这些结果来运行和调整，再针对性解决出现的问题 – 把使用和学习 R 语言放在解决实际问题的过程中才是最高效的。而类似 R 语言教程可以作为词典辅助使用。
建议加入我们的互助社群，可以帮助你进行临床医学研究设计，论文预审，R 语言安装和运行答疑；还可以获取课程的 RMD/MD 文件，以便直接复制到自己电脑上操作学习。
我们会有更多“数据-统计-临床”三位一体的讲解，结合最前沿的权威知识，甚至英文和 R 语言都呈现到位，敬请关注。Keep tuned!!
hellomedstat 所有平台同名。
微信：Positivism_y
QQ 群：760447631
邮箱：hellomedstat@gmail.com
