题目
image.png
背景
OLS:线性回归
案例
image.png
##logistic回归与判别分析
library(MASS)
data(biopsy)#加载数据包
#View(biopsy)
str(biopsy)#查看数据格式
biopsy$ID = NULL##删掉ID项
names(biopsy) = c("thick", "u.size", "u.shape", "adhsn",
"s.size", "nucl", "chrom", "n.nuc", "mit", "class")#重命名，帮助文档有
names(biopsy)
biopsy.v2 <- na.omit(biopsy)#去除缺失数据，保存为另一数据，这样就未动原数据
y <- ifelse(biopsy.v2$class == "malignant", 1, 0)
#判别良恶性01代替；可以加载到数据框中，作为一列
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)
#开始创建训练集和验证集7：3
set.seed(123) #random number generator#设定种子
ind <- sample(2, nrow(biopsy.v2), replace = TRUE, prob = c(0.7, 0.3))
View(biopsy.v2)
#抽样函数，分成两个水平1和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)
#逻辑回归
full.fit <- glm(class ~ ., family = binomial, data = train)#class后面的点，表示接下来所有
summary(full.fit)#总结
confint(full.fit)#置信区间
exp(coef(full.fit))#OR值
#优势比可以解释为特征中1个单位的变化导致结果发生比的变化，如果系数大于1，说明当特征值增加时，结果的特征比会增加，反之。
#在进行数据探索时发现多重共线性的问题，同线性回归一样，在logistic回归当中，也可以计算出VIF统计量
library(car)
vif(full.fit)
#如果VIF的值不大于5的，说明没有共线性问题；下面可以进行特征选择
#先看看模型在训练集和验证集中表现如何
train.probs <- predict(full.fit, type = "response")
#计算每个观测的预测概率，即，良恶性的概率，对象就是我们构建回归方程的对象
#predict中不写response,其默认为线性预测值(有正有负),这样可以计算每个观测的预测概率
train.probs[1:5] #inspect the first 5 predicted probabilities#前5行预测概率
contrasts(train$class)#设置参照
###评价模型在训练集上的效果，之后，评价其在测试集上的拟合程度，快速实现评价的方法是生存一个混淆矩阵
#我们使用的混淆矩阵由caret包实现的，InformationValue也可以实现混淆矩阵；这时我们用0，1表示结果；
#而函数区分良性和恶性的默认结果是0.5，即当概率大于0.5，被认为是恶性的。
#install.packages("InformationValue")
library(InformationValue)
trainY <- y[ind==1]
testY <- y[ind==2]
confusionMatrix(trainY, train.probs)
#矩阵的行表示预测值，列表示实际值；对角线上的元素是预测正确的分类，7表示误为恶性的预测值；
#8表示误为良性的预测值；可以查看错误预测的比例，如下所示
# optimalCutoff(trainY, train.probs)
misClassError(trainY, train.probs)
#表示在训练集上的错误率，必须正确预测测试集，在测试集上建立建立混淆矩阵的方法和训练集上一样。
confusionMatrix(trainY, train.probs)
###测试集的表现
test.probs <- predict(full.fit, newdata = test, type = "response")#test表示测试集
misClassError(testY, test.probs)
confusionMatrix(testY, test.probs)
####
####虽然本次预测准确率高，接下来看看还有没有改进的空间；
##交叉验证的目的，就是提高测试集上的预测准确率，以及避免过拟合；
##K折交叉验证的做法是将数据集拆分成k个子集，每个等分称为一个k子集（k-set）
#算法每次留一个子集，使用k-1个子集拟合模型，然后再留出的那个子集上做预测；
#将上面K次验证的结果进行平均，使得误差最小化，并且获得合适的特征选择
#也可以使用留一交叉验证方法，这里的k等于N，模拟表明，LOOCV可以获得近乎无偏的估计，但方差很很高；
#所以大多数机器学习建议k的值为5或10
#bestglm包可以自动进行交叉验证，该包的使用依赖于我们在线性回归中使用过的leaps
#交叉验证的语法和数据格式存在注意事项，所以我们按部就班的就行就可以。
#install.packages("bestglm")
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是指已经格式化的数据框；CV是指交叉验证；CVArgs是指交叉认证的参数；
#HTF就是K折交叉验证法；10就是拆分的份数；
#REP = 1表示随机使用等分，并且只迭代一次；
#像glm()函数一样，需要指定参数family=binomial，
#如果指定参数family=gaussian,就可以使用bestglm就行线性回归；
#完成上面的交叉验证之后，会得到下面的结果：best model有三个特征；thick,u.size和nucl;
#我们可以把特征放到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)
#精简了特征的模型和全特征模型相比，精确度略有下降，我们可以使用bestglmb奥再测试一遍，使用信息准则为BIC最优子集方法。
bestglm(Xy = Xy, IC = "BIC", family = binomial)
#试过所有的可能的子集之后，以上4个特征提供了最小的BIC评分，该模型在测试集的效果如何？下面
#BIC越小，表示模型越优；
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)
#我们得到5个错误的预测，和全特征模型一样。所以，哪一个模型更好？在任何情况下，如果拥有相同的泛化效果；
#经验法则会告诉我们选择最简单的或解释性最好的模型；
#下面讨论判别分析的方法，在最后的模型推荐中，我们有可能使用这种方法建立模型
#总结，验证集验证讲了两种方法，K折验证和BIC验证，均为最优子集法。
判别分析
image.png
image.png
image.png
image.png
#4）判别分析
####z载入数据
library(MASS)
data(biopsy)#加载数据包
#View(biopsy)
str(biopsy)#查看数据格式
####
lda.fit <- lda(class ~ ., data = train)
lda.fit
#在分组先验概率中，良性大约为64%，恶性大约为36%。
#还有分类特征的均值；线性判别系数是标准线性组合，用来确定观测的判别评分特征；评分越高，越可能是恶性
#对LDA模型使用PLOT()函数，可以画出判别分数和直方图和密度图。
plot(lda.fit, type="both")
#图中可以看到，部分有重叠，表明观测被错误分类；
#LDA模型可用predict()函数得到3种元素；
#class(对良恶性的预测)/posterior(值为x的评分可能属于某个类别的概率)/x（为判别评分）的列表；
#通过下面的函数，我们仅提取恶性观测的概率
train.lda.probs <- predict(lda.fit)$posterior[, 2]#
misClassError(trainY, train.lda.probs)#错分率
confusionMatrix(trainY, train.lda.probs)#建立混淆矩阵
#很不幸，我们的LDA模型在在训练集上表现的比logistic回归差。
#LDA模型在测试集上的表现如下：
test.lda.probs <- predict(lda.fit, newdata = test)$posterior[, 2]
misClassError(testY, test.lda.probs)#错分率
confusionMatrix(testY, test.lda.probs)#混淆矩阵
#LDA模型在训练集上的表现差，在测试集的表现好。
#从正确分类率的角度看，其总体表现(96%)不如logistic回归模型(98%)好。
#下面用二次判别分析(QDA)模型来拟合数据，在R中，QDA也是MASS的一部分，函数为QDA()
#我们将模型存储在一个名为qda.fit的对象中，如下所示：
####
qda.fit <- qda(class ~ ., data =train)#仍然是训练集建模
qda.fit
#结果中有分组均值，和LDA一样，但是没有系数，因为这是二次函数，我们前面讨论过了。
#使用与LDA相同的代码，可以得到QDA模型在训练集和测试集上的结果。
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个预测错误，而误为恶性的比例尤高
#
#5）多元自适应回归样条法
#同样也是一种分类方法。
#如果以下有一种分类方法有以下特点，会不会喜出望外
#1）对于回归和分类问题，都可以灵活选择线性模型和非线性模型；
#2）简单易懂，易于解释；
#3）几乎不需要数据预处理；
#4) 可以处理所有类型的数据，数值型，因子型等
#5）可以很好的预测为止数据，在数据偏差方差权衡方面做的很好。
## MARS###########################################
#先根据一些数据建立一个标准线性回归模型；它的响应变量内部编码为0和1；
#在对特征进行编码之后，生存最终模型，生成一个GLM模型，所以不用考虑R方的值。
install.packages("earth")
library(earth)#多元自适性函数
set.seed(1)#设定种子，有交叉验证，结果可重复
earth.fit <- earth(class ~ ., data = train,
pmethod = "cv",#交叉验证
nfold = 5,#5折
ncross = 3,#重复三次
degree = 1,
minspan = -1,
glm=list(family=binomial)
)
summary(earth.fit)
#模型有8项，包括7个预测变量，其中两个变量有铰链函数，这就是浓度和染色质变量。
plotmo(earth.fit)#展示了保持其它变量不变，
#某个预测变量发生改变时，可以清楚的看到铰链函数对浓度所起的作用
plotd(earth.fit)
#可以生成按类别标签分类的概率密度图
evimp(earth.fit)
#变量之间的相对重要性，nsubset是精简过程完成过后包含这个变量的模型的个数；
#对于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)
#这个完全可以和logistic回归模型相比较，下面介绍各个模型，看看哪个模型时最好的选择。
#6）模型选择
#我们从模型中计算出的混淆矩阵和错分率就是给模型选择提供一个依据；
#对于分类模型的比较，受试者工作曲线ROC曲线是一个有效的工具。
#ROC基于分类器的性能对其进行可视化、组织和选择。
#在ROC曲线中，y轴是真阳性率(TPR)，x轴是假阳性率(FPR)。
#TPR:真阳/所有阳性样本
#FPR：假阳/所有阳性样本
##################################################
#下面将展示不同的ROC图:全特征模型；基于BIC选择特征的简化模型；MARS模型和一个糟糕迷行。
#该糟糕模型只有一个预测特征，会和其它模型形成对比。
#下面加载ROCR包，使用thick特征建立一个糟糕模型，并命名为bad.fit
library(ROCR)
#首先建立模型，该模型只纳入一个变量thick，该变量对整个模型的贡献度。
#在evimp(earth.fit)结果中看到。
bad.fit <- glm(class ~ thick, family = binomial, data = train)
#接着，验证该模型在验证集中的每个观测的概率
test.bad.probs <- predict(bad.fit, newdata = test, type = "response") #save probabilities
#下面就使用测试集数据绘制ROC图，每个模型3行代码。
#首先建立一个对象，保存对实际分类的预测概率。
#然后使用该对象建立一个带FPR和TPR的对象
#之后，使用plot()函数绘制ROC图。
#从全特征模型开始
#生成不同cut.off值下面的敏感度和特异度
pred.full <- prediction(test.probs, test$class)
#建立带TPR和FPR的 performance对象
perf.full <- performance(pred.full, "tpr", "fpr")
#下面用tpr和fpr进行plot画图,题目为ROC，字体颜色为黑色
plot(perf.full, main = "ROC", col = 1)
abline(0,1,col=5,lty=2)
#BIC建模
pred.bic <- prediction(test.bic.probs, test$class)
perf.bic <- performance(pred.bic, "tpr", "fpr")
plot(perf.full, main = "ROC", col = 1)
plot(perf.bic, col = 2, add = TRUE)#add = TRUE,加入现有图中
abline(0,1,col=5,lty=2)
#最后加入MARS模型和表现糟糕的模型
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)#perf.earth多元自适应回归样条法
legend(0.6, 0.6, c("FULL", "BIC", "BAD", "EARTH"), 1:4)
#可以看到全特征模型、BIC模型和MARS基本重叠一起。糟糕模型表现很差。
#接下来，计算AUC曲线：通过ROCR包建立performance对象来实现，将tpr和fpr换成auc
#代码如下
performance(pred.full, "auc")@y.values
performance(pred.bic, "auc")@y.values
performance(pred.bad, "auc")@y.values
performance(pred.earth, "auc")@y.values
#结果 pred.full为0.9972672 pred.bic为0.9944293 pred.bad为0.8962056 pred.earth为0.9967416
#我们发现，各个模型的预测能力没有明显的差别，
#解决方案就是我们可以重新将训练集合验证集重新随机化分组，将各种建模方法重新做一遍；
#统计学加往往选择最简约模型，其他人选择全特征模型，这就涉及到模型准确性和解释性的问题，简约型和拓展性的问题。
#仅通过GLM模型和判别分析不可能每次都有很好的预测效果，其它章节有更好的技术。
##end
#总结，本章研究了如何使用基于概率的线性模型预测定性的响应变量
#介绍了三种方法：logistic回归、判别分析和MARS。
