[R 语言]一篇入门逻辑回归建模与可视化(克利夫兰心血管病临床数据)
学习一例克利夫兰心血管病临床数据的逻辑回归建模与可视化
阅读本文的前提是你已经大概理解了逻辑回归，不懂也没关系，既往文章的前言 （大约3分钟阅读）和我们的生动视频 （需10分钟观看）让你真正搞懂逻辑回归的原理。
学习材料：
这一课时将利用克利夫兰医学的公开数据 “Heart Disease” ，存放在 UCI 机器学习数据库（著名的 IRIS 数据集就来自于此）。
R 代码取自 Josh Starmer 的 GitHub 公开代码 。hellomedstat 对每一步做了详细的统计学和医学解释，并做了拓展延伸。
学习目标：
-
有多种变量类型的数据集清理，以用于 R 语言统计分析
-
多变量的单因素分析和参数理解
-
逻辑回归建模和结果可视化
下面开始正文:
1. 数据理解
先装好需要用到的 R 包 (package)。这两个包会在后面用来将数据作图（即“可视化”）：
library(ggplot2)
library(cowplot)
把网页的 csv 格式数据集传给 R。原数据集没有标题 (column names)，所以先另 header=False，否则第一行数据会被当成标题：
url <- "https://raw.githubusercontent.com/StatQuest/logistic_regression_demo/master/processed.cleveland.data"
data <- read.csv(url, header=FALSE)
用 head() 查看列名，head() 也会返回前几行数据，列名是 V1-V14：
head(data)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12
## 1 63 1 1 145 233 1 2 150 0 2.3 3 0.0
## 2 67 1 4 160 286 0 2 108 1 1.5 2 3.0
## 3 67 1 4 120 229 0 2 129 1 2.6 2 2.0
## 4 37 1 3 130 250 0 0 187 0 3.5 3 0.0
## 5 41 0 2 130 204 0 2 172 0 1.4 1 0.0
## 6 56 1 2 120 236 0 0 178 0 0.8 1 0.0
## V13 V14
## 1 6.0 0
## 2 3.0 2
## 3 7.0 1
## 4 3.0 0
## 5 3.0 0
## 6 3.0 0
为了方便数据分析，要定义符合数据含义的列名。查看 UCI 对数据集 “Heart Disease” 的介绍后，我们赋予以下列名：
colnames(data) <- c("age","sex","cp","trestbps","chol","fbs","restecg","thalach","exang","oldpeak","slope","ca","thal","hd")
根据变量的意义，明确属于哪种变量类型：
-
sex -> 0 = female, 1 = male.
性别：分类变量 -
cp -> chest pain; 1 = typical angina, 2 = atypical angina, 3 = non-anginal pain, 4 = asymptomatic.
胸痛分级：分类变量 -
trestbps -> resting blood pressure (in mm Hg).
静息状态血压：数值变量 -
chol -> serum cholestoral in mg/dl.
血清胆固醇：数值变量 -
fbs -> fasting blood sugar if less than 120 mg/dl, 1 = TRUE, 0 = FALSE.
空腹血糖是否小于 120 mg/dl：分类变量 -
restecg -> resting electrocardiographic results; 1 = normal, 2 = having ST-T wave, 3 = showing probable or definite left ventricular hypertrophy.
静息状态心电图类型：分类变量 -
thalach -> maximum heart rate achieved.
最大心率：数值变量 -
exang -> exercise induced angina, 1 = yes, 0 = no.
是否运动导致心绞痛：分类变量 -
oldpeak -> ST depression induced by exercise relative to rest.
心电图 ST 段参数：数值变量 -
slope -> the slope of the peak exercise ST segment; 1 = upsloping, 2 = flat, 3 = downsloping.
心电图 ST 段斜率类型：分类变量 -
ca -> number of major vessels (0-3) colored by fluoroscope.
主血管荧光显色数量：分类变量（数值为整数且取值较少的可以作为分类变量） -
thal -> this is short of thallium heart scan; 3 = normal (no cold spots), 6 = fixed defect (cold spots during rest and exercise), 7 = reversible defect (when cold spots only appear during exercise).
心肌灌注扫描情况：分类变量 -
hd (the predicted attribute) -> diagnosis of heart disease; 0 if less than or equal to 50% diameter narrowing, 1 if greater than 50% diameter narrowing.
有无心脏病（模型需要预测的结局）：分类变量
查看是否成功命名：
head(data)
## age sex cp trestbps chol fbs restecg
## 1 63 1 1 145 233 1 2
## 2 67 1 4 160 286 0 2
## 3 67 1 4 120 229 0 2
## 4 37 1 3 130 250 0 0
## 5 41 0 2 130 204 0 2
## 6 56 1 2 120 236 0 0
## thalach exang oldpeak slope ca thal hd
## 1 150 0 2.3 3 0.0 6.0 0
## 2 108 1 1.5 2 3.0 3.0 2
## 3 129 1 2.6 2 2.0 7.0 1
## 4 187 0 3.5 3 0.0 3.0 0
## 5 172 0 1.4 1 0.0 3.0 0
## 6 178 0 0.8 1 0.0 3.0 0
查看各变量数据类型是否正确：
str(data)
## 'data.frame': 303 obs. of 14 variables:
## $ age : num 63 67 67 37 41 56 62 57 63 53 ...
## $ sex : num 1 1 1 1 0 1 0 0 1 1 ...
## $ cp : num 1 4 4 3 2 2 4 4 4 4 ...
## $ trestbps: num 145 160 120 130 130 120 140 120 130 140 ...
## $ chol : num 233 286 229 250 204 236 268 354 254 203 ...
## $ fbs : num 1 0 0 0 0 0 0 0 0 1 ...
## $ restecg : num 2 2 2 0 2 0 2 0 2 2 ...
## $ thalach : num 150 108 129 187 172 178 160 163 147 155 ...
## $ exang : num 0 1 1 0 0 0 0 1 0 1 ...
## $ oldpeak : num 2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...
## $ slope : num 3 2 2 3 1 1 3 1 2 3 ...
## $ ca : chr "0.0" "3.0" "2.0" "0.0" ...
## $ thal : chr "6.0" "3.0" "7.0" "3.0" ...
## $ hd : int 0 2 1 0 0 0 3 0 2 1 ...
从 str() 返回结果发现有两个问题会影响数据分析 —— “ca” 和 “thal” 是 chr 的变量类型，说明该列错误包含了字符格式数据；所有的分类变量都被识别为了数值变量 “num”。进一步用 View(data) 查看原数据，发现是因为有“？”存在于“ca”和“thal”。所以，我们要进行一些数据清理，在这之前我们先小结一下...
小结：
命名时尽量避免大写字母，可以省去以后潜在的麻烦。查看是否命名成功，用 head()。
快速查看数据是否“大致正常”，用 str 函数。发现在 “ca” 和 “thal” 这两列变量是异常的 “character” 类型，查看原数据（用view(data)功能）发现是有 “？” 存在而导致。
以上数据量很小，我们很容易通过查看原数据来发现“？”。当数据量大时这种方法则不可行。最根本的解决办法是，在原始数据记录的时候，不要用“？”或其他符号代表缺失数据，而可以用类似999，或9999这种不可能与真实数据重复的数值。 在相应数据处理时，可以设置条件来过滤掉这些数值。
2. 数据清理
"NA" 代表 not available。它能让 R 把缺失的数据忽略，而不要形成干扰。所以，现在需要将数据中的问号更换为 “NA”：
data[data == "?"] <- NA
处理分类变量时要用到 “因子型变量 (factor)” ，它就是用来容纳分类变量的。根据临床意义，用 as.factor 将现实的分类数据（性别、胸痛程度等）从原来的数值型变量（numeric） 转换成分类变量（categorical）：
# 在 R 语言中，“==”表示“等于”，“=” 则用来赋值。
data[data$sex == 0,]$sex <- "F" # 当 sex 为 0 时，另变量取值为 F
data[data$sex == 1,]$sex <- "M"
data$sex <- as.factor(data$sex)
data$cp <- as.factor(data$cp)
data$fbs <- as.factor(data$fbs)
data$restecg <- as.factor(data$restecg)
data$exang <- as.factor(data$exang)
data$slope <- as.factor(data$slope)
注意，原始数据中的 “?” 使得相应变量类型为 character，被替换成 NA 后剩下数值仍为 character 类型（一般会在再次检查变量类型时发现，运行代码就是要反复确认一些很基本的东西）。
先将 character 类型转换为 integer 类型，再转换为 factor 类型，从而让 R 知道它们是分类变量:
data$ca <- as.integer(data$ca)
data$ca <- as.factor(data$ca)
data$thal <- as.integer(data$thal)
data$thal <- as.factor(data$thal)
对于 "hd"，它是因变量（dependent variable），是我们要预测的结局。有必要让该变量中的 0 显示为 “Healthy”，非 0 显示为 “Unhealthy”，可以用下列方式。最后用 str() 查看所有的替换结果：
data$hd <- ifelse(test=data$hd == 0, yes="Healthy", no="Unhealthy")
data$hd <- as.factor(data$hd)
str(data)
## 'data.frame': 303 obs. of 14 variables:
## $ age : num 63 67 67 37 41 56 62 57 63 53 ...
## $ sex : Factor w/ 2 levels "F","M": 2 2 2 2 1 2 1 1 2 2 ...
## $ cp : Factor w/ 4 levels "1","2","3","4": 1 4 4 3 2 2 4 4 4 4 ...
## $ trestbps: num 145 160 120 130 130 120 140 120 130 140 ...
## $ chol : num 233 286 229 250 204 236 268 354 254 203 ...
## $ fbs : Factor w/ 2 levels "0","1": 2 1 1 1 1 1 1 1 1 2 ...
## $ restecg : Factor w/ 3 levels "0","1","2": 3 3 3 1 3 1 3 1 3 3 ...
## $ thalach : num 150 108 129 187 172 178 160 163 147 155 ...
## $ exang : Factor w/ 2 levels "0","1": 1 2 2 1 1 1 1 2 1 2 ...
## $ oldpeak : num 2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...
## $ slope : Factor w/ 3 levels "1","2","3": 3 2 2 3 1 1 3 1 2 3 ...
## $ ca : Factor w/ 4 levels "0","1","2","3": 1 4 3 1 1 1 3 1 2 1 ...
## $ thal : Factor w/ 3 levels "3","6","7": 2 1 3 1 1 1 1 1 3 3 ...
## $ hd : Factor w/ 2 levels "Healthy","Unhealthy": 1 2 2 1 1 1 2 1 2 2 ...
注意！！上述部分分类变量中具备顺序等级（比如胸痛程度有等级之分，所以实际上是有序型变量 (ordinal variable)，而性别不具有等级，实际是名义型变量 (nominal variable)。用 as.factor 转换并不会自动生成顺序。
比如，即便胸痛 cp 是被 1～4 编码的，当用 is.ordered 查看这时 cp 是否为有序型（ordinal）变量，答案为FALSE，即胸痛没有被识别为有序型变量：
is.ordered(data$cp)
## [1] FALSE
现在查看有多少观测值含有缺失值（NA），如果少的话（不影响所需要的样本量；一般不超过样本的5%-10%），可以直接去除。本例只有 6 行（2%）缺失值，可以去除。否则，需要我们将其补充 （imputation）：
nrow(data[is.na(data$ca) | is.na(data$thal),])
## [1] 6
data[is.na(data$ca) | is.na(data$thal),]
## age sex cp trestbps chol fbs restecg
## 88 53 F 3 128 216 0 2
## 167 52 M 3 138 223 0 0
## 193 43 M 4 132 247 1 2
## 267 52 M 4 128 204 1 0
## 288 58 M 2 125 220 0 0
## 303 38 M 3 138 175 0 0
## thalach exang oldpeak slope ca thal
## 88 115 0 0.0 1 0 <NA>
## 167 169 0 0.0 1 <NA> 3
## 193 143 1 0.1 2 <NA> 7
## 267 156 1 1.0 2 0 <NA>
## 288 144 0 0.4 2 <NA> 7
## 303 173 0 0.0 1 <NA> 3
## hd
## 88 Healthy
## 167 Healthy
## 193 Unhealthy
## 267 Unhealthy
## 288 Healthy
## 303 Healthy
其实以上代码只查看了“ca”和“thal”两列。如果其他列（包括那些数值类型的列）也存在 NA，我们也应该提取，但是我们人工查看过原数据，只有“ca”和“thal”两列包含缺失值。
其实还有其他奇奇怪怪的空缺值定义，如何识别？不用担心这种不重要的问题。对于任何 R 语言的细节问题，等真正遇到的时候直接问会的人或者网上搜就可以。况且如果不会用 R 语言处理，上述问题基本都可以在 Excel 中解决，只要能达到目的就可以。
我们只要知道这些问题是存在的，就不会在实操中忽略了。
查看原来有多少个观测值：
nrow(data)
## [1] 303
去掉含有 NA 的观测值，确认剩下 297 个观测值：
data <- data[!(is.na(data$ca) | is.na(data$thal)),]
nrow(data)
## [1] 297
好了，看来即便来源很好的数据也会给我们挖很多坑，但是都会有办法识别和处理，这就够了。休息一下...
3. 探索性数据分析 EDA （单因素分析）
结束了原始数据 (raw data) 层面的“数据清洗”，下面要开始统计学的操作，先从探索性数据分析EDA（exploratory data analysis）开始。
为了进行有效的统计分析，需要先进行单因素分析（univariate analysis），比如检查数据分布情况，否则分析结果的可靠性会产生问题。针对本例：
- 查看是否每个“分类变量”下的各类别（如性别男和女）分别对应一定数量的两种结局（Healthy 和 Unhealthy）
- 如果上述某一分类对应的样本量太少，也会导致模型的不稳定，比如当换个样本进行模型验证时，该分类对应的样本较小的变化就会引起较大的统计值改变
- 1 和 2 这两个需求在本例用列联表（contingency table）实现。
- 用 xtabs() 生成列联表是一种很简便的方法。
本例使用 xtabs 查看各分类变量相对于结局的列联表（不适用于数值变量）：
xtabs(~ hd + sex, data=data)
## sex
## hd F M
## Healthy 71 89
## Unhealthy 25 112
xtabs(~ hd + cp, data=data)
## cp
## hd 1 2 3 4
## Healthy 16 40 65 39
## Unhealthy 7 9 18 103
xtabs(~ hd + fbs, data=data)
## fbs
## hd 0 1
## Healthy 137 23
## Unhealthy 117 20
xtabs(~ hd + restecg, data=data)
## restecg
## hd 0 1 2
## Healthy 92 1 67
## Unhealthy 55 3 79
xtabs(~ hd + exang, data=data)
## exang
## hd 0 1
## Healthy 137 23
## Unhealthy 63 74
xtabs(~ hd + slope, data=data)
## slope
## hd 1 2 3
## Healthy 103 48 9
## Unhealthy 36 89 12
xtabs(~ hd + ca, data=data)
## ca
## hd 0 1 2 3
## Healthy 129 21 7 3
## Unhealthy 45 44 31 17
xtabs(~ hd + thal, data=data)
## thal
## hd 3 6 7
## Healthy 127 6 27
## Unhealthy 37 12 88
另外常用的是 epiDisplay 包里的 tabpct() 函数，可以一下子实现频数分布，对应占比，和可视化结果。可以用以下代码运行函数：tabpct(data$hd, data$sex)
小结：
以上结果显示两种结局分布合理，现在大致确认了数据可以被用来进行有效的统计学分析，我们终于可以“建模”了。
4. 单因素逻辑回归模型
从列联表看出性别似乎与疾病结果有明显的相关性（association）。所以，我们先将性别作为自变量进行建模，叫做单因素逻辑回归（single-variate logistic regression）。用 R 语言自带的 glm 函数 即可。 glm 是指 generalized linear models，线性回归和逻辑回归都隶属其中。
需要在 glm 函数中用 “binomial” （二分类）来指定逻辑回归，因为预测的结局是二分类的 (Healthy & Unhealthy)。建完之后用 summary() 查看模型参数：
logistic <- glm(hd ~ sex, data=data, family="binomial")
summary(logistic)
##
## Call:
## glm(formula = hd ~ sex, family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -1.0438 0.2326 -4.488
## sexM 1.2737 0.2725 4.674
## Pr(>|z|)
## (Intercept) 7.18e-06 ***
## sexM 2.95e-06 ***
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 409.95 on 296 degrees of freedom
## Residual deviance: 386.12 on 295 degrees of freedom
## AIC: 390.12
##
## Number of Fisher Scoring iterations: 4
注意，相关性不代表因果关系（causal relationship）。此时我们没有考虑临床意义和混杂因素（confounder）。
比如，我们先不管性别在医学上是否真的影响心脏疾病发生，也不管是不是有其他因素影响了疾病发生，而伪装成性别因素（比如吸烟的男性显著多于女性，吸烟因素就可能“伪装”为性别因素）。如果想要进一步查看混杂因素，可以用逻辑回归实现。如果要明确因果关系，一般需要对照、前瞻性临床实验 (controlled, prospective clinical trials)，有时候观察性研究 (observational study) 也可以。
5. 理解逻辑回归参数
截距 Intercept
一般教材会说 Intercept 是所有变量取值为 0 时的结果。但是从实际意义来看，不一定存在所有变量 = 0 的现实情况，且有时候有的变量取值不可能为 0。
仅就本例来说，因为 R 语言识别字母顺序， 女性 Female 被当作性别变量的参考水平 (reference level)，因而被设为 0，但这不代表性别是有序型变量，这个 0 也不是来自于原始数据 sex 中的赋值 “0”。所以 intercept (-1.0438) 表示当变量取值为“女性”时，结局的发生几率取 log 对数之后的数值。在建模的学术语言中，“结局 (Event) ” 特指感兴趣的结果，这里指 Unhealthy。
引用一段文献 中的话帮助理解 Intercept:
The statistical program first calculates the baseline odds of having the outcome versus not having the outcome without using any predictor. This gives us the constant (also known as the intercept).
这句话是说，应用程序首先计算在不使用任何预测因子 (predictor) 的情况下，出现结局与不出现结局的基线几率。这样就得到了常数 Intercept（截距）。
回归系数 Coefficient
1.2737 是性别因素的回归系数。当作为自然常数 e 的指数，结果为 exp(1.2737) = 3.57，也就是性别的优势比 OR，可以理解为性别给结果带来的影响大小：男性和女性（参考水平）比，被预测为不健康的几率是 3.57 倍，英文表述为 "the odds of ...(the event) is 3.57 times than female"。实际上也可以通过列联表的数值算出：(112 / 89) / (25/71) = 3.57。
回归系数更本质的含义（尤其在纳入多个自变量的时候）可参考我们的相应视频课程来理解。当纳入模型的自变量为 3 个及以上时，通过列联表计算则不再适用，因为每个列联表只单独考虑了一个自变量，没有考虑其他自变量的影响。需要用逻辑回归来计算调整后优势比 (adjusted OR, aOR)，在 R 语言中就是 "glm()"。
用女性的发病几率计算，验证上述对截距的解释：
female.log.odds <- log(25 / 71)
female.log.odds
## [1] -1.043804
用是否暴露于风险因素（此处“风险因素”为“男性”），验证上述对性别的回归系数的解释：
male.log.odds.ratio <- log((112 / 89) / (25/71))
male.log.odds.ratio
## [1] 1.273667
得到模型参数截距和回归系数也就确定了模型。接下来要评估模型表现，先计算模型的 伪 R 方 (McFadden's Pseudo R-squared)。该数值可以理解为模型能够解释结果的程度。在我们的另一课 中有详细讲解。
评估模型表现的伪 R 方计算公式为：
R squared = 1 - (log_likelihood(proposed_model_deviance)/log_likelihood(null_model_deviance))
"proposed_model_deviance/null_model_deviance" 即被衡量的模型的偏差程度除以零模型 (null model) 的偏差程度。零模型可以理解为完全不动脑筋的情况下给出的一个最最简单的模型，是一个常数。伪 R 方取值区间为 [0,1], 1 代表 “完美” 模型，即所有的点都落在拟合曲线上 (perfectly fitted)。
统计学上，上面公式里的 “deviance” 被定义为 −2 乘以 log(likelihood)，所以：
null_model_devaince = 2*(0 - LogLikelihood(null_model)) = -2*LogLikihood(null_model)
proposed_model_deviacne = 2*(0 - LogLikelihood(proposed_model)) = -2*LogLikelihood(proposed_model)
上述伪 R 方的计算用 R 代码实现就是：
ll.null <- logistic$null.deviance/-2 # "ll" 表示 log(likelihood)
ll.proposed <- logistic$deviance/-2
(ll.null - ll.proposed) / ll.null
# 这三行缩减为一行代码就是 r_squared <- 1 - as.numeric(logLik(proposed_model)/logLik(null_model))，
## [1] 0.05812569
ll.null - ll.proposed 实际就是似然比检验 的统计量，它的结果是符合卡方分布的，所以它的数值可以被用来做统计检验：
ll.null <- logistic$null.deviance/-2
ll.proposed <- logistic$deviance/-2
1 - pchisq(2*(ll.proposed - ll.null), df=1)
# 或者也可以只用一行代码 1 - pchisq((logistic$null.deviance - logistic$deviance), df=1)
## [1] 1.053157e-06
小结：
似然比检验得到了一个非常小的 p 值，说明当前模型比零模型更好。但是结合前面非常小的伪 R 方 （伪 R 方越大越好），又说明当前模型还不够好。假设检验就是统计学的精髓，但统计学的结论不一定等于现实意义上的结论，所以这个模型到底好不好，目前还不确定。
似然比检验也可以直接用函数 anova(model.1, model.2, test=“Chisq”)。这里 model.1 和 model.2 就是做比较的两个模型，如果另其中一个为 null.model，那么就是上述的似然比检验；如果两个模型都不是零模型，那么计算的就是一个复杂一点的模型是否比一个简单一点的模型更优，其中，简单模型纳入的自变量是复杂模型自变量的子集 (subset)。详见另一课时似然比检验 。
6. 作图/可视化
为什么说当前模型还不够好？光看伪 R 方和 p 值都不够直观，我们需要将数据可视化。本例可视化是由 ggplot 函数实现（R 包：ggplot2）。这个包非常实用和主流，可以多了解。
我们将结局预测值（纵坐标）和性别信息（横坐标）赋值给一个数据框 (data frame)，并用 $fitted.values 提取预测值：
predicted.data <- data.frame(
probability.of.hd=logistic$fitted.values,
sex=data$sex)
ggplot(data=predicted.data, aes(x=sex, y=probability.of.hd)) +
geom_point(aes(color=sex), size=5) +
xlab("Sex") +
ylab("Predicted probability of getting heart disease")
上面的图只得到两个点（其实是所有的点都重合在这两个点上）。用代码 View(predicted.data)）会看到，男性和女性对应的预测值分别是 0.5572139 和 0.2604167。这显然是不符合现实情况的，因为不可能所有男性都是同一水平的高风险，所有女性都是同一水平的低风险。
由于这个例子比较夸张，所以如果用 xtabs 函数做一个 2x2 列联表来看频数分布，也很容易体现该模型与结局明显不符：
xtabs(~ probability.of.hd + sex, data=predicted.data)
## sex
## probability.of.hd F M
## 0.260416666667239 96 0
## 0.55721393034826 0 201
7. 多元逻辑回归
现在已经尝试并理解了简单的建模过程，现在用相同的思路将所有自变量纳入，再建一个多元逻辑回归模型 (multiple logistic regression)。
用 glm() 建模，“.” 代表除了 hd 以外的所有变量：
logistic <- glm(hd ~ ., data=data, family="binomial")
summary(logistic)
##
## Call:
## glm(formula = hd ~ ., family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -6.253978 2.960399 -2.113
## age -0.023508 0.025122 -0.936
## sexM 1.670152 0.552486 3.023
## cp2 1.448396 0.809136 1.790
## cp3 0.393353 0.700338 0.562
## cp4 2.373287 0.709094 3.347
## trestbps 0.027720 0.011748 2.359
## chol 0.004445 0.004091 1.087
## fbs1 -0.574079 0.592539 -0.969
## restecg1 1.000887 2.638393 0.379
## restecg2 0.486408 0.396327 1.227
## thalach -0.019695 0.011717 -1.681
## exang1 0.653306 0.447445 1.460
## oldpeak 0.390679 0.239173 1.633
## slope2 1.302289 0.486197 2.679
## slope3 0.606760 0.939324 0.646
## ca1 2.237444 0.514770 4.346
## ca2 3.271852 0.785123 4.167
## ca3 2.188715 0.928644 2.357
## thal6 -0.168439 0.810310 -0.208
## thal7 1.433319 0.440567 3.253
## Pr(>|z|)
## (Intercept) 0.034640 *
## age 0.349402
## sexM 0.002503 **
## cp2 0.073446 .
## cp3 0.574347
## cp4 0.000817 ***
## trestbps 0.018300 *
## chol 0.277253
## fbs1 0.332622
## restecg1 0.704424
## restecg2 0.219713
## thalach 0.092781 .
## exang1 0.144267
## oldpeak 0.102373
## slope2 0.007395 **
## slope3 0.518309
## ca1 1.38e-05 ***
## ca2 3.08e-05 ***
## ca3 0.018428 *
## thal6 0.835331
## thal7 0.001141 **
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 409.95 on 296 degrees of freedom
## Residual deviance: 183.10 on 276 degrees of freedom
## AIC: 225.1
##
## Number of Fisher Scoring iterations: 6
计算伪 R 方和似然比检验的 p 值。同样地，McFadden's Pseudo R^2 = [ LL(Null) - LL(Proposed) ] / LL(Null)：
ll.null <- logistic$null.deviance/-2
ll.proposed <- logistic$deviance/-2
(ll.null - ll.proposed) / ll.null
这次的数值好很多，达到 55.34%，可以认为取值 50% 以上属于拟合度较好的模型（但不绝对）：
## [1] 0.5533531
计算 p 值：
1 - pchisq(2*(ll.proposed - ll.null), df=(length(logistic$coefficients)-1))
## [1] 0
该 p 值接近于 0，说明当前模型的拟合度很不错，有必要作图直观地看一下。先将需要作图的数据打包进一个 data.frame 里：
predicted.data <- data.frame(
probability.of.hd=logistic$fitted.values,
hd=data$hd)
这次纵坐标依然是各个预测值，即 "probability.of.hd"，但横坐标是什么呢？我们在视频课里讲过 “Z” 值，Z = β0 + β1X1 + β2X2 +...，其中 β0 就是截距 Intercept，其余 β 对应每个变量 X 的回归系数，所以 Z 也就是“所有影响因素的作用总和”。但是用 Z 作为横坐标比较麻烦，我们可以将每个样本的 Z 值大小排序，用序号代替其具体数值来做横坐标。这一排序实际上就是将预测值排序，因为 Z 值与预测值是成正比的。
将这个数据框的预测值（probability.of.hd）按升序排序：
predicted.data <- predicted.data[
order(predicted.data$probability.of.hd, decreasing=FALSE),]
再将排序后的每一行序号作为预测值的秩次 (rank)，这在之后会被作为横坐标：
predicted.data$rank <- 1:nrow(predicted.data)
最后，我们将秩次作为横坐标（1～297），预测值作为纵坐标：
ggplot(data=predicted.data, aes(x=rank, y=probability.of.hd)) +
geom_point(aes(color=hd), alpha=1, shape=4, stroke=2) +
xlab("Index") +
ylab("Predicted probability of getting heart disease")
除了横纵坐标，ggplot 帮我们展示了第三个重要的信息，就是每个点对应的实际结局，用红色和蓝色形成鲜明对比：
养成保存图片的习惯，用 ggsave(file_name.format)：
ggsave("heart_disease_probabilities.pdf")
## Saving 7 x 7 in image
注意，这条曲线并不是该逻辑回归模型本来的曲线，因为它的横坐标只是 Z 的秩次，这的确会影响曲线的形状。但还是能呈现出两种结局（红和蓝）的分布具有显著趋势，所以用秩次做横坐标并不影响判断。
查看实际结果和预测结果的一致性还可以用一个常规思路，就是直接对比预测值和实际值，看散点图(scatter plot) 的分布情况，可以用以下代码实现：
probability.of.hd=logistic$fitted.values
plot_logistic <- data.frame(Actual = as.numeric(data$hd), PredictedProb = probability.of.hd)
ggplot(plot_logistic, aes(x = Actual, y = PredictedProb)) +
geom_jitter(alpha = 0.5) +
labs(x = "Actual Outcome", y = "Predicted Probability")
可以看出点倾向于集中在左下和右上：
注释：因为“Actual”真实的取值只有两个，会导致所有的点都集中在一起，所以上述代码用 jitter 功能 打散这些点，使得它们分布在以实际取值为中心，左右 0.5 的坐标范围内。
8. 模型应用之阈值
在临床应用中，我们可能不仅用模型给出疾病风险预测值，还要进一步判定其对应的二元结局（binary outcomes），这个结果只能为 “Healthy“ 或 ”Unhealthy”，而不是结局 （Unhealthy) 的概率。这时，需要选择一个临界值 (cutoff value) 作为阈值。如果选择 0.50，那么预测值只要大于 0.50，我们就给出 Unhealthy 的诊断。这时，我们选择阈值可以依据前面的曲线或是上述补充的散点图，从视觉上找一个满意的界限，作为粗略估计。
选择阈值要有明确定义的模型使用场景和目的。可以是为了最准确地区分两个结局 (discrimination)，或是为了最大程度避免漏诊 (high sensitivity)，等等。不同的诉求选择的阈值是不一样的。当然，有很多方法和衡量指标帮助我们实现阈值选择。
结语
恭喜你，已经完成了一次逻辑回归模型建立，还可以思考下面这些不同学科角度的问题：
- 有关数据处理的问题：
- 如果将分类变量赋予该有的等级，上述代码和结果会有哪些变化？
- 关于缺失值填充，本例认为缺失比例不足 2%，也有学者认为不该用缺失值比例来决定 。
- 有关统计分析的问题：
为什么叫 “McFadden's pseudo R-squared” - 伪 R 方？有兴趣了解的可以去看Josh Starmer的生动解释："Logistic Regression Details Pt 3: R-squared and p-value" (YouTube)，也可以看我们的另一课 讲解。
知道 null model, full model, saturated model 各自的含义吗？
常规的逻辑回归统计分析是先进行单因素分析，单独计算每个因素/变量和结果的关联，然后一般会将“表现较好”（如 p < 0.2）的因素筛选出来纳入后面的多因素分析（multivariable analysis），才会得到调整后的优势比。所谓调整后的优势比，就是在纳入所有想要纳入的因素之后，各因素与结果的关联程度。上述这种变量筛选方式只是统计学角度的方法论，在研究的全局下不一定合适，我们的视频课有讲述。
- 有关医学应用的问题：
关于阈值选择，需要结合模型应用端的需求，一定是 0.50 吗？
ROC 曲线也可以帮助选择阈值，因为它描绘了每个阈值情形下对应的模型敏感度和特意度，但 ROC 是最好的模型选择指标吗？
发散一下，针对本案例，我们还可以选择两个阈值，将结果分为三类，中间一类的人群需要用进一步措施来明确诊断或预后。这样的案例（分析自一篇JAMA文献）在我们的视频课中有生动讲述。
硬广
hellomedstat 用 "R 语言-统计-医学" 三位一体的讲述，希望让不同的学科背景的读者都能够 “知其然，知其所以然”，真正明白每一步“为什么那么做”、”下次我自己会如何实现“。
建议加入我们的互助社群(文末二维码)。有群主为核心的精锐团队近距离回答临床医学研究的问题，并帮助预审论文，更有技术成员进行 R 语言答疑；入群还可以获取课程的 RMD/MD 文件，以便直接复制到自己电脑上操作学习。
关于 R 语言学习，大家不必焦虑，更不要啃下一本教程 - Preparation isn't readiness!
无论学习还是使用 R，应该在遇到实际问题时用高效的人工智能辅助，可以广泛利用开发好的资源，比如找代码模板，甚至 “开口要” —— 去问 GPT 描述需求获取代码框架已经成为可能。总之，如果你是医学或流统背景，没有必要从零开始学 R 语言！
hellomedstat 全平台同名，欢迎小伙伴一起讨论学习！
微信：Positivism_y
QQ 群：760447631
邮箱：hellomedstat@gmail.com