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)

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,884评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,755评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 158,369评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,799评论 1 285
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,910评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,096评论 1 291
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,159评论 3 411
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,917评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,360评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,673评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,814评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,509评论 4 334
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,156评论 3 317
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,882评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,123评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,641评论 2 362
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,728评论 2 351

推荐阅读更多精彩内容