Logistic回归 R语言

 

对于二分类结果变量为了便于计算，结果通常用0和1表示。对于体重或身高这样的连续型变量，总体或祥本的代表性指标是均值或中位数。而对于二分类资料，其代表性指标是结果变量中某种结果所占的比例。这些比例可以作为概率的估计值。
   概率虽然易于理解，但是在Logistic回归模型中，用优势的对数值即更方便。假设P代表患病的概率，1-P 就是不患病的概率，优势的定义为P/（1-P）。

1. x1-xm表示各危险因素、混杂因素或它们之间的交互项。(在这些自变量中，一部分是研究的目标变量，其他的变量代表潜在混杂因素，也称之为协变量）
2. p/(1-p)为发病与不发病之比，称为比值或比数(odds) 。
3. β0为常数项，表示所有自变量都不存在时正常人群中该病的基准发病率。
4. β1-βm为需要估计的各自变量的偏回归系数，反映危险因素、混杂因素及交互项的效应。

 

偏回归系数β表示自变量X每变化一个单位，所引起的OR值自然对数改变量。

代码

步骤一：导入数据，处理0-1变量为factor

install.packages("rJava")
library(rJava)
install.packages("xlsxjars")
library(xlsxjars)
install.packages("xlsx")
library(xlsx)
a<-read.xlsx("C:/Users/Administrator/Desktop/r.xlsx",1)


a$y<-factor(a$y,levels=c(0,1),labels=c("挂了","没挂"))
#将二分裂变量处理成factor

summary(a)
#查看各变量的基本统计信息

步骤二：查看数据

str(a)
#查看变量类型

步骤三：逻辑回归

model<-glm(y~.,data=a,family=binomial())
#将所有变量纳入模型中
#binomial表示二分类
summary(model)

#结果解读：

Deviance Residuals——表示偏差残差统计量。在理想情况下应服从正态分布，均值应为0。

Estimate——是回归系数和截距，β值（这里就是Estimate）

Std. Error——表示回归系数的标准误

z value——是统计量值（z的平方就是Wald值），Wald是一个卡方值，等于β除以它的标准误（这里是Std. Error），然后取平方（也就是z值的平方）。Wald用于对β值进行检验，考察β值是否等于0。若β值等于0，其对应的OR值，也就是Exp(β)为1，表明两组没有显著差异。OR等于β值的反自然对数。Wald值越大，β值越不可能等于0。

Pr(>|z|)——是P值

Null deviance——无效偏差，只包含常数项的离差平方和

Residual deviance——残差偏差，当前模型的离差平方和。

无效偏差和残差偏差之间的差异越大越好，用来评价模型与实际数据的吻合情况

AIC——赤池信息准则，表示模型拟合程度的好坏，AIC越低表示模型拟合越好。

 

步骤四：卡方分布检验

（下图来自百度）

pchisq(885.79-872.95,df=855-851,lower.tail = FALSE)
#“885.79-872.95”表示卡方分布的值，DF表示自由度，lower.tail表示是否计算小于“885.79-872.95”的概率

#得出p=0.01208463,小于0.05在95%置信区间，两个模型的差异具有统计学意义

步骤五：查看OR值和置信区间

exp(cbind("OR"=coef(model),confint(model)))
#计算OR值及其可信区间

步骤六：方差分析，筛选有统计学意义的变量

anova(model,test="Chisq")
#方差分析

#NULL表示只有截距

#每增加一个变量残差就会降低

#本例中a4的p值显著，对模型方差的改进程度=884.55-872.95

步骤七：基于筛选后的变量，向后法逐步回归

form<-as.formula(y~a3+a4)
model.new<-glm(form,data=a,family=binomial())
summary(model.new)$coefficients
#向后法逐步回归

coef(model.new)
#查看回归函数的系数

confint(model.new)
#查看参数置信区间

exp(coef(model.new))
#查看OR值

第二种方法——调用colgit函数

x<-read.xlsx("C:/Users/Administrator/Desktop/map.xlsx",1)
#分类结果按0-1读取，不转为factor
install.packages("survival")
library(survival)
clogitl<-clogit(y~a4,data=a)
summary(clogitl)

 绘制ROC曲线

data=a
data$predvalue<-predict(model,type="response")
install.packages("pROC")
library(pROC)
ROC<-roc(data$y,data$predvalue)
auc(ROC)#查看曲线下面积
ci(auc(ROC))#查看置信区间
plot(ROC,print.auc=T,auc.polygon=F,lwd=3,max.auc.polygon=F,print.thres=F)

显示ROC的详细信息

ROC.reuslts<-coords(ROC,"best",ret="all",transpose=FALSE)
as.matrix(ROC.reuslts)

#其中，threshold代表截断值，specificity代表灵敏度，sensitivity代表特异度