Xianxiong Ma
2020年3月14日
案例：威斯康辛大学威廉.沃尔伯格博士1990年发布了一个乳腺数据集，为明确肿瘤活检结果是良性还是恶性，研究者使用细针穿刺（FNA）技术收集样本，病理学家对活体组织进行检查，确定诊断结果(良性或者恶性)。良性的乳腺肿瘤并不危险，恶性肿瘤必须进行干预。误诊为恶性会导致昂贵但不必要的治疗费用，还会使患者背负巨大的心理和生理压力。误诊为良性（false negative），会使患者得不到应有的治疗，造成癌细胞扩散，引起过早死亡。乳腺癌患者的早期干预，可大大提高成活率，我们的任务就是开发尽可能准确地积极学习诊断算法，帮助医疗团队确定肿瘤性质。 数据集包含699个患者的组织样本，保存在有11个变量的数据框中，如下所示：
ID : 样本编码
V1 : 细胞浓度
V2 : 细胞大小均匀度
V3 : 细胞形状均匀度
V4 : 边缘黏着度
V5 : 单上皮细胞大小
V6 : 裸细胞核（16个观测值确实）
V7 : 平和染色质
V8 : 正常核仁
V9 : 有丝分裂状态
class : 肿瘤诊断结果，良性或恶性；这就是我们要预测的结果变量。
对9个特征进行了评分和编码，评分的范围是1~10。可以在R的MASS包中找到该数据框，名为biopsy。我们加载这个数据框，确认数据结构，将变量重命名为有意义的名称，还有删除有缺失项的观测，然后即可开始对数据进行可视化探索。
library(MASS)
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 ...
biopsy$ID = NULL
names(biopsy) = c("thick", "u.size", "u.shape", "adhsn", "s.size", "nucl", "chrom",
"n.nuc", "mit", "class")
names(biopsy)
## [1] "thick" "u.size" "u.shape" "adhsn" "s.size" "nucl" "chrom"
## [8] "n.nuc" "mit" "class"
biopsy.v2 <- na.omit(biopsy)
y <- ifelse(biopsy.v2$class == "malignant", 1, 0)
library(reshape2)
library(ggplot2)
biop.m <- melt(biopsy.v2, id.var = "class")
ggplot(data = biop.m, aes(x = class, y = value)) + geom_boxplot() + facet_wrap(~variable,
ncol = 3)
library(corrplot)
bc <- cor(biopsy.v2[, 1:9]) #create an object of the features
corrplot.mixed(bc)
set.seed(123) #random number generator
ind <- sample(2, nrow(biopsy.v2), replace = TRUE, prob = c(0.7, 0.3))
train <- biopsy.v2[ind == 1, ] #the training data set
test <- biopsy.v2[ind == 2, ] #the test data set
str(test) #confirm it worked
## '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 ...
## $ adhsn : 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
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
## adhsn 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
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
## adhsn -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
exp(coef(full.fit))
## (Intercept) thick u.size u.shape adhsn 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
library(car)
vif(full.fit)
## thick u.size u.shape adhsn s.size nucl chrom n.nuc
## 1.235204 3.248811 2.830353 1.302178 1.635668 1.372931 1.523493 1.343145
## mit
## 1.059707
没有一个VIF值大于5，根据VIF经验法则，共线性看来不成为一个问题，下面的任务就是进行特征选择，先别急，我们先编写一些代码，看看模型在训练集和测试集表现如何。
train.probs <- predict(full.fit, type = "response")
train.probs[1:5] #inspect the first 5 predicted probabilities
## 1 3 6 7 9
## 0.02052820 0.01087838 0.99992668 0.08987453 0.01379266
contrasts(train$class)
## malignant
## benign 0
## malignant 1
下一步需要评价模型在训练集上执行的效果，然后再评价它在测试集上的拟合程度，快速实现评价的方法是生成一个混淆矩阵。在后面的章节中，我们使用的混淆矩阵是由caret包实现的，InformationValue包也可以实现混淆矩阵，这时，我们需要用0和1来表示结果，函数区别良性结果和恶性结果使用的默认值是0.5，也就是说当概率大于0.5时，就认为这个结果是恶性的。
library(InformationValue)
trainY <- y[ind == 1]
testY <- y[ind == 2]
confusionMatrix(trainY, train.probs)
## 0 1
## 0 294 7
## 1 8 165
# optimalCutoff(trainY, train.probs)
misClassError(trainY, train.probs)
## [1] 0.0316
confusionMatrix(trainY, train.probs)
## 0 1
## 0 294 7
## 1 8 165
看上去我们干得相当不错，训练集上只有3.16%的预测错误率，如前所述，我们必须正确预测未知数据，换句话说，我们必须正确预测测试集，在测试集上建立混淆矩阵的方法和训练集一样。如下所示：
test.probs <- predict(full.fit, newdata = test, type = "response")
misClassError(testY, test.probs)
## [1] 0.0239
confusionMatrix(testY, test.probs)
## 0 1
## 0 139 2
## 1 3 65
看上去我们使用全部特征建立的模型效果非常好，差不多98%的预测正确率真是让人感动，但是我们还是要看看是否有改进的空间。想象一下。假如你和你的亲人是被误诊的患者会怎么样，就像前面提到的，后果很严重，出于这种考虑，还有可能更好的方式建立分类算法吗？下面我们来看一下。
交叉验证的目的是提高测试集的预测正确率，以及尽可能避免过拟合，K折交叉验证的做法就是将数据分成K个相等的等份，每个等份称为一个K子集（K-set）。算法每次留出一个子集，使得其余K-1个子集拟合模型，然后用模型在留出来的那个子集上做预测，将上面的k次验证的结果进行平均。可以使误差最小化，并且获得合适的特征选择，你也可以使用留一交叉验证方法，这里K等于N。模型表明。LOOCV可以获得近乎无偏的估计，但是会有很高的方差，所以大多数机器学习专家都建议将K的值定为5或者10.
bestglm包可以自动进行交叉验证，这个包依赖于我们在线性回归中使用的leaps包，交叉验证的语法和数据格式存在注意事项，所以我们按部就班的进行。
加载程序包之后，需要将结果编码为0或者1，如果结果仍为因子，程序将不起作用。使用这个程序包的另外一个要求是结果变量(或y)必须是最后一列，而且要删除所有没有用的列，通过以下代码新建一个数据框即可快速实现上述要求。
library(bestglm)
X <- train[, 1:9]
Xy <- data.frame(cbind(X, trainY))
bestglm(Xy = Xy, 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
在上面的代码中，Xy=Xy指的是我们已经格式化的数据框，IC=“CV"告诉程序使用的信息准则为交叉验证，CVArgs是我们要使用的交叉验证参数。HIF方法就是K折交叉验证，后面的数字K=10指定了均分的份数，REP=1告诉我们程序随机使用等分并且只迭代一次，像glm()函数一样，我们需要指定参数family=binomial，顺便说一句，如果指定参数family=gaussian，你就可以使用bestglm()函数进行线性回归。完成上面的交叉验证分析之后，会得到下面的结果：Best Model有三个特征，他们是thick, u.size和nucl。Morgan-Tatar搜索的结果表明，程序对所有可能的子集进行了简单而彻底的探索。
我们可以把这些特征放到glm()函数中，看看模型在测试集上表型如何，predict()函数不能用于bestglm生成的模型，所以下面的步骤是必须的：
reduce.fit <- glm(class ~ thick + u.size + nucl, family = binomial, data = train)
test.cv.probs = predict(reduce.fit, newdata = test, type = "response")
misClassError(testY, test.cv.probs)
## [1] 0.0383
confusionMatrix(testY, test.cv.probs)
## 0 1
## 0 139 5
## 1 3 62
精简了特征的模型与全特征模型相比，其精确度略有下降，但这并不意味着世界末日，我们可以用bestglm包来再试一次，这次使用信息准则为BIC的最优子集方法
bestglm(Xy = Xy, 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
## adhsn 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
bic.fit <- glm(class ~ thick + adhsn + nucl + n.nuc, family = binomial, data = train)
test.bic.probs = predict(bic.fit, newdata = test, type = "response")
misClassError(testY, test.bic.probs)
## [1] 0.0239
confusionMatrix(testY, test.bic.probs)
## 0 1
## 0 138 1
## 1 4 66
我们得到5个错误的预测，和全特征模型一样，那么问题来了，哪一个模型更好？在任何正常情况下，如果具有相同的泛化效果，经验法则会选择最简单或者解释性最好的模型，我们当然可以开启一个全面的分析，使用新的随机数以及新的训练集和测试集划分比例。
判别分析也是一项常用分类技术。当分类很确定时，判别分析可有效替代logistics回归。当分类结果很确定时，logistics回归的估计结果可能是不稳定的，即置信区间很宽，不同样本之间的估计值会有很大变化。判别分析会比logistic做的很好，泛化能力很强。反之，如果特征和结果变量之间具有错综复杂的关系，判别分析在分类任务上的表现就会非常差。
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 adhsn 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
## adhsn 0.04752757
## s.size 0.10678521
## nucl 0.26196145
## chrom 0.08102965
## n.nuc 0.11691054
## mit -0.01665454
plot(lda.fit, type = "both")
可以看出，组间有些重合，这表明有些观测被错误分类，LDA(线性判别分析)模型可以用predict()函数得到3种元素（class，posterior和x）。class元素是对良性或者恶性的观测，posterior是值为x的评分可能属于某个类别的概率，x是线性判别评分。通过下面的函数，我们仅提取恶性观测的概率：
train.lda.probs <- predict(lda.fit)$posterior[, 2]
misClassError(trainY, train.lda.probs)
## [1] 0.0401
confusionMatrix(trainY, train.lda.probs) #LDA模型训练集上不如逻辑斯蒂回归模型
## 0 1
## 0 296 13
## 1 6 159
test.lda.probs <- predict(lda.fit, newdata = test)$posterior[, 2]
misClassError(testY, test.lda.probs)
## [1] 0.0383
confusionMatrix(testY, test.lda.probs)
## 0 1
## 0 140 6
## 1 2 61
基于LDA模型在训练集上的糟糕表现，它在测试集上的表现比我预期的要好。从正确分类率的角度看，LDA模型表现得依然不如逻辑斯蒂回归模型（LDA:96%; 逻辑斯蒂回归：98%）。下面用二次判别分析（QDA）模型来拟合数据。在R中，QDA也是MASS包的一部分，函数qda(),建模过程直截了当，我们将模型存储在一个名为qda.fit的对象中。
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 adhsn 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
train.qda.probs <- predict(qda.fit)$posterior[, 2]
misClassError(trainY, train.qda.probs)
## [1] 0.0422
confusionMatrix(trainY, train.qda.probs)
## 0 1
## 0 287 5
## 1 15 167
test.qda.probs <- predict(qda.fit, newdata = test)$posterior[, 2]
misClassError(testY, test.qda.probs)
## [1] 0.0526
confusionMatrix(testY, test.qda.probs)
## 0 1
## 0 132 1
## 1 10 66
R中可以调整的参数非常少。在下面例子中，演示一种简单而有效MARS方法。如果你的需求很多，可以学习stephen Milborrow非常棒在线教程Notes on the earth package,通过以下链接可以学习MARS更多灵活用法: http://www.milbo.org/doc/earth-notes.pdf
介绍完MARS以后，我们开始学习其用法。你可以使用MDA包。但因为我是通过earth包学习的，所以我还是使用这个包。代码和前面的例子很相似，在那里我们使用了g1m()函数。但需要注意，一定要设定如何精简模型，并且使响应变量服从二元分布。我的设定是使用5折交叉验证来选择模型(pmethod="cv”. nfold=5) 。重复3次(ncross=3): 使用没有交互项的加法模型(degree=1)。每个输入特征只使用一个较链函数(minspan=l ) .在我使用的数据中，交互项和多个较链函数都会导致过拟合。当然，你的情况会有所不同。代码如下所示:
library(earth)
set.seed(1)
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.2980613
## u.size 0.1183107
## adhsn 0.3022766
## s.size 0.2783728
## nucl 0.4133847
## n.nuc 0.2208604
## h(3-thick) -0.3658432
## h(thick-3) 0.6670414
## h(3-chrom) -0.5960837
## h(chrom-3) 0.2151845
##
## GLM (family binomial, link logit):
## nulldev df dev df devratio AIC iters converged
## 620.989 473 80.9185 464 0.87 100.9 8 1
##
## Earth selected 10 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, adhsn, ...
## Number of terms at each degree of interaction: 1 9 (additive model)
## Earth GRSq 0.8327277 RSq 0.8452165 mean.oof.RSq 0.8390962 (sd 0.0403)
##
## pmethod="backward" would have selected:
## 8 terms 7 preds, GRSq 0.8354593 RSq 0.8450554 mean.oof.RSq 0.8325594
模型有8项，包括截距和7个预测变量。其中两个预测变量有较链函数，这就是浓度和染色质变量。如果浓度大于3.就会用系数0.7019乘以较链函数的值。否则这一项就是0。对于染色质，如果它的值小于3.那么就用系数乘以较链函数值，否则这一项就是0。还可以做出统计图。第一张图是用p1otmo()函数生成的，它展示了保持其他预测变量不变，某个预测变量发生变化时，响应变量发生的改变。你可以清楚地看到较链函数对浓度所起的作用:
plotmo(earth.fit)
## plotmo grid: thick u.size u.shape adhsn s.size nucl chrom n.nuc mit
## 4 1 2 1 2 1 3 1 1
plotd(earth.fit)
下面看看变量之间的相对重要性。先解释一下，nsubsets是精简过程完成之后包含这个变量的模型的个数。对于gcv和rss列，其中的值表示这个变量贡献的gcv和rss值的减少量(gcv 和Irss的范围都是0-100)
evimp(earth.fit)
## nsubsets gcv rss
## nucl 9 100.0 100.0
## u.size 8 43.9 44.9
## thick 7 23.1 25.2
## n.nuc 6 14.0 16.9
## chrom 5 6.1 10.8
## s.size 4 1.7 8.2
## adhsn 3 -5.2 4.8
test.earth.probs <- predict(earth.fit, newdata = test, type = "response")
misClassError(testY, test.earth.probs)
## [1] 0.0191
confusionMatrix(testY, test.earth.probs)
## 0 1
## 0 139 1
## 1 3 66
library(ROCR)
bad.fit <- glm(class ~ thick, family = binomial, data = train)
test.bad.probs <- predict(bad.fit, newdata = test, type = "response") #save probabilities
pred.full <- prediction(test.probs, test$class)
perf.full <- performance(pred.full, "tpr", "fpr")
plot(perf.full, main = "ROC", col = 1)
pred.bic <- prediction(test.bic.probs, test$class)
perf.bic <- performance(pred.bic, "tpr", "fpr")
plot(perf.bic, col = 2, add = TRUE)
pred.bad <- prediction(test.bad.probs, test$class)
perf.bad <- performance(pred.bad, "tpr", "fpr")
plot(perf.bad, col = 3, add = TRUE)
pred.earth <- prediction(test.earth.probs, 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)
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.9967416