R笔记:二分类资料的logistic回归
转自个人微信公众号【Memo_Cleon】的统计学习笔记:R笔记:二分类资料的logistic回归。
写在前面:最近发现了一个非常“watch out”的事情,以前搜寻程序包的R Site Search用不了了。不过好在发现程序包pkgsearch,如果你想通过R实现某项统计分析但不知道使用哪个函数时,不妨试试用这个包里的函数pkg_search()或ps()来找找。本次笔记会涉及用AUC评估模型拟合效果,但我不知道绘制ROC曲线的R函数有哪些,用这个函数搜一下:
library(pkgsearch)
ps(query="ROC",size=10) ##寻找ROC相关程序包,显示10个
结果显示有92个相关程序包,其中pROC比较受欢迎,我们后面就使用了这个函数。
二分类资料其实还有很多其他分析模型,比如对数线性模型,但其应用远不及logistic回归普遍,除了对数线性模型相对复杂外,还有一个重要的原因logistic在对结果的解释上直接与OR相关联,简单易懂。
当控制其他变量,变量Xm增加一个单位(连续变量为增加一个单位,分类变量由取值水平发生改变)时,模型变为:
两个方程两侧相减,可得:
即OR=exp(βm),当控制其他变量,变量Xm增加一个单位时,因变量事件发生率大约是原来的OR倍左右(注意OR的解释)。
data(SAheart,package="bestglm") ##导入bestglm包中的SAheart数据
head(SAheart) ##查看数据结构
by(SAheart,SAheart$chd,summary) ##分组展示两个组(患冠心病与否)的描述性统计信息
lgrfit1<-glm(chd~age+typea+famhist+tobacco+alcohol+sbp+ldl+adiposity, family=binomial(link = logit), data=SAheart)
##二分类logistic回归。Family to use: binomial(link = "logit"),gaussian(link = "identity"),Gamma(link = "inverse"),inverse.gaussian(link = "1/mu^2"),poisson(link = "log"),quasi(link = "identity", variance = "constant"),quasibinomial(link = "logit"),quasipoisson(link = "log")
summary(lgrfit1)##展示模型拟合的细节信息。单独展示模型的参数可采用函数coefficients(lgrfit1)或coef(lgrfit1)
结果中首先给出了模型的表达式,其后是偏差残差的最小值、四分位数和最大值,接着是模型中各个自变量的系数、系数标准误及统计学检验结果,在logistic回归中我们更关注的是由系数转换成的OR,其解读见后面。本部分结果的后面还出现了两个偏差(Deviance):Null deviance和Residual deviance,Deviance就是-2lnL,Null deviance指的是空模型(只含截距的模型)的-2lnL,而Residual deviance指的是当前模型的-2lnL,当前模型的Deviance也可以通过函数deviance(lgrfit1)直接获得,其值越小模型拟合越好。两者的差值即为卡方值,通过卡方值(596.11-474.23)和自由度(461-453)即可获得P值,可以说明新引进的变量是否具有统计学意义,其实这就是所谓的似然比检验(Likelihood Ratio Test)。最后还给出了赤池信息准则(AIC),值越小代表模型拟合越好。
confint(lgrfit1) ##给出模型参数的置信区间。confint(object, parm, level = 0.95),parm可指定想显示置信区间的参数,level指具体定置信区间,不指定默认为95%。如confint(lgrfit1, “age”, level = 0.99)显示age的99%的置信区间
exp(cbind(OR=coef(lgrfit1),confint(lgrfit1))) ##给出模型参数的OR和置信区间
以年龄和家族史为例对OR值做解释:在固定其他因素之后,年龄每增加一岁,患者发生冠心病的风险大约是原来的1.05倍左右(注意OR以这种容易理解的形式解读,必须带上大约/左右此类的词);在固定其他因素之后,有家族史的患者发生冠心病的风险大约没有家族史患者的2.50倍左右。
residuals(lgrfit1) ##生成模型的残差。跟SPSS相比,此处的残差应该是偏差残差(Deviance residual)
predict(lgrfit1) ##函数predict(object, newdata = NULL, type = c("link", "response", "terms"))可实现用拟合模型对新数据集进行预测。newdata缺失,默认使用模型对建模数据进行预测【optionally, a data frame in which to look for variables with which to predict. If omitted, the fitted linear predictors are used】。tpye缺省默认生成线性模型的预测结果即logitP,设定type = "response"则生成预测概率,等同于fitted(lgrfit1)【for a default binomial model the default predictions are of log-odds (probabilities on logit scale) and type = "response" gives the predicted probabilities】
deviance(lgrfit1)##模型的偏差(即-2lnL)
plot(lgrfit1) ##生成评价拟合模型的诊断图
library(pROC)
roc(SAheart$chd,fitted(lgrfit1),plot=TRUE)
结果显示AUC=0.7942。一般来说0.5-0.7表示诊断价值较低,0.7-0.9表示诊断价值中等,0.9以上表示诊断价值较高。
Hosmer-Lemeshow检验:可以使用函数HosmerLemeshowTest {DescTools}、
hoslem.test{ResourceSelection}等。
library(ResourceSelection)
hoslem.test(SAheart$chd,fitted(lgrfit1),g=10)
Hosmer and Lemeshow goodness of fit (GOF) test
data: SAheart$chd, fitted(lgrfit1)
X-squared = 10.049, df = 8, p-value = 0.2616
结果显示当前模型与饱和模型没有统计学意义(P=0.2616),表示当前模型已经对数据信息充分提取。
lgrfit2<-glm(chd~age+typea+famhist+tobacco+ldl, family=binomial(link = logit), data=SAheart)
summary(lgrfit2)
新模型lgrfit2的每个回归系数都具有统计学意义(p<0.05)。
anova(lgrfit1,lgrfit2, test="Chisq")
Analysis of Deviance Table
Model 1: chd ~ age + typea + famhist + tobacco + alcohol + sbp + ldl +
adiposity
Model 2: chd ~ age + typea + famhist + tobacco + ldl
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 453 474.23
2 456 475.69 -3 -1.4523 0.6933
AIC(lgrfit1,lgrfit2)
df AIC
lgrfit1 9 492.2332
lgrfit2 6 487.6856
AIC也显示lgrfit2优于lgrfit1。
本例只有一个typea(A型行为)这一个变量为分类变量,二分类变量在统计软件中会按照0和1进行赋值处理,无论是按连续变量还是哑变量纳入结果都是一样的。哑变量设置参见《R笔记:哑变量》。在R语言中大部分函数会把字符变量和因子变量直接按哑变量处理,所以在遇到多分类变量时我们首先需要将其设为因子,因子设置参见《R笔记:因子》。
SAheart$sbpG<-ifelse (SAheart$sbp<120,0,SAheart$sbp)
SAheart$sbpG<-ifelse (SAheart$sbp>=120 & SAheart$sbp<140,1,SAheart$sbpG)
SAheart$sbpG<-ifelse (SAheart$sbp>=140 & SAheart$sbp<160,2,SAheart$sbpG)
SAheart$sbpG<-ifelse (SAheart$sbp>=160 & SAheart$sbp<180,3,SAheart$sbpG)
SAheart$sbpG<-ifelse (SAheart$sbp>=180,4,SAheart$sbpG)
library(rms);library(Hmisc);library(lattice);library(survival);library(Formula);library(ggplot2);library(SparseM)
lgrfit3<-lrm(chd~age+typea+famhist+tobacco+alcohol+sbpG+ldl+adiposity, data=SAheart)
lgrfit3
cbind(coef=coef(lgrfit3),OR=exp(coef(lgrfit3)))
lrm{rms}拟合的模型展示了更多的模型信息:(1)拟合模型的观测数量,本例总共462例,其中冠心病160例,无冠心病患者302例;(2)模型似然比检验结果:当前模型与空模型相比具有统计学意义(chi2=124.55,P<0.001);(3)区分度指标包括Nagelkerke R2(0.326)、C-指数(0.796,C指数就是AUC)等;(4)各自变量的参数估计结果。这个前面用glm拟合的结果解释是类似的,只是这里面有一个多分类哑变量的解释,哑变量在模型中要遵循同进同出,R中参照水平是赋值最低的水平。高值血压患者(sbpG=1)发生冠心病的风险是正常血压的患者0.74倍,或者说正常血压的患者发生冠心病的风险是高值血压患者(sbpG=1)的1/0.743=1.35倍,但没有统计学意义(P=0.4122),其他几个哑变量解释与此类似,不再累述。
转自个人微信公众号【Memo_Cleon】的统计学习笔记:R笔记:二分类资料的logistic回归。
END