41-R语言机器学习：逻辑回归与判别分析
2020-03-06 本文已影响0人
wonphen
《精通机器学习：基于R 第二版》学习笔记
1、数据理解与准备
> library(pacman)
> p_load(MASS, dplyr, ggplot2)
> data("biopsy")
> str(biopsy)
## 'data.frame': 699 obs. of 11 variables:
## $ ID : chr "1000025" "1002945" "1015425" "1016277" ...
## $ V1 : int 5 5 3 6 4 8 1 2 2 4 ...
## $ V2 : int 1 4 1 8 1 10 1 1 1 2 ...
## $ V3 : int 1 4 1 8 1 10 1 2 1 1 ...
## $ V4 : int 1 5 1 1 3 8 1 1 1 1 ...
## $ V5 : int 2 7 2 3 2 7 2 2 2 2 ...
## $ V6 : int 1 10 2 4 1 10 10 1 1 1 ...
## $ V7 : int 3 3 3 3 3 9 3 3 1 2 ...
## $ V8 : int 1 2 1 7 1 7 1 1 1 1 ...
## $ V9 : int 1 1 1 1 1 1 1 1 5 1 ...
## $ class: Factor w/ 2 levels "benign","malignant": 1 1 1 1 1 2 1 1 1 1 ...
ID ：样本编码
V1 ：细胞浓度
V2 ：细胞大小均匀度
V3 ：细胞形状均匀度
V4 ：边缘黏着度
V5 ：单上皮细胞大小
V6 ：裸细胞核（16个观测值缺失）
V7 ：平和染色质
V8 ：正常核仁
V9 ：有丝分裂状态
class ：肿瘤诊断结果，良性或恶性；这就是我们要预测的结果变量
重命名变量：
> # 删除ID列
> biopsy$ID <- NULL
>
> # 换成有意义的列名
> names(biopsy) <- c("thick", "u.size", "u.shape", "ashsn", "s.size", "nucl", "chrom", "n.nuc", "mit", "class")
>
> # 删除缺失值
> df <- na.omit(biopsy)
>
> # 用箱线图检查各个特征
> df2 <- reshape2::melt(df[, -11], value.name = "value")
> head(df2)
## class variable value
## 1 benign thick 5
## 2 benign thick 5
## 3 benign thick 3
## 4 benign thick 6
## 5 benign thick 4
## 6 malignant thick 8
> ggplot(df2, aes(class, value)) + geom_boxplot(outlier.color = "red") +
labs(y = "", x = "") + facet_wrap(~variable, ncol = 3)
箱线图检查各变量
从中位数的间隔距离和相关分布来看，我觉得可以有把握地认为nucl是一个重要特征。与之相反，不同class组的mit特征的中位数几乎没有区别，这说明它很可能是个无关特征。
1.1 检查相关性
> df[, 1:9] %>% cor %>% corrplot::corrplot.mixed()
相关性检查
从相关系数可以看出，我们会遇到共线性问题，特别是细胞大小均匀度(u.size)和细胞形状均匀度(u.shape)表现出非常明显的共线性。
1.2 划分训练集和测试集
> # 按70/30比例划分
> set.seed(123)
> ind <- sample(2, nrow(df), replace = T, prob = c(0.7, 0.3))
> train <- df[ind == 1, ]
> test <- df[ind == 2, ]
> str(test)
## 'data.frame': 209 obs. of 10 variables:
## $ thick : int 5 6 4 2 1 7 6 7 1 3 ...
## $ u.size : int 4 8 1 1 1 4 1 3 1 2 ...
## $ u.shape: int 4 8 1 2 1 6 1 2 1 1 ...
## $ ashsn : int 5 1 3 1 1 4 1 10 1 1 ...
## $ s.size : int 7 3 2 2 1 6 2 5 2 1 ...
## $ nucl : int 10 4 1 1 1 1 1 10 1 1 ...
## $ chrom : int 3 3 3 3 3 4 3 5 3 2 ...
## $ n.nuc : int 2 7 1 1 1 3 1 4 1 1 ...
## $ mit : int 1 1 1 1 1 1 1 4 1 1 ...
## $ class : Factor w/ 2 levels "benign","malignant": 1 1 1 1 1 2 1 2 1 1 ...
## - attr(*, "na.action")= 'omit' Named int 24 41 140 146 159 165 236 250 276 293 ...
## ..- attr(*, "names")= chr "24" "41" "140" "146" ...
> # 检查拆分后是否均衡
> table(train$class)
##
## benign malignant
## 302 172
> table(test$class)
##
## benign malignant
## 142 67
2、模型构建与评价
> full.fit <- glm(class ~ ., family = binomial, data = train)
> summary(full.fit)
##
## Call:
## glm(formula = class ~ ., family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.3397 -0.1387 -0.0716 0.0321 2.3559
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.4293 1.2273 -7.683 1.55e-14 ***
## thick 0.5252 0.1601 3.280 0.001039 **
## u.size -0.1045 0.2446 -0.427 0.669165
## u.shape 0.2798 0.2526 1.108 0.268044
## ashsn 0.3086 0.1738 1.776 0.075722 .
## s.size 0.2866 0.2074 1.382 0.167021
## nucl 0.4057 0.1213 3.344 0.000826 ***
## chrom 0.2737 0.2174 1.259 0.208006
## n.nuc 0.2244 0.1373 1.635 0.102126
## mit 0.4296 0.3393 1.266 0.205402
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 620.989 on 473 degrees of freedom
## Residual deviance: 78.373 on 464 degrees of freedom
## AIC: 98.373
##
## Number of Fisher Scoring iterations: 8
可以看到，只有两个特征的p值小于0.05（thick和nucl）。使用 confint() 函数可以对模型进行95%置信区间的检验。如下所示：
> confint(full.fit)
## 2.5 % 97.5 %
## (Intercept) -12.23786660 -7.3421509
## thick 0.23250518 0.8712407
## u.size -0.56108960 0.4212527
## u.shape -0.24551513 0.7725505
## ashsn -0.02257952 0.6760586
## s.size -0.11769714 0.7024139
## nucl 0.17687420 0.6582354
## chrom -0.13992177 0.7232904
## n.nuc -0.03813490 0.5110293
## mit -0.14099177 1.0142786
2.1 计算优势比
> exp(coef(full.fit))
## (Intercept) thick u.size u.shape ashsn s.size
## 8.033466e-05 1.690879e+00 9.007478e-01 1.322844e+00 1.361533e+00 1.331940e+00
## nucl chrom n.nuc mit
## 1.500309e+00 1.314783e+00 1.251551e+00 1.536709e+00
优势比可以解释为特征中1个单位的变化导致的结果发生比的变化。如果系数大于1，就说明当特征的值增加时，结果的发生比会增加。反之，系数小于1就说明，当特征的值增加时，结果的发生比会减小。在本例中，除u.size之外的所有特征都会增加对数发生比。
2.2 检查潜在的多重共线性问题
> car::vif(full.fit)
## thick u.size u.shape ashsn s.size nucl chrom n.nuc mit
## 1.235204 3.248811 2.830353 1.302178 1.635668 1.372931 1.523493 1.343145 1.059707
没有一个VIF值大于5，根据VIF经验法则，共线性看来不成为一个问题。
2.3 检查模型的准确性
> train.prob <- predict(full.fit, type = "response")
> train.p <- ifelse(train.prob < 0.5, 0L, 1L)
> y <- ifelse(df$class == "benign", 0L, 1L)
> train.y <- y[ind == 1]
> test.y <- y[ind == 2]
> caret::confusionMatrix(as.factor(train.y),as.factor(train.p))
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 294 8
## 1 7 165
##
## Accuracy : 0.9684
## 95% CI : (0.9483, 0.9822)
## No Information Rate : 0.635
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.9316
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9767
## Specificity : 0.9538
## Pos Pred Value : 0.9735
## Neg Pred Value : 0.9593
## Prevalence : 0.6350
## Detection Rate : 0.6203
## Detection Prevalence : 0.6371
## Balanced Accuracy : 0.9653
##
## 'Positive' Class : 0
> p_load(InformationValue)
> # 绘制ROC曲线
> plotROC(actuals = train.y, predictedScores = train.p, Show.labels = F, returnSensitivityMat = F)
ROC曲线
> # 查看错误率
> misClassError(actuals = train.y, predictedScores = train.p)
## [1] 0.0316
从ROC曲线和错误率上看我们干得相当不错，训练集上只有3.16%的预测错误率。如前所述，我们必须正确预测未知数据，换句话说，必须正确预测测试集。
> test.prob <- predict(full.fit, newdata = test, type = "response")
> test.p <- ifelse(test.prob < 0.5, 0, 1)
> misClassError(test.y, test.p)
## [1] 0.0239
看上去，我们使用全部特征建立的模型效果非常好，差不多98%的预测正确率，那么还有没有更好的方式来建立模型呢？
3、使用交叉验证的逻辑回归
模拟表明，LOOCV可以获得近乎无偏的估计，但是会有很高的方差。所以，大多数机器学习专家都建议将K的值定为5或10。
> p_load(bestglm)
> # bestglm包需要将因子变量编码为0或1，并且要求结果变量必须是最后一列，删除所有没有用的列
> df3 <- df %>%
+ mutate(y = case_when(class == "benign" ~ 0L, class == "malignant" ~ 1L)) %>%
+ select(-"class")
> head(df3)
## thick u.size u.shape ashsn s.size nucl chrom n.nuc mit y
## 1 5 1 1 1 2 1 3 1 1 0
## 2 5 4 4 5 7 10 3 2 1 0
## 3 3 1 1 1 2 2 3 1 1 0
## 4 6 8 8 1 3 4 3 7 1 0
## 5 4 1 1 3 2 1 3 1 1 0
## 6 8 10 10 8 7 10 9 7 1 1
> train.2 <- df3[ind == 1, ]
>
> # Xy = train.2 指的是我们已经格式化的数据框
> # IC = 'CV' 告诉程序使用的信息准则为交叉验证
> # CVArgs 是我们要使用的交叉验证参数
> # HTF 方法就是K折交叉验证，后面的数字 K = 10 指定了均分的份数
> # REP = 1 告诉程序随机使用等份并且只迭代一次
> # family = binomial，表示逻辑回归，如果family =gaussian 表示线性回归
> # bestglm不支持tibble数据框
> bestglm(Xy = train.2, IC = "CV", CVArgs = list(Method = "HTF", K = 10, REP = 1),
+ family = binomial)
## CV(K = 10, REP = 1)
## BICq equivalent for q in (7.16797006619085e-05, 0.273173435514231)
## Best Model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.8147191 0.90996494 -8.587934 8.854687e-18
## thick 0.6188466 0.14713075 4.206100 2.598159e-05
## u.size 0.6582015 0.15295415 4.303260 1.683031e-05
## nucl 0.5725902 0.09922549 5.770596 7.899178e-09
将thick,u.size,nucl特征放到 glm() 函数中，看看模型在测试集上表现如何。 因为predict() 函数不能用于bestglm生成的模型，所以下面的步骤是必需的：
> reduce.fit <- glm(class ~ thick + u.size + nucl, family = binomial, data = train)
> test.prob.2 <- predict(reduce.fit, newdata = test, type = "response")
> test.p.2 <- ifelse(test.prob.2 < 0.5, 0L, 1L)
> misClassError(test.y, test.prob.2)
## [1] 0.0383
精简了特征的模型和全特征模型相比，精确度略有下降。我们使用bestglm包信息准则为BIC的最优子集再试一次：
> bestglm(Xy = train.2, IC = "BIC", family = binomial)
## BIC
## BICq equivalent for q in (0.273173435514231, 0.577036596263757)
## Best Model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.6169613 1.03155250 -8.353391 6.633065e-17
## thick 0.7113613 0.14751510 4.822295 1.419160e-06
## ashsn 0.4537948 0.15034294 3.018398 2.541153e-03
## nucl 0.5579922 0.09848156 5.665956 1.462068e-08
## n.nuc 0.4290854 0.11845720 3.622282 2.920152e-04
thick,ashsn,nucl,n.nuc4个变量提供了最小的BIC评分，检测下预测效果：
> bic.fit <- glm(class ~ thick + ashsn + nucl + n.nuc, family = binomial, data = train)
> test.prob.3 <- predict(bic.fit, newdata = test, type = "response")
> test.p.3 <- ifelse(test.prob.3 < 0.5, 0L, 1L)
> misClassError(test.y, test.p.3)
## [1] 0.0239
错误率为2.39%，有所下降。那么问题来了：哪一个模型更好？在任何正常情况下，如果具有相同的泛化效果，经验法则会选择最简单的或解释性最好的模型。
4、判别分析
当遇到分类结果很确定的分类问题时，逻辑斯蒂回归的估计结果可能是不稳定的，即置信区间很宽，不同样本之间的估计值会有很大变化。判别分析不会受到这个问题的困扰，实际上，它会比逻辑回归做得更好，泛化能力更强。反之，如果特征和结果变量之间具有错综复杂的关系，判别分析在分类任务上的表现就会非常差。
4.1 线性判别分析
线性判别分析技术比逻辑斯蒂回归更具灵活性，同时还要时刻牢记偏差—方差权衡的问题。使用更有灵活性的技术可以得到偏差更小的结果，但很可能具有更高的方差。和很多灵活的技术一样，需要一个高鲁棒性的训练数据集来降低高分类方差。
> lda.fit <- lda(class ~ ., data = train)
> lda.fit
## Call:
## lda(class ~ ., data = train)
##
## Prior probabilities of groups:
## benign malignant
## 0.6371308 0.3628692
##
## Group means:
## thick u.size u.shape ashsn s.size nucl chrom
## benign 2.92053 1.304636 1.413907 1.324503 2.115894 1.397351 2.082781
## malignant 7.19186 6.697674 6.686047 5.668605 5.500000 7.674419 5.959302
## n.nuc mit
## benign 1.225166 1.092715
## malignant 5.906977 2.639535
##
## Coefficients of linear discriminants:
## LD1
## thick 0.19557291
## u.size 0.10555201
## u.shape 0.06327200
## ashsn 0.04752757
## s.size 0.10678521
## nucl 0.26196145
## chrom 0.08102965
## n.nuc 0.11691054
## mit -0.01665454
从结果可以看出，在分组先验概率中，良性概率大约为64%，恶性概率大约为36%。下面再看看分组均值，这是按类别分组的每个特征的均值。线性判别系数是标准线性组合，用来确定观测的判别评分的特征。评分越高，越可能被分入恶性组。
判别评分的直方图和密度图：
> plot(lda.fit, type = "both")
判别评分的直方图和密度图
可以看出，组间有些重合，这表明有些观测被错误分类。
查看LDA模型的准确率：
> train.lda.class <- predict(lda.fit)$class
>
> # 创建一个计算错误率的函数
> test_err <- function(df, pre) {
+ n <- 0
+ for (i in 1:nrow(df)) {
+ if (df$class[i] != pre[i])
+ n = n + 1
+ }
+ error <- n/nrow(df)
+ return(error)
+ }
> train.err <- test_err(train, train.lda.class)
> print(train.err)
## [1] 0.04008439
LDA模型的错误率为4.01%，比逻辑斯蒂回归模型差多了。看看在测试集上的表现：
> test.lda <- predict(lda.fit, newdata = test)$class
> test.err <- test_err(test, test.lda)
> print(test.err)
## [1] 0.03827751
从正确分类率的角度看，LDA模型表现得依然不如逻辑斯蒂回归模型（LDA模型：96%，逻辑斯蒂回归模型：98%）。
4.2 二次判别分析
> qda.fit <- qda(class ~ ., data = train)
> qda.fit
## Call:
## qda(class ~ ., data = train)
##
## Prior probabilities of groups:
## benign malignant
## 0.6371308 0.3628692
##
## Group means:
## thick u.size u.shape ashsn s.size nucl chrom
## benign 2.92053 1.304636 1.413907 1.324503 2.115894 1.397351 2.082781
## malignant 7.19186 6.697674 6.686047 5.668605 5.500000 7.674419 5.959302
## n.nuc mit
## benign 1.225166 1.092715
## malignant 5.906977 2.639535
> # 查看QDA模型在训练集和测试集上的准确率
> train.qda <- predict(qda.fit)$class
> train.qda.err <- test_err(train, train.qda)
> print(train.qda.err)
## [1] 0.04219409
> test.qda <- predict(qda.fit, newdata = test)$class
> test.qda.err <- test_err(test, test.qda)
> print(test.qda.err)
## [1] 0.05263158
根据错误率可以立即断定，QDA模型在训练集和测试集上表现得最差。
4.3 多元自适应回归样条方法
> p_load(earth)
> set.seed(111)
> # nfold = 5：5折交叉验证
> # ncross = 3：重复3次
> # degree = 1：使用没有交互项的加法
> # minspan = 1：每个输入特征只使用一个铰链函数。负值表示每个预测器的最大节数。这些是等距的。
> # 例如，minspan=-3允许每个预测器有三个均匀间隔的节点
> earth.fit <- earth(class ~ ., data = train, pmethod = "cv", nfold = 5, ncross = 3,
+ degree = 1, minspan = -1, glm = list(family = binomial))
> summary(earth.fit)
## Call: earth(formula=class~., data=train, pmethod="cv",
## glm=list(family=binomial), degree=1, nfold=5, ncross=3,
## minspan=-1)
##
## GLM coefficients
## malignant
## (Intercept) -6.5746417
## u.size 0.1502747
## ashsn 0.3058496
## s.size 0.3188098
## nucl 0.4426061
## n.nuc 0.2307595
## h(thick-3) 0.7019053
## h(3-chrom) -0.6927319
##
## GLM (family binomial, link logit):
## nulldev df dev df devratio AIC iters converged
## 620.989 473 81.9098 466 0.868 97.91 8 1
##
## Earth selected 8 of 10 terms, and 7 of 9 predictors (pmethod="cv")
## Termination condition: RSq changed by less than 0.001 at 10 terms
## Importance: nucl, u.size, thick, n.nuc, chrom, s.size, ashsn, ...
## Number of terms at each degree of interaction: 1 7 (additive model)
## Earth GRSq 0.8354593 RSq 0.8450554 mean.oof.RSq 0.8330208 (sd 0.0468)
##
## pmethod="backward" would have selected the same model:
## 8 terms 7 preds, GRSq 0.8354593 RSq 0.8450554 mean.oof.RSq 0.8330208
模型有8项，包括截距和7个预测变量。其中两个预测变量有铰链函数，这就是浓度（thick）和染色质（chrom）变量。如果浓度大于3，就会用系数0.7019乘以铰链函数的值，否则这一项就是0。对于染色质，如果它的值小于3，那么就用系数乘以铰链函数值，否则这一项就是0。
统计图展示了保持其他预测变量不变，某个预测变量发生变化时，相应变量发生的改变。可以清楚地看到铰链函数对浓度（thick）所起的作用。
> plotmo(earth.fit)
## plotmo grid: thick u.size u.shape ashsn s.size nucl chrom n.nuc mit
## 4 1 2 1 2 1 3 1 1
统计图
生成按类别标签分类的预测概率密度图：
> plotd(earth.fit)
预测概率密度图
查看变量的相对重要性：
> # nsubsets表示精简过程完成之后包含这个变量的模型的个数
> # gcv和 rss的范围都是0~100，值表示这个变量贡献的 gcv 和 rss 值的减少量
> evimp(earth.fit)
## nsubsets gcv rss
## nucl 7 100.0 100.0
## u.size 6 44.2 44.8
## thick 5 23.8 25.1
## n.nuc 4 15.1 16.8
## chrom 3 8.3 10.7
## s.size 2 6.0 8.1
## ashsn 1 2.3 4.6
模型在测试集上的错误率：
> test.earth.prob <- predict(earth.fit, newdata = test, type = "class")
> test.earth.err <- test_err(test, test.earth.prob)
> print(test.earth.err)
## [1] 0.02870813
错误率为2.875%，跟逻辑斯蒂回归模型差不多。
5、模型选择
> p_load(ROCR)
>
> # 最差bad
> bad.fit <- glm(class ~ thick, family = binomial, data = test)
> test.bad.prob <- predict(bad.fit, type = "response")
>
> pred.full <- prediction(test.prob, test$class)
> perf.full <- performance(pred.full, "tpr", "fpr")
> plot(perf.full, main = "ROC", col = 1)
>
> pred.bic <- prediction(test.prob.3, test$class)
> perf.bic <- performance(pred.bic, "tpr", "fpr")
> plot(perf.bic, col = 2, add = TRUE)
>
> # bad
> pred.bad <- prediction(test.bad.prob, test$class)
> perf.bad <- performance(pred.bad, "tpr", "fpr")
> plot(perf.bad, col = 3, add = TRUE)
>
> # earth
> test.earth.p <- predict(earth.fit, newdata = test, type = "response")
> pred.earth <- prediction(test.earth.p, test$class)
> perf.earth <- performance(pred.earth, "tpr", "fpr")
> plot(perf.earth, col = 4, add = TRUE)
>
> legend(0.6, 0.6, c("FULL", "BIC", "BAD", "EARTH"), 1:4)
ROC
可以看到，全特征模型、BIC模型和MARS模型基本上重叠在一起。显而易见， 糟糕的模型表现得和我们预想的一样差。
计算AUC：
> performance(pred.full, "auc")@y.values
## [[1]]
## [1] 0.9972672
> performance(pred.bic, "auc")@y.values
## [[1]]
## [1] 0.9944293
> performance(pred.bad, "auc")@y.values
## [[1]]
## [1] 0.8962056
> performance(pred.earth, "auc")@y.values
## [[1]]
## [1] 0.9952701
最高的AUC值是全模型的99.7%，BIC模型是99.4%，糟糕模型是89.6%，MARS模型是99.5%。所以，从各个方面来看，排除掉糟糕模型，其他几个模型在预测能力方面没有什么区别。我们该怎么做？一个简单的解决方案就是，将训练集和测试集重新随机化，把各种分析再做一遍，比如使用60/40的划分比例和一个不同的随机数种子。但是，如果我们还是得到相同的结果，又该怎么办？我认为，一个纯粹的统计学家会建议选择最简约的模型，但其他人可能更倾向于全特征模型。这就又归结到各种因素的权衡问题了，比如模型准确性与解释性、简约性与扩展性之间的权衡。在本章的例子中，我们完全可以选择更简单的模型，它具有和全模型一样的正确性。
