【原】R笔记：二分类资料的logistic回归

Memo_Cleon 2021-07-31

展开全文

转自个人微信公众号【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回归我们曾在<<因变量二分类资料的logistic回归>>中做出过详细的操作介绍及结果解读。

对结局为二分类的资料而言，因变量是事件发生率P，我们当然希望建立这样一个模型：P=β₀+β₁X₁+…β_mX_m，但是由于P取值在[0,1]，而β₀+β₁X₁+…β_mX_m的取值却在[-∞,+∞]。Logit变换解决了这个问题，P经过logit变换之后取值范围就变成了[-∞,+∞]，

二分类logistic回归自然要求因变量是二分类，比如患病与否、是否有效。从上面模型可以看出，Logistic回归仍然属于线性回归的范畴，但在构建方程时，模型因变量是对原来的因变量P进行了变换，变成了logit(P)或者说是ln(P/(1-P))，而不再是原来的P。

既然属于线性回归，其适用条件依然是独立、线性、正态、方差齐、自变量间不存在多重共线。只是需要注意的是模型方程中的因变量是logit(P)。线性模型的适应条件考察我们在以往的大量文中都有涉及：《线性回归中的线性考察》、《线性回归中的正态分布》、《线性回归中的方差齐性探察》、《R笔记：多重线性回归适用条件考察》、《多重线性回归分析》、《偏回归图与偏残差图》。但很少有logistic回归模型适应条件考察和模型的教程。大概是因为这些评估大部分需要用到“模型残差”，注意这里的“模型残差”应该是因变量为logitP的模型残差，而不是我们关注的P的残差。残差是实际观测值与模型拟合值之差，虽然我们可以通过模型获得logitP的预测值，但对某个个案而言，实际观察结果要么是0，要么是1，即P要么0%，要么100%，实测值的logitP=ln(P/(1-P))就无法计算，从而logitP残差无法获得。当然如果自变量是分类变量，倒是可以通过每个单元格中的事件发生与否来确定每个单元格的logitP的实际观测值，但这样的计算也无法针对每个个案，如果自变量是连续变量就更无能为力了。虽然在SPSS、R等软件中可以计算未标准化残差、logit残差、学生化残差、标准化残差和偏差残差（Deviance）等各种残差，但很遗憾这些都是针对因变量P的残差，而非参与建模的logitP的残差。虽然软件中可能也提供了适用条件/模型诊断的一些分析，但涉及残差的方法如果直接用P的残差来代替logitP的残差进行个人觉得并不妥当。因为二分类结局有0和1两个值，可以预见的是诊断图中会有两条分开的线。

二分类资料其实还有很多其他分析模型，比如对数线性模型，但其应用远不及logistic回归普遍，除了对数线性模型相对复杂外，还有一个重要的原因logistic在对结果的解释上直接与OR相关联，简单易懂。

当控制其他变量，变量Xm增加一个单位（连续变量为增加一个单位，分类变量由取值水平发生改变）时，模型变为：

两个方程两侧相减，可得：

即OR=exp(β_m)，当控制其他变量，变量Xm增加一个单位时，因变量事件发生率大约是原来的OR倍左右（注意OR的解释）。

二分类logistic回归R示例

示例数据来源：R语言bestglm包中的数据，数据来自一项南非西开普省心脏病高危地区男性的回顾性样本。响应变量chd表示冠心病发生与否，解释变量有sbp（收缩压）、tobacco（累积烟草量(kg)）、ldl（低密度脂蛋白胆固醇）、adiposity（肥胖）、famhist（家族史）、typea（A型行为）、obesity（肥胖）、alcohol（当前饮酒量）、age（发病年龄）

数据中只有famhist是二分类变量，其余变量均为连续变量。为更好的演示实际中可能遇到的情况，后面我们将sbp转化为有序多分类变量后利用不同的函数重新进行二次分析。

此次笔记建立二分类logistic回归模型采用函数glm {stats}和lrm{rms}。软件是R4.03。

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) ##生成评价拟合模型的诊断图

模型建立后，模型拟合效果如何，或者说模型区分结局的预测能力怎么样（按模型预测能否区分结局）？当前数据的信息是否被充分的提取？这就涉及到两个检验：模型拟合效果检验和拟合优度检验。

模型拟合效果代表的是模型区分试验结局的能力，现有的模型能否将因变量的各个分类区分开来，常用的指标有-2倍的对数似然值（-2lnL，即偏差(Deviance)，其值越小越接近0说明模型拟合越好）以及通过似然值计算出的伪决定系数Cox & Snell R²和Negelkerke R²，模型预测的正确率（对因变量结局预测的准确程度），AUC（ROC曲线下面积）。

模型拟合优度代表的是当前模型是否已经充分提取了数据信息，或者说当前模型是否还有进一步改善的空间，一般通过检验当前模型与饱和模型的预测效果是否具有统计学意义来实现，所谓饱和模型指的是纳入了所有自变量的主效应和各介交互作用的模型，饱和模型的预测值就是样本实测值。常用模型拟合优度检验方法有Pearson拟合优度检验和Deviance拟合优度检验、似然比检验(Likelihood Ratio Test)以及Hosmer-Lemeshow检验。Pearson拟合优度检验和Deviance拟合优度检验都是基于当前模型预测值与样本实测值间的比较，似然比检验是考察当前模型与饱和模型的-2lnL值之差是否具有统计学意义，Hosmer-Lemeshow检验是根据模型的预测概率从小到大排列将观察单位十等份，根据每一组中的实测值与理论值进行Pearson卡方分析，如无统计学意义则认为拟合模型与饱和模型（实际模型）基本没有偏差。

ROC曲线函数我们在本文一开始介绍了，此次笔记使用pROC包的roc函数来绘制ROC曲线和AUC。

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），表示当前模型已经对数据信息充分提取。

似然比检验用于嵌套模型间的比较，比如比较当前模型和饱和模型可以使用函数lrtest{rms}，还有一个函数anova()生成两个拟合模型的方差分析表，也可以比较两个嵌套模型的拟合优度。

在上面这个模型中，饮酒量、SBP和肥胖是没有统计学意义的，去除这些变量重新拟合模型（只是为了模型比较，直接这样删除方式有待商榷）：

lgrfit2<-glm(chd~age+typea+famhist+tobacco+ldl, family=binomial(link = logit), data=SAheart)

summary(lgrfit2)

新模型lgrfit2的每个回归系数都具有统计学意义（p<0.05）。

模型lgrfit2嵌套于模型lgrfit1中（模型lgrfit2是模型lgrfit1的一个子集），可以使用anova()函数来比较两个模型是否具有统计学意义。

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

结果的两个模型没有统计学意义（p=0.6933），表明用当前5个预测变量的模型lgrfIt2与8个预测变量的模型lgrfit1拟合程度相同。我们可以采用更简单的模型进行解释，OR转换及解释同上，不再累述。

模型间的比较还可是使用似然比检验和信息准则来进行，似然比检验前面已经介绍。这里简单说下AIC（赤池信息准则），信息准则无需考虑是否嵌套模型。信息准则考虑了模型的统计拟合度以及拟合的参数数量，较小信息准则的模型说明模型用较少的参数获得了足够的拟合度。

AIC(lgrfit1,lgrfit2)

df AIC

lgrfit1 9 492.2332

lgrfit2 6 487.6856

AIC也显示lgrfit2优于lgrfit1。

当然变量的筛选单纯从统计学上应考虑逐步回归（向前、向后、向前向后法）或者全子集回归，这些不是本次笔记的内容，全子集回归我们在《R笔记：全子集回归 | 最优子集筛选》中已有过介绍，逐步回归以后单独成文吧。

本例只有一个typea（A型行为）这一个变量为分类变量，二分类变量在统计软件中会按照0和1进行赋值处理，无论是按连续变量还是哑变量纳入结果都是一样的。哑变量设置参见《R笔记：哑变量》。在R语言中大部分函数会把字符变量和因子变量直接按哑变量处理，所以在遇到多分类变量时我们首先需要将其设为因子，因子设置参见《R笔记：因子》。

为演示多分类变量的纳入模型的形式，根据临床分级[sbp按高血压诊断标准分为0=正常、1=正常高值（120-140mmHg），2=1级高血压（140-160mmHg）、3=2级高血压（160-180mmHg）、4=3级高血压（≥180mmHg）]将sbp转换为多分类变量sbpG，对应的是高血压分级。

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)

SAheart$sbpG<-factor(SAheart$sbpG) ##将sbpG设成为因子，其排序方式刚好与逻辑顺序相一致，不再使用levels属性进行制定，如想设置可参见《R笔记：因子》

其后使用前面的函数glm {stats}进行拟合就可以了。接下来我们换用另外一个函数lrm{rms}进行拟合。

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 R²（0.326）、C-指数（0.796，C指数就是AUC）等；（4）各自变量的参数估计结果。这个前面用glm拟合的结果解释是类似的，只是这里面有一个多分类哑变量的解释，哑变量在模型中要遵循同进同出，R中参照水平是赋值最低的水平。高值血压患者(sbpG=1)发生冠心病的风险是正常血压的患者0.74倍，或者说正常血压的患者发生冠心病的风险是高值血压患者(sbpG=1)的1/0.743=1.35倍，但没有统计学意义（P=0.4122），其他几个哑变量解释与此类似，不再累述。