[R 语言]以白血病结局预测为例说明逻辑回归建模（上）
一例白血病临床数据的逻辑回归建模和优化（上篇）
前言
什么叫逻辑回归（logistic regression）？临床场景中，很多诊断和预后的结局是二分类的（是或非），比如阳性结果的“发生”和“未发生”。医生需要根据病人身上得到的各类有价值的线索来综合预判最有可能的结果，这就是逻辑回归最常应用的场景。本课时用一组临床的白血病数据集为范例，目标是建立逻辑回归预测模型，预测疾病复发 (remission)。
复发与否就是一个二分类的结局，我们要从数据中找规律，即各个因素（“自变量”）与结局（“因变量”）之间是否存在联系，和联系的强弱。比如，某一个自变量在大部分未复发的病例中都等于1，在大部分复发的病例中都等于5，那这就呈现了一定的关联。这就是逻辑回归通过处理样本数据而帮我们归纳的信息。如果存在联系，逻辑回归就认为这个自变量应该成为预测因子，并根据整个数据集，赋予这个自变量一个系数来表示联系的强弱，那么预测因子的取值乘以相应系数就是对结局产生的作用。对于数据中提供的每个病例，其所有预测因子的作用总和会产生一个最终的预测值（即结局发生的概率，在 0 到 1 之间），逻辑回归尝试最大程度地让预测值符合实际情况，而实际情况也就是现有数据告诉我们的复发和未复发，我们需要将预测值与实际值对比，看模型表现。如果把这个模型拿到当前数据之外的其他数据集，也可以查看其预测表现，这叫做验证。
可以用我们的视频 来理解上述内容，会更加浅显易懂。
回到这节课的任务。当拿到一个现实世界来源的数据，第一步要先做一些清理，比如查找缺失和重复，发现不合理的数值，规范格式等基本操作，才能使数据提供分析价值。然后，根据数据内容做一些浅显的分析（如单因素和多因素分析），进一步查看数据合理性与数据特征，这也叫做探索性数据分析。在这之后，分析人员心中才会对接下来的模型建立有了初步的想法，于是就开始了建模的过程。建好的逻辑回归模型会对，会呈现一系列参数来说明它的表现。接下来讲解这些内容。
共分三篇：
-
上篇：数据导入、数据清理；探索性数据分析包括单因素、多因素分析（本篇）
-
还有贯穿这些内容的主流 R 语言代码；对每一步方法和结果的辩证判读。
学习目标：
a. 根据原始数据，以 “REMISS” （复发）为结局，用其他变量预测结局，优化模型
b. 批判地理解并尝试逐步回归选择法、LASSO 回归
c. 评价模型的主流方法，熟悉模型诊断 (model diagnostics)
d. 理解并比较几种拟合优度检验 (Goodness-of-Fit tests)
本课时相关资源：
高校资源
经典书籍
人工智能 R 语言工具
其他见文中附超链接
以上资源并非生硬搬运，而是以一种注重读者理解的方式融合、发散性地呈现。
（阅读过程中可跳过斜体英文，下方会有相应中文，不影响理解）
下面开始正文：
1. 数据获取
本地（自己电脑上的）数据可以用以下方式直接打开：
RStudio -> File –> Open File
或是用读取指令：
- 对 csv 文件 -> read.csv(path.csv) 读取一个 csv 文件的路径；
- 对Excel 文件 -> 先加载 readxl 包，再用 read_excel (path.xlsx)。
网络来源的数据可以先下载到本地，再用上述方法。也可以直接读取：把网络链接加引号放在 read.table() 里。接着，如果保持使用原来表格的标题，则 header=TRUE。
一般只要写完路径和 header，代码就结束了，但这份数据需要纠正，另 fileEncoding = UCS-2 才能读取：
#
leukemia <- read.table
('https://online.stat.psu.edu/stat462/sites/onlinecourses.science.psu.edu.stat462/files/data/leukemia_remission/index.txt',
header=T,
fileEncoding = "UCS-2")
用 attach() 可以让 R 在评估变量时会搜索数据库，因此只需提供数据库中的对象名称即可访问这些对象:
attach(leukemia)
接下来，对手中的数据要按顺序做三件事：查看数据内容（有无异常、类型是否正确），数据分布（是否在正常范围内合理分布），初步探索内在联系（如两两相关性）...
2.原始数据处理
这一步主要是与纯数据打交道，查看、清理数据，还有关键的一步是理解数据内容。
可以先从查看数据结构开始，str 表示structure：
str(leukemia)
## 'data.frame': 27 obs. of 7 variables:
## $ REMISS: int 1 1 0 0 1 0 1 0 0 0 ...
## $ CELL : num 0.8 0.9 0.8 1 0.9 1 0.95 0.95 1 0.95 ...
## $ SMEAR : num 0.83 0.36 0.88 0.87 0.75 0.65 0.97 0.87 0.45 0.36 ...
## $ INFIL : num 0.66 0.32 0.7 0.87 0.68 0.65 0.92 0.83 0.45 0.34 ...
## $ LI : num 1.9 1.4 0.8 0.7 1.3 0.6 1 1.9 0.8 0.5 ...
## $ BLAST : num 1.1 0.74 0.18 1.05 0.52 0.52 1.23 1.35 0.32 0 ...
## $ TEMP : num 1 0.99 0.98 0.99 0.98 0.98 0.99 1.02 1 1.04 ...
显示都是 integer (int) 和 number (num) 格式，说明原始数据不包含数字以外的内容，否则要查看是否有 “？” 等代替缺失值的符号，同时也能看到只有 27 个观测值（observations）。另外，也没有需要转换成分类变量的自变量。
再用 View() 查看原始数据的全貌（还可以用 glimpse 函数查看一部分）：
View(leukemia)
如果随机抽取部分数据查看，用sample_n(): 这是dplyr包中的函数，用于从数据集进行随机抽样。它可以按照指定的数量来随机抽取观测。如果想再不重复地抽取一次，只要再运行一次代码：
library(dplyr)
sample_n(leukemia, 4)
## REMISS CELL SMEAR INFIL LI BLAST TEMP
## 1 0 0.50 0.44 0.22 0.6 0.11 0.99
## 2 0 0.70 0.76 0.53 1.2 0.15 0.98
## 3 0 1.00 0.73 0.73 0.7 0.40 0.99
## 4 0 0.95 0.87 0.83 1.9 1.35 1.02
注意 sample_n 与 sample 函数的区别：sample 函数是随机选取几个列。
除了单纯的查看，对于较大的数据集常常要进行清理 (data cleaning)，处理**缺失数据 (missing data) **等。这在心血管病那一例课时中有详细操作。本例只有27个样本，逐个查看之后不需要上述处理。
接下来我们要对数据做一些理解。原资料对变量名的解释：
（翻译）自变量包括：骨髓凝块切片的细胞度（CELL）、涂片差异囊泡百分比（SMEAR）、骨髓白血病细胞绝对浸润百分比（INFIL）、骨髓白血病细胞百分比标记指数（LI）、外周血中囊泡的绝对数量（BLAST）和开始治疗前的最高温度（TEMP）。
所以，原始数据的变量与疾病复发诊断存在不同意义上的关联。这些名称已经具备较好的自释性（self-explanatory），足以代表各自含义，不需要修改。
3. 探索性数据分析 EDA
经过上述过程有了真正准备好的数据，才可以进行一些简单的统计学处理。处理方法通常是单因素和多因素分析，以查看数据特征。EDA (exploratory data analysis) 是用来分析和研究数据集并总结其主要特征的方法，通常还采用数据可视化方法。
3.1 查看数据分布
初步带着问题查看数据：有多少病例（观测值）、想要考虑多少变量、结果的性质（连续、二元结果，生存数据等）、以及对目前临床实践的了解程度。在进行任何分析之前，先有一个初步假设或预期，避免盲目分析和无意义的探究。
以下方法可以看出大致频数分布，离群值 (Outliers) ，和数据重复情况：
1. histogram
2. qq-plot, 比 histogram 更好
3. 茎叶图 stem(leukemia$variable) 默认scale = 1；
例如用代码：stem(leukemia$REMISS)查看REMISS分布
另外在进一步数据处理中，很多常见的统计学方法都要求数据满足正态性，如常见的 t 检验、单因素方差分析等。有的数据进行 log 变换后会呈现正态性。本例使用的逻辑回归无需事先假设数据分布，所以不进行参数的分布类型检验。
3.2 数据内在联系和单因素分析
此时，我习惯先查看各变量之间的两两关联，可以暴露一些问题并提供一些思路。这一步在向数据分析过渡。举例方法一 ：用 R 自带的 pairs() 就可以实现。方法二：用 GGally 包中的 ggpairs 实现。
在数据处理初始阶段，我们对数据分析和呈现要求不高，所以用 pairs() 就够了：
pairs(leukemia)
这张图将所有变量进行两两配对，就是将每个变量上的取值画在坐标系中。
对比 ggpairs：
library(GGally)
ggpairs(leukemia)
ggpairs 除了进行两两配对，还计算了相关性参数
我们从上述两个散点图都看到 SMEAR 和 INFIL 之间有显著的线性关联。两个自变量之间存在线性相关性叫做共线性 (collinearity/multi-collinearity)。共线性在我们的视频课 "Feature Selection for Logistic Regression (逻辑回归的变量筛选)" 中有讲解。简单地理解共线性，就是说两个变量的作用可能是相同的，我们或许可以去掉一个。
3.3 单因素 多因素分析 优势比 模型参数
对于多个变量 (multiple variables) 的优势比 (Odds Ratio, OR)，一般先进行单因素回归分析，也就是每次只查看一个变量与结果的关联，通过简单的列联表就可以计算有该因素和没有该因素分别对应多少事件发生，从而生成单因素分析下的 OR 等。这样的单因素回归分析在 Hosmer 和 Lemeshow 的书 Applied Logistic Regression 中也被描述为 “purposeful variable selection”，与后续的建模相关。
但是要警惕完全依赖单因素回归分析的结果。比如，一个在单因素分析中不呈现关联的变量可能在多因素分析中反而具有重要性。
本例的单因素分析其实前面的散点图就有呈现，其进一步的 OR 计算水到渠成，具体的计算在我们的心血管病的相应课时中有介绍。
接下来直接讲多因素分析，也就是将所有变量纳入到同一个计算系统，这样他们之间不再孤立，计算出的 OR 叫做调整后优势比 (adjusted OR, aOR)，比如逻辑回归就是这样一个“计算系统”。所以，逻辑回归模型不仅是建模工具，也是计算 aOR 的工具。
用 glm() 拟合逻辑回归曲线，将这个包含所有变量的模型取名叫 full_model，binomial 表示结局为二分类（复发 1 和未复发 0）:
full_model <- glm(REMISS ~ CELL + SMEAR + INFIL + LI + BLAST + TEMP, data = leukemia, family = binomial)
summary(full_model)
##
## Call:
## glm(formula = REMISS ~ CELL + SMEAR + INFIL + LI + BLAST + TEMP,
## family = binomial, data = leukemia)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 64.25808 74.96480 0.857
## CELL 30.83006 52.13520 0.591
## SMEAR 24.68632 61.52601 0.401
## INFIL -24.97447 65.28088 -0.383
## LI 4.36045 2.65798 1.641
## BLAST -0.01153 2.26634 -0.005
## TEMP -100.17340 77.75289 -1.288
## Pr(>|z|)
## (Intercept) 0.391
## CELL 0.554
## SMEAR 0.688
## INFIL 0.702
## LI 0.101
## BLAST 0.996
## TEMP 0.198
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 34.372 on 26 degrees of freedom
## Residual deviance: 21.594 on 20 degrees of freedom
## AIC: 35.594
##
## Number of Fisher Scoring iterations: 8
summary(full_model) 给出了模型的各变量回归系数和相应的检验 (Wald test) ，这一套结果只适用于当前模型，如果增减变量或改变数据集，结果就会发生变化。
目前来看变量 LI 有最小的 p 值。从统计学角度，p 值大小与该变量的预测能力相关。如果用 p 值筛选变量，应比 0.5 更宽松（因为并不是进行研究问题的假设检验），一般选 p = 0.10 或 0.20 (cut-off value)。
注意，要在上述变量 p 值计算出来之前设定好这个选择界限，而不是算出来之后再设定，这在我们的视频课中有提及。虽然用 p 值筛选是一个“常规牌”，但是要知道，一些主流学术机构已经不用 p 值筛选变量，而是鼓励用更有现实意义的标准，比如 DAG (Directed acyclic graph)。
对 AIC (Akaike information criterion) 的理解：单个 AIC 数值没有意义。但是，如果拟合了多个回归模型，AIC 可以用来比较模型的表现。AIC 值最小的模型拟合效果最好。因此，从一组可能的模型中选出的最佳模型就是 AIC 值最小的模型。这是一种 AIC 法（AIC-based approach），因为它没有穷尽所有候选模型，这类方法计算速度相当快，而且可以处理复杂的模型公式，包括由两个自变量叠加成的新变量交互项 (interaction terms)。
3.4 回归系数的置信区间
要进一步计算每个回归系数的置信区间 CI (confidence intervals)：
置信区间有助于确定系数的统计意义。如果置信区间不包括零值，则表明参数在特定的显著性水平（通常为 95%）上具有统计意义。
基于渐进正态性（asymptotic normality）和最小二乘估计（least-squares estimation，LSE）两种方法计算参数的 95% CI。运行后可以看到，只有 LI 的回归系数 CI 区间不包括 0，说明只有 LI 与结果的关系方向是恒定的（有95%可能性）。还可以看到一些很大的置信区间，这可能与“共线性” (multi-collinearity) 相关，本案例的中篇会讲。
置信区间代码和结果见下：
confint.default(model.1) #基于渐进正态性
## 2.5 % 97.5 %
## (Intercept) -82.6702241 211.186392
## CELL -71.3530610 133.013183
## SMEAR -95.9024389 145.275071
## INFIL -152.9226458 102.973698
## LI -0.8490839 9.569992
## BLAST -4.4534852 4.430423
## TEMP -252.5662623 52.219462
confint(model.1) #基于最小二乘估计
## 2.5 % 97.5 %
## (Intercept) -68.6041096 241.051733
## CELL -29.2192655 153.577371
## SMEAR -62.1384936 163.961186
## INFIL -171.6480455 70.015471
## LI 0.3514287 10.958965
## BLAST -4.6844284 4.511037
## TEMP -284.3010127 23.356799
以上为上篇的全部。至此已经学习很多内容，应该奖励自己！！之后的中篇会讲到建模和模型诊断（即衡量模型表现用于优化）和介绍 LASSO 回归，下篇会专门讲几种拟合优度的参数。
我们详细的项目解析，除了融合最实用的 R 语言代码和最前沿的统计学解释，更强调相应的医学背景！！我们希望让不同的学科背景的读者都能够“知其然，知其所以然”，真正明白每一步“为什么那么做”、”下次我自己会如何实现“。
同时为了大家能够更加习惯学术英文的使用，确保关键词汇遵循学术使用习惯，并且能直接用在学术写作中，我们在一些关键部分使用了中英双语注解。
建议加入我们的互助社群，可以帮助你进行临床医学研究设计，论文预审，R 语言安装和运行答疑；还可以获取课程的 RMD/MD 文件，以便直接复制到自己电脑上操作学习。
关于 R 语言学习，大家不必焦虑，更不要啃下一本教程 - Preparation isn't readiness!
无论学习还是使用 R，应该在遇到实际问题时用高效的人工智能辅助，可以广泛利用开发好的资源，比如找代码模板，甚至 “开口要” —— 去问 GPT 描述需求获取代码框架已经成为可能。
又比如公益在线工具 RTutor：将数据上传至 RTutor，用自然语言（中英文皆可）告诉它想要做什么数据处理，它会返回代码和结果，并且可以根据结果追问。以此为基础再针对性解决出现的问题。把使用和学习 R 语言放在解决实际问题的过程中才是最高效的。而类似 R 语言教程可以作为词典辅助使用。
hellomedstat 所有平台同名。
微信：Positivism_y
QQ 群：760447631
邮箱：hellomedstat@gmail.com
