温馨提示：代码运行过程中，有任何问题可留言，定答复！

#有一种方法可以看出逻辑斯蒂回归和线性回归之间的联系， #逻辑斯蒂回归就是以对数发生比 为响应变量进行线性拟合，即log(P(Y)/1  P(Y)) =B0+B1x。 #这里的系数是通过极大似然估计得到 的，而不是通过OLS。 #极大似然的直观意义就是，我们要找到一对B0和B1的估计值，使它们产生 的对观测的预测概率尽可能接近Y的实际观测结果，这就是所谓的似然性。 #请记住，逻辑斯蒂回归是一项强大的技术，可以预测分类问题，通常也是对分类问题进行建 模的起点。 #因此，对于本章中将要面临的商业问题，我们使用逻辑斯蒂回归作为解决问题的首选 方法。 #为进行数据准备，我们加载这个数据框， 确认数据结构，将变量重命名为有意义的名称，还要删除有缺失项的观测，然后即可开始对数据 进行可视化探索。 library(MASS) data(biopsy) str(biopsy) #删除ID列 biopsy$ID = NULL #重命名变量，并确认操作结果 names(biopsy) = c("thick", "u.size", "u.shape", "adhsn", "s.size", "nucl", "chrom", "n.nuc", "mit", "class") names(biopsy) #删除那些有缺失数据的观测。 #既然只有16个观测有缺失数据，所以删除它们不 会有什么危险，因为它们只占所有观测的2%。 biopsy.v2 <- na.omit(biopsy) #根据分析数据时使用的R包的不同，结果变量也可能要求是数值型，比如0或者1。 #为了满 足这个要求，可以创建一个变量y，用0表示良性，用1表示恶性 y <- ifelse(biopsy.v2$class == "malignant", 1, 0) library(reshape2) library(ggplot2) #melt()函数建立一个数据框，这样做是为了在特征融合后就建立一个箱线图矩阵，使我们可以 轻松地进行下面的可视化检查： biop.m <- melt(biopsy.v2, id.var = "class") ggplot(data = biop.m, aes(x = class, y = value)) + geom_boxplot() + facet_wrap(~ variable, ncol = 3) #通过对上面箱线图的研究和判断，我们发现，很难确定哪个特征在分类算法中是重要的。 #但 是，从中位数的间隔距离和相关分布来看，我觉得可以有把握地认为nuclei是一个重要特征。 #与 之相反，不同class组的mitosis特征的中位数几乎没有区别，这说明它很可能是个无关特征。 #逻辑斯蒂回归中的共线性也会使估计发生偏离,加载corrplot包，检查特征之间的 相关性 library(corrplot) bc <- cor(biopsy.v2[ ,1:9]) #create an object of the features corrplot.mixed(bc) #有多种方式可以将数据恰当地划分成训练集和测试集：50/50、60/40、70/30、80/20，诸如此 类。 #你应该根据自己的经验和判断选择数据划分方式。在本例中，我选择按照70/30的比例划分 数据。 set.seed(123) #random number generator ind <- sample(2, nrow(biopsy.v2), replace = TRUE, prob = c(0.7, 0.3)) sss<-ind[ind==2] train <- biopsy.v2[ind==1, ] #the training data set test <- biopsy.v2[ind==2, ] #the test data set str(test) #confirm it worked #查看两个数据集的结果变量是否均衡 table(train$class) table(test$class) #我们必须在函数中使用family = binomial这个参 数，告诉R运行逻辑斯蒂回归，而不是其他版本的广义线性模型。 #首先在训练数据集上使 用所有特征建立一个模型，看看这个模型在测试数据集上运行的效果。 full.fit <- glm(class ~ ., family = binomial, data = train) summary(full.fit) #使用confint()函数可以对模型系数进行95%置信区间的检验 confint(full.fit) #请注意，两个显著特征的置信区间不包括0。 #对于逻辑斯蒂回归模型的系数，你不能解释为“当 X改变1个单位时Y会改变多少”。 #这时，优势比的作用就体现出来了。对数函数log(P(Y)/1-P(Y)) = B0 + B1x的B系数可以通过exponent(beta)指数函数转化为优势比。 exp(coef(full.fit)) #优势比可以解释为特征中1个单位的变化导致的结果发生比的变化。如果系数大于1，就说明 当特征的值增加时，结果的发生比会增加。 #反之，系数小于1就说明，当特征的值增加时，结果 的发生比会减小。 多重共线性 library(car) vif(full.fit) #没有一个VIF值大于5，根据VIF经验法则，共线性看来不成为一个问题。 train.probs <- predict(full.fit, type = "response") train.probs[1:5] #inspect the first 5 predicted probabilities contrasts(train$class) library(InformationValue) trainY <- y[ind==1] testY <- y[ind==2] confusionMatrix(trainY, train.probs) # optimalCutoff(trainY, train.probs) misClassError(trainY, train.probs) 在测试集进行预测 test.probs <- predict(full.fit, newdata = test, type = "response") misClassError(testY, test.probs) confusionMatrix(testY, test.probs) 使用交叉验证的逻辑斯蒂回归 #交叉验证的目的是提高测试集上的预测正确率，以及尽可能避免过拟合。 #K折交叉验证的做 法是将数据集分成K个相等的等份，每个等份称为一个K子集（K-set）。 #算法每次留出一个子集， 使用其余K 1个子集拟合模型，然后用模型在留出的那个子集上做预测。 #将上面K次验证的结果 进行平均，可以使误差小化，并且获得合适的特征选择。 #你也可以使用留一交叉验证方法，这 里的K等于N。 #模拟表明，LOOCV可以获得近乎无偏的估计，但是会有很高的方差。 #所以，大多 数机器学习专家都建议将K的值定为5或10。 #bestglm包可以自动进行交叉验证，这个包依赖于我们在线性回归中使用过的leaps包。 library(leaps) library(shape) library(bestglm) #加载程序包之后，需要将结果编码成0或1。如果结果仍为因子，程序将不起作用。 #使用这个 程序包的另外一个要求是结果变量（或y）必须是后一列，而且要删除所有没有用的列。 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) #在上面的代码中，Xy = Xy指的是我们已经格式化的数据框， #IC = "CV"告诉程序使用的信 息准则为交叉验证， #CVArgs是我们要使用的交叉验证参数 #。HTF方法就是K折交叉验证，后面的 数字K = 10指定了均分的份数，参数REP = 1告诉程序随机使用等份并且只迭代一次。 #像glm() 函数一样，我们需要指定参数family = binomial。 #顺便说一句，如果指定参数family = gaussian，你就可以使用bestglm()函数进行线性回归。 #我们可以把这些特征放到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) confusionMatrix(testY, test.cv.probs) #精简了特征的模型和全特征模型相比，精确度略有下降，但这并不是“世界末日”。 #我们可 以用bestglm包再试一次，这次使用信息准则为BIC的优子集方法 bestglm(Xy = Xy, IC = "BIC", family = binomial) 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) confusionMatrix(testY, test.bic.probs) #在任何正 常情况下，如果具有相同的泛化效果，经验法则会选择简单的或解释性好的模型。 #判别分析又称费舍尔判别分析，也是一项常用的分类技术。 #当分类很确定时，判别分析可以 有效替代逻辑斯蒂回归。 #当你遇到分类结果很确定的分类问题时，逻辑斯蒂回归的估计结果可能 是不稳定的，即置信区间很宽，不同样本之间的估计值会有很大变化（James，2013）。 #判别分析 不会受到这个问题的困扰，实际上，它会比逻辑回归做得更好，泛化能力更强。 #反之，如果特征 和结果变量之间具有错综复杂的关系，判别分析在分类任务上的表现就会非常差。 判别分析 线性判别分析（LDA） library(MASS) lda.fit <- lda(class ~ ., data = train) lda.fit #从结果可以看出，在分组先验概率中，良性概率大约为64%，恶性概率大约为36%。 #下面再 看看分组均值，这是按类别分组的每个特征的均值。 #线性判别系数是标准线性组合，用来确定观 测的判别评分的特征。评分越高，越可能被分入恶性组。 #对LDA模型使用plot()函数，可以画出判别评分的直方图和密度图 plot(lda.fit, type="both") #可以看出，组间有些重合，这表明有些观测被错误分类。 train.lda.probs <- predict(lda.fit)$posterior[, 2] misClassError(trainY, train.lda.probs) confusionMatrix(trainY, train.lda.probs) #predict(lda.fit, newdata = test) test.lda.probs <- predict(lda.fit, newdata = test)$posterior[, 2] misClassError(testY, test.lda.probs) confusionMatrix(testY, test.lda.probs) #基于LDA模型在训练集上的糟糕表现，它在测试集上的表现比我预期的要好。 #从正确分类率 的角度看，LDA模型表现得依然不如逻辑斯蒂回归模型（LDA模型：96%，逻辑斯蒂回归模型：98%）。 二次判别分析（QDA） qda.fit <- qda(class ~ ., data =train) qda.fit train.qda.probs <- predict(qda.fit)$posterior[, 2] misClassError(trainY, train.qda.probs) confusionMatrix(trainY, train.qda.probs) test.qda.probs <- predict(qda.fit, newdata = test)$posterior[, 2] misClassError(testY, test.qda.probs) confusionMatrix(testY, test.qda.probs) #根据混淆矩阵可以立即断定，QDA模型在训练数据集上表现得差。 #QDA在测试集上的分 类效果也很差，有11个预测错误，而误为恶性的比例尤其高。 #可以学习Stephen Milborrow制作的非常棒的在线教程Notes on the earth package，通过以下链接可以学习MARS的更多灵活用法： #http://www.milbo.org/doc/earth-notes.pdf 多元自适应回归样条 MARS# 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) ) #一定要设定如何精简模型，并且使响应变量服从二元分布: #此处，使用5折交叉验证 来选择模型（pmethod="cv"，nfold = 5），重复3次（ncross = 3）； #使用没有交互项的加法 模型（degree = 1），每个输入特征只使用一个铰链函数（minspan = 1）。 #在我使用的数据 中，交互项和多个铰链函数都会导致过拟合。当然，你的情况会有所不同。 summary(earth.fit) #。MARS模型的生成过 程是： #先根据数据建立一个标准线性回归模型，它的响应变量在内部编码为0和1； #对特征（变量） 进行精简之后，生成终模型，并将其转换为一个GLM模型。 #MARS 算法采用形式为max(0,x-t)和max(0,t-x)的两组函数作为基对输入信号进行展开 plotmo(earth.fit) #展示了保持其他预测变量不变， 某个预测变量发生变化时，响应变量发生的改变。 plotd(earth.fit) #可以生成按类别标签分类的预测概率密度图 evimp(earth.fit) #看看变量之间的相对重要性。先解释一下，nsubsets是精简过程完成之后包含这个变 量的模型的个数。 #对于gcv和rss列，其中的值表示这个变量贡献的gcv和rss值的减少量（gcv 和rss的范围都是0~100） test.earth.probs <- predict(earth.fit, newdata = test, type = "response") misClassError(testY, test.earth.probs) confusionMatrix(testY, test.earth.probs) #模型选择 #。对于分类模型的比较，受试者工作特征（ROC） 图是一个很有效的工具。 #简言之，ROC基于分类器的性能对其进行可视化、组织和选择（Fawcett， 2006）。 #在ROC图中，Y轴是真阳性率（TPR），X轴是假阳性率（FPR）。 #TPR = 正确分类的阳性样本数/所有阳性样本数 #FPR = 错误分类的阴性样本数/所有阴性样本数 #AUC等于一个观测者在面对一对随机选出的案例（其中一个 是阳性案例，一个是阴性案例）时，能正确识别出阳性案例的概率 #用3行代码就可以画出ROC图。这个程序包还配有一个网站，上面有示例和演示，你可以通过以 下网址找到该网站： #http://rocr.bioinf.mpi-sb.mpg.de #展示4个不同的ROC图：全特征模型、使用BIC选择特征的简化模型、MARS模型和 一个糟糕模型。 # 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) #然后使用这个对象建立一个带有TPR和FPR的对象 perf.full <- performance(pred.full, "tpr", "fpr") #。在此之后，使用plot() 函数绘制ROC图 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 performance(pred.bic, "auc")@y.values performance(pred.bad, "auc")@y.values performance(pred.earth, "auc")@y.values #我认为，一个纯粹的统计学家会建议选择简约的模型，但其他人可能更倾向于全特征模 型。 #这就又归结到各种因素的权衡问题了，比如模型准确性与解释性、简约性与扩展性之间的权 衡。