logistic回归系列

1、常用的总结性函数

summary()展示拟合模型细节
coefficients()
coef()列出拟合模型的参数（截距项和斜率）
confint()给出模型参数的置信区间（默认95%）
residuals()列出拟合模型的残差值
anova()生成两个拟合模型方差分析表
plot()生成评价拟合模型的诊断图
predict()用拟合模型对新数据集进行预测
deviance()拟合模型的偏差
df.residual()拟合模型的残差自由度
MASS包做判别分析
earth包做MARS分析，多元自适应回归样条
InformationValue包生成混淆矩阵，评估模型错误预测率
ROCR包话ROC曲线
ggDCA决策曲线

2、变量筛选方法：

1、全部变量全部纳入，不做筛选；
2、用逐步回归筛选自变量（step()函数）；此步没讲解，我也不懂
3、先做单因素Logistic回归，P<0.1纳入最后的回归方程；
4、本教程先后用全部变量建模，bestglm包筛选最优建模变量，以及多元自适应回归样条（MARS）三种方法建模，最后画ROC曲线评估。
本文利用bestglm函数。

3、数据练习,

3.1、加载需要的R包+读取数据

library(foreign)
library(nnet)
library(ggplot2)
library(reshape2)
library(car)
library(MASS)
library(ROCR)
library(InformationValue)
library(ggDCA)
library(earth)
rm(list = ls())
setwd("D:\\雷爷\\吻合口瘘数据")
mydata=read.table('logtttt.txt',sep = "\t",
                      comment.char = "#",stringsAsFactors = F,
                      header = T,fill = T,quote = "")

na.omit函数去掉NA值,
str()函数看数据结构,
并且把理论上是数值的变为数值用as.numeric()函数；或者分类归分类

mydata$distance=as.numeric(mydata$distance)
mydata$Times=as.numeric(mydata$Times)
mydata1=na.omit(mydata)

3.2、卡方检验，与fish精确检验

chisq.test(mydata1$leakage, y = mydata1$Radiotherapy, correct = TRUE,
           p = rep(1/length(x), length(x)), rescale.p = FALSE,
           simulate.p.value = FALSE, B = 2000)
fisher.test(mydata1$Radiotherapy, y = mydata1$leakage, workspace = 200000, hybrid = FALSE,
            hybridPars = c(expect = 5, percent = 80, Emin = 1),
            control = list(), or = 1, alternative = "two.sided",
            conf.int = TRUE, conf.level = 0.95,
            simulate.p.value = FALSE, B = 2000)

3.3、设置验证集以及训练集(70%训练集，30%为测试集)

set.seed(123)
ind=sample(2,nrow(mydata1),replace = TRUE,prob = c(0.7,0.3))
train=mydata1[ind==1,2:32]
test=mydata1[ind==2,2:32]

利用table()函数可看两个子集的构成，似乎不是每次都会均衡哦；

3.3.2其他方法

另一种设置训练集以及验证的方法是基于caret包的createDataPartition函数，以下根据患者ID分验证集训练集

library(caret)
datad=createDataPartition(y=mydata1$ID,p=0.70,list = F)
train=mydata1[datad,]
test=mydata1[-datad,]
table(train$leakage)
table(test$leakage)

3.4.1全部变量（特征）一起建模(.表示全选)，看看结果如何

logistic回归,全部一起分析

model_full=glm(leakage~.,data = train,family = binomial())
summary(model_full)

用exp()，confint()看95%置信区间

confint(model_full)
exp(coef(model_full))

用VIF值（方差影响因子？）检测每个特征之间的多重共线性
一般认为vif<5可排除共线性
方差膨胀系数(variance inflation factor，VIF)是衡量多元线性回归模型中复 (多重)共线性严重程度的一种度量。它表示回归系数估计量的方差与假设自变量间不线性相关时方差相比的比值。
VIF的取值大于1。VIF值越接近于1，多重共线性越轻，反之越重；
通常以10作为判断边界。当VIF<10,不存在多重共线性；当10<=VIF<100,存在较强的多重共线性；当VIF>=100, 存在严重多重共线性。

vif(model_full)

计算预测概率，predict()函数实现

train_p=predict(allmodel,type = "response")

评价模型在训练集上的执行效果，使用R包InformationValue，
先实现混淆矩阵

trainY=mydata1$lavage[ind==1]
testY=mydata1$leakage[ind==2]

以上得出结果，行表示预测值，列表示实际值，对角上的元素表述预测正确的
分类，右上角为误认为瘘的预测数，左下角是误认为没有瘘的预测数。
misClassError()查看错误预测的比例！

misClassError(trainY,train_p)

将近8%,好高啊！接下来在验证集上预测，计算

test_full_p=predict(model_full,type = "response",newdata = test)
confusionMatrix(testY,test_full_p)
misClassError(testY,test_full_p)

嗯。。。2%差强人意！

3.4.2、bestglm包筛选最优建模变量（BIC\CV等方法）

接下来，bestglm()包自带交叉验证方法来筛选最佳建模变量
IC = "CV"，可以改其他方法，例如最有子集BIC，不懂可以利用help函数查看帮助。

library(bestglm)
x=train[,1:30]
y=train[,31]
Xy=data.frame(cbind(x,y))
bestglm(Xy, family = gaussian, IC = "BIC", t = "default", 
        CVArgs = "default", qLevel = 0.99, TopModels = 5, 
        method = "exhaustive", intercept = TRUE, weights = NULL, 
        nvmax = "default", RequireFullEnumerationQ = FALSE)

接下来用挑选出来的变量建模

model_BIC=glm(leakage~Defecation+drainage_tube,
              family = binomial(),data = train)
summary(model_BIC)
exp(coef(model_BIC))

验证集验证预测准确率

test_BIC_p=predict(model_BIC,type = "response",newdata = test)
confusionMatrix(testY,test_BIC_p)
misClassError(testY,test_BIC_p)

画图，当评价模型适用性，绘制初始响应变量的预测值余残差的图形

plot(predict(model_BIC,type="response"),residuals(model_BIC,type="deviance"))

3.4.3、判别分析

3.4.3.1、线性判别分析
model_lda=lda(leakage~.,data = train)
model_lda
plot(model_lda,type="both")

两类结果的百分比，两类结果下的各个变量均值，线性判别回归系数，plot可看图，两类结果的得分情况有些变量被错误分类。

train_lda_p=predict(model_lda)$posterior[,2]
misClassError(trainY,train_lda_p)
confusionMatrix(trainY,train_lda_p)

看看看测试集表表现

test_lda_p=predict(model_lda,newdata = test)$posterior[,2]
misClassError(testY,test_lda_p)
confusionMatrix(testY,test_lda_p)

用二次判别分析（QDA）

model_qda=qda(leakage~.,data = train)
model_qda

嗯。。。。数据太少（瘘的发生率）。好吧！后续一样计算错误预测率

3。4.4、多元自适应回归样条（）_MARS模型——earth包

set.seed(1)
model_earth=earth(leakage~.,data=train,
                  pmethod="cv",
                  nfold=5,
                  ncross=3,
                  degree=1,
                  minspan=-1,
                  glm=list(family=binomial)
                  )
summary(model_earth)

画图

plotmo(model_earth)
plotd(model_earth)
plot(model_earth)

a

evimp(model_earth)

看看验证集的表现

test_earth_p=predict(model_earth,newdata = test,type = "response")
misClassError(testY,test_earth_p)
confusionMatrix(testY,test_earth_p)

3.5.ROC曲线-评价模型

ROC曲线是好的选择！不同模型画图即可。三个模型一起出图舒服！

pred.full=prediction(test_full_p,test$leakage)
perf.full=performance(pred.full,"tpr","fpr")
pred.BIC=prediction(test_BIC_p,test$leakage)
perf.BIC=performance(pred.BIC,"tpr","fpr")
pred.earth=prediction(test_earth_p,test$leakage)
perf.earth=performance(pred.earth,"tpr","fpr")
plot(perf.full,main="ROC",col=1)
plot(perf.BIC,col=2,add=TRUE)
plot(perf.earth,col=3,add=TRUE)
abline(0,1,col=5,lty=2)
legend(0.8,0.4,c("FULL","BIC","EARTH"),1:4)

计算AUC

performance(pred.full,"auc")@y.values
performance(pred.BIC,"auc")@y.values
performance(pred.earth,"auc")@y.values

3.6、策曲线分析（decision curve analysis, DCA）

DCA也可用于cox回归分析。
ggDCA包，用于logistic回归中,可以同时画多个模型(不包括earth选出来的模型)，
那就把m1，m2和m3都赋值给dca()即可。

dca_full=dca(model_BIC,model.names='BIC')
ggplot(dca_full)

3.6.2、验证集画DCA

当验证1个模型的时候，我们会讲原始模型(m3)和验证模型(validate)画在一张图里。
但是当验证的模型在2个及以上的时候，所画出的图形仅仅包含了验证结果，不包含原始模型结果。

dca_test=dca(model_BIC,new.data=test,model.names='BIC')
ggplot(dca_test)

最后的最后，简单的线性回归模型——lm()

4.1

linemodel=lm(age~BMI+CEA,data = mydata)
summary(linemodel)
©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成，浏览时请结合常识与多方信息审慎甄别。
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

相关阅读更多精彩内容

友情链接更多精彩内容