第3章
?Logistic回归分析
在
实际的临床研究中，有多种类型的变量，如连续型变量或分类变量。响应变量除可能有连续的取值外，还可能会是二分类变量，即只有两种状态，例如某诊断结果是“阳性”或“阴性”，某结局事件是“生存”或“死亡”，某药物治疗效果是“有效”或“无效”等。此外，响应变量还可能会是多分类变量，即有多种状态，此时要综合考虑响应变量与其他变量之间的关系时，多元线性回归已不再适用，应采用Logistic回归来解决这类问题。
3.1 Logistic回归分析的基本概念
在医学研究中，Logistic回归是分析疾病与致病因子间联系的重要统计方法。它是以疾病发生概率为因变量，影响疾病发生的因子为自变量的一种回归方法。医学研究中的因变量有时并不是呈正态分布的连续型随机变量，其取值可能只有两个，如发病与未发病、阳性与阴性、暴露与未暴露等。此时，线性回归不再适用，而Logistic回归模型成功地解决了这一问题。根据响应变量（因变量）的类型，Logistic回归可分为二分类响应变量的Logistic回归和多分类响应变量的Logistic回归。
3.2 Logistic回归的模型结构
在多元线性回归中，Y可以是任意的数值变量。若Y表示的是疾病发生的概率P，则可以把P作为因变量并建立与各自变量xi的回归方程。经过研究，如果把P转化为ln，则会使回归方程的统计性能更好。此变换被称为P的logit转换。即：
logit(P)=ln
Logistic回归的模型结构：设二分类因变量Y（Y=1或Y=0），令Y=1的概率为，则Y=0的概率为1?π。令，则以为因变量建立回归方程：
进一步可推导出概率预报模型：
实际上，Logistic回归模型可以看作是多元线性回归模型的推广。
Logistic回归系数通常采用最大似然估计法，此方法适用范围广，不需要自变量呈多元正态分布，但计算较为繁琐，不过有了计算机，已经变得很方便了。在Logistic回归模型中，通常采用wald 检验对参数估计值进行假设检验，可检验参数是否为0。目前大多数软件都采用这种检验方法。将第j个自变量对应的优势比定义为：。其中，优势比表示当其他自变量保持不变时，该自变量每增加一个单位，所引起的变化量。每个自变量对应的优势比ORj的95%的可信区间为：。当自变量为连续型变量时，OR值是指自变量每增加一个单位，其优势的变化量。当OR>1，说明该因素是危险因素；而当OR<1，说明该因素是保护因素。
3.3 应用实例1：一般资料的Logistic回归
当前一些研究表明糖尿病与促甲状腺激素（TSH）血清水平有一定的潜在关系。加强对糖尿病患者甲状腺功能指标的检测可以早期诊断防治糖尿患者中无症状的甲状腺功能 异常。
本例中将126名65岁以上糖尿病患者分为两组：TSH<4组编码为0，表示TSH水平正常；而TSH>4组编码为1，表示TSH水平异常。将年龄、病程、糖化血红蛋白、空腹血糖、右眼底病变（分为8个等级，其中0表示正常，1～7表示病变等级，等级越高表明病变越严重）作为自变量，采用Logistic回归分析TSH血清水平的影响因素。
这里采用SPSS软件进行分析。首先进行数据录入，如图3.1所示。
选择菜单Analyze→Regression→Binary Logistic，如图3.2所示。
在Logistic回归对话框中，因变量“Dependent”放入“tsh分组”，自变量“Covariates”放入“年龄、病程、血红蛋白、空腹血糖和右眼病变”。“Method”的主要方法如下。
（1）Enter：所有变量一次全部进入方程；
（2）Forward LR：逐步向前法，自变量根据似然比检验结果依次进入方程；
（3）Backward LR：后退法，自变量根据似然比检验结果依次移出方程。
图3.1 糖尿病与血清水平的数据录入界面（部分数据）
本例采用“Enter”方法来进行分析（图3.3）。
图3.2 Logistic回归分析菜单
图3.3 Logistic回归分析主对话框
单击【Options…】按钮，打开如图3.4所示的对话框，勾选“CI for exp（B）”，可以计算出OR值的95%置信区间。
主要输出结果如图3.5～图3.8所示。
如图3.5所示的输出表输出的是模型总的全局检验，即似然比卡方检验，3个结果分别为：Step统计量为每一步与前一步相比的似然比检验结果；Block统计量是指将block1与block0相比的似然比检验结果；而Model统计量则是上一个模型与现在方程中变量有变化后模型的似然比检验结果。
图3.4 Logistic回归Options子对话框
图3.5 Logistic回归分析结果（全局检验）
如果采用Enter法，得到的3个统计量及假设检验结果是一致的，即=11.804，P=0.038，说明模型有统计学意义。
如图3.6所示的输出表输出的“?2 Log likelihood”可用于拟合优度检验结果。而“Cox&Snell R Square”和“Nagelkerke R Square”类似于线性回归中的决定系数。
图3.6 Logistic回归分析结果（Model Summary）
如图3.7所示的输出表是Logistic回归分析中最重要的表。表中的B为回归系数，S.E.为标准误，Wald为Wald 值，df为自由度，Sig.为P值，Exp(B)为OR值，95.0% C.I.for EXP(B)为OR值的95%置信区间，其中Lower为置信区间的下限，Upper为置信区间的上限。该回归分析结果说明，只有年龄是TSH水平的影响因素（P=0.021），而其他因素的回归系数均无统计学意义（P>0.05）。年龄的OR值（OR=1.106，95%CI：1.015～1.206）可解释为年龄每增加一岁，TSH异常的风险增加到原来的1.106倍。Logistic回归方程可以表示如下：
图3.7 Logistic回归分析结果（Variables in the Equation）
?6.977+0.101*年龄
如果采用逐步向前法（Forward LR），最终的输出结果如图3.8所示。该结果与采用Enter方法获得的结果基本相似，年龄仍然是唯一有统计学意义的变量（P=0.016），年龄的OR值（OR=1.106，95%CI：1.019～1.202）。
图3.8 Logistic回归分析结果（Variables in the Equation）
3.4 应用实例2：列联表资料的Logistic回归
列联表资料的数据除了可以用卡方检验进行比较之外，同样也可以采用Logistic回归进行分析。来看下面的例子，如表3.1所示。
表3.1 COPD患者与非患者的吸烟情况资料
有吸烟史
无吸烟史
合计
患者
231
125
356
非患者
183
296
479
合计
414
421
835
下面建立Logistic回归模型方程分析病例与吸烟的关系。令Y＝1表示患者，Y＝0表示非患者；X＝1表示吸烟，X＝0表示不吸烟。由此，可建立Logistic回归方程为：
其中，π表示COPD发病率。下面用SPSS软件完成Logistic回归分析。首先进行数据录入：病例否变量中1表示COPD患者，0表示非患者；吸烟否变量中，1表示吸烟，0表示不吸烟（图3.9）。
然后对人数进行加权，选择Data→Weight Cases（图3.10）。
图3.9 列联表资料的Logistic回归数据录入界面
图3.10 列联表资料的加权菜单
弹出对话框后，对人数进行加权（图3.11）。
图3.11 列联表资料的加权对话框
对人数进行加权后，选择菜单Analyze→Regression→Binary Logistic，在弹出的对话框中选择因变量“Dependent”为“病例否”，协变量“Covariates”为“吸烟否”（图3.12）。
图3.12 列联表资料的Logistic回归主对话框
主要的输出结果如图3.13所示。
图3.13 Logistic回归分析结果（Variables in the Equation）
如图3.13所示的输出结果说明，吸烟是罹患COPD的影响因素（P<0.001，OR=2.989，95%CI：2.247～3.976）。如果写成方程，则可以写成：吸烟，说明吸烟对罹患COPD有影响，吸烟者患病的优势是不吸烟的近3倍。
3.5 应用实例3：多项Logistic回归分析
如果响应变量是多分类的，此时就需要采用多项Logistic回归分析。与二项Logistic回归不同，它是通过拟合广义Logit模型方法进行的。若响应变量有K个分类，则有一个参照类，其余每一类与参照类比较，拟合K-1个Logit模型。例如，因变量取3个值a、b、c，如果以a为参照水平，则得到2个Logistic函数，一个是b与a相比，另一个是c与a相比。下面通过一个案例进行分析。
这里选择R软的流行病学软件包epicalc中的数据Ectopic。先来看一下数据，在R窗口中输入语句：
install.packages(pkgs="epicalc")
（安装R软件流行病学epicalc软件包）
library(epicalc)
（加载R软件流行病学epicalc软件包）
data (Ectopic)
（加载Ectopic数据集）
Ectopic
（输出Ectopic数据）
输出的结果如图3.14所示。
图3.14 epicalc软件包中的数据Ectopic
这里仅列出了部分数据。该数据是检验先前的人工流产是否为当前宫外孕的一个危险因素的病例对照研究。总样本数为723人。研究的患者（变量outc）分为3组，分别为宫外孕患者（EP）、来做人工流产的孕妇（IA）和来分娩的孕妇（Deli）。变量hia表示患者的人工流产史：分为从未做过人工流产组（never IA）和曾经做过人工流产组（ever IA）。变量gravi表示患者的怀孕次数：分为怀孕1～2次（1～2），怀孕3～4次（3～4）和怀孕>4次（>4）。下面以outc变量作为分类响应变量，hia和gravi作为影响因素进行多项logistc回归分析。
由于SPSS软件做多项Logistic回归有一些不足，因此这里采用R软件的nnet软件包分析，简单又便于解释。在加载epicalc软件包时，会自动加载上nnet软件包。
由于在多项分类的Logistic回归中，响应变量都是高于2个水平。通常系统会以默认的第一个水平作为参照水平。但有时如果以第一个水平作为参照，可能不便于结果的解释，因此可以通过自行设定参照水平进行分析。
例如，想以Deli作为参照水平，可以在R窗口中输入如下语句：
outc<-Ectopic[,2]
（从数据集中选择分类响应变量outc）
EP<-outc=="EP"
（筛选出outc中的EP样本）
IA<-outc=="IA"
（筛选出outc中的IA样本）
Deli<-outc=="Deli"
（筛选出outc中的Deli样本）
multi<-multinom(cbind(Deli,EP,IA)~hia+gravi,data=Ectopic)
（进行多项logistic回归分析，其中cbind表示将outc的三个水平合并，列在第一个的Deli被作为参照水平，如果想让EP作为参照水平，只需要将EP列在第一位即可）
mlogit.display(multi)
（显示多项logistic回归分析结果）
输出的结果如图3.15所示。
图3.15 多项Logistic回归分析结果
上述结果解释为：把分娩孕妇（Deli）作为参照水平，对有人工流产史的妇女（ever IA），在此次入院中出现宫外孕的风险增加为没有人工流产史妇女（never IA）的4.44倍（P<0.001，具有高度显著性）。把分娩孕妇（Deli）作为参照水平，对有人工流产史的妇女（ever IA），在此次入院中再次进行人工流产的风险是没有人工流产史妇女（never IA）的1.47倍（P>0.1，无统计学显著性）。
对于怀孕次数来说，把分娩孕妇（Deli）作为参照水平，怀孕3～4次的妇女（gravi3-4），在此次入院中出现宫外孕的风险增加为怀孕1～2次妇女（gravi1-2）的1.59倍（P>0.1，无统计学显著性）。怀孕超过4次（gravi>4）的妇女，在此次入院中出现宫外孕的风险增加为怀孕1～2次妇女（gravi1-2）的2倍（P>0.1，无统计学显著性）。同样，对于怀孕次数来说，把分娩孕妇（Deli）作为参照水平，怀孕3～4次的妇女（gravi3-4），在此次入院中再次人工流产的风险增加为怀孕1～2次妇女（gravi1-2）的2.35倍（P<0.001，具有高度显著性）。怀孕超过4次的妇女（gravi>4），在此次入院中再次人工流产的风险增加为怀孕1～2次妇女（gravi1-2）的3.2倍（P=0.001，具有高度显著性）。
研究者在进行疾病的风险因素分析时，常常会受到混杂因素的干扰，如年龄、性别、病程长短和病情轻重等。如果不对混杂因素加以控制，会使研究结果产生偏倚。Logistic回归方法能够充分利用数据信息，有效地控制混杂因素，得到校正后风险因素的OR估计值和可信区间。
Logistic回归与多元线性回归的区别在于多元线性回归模型的响应变量值是具体数值，因此模型的预测值具有实际意义。而Logistic回归模型的响应变量是分类变量，因此模型的预测更关注于样本的分类标识。Logistic回归中的自变量既可以是连续变量，也可以是分类变量。但应用Logistic回归时应注意样本量的问题。一般样本量应是自变量个数的20～30倍。在样本量足够大，且自变量之间无相关性的情况下，可以将全部自变量进入回归方程中，并采用逐步回归方法筛选有统计学意义的变量。在样本含量不高的情况下，可以先采用单因素分析筛选出有统计学意义的变量，将这些变量进入到回归方程中再进行筛选。在作多项Logistic回归时，由于变量关系较为复杂，且涉及到参照组的选择，因此在作结果解释和下结论的时候要谨慎。
3.6 Logistic回归模型的Nomogram图展示
诺莫（Nomogram）图，又称为列线图，是一种综合分析多个定量变量和定性变量以预测某特定事件发生的绘图法预测模型，它可以用一种直观的绘图来对个体患者进行风险评估。该模型可以基于Logistic回归模型和Cox回归模型，将其结果进行可视化的呈现。它根据模型回归系数的大小来制定评分标准，给每个自变量的每种取值赋值一个评分，对每个患者，均可计算得到一个总分，再通过得分与结局发生概率之间的转换函数来计算每个患者的结局时间发生的概率，其轴结构和风险点反映了各个变量对预测结果的影响和重 要性。
目前在一些国家和地区，这种预测模型已经受到广大患者和临床医师的认可，患者和临床医师都被鼓励使用这种模型进行预后风险评估工作。下面通过一个案例讲解如何绘制Logistic回归模型的诺莫图。
首先安装R软件的rms软件包。在R窗口中输入语句：
install.packages("rms")
library(rms)
（安装R软件的rms软件包）
（加载R软件的rms软件包）
这里提供一个案例数据，其中gene1和gene2中的0表示基因低表达，1表示基因高表达；stage中的0表示肿瘤早的分期，1表示肿瘤晚的分期；acid表示某项临床指标；outcome中的0表示存活，1表示死亡。图3.16列出了数据集中的部分信息。数据保存为nuomo.csv文件，并存储于D盘中。
图3.16 绘制诺莫图的案例数据（部分数据）
然后读入数据并构建数据列表。
read.table("D:\\nuomo.csv",header=TRUE,sep=",")->data
ddist<-datadist(data)
options(datadist="ddist")
（读入数据）
（构建数据列表）
以outcome为因变量，以gene1、gene2、age和stage为自变量构建Logistic回归模型，在R窗口中输入如下语句：
model<-lrm(outcome~gene1+gene2+age+stage,data)
（构建Logistic回归模型）
输出结果如图3.17所示。
其中，在“Rank Discrim.Indexes”中的C为模型区分度评价的统计量，也就是AUC=0.813，说明模型的区分度较好。gene1(P=0.0036)和stage(P=0.0356)是outcome的影响因素。
下面构建诺莫图的绘图函数。在R窗口中输入如下语句：
nom<-nomogram(model,lp.at=seq(-2,4,by=0.5),
fun=function(x)1/(1+exp(-x)),
funlabel="Risk of Death",
fun.at=c(0.05,seq(0.1,0.9,by=0.1),0.95),
conf.int=c(0.1,0.3))
nom
（设置坐标轴的范围，从-2到4，步长为0.5）
（转化为风险概率）
（给新的坐标轴设置名称）
（给新的坐标轴设置范围）
（显示置信区间）
（输出nom函数结果）
图3.17 构建Logistic回归模型
nom函数的输出结果如图3.18所示。
图3.18 nom函数的输出结果
该输出结果中含有给每个变量的打分及总分（Total Points），并根据总分计算出相应的死亡风险概率。例如，如果gene1低表达就赋予0分，高表达则赋予100分。下面根据输出的nom函数结果进行绘图。在R窗口中输入语句：
plot(nom,lplabel="linear Predictor",
fun.side=c(3,1,1,1,3,1,3,1,1,1,3),
label.every=2,
col.conf=c("red","green"),
conf.space=c(0.1,0.3),
col.grid=gray(c(0.8,0.95)))
（绘制线性预测坐标轴）
（设置风险概率坐标轴刻度；1，3表示刻度出现在轴的上下方）
（设置变量之间的间隔）
（设置置信区间颜色）
（显示置信区间）
（设置图片沿格线的灰度颜色）
绘制的诺莫图如图3.19所示。
图3.19 Logistic回归模型的诺莫图
下面对诺莫图进行解释。在图3.19中，患者的取值都位于每个变量轴上，向上绘制一条直线以确定每个变量值对应的分数（Points）。所有变量获得的分数总和为Total Points。将Total Points向下垂直延伸到Risk of Death概率轴，就可以获得患者发生死亡的风险概率了。例如，某个患者gene1=1（高表达，Points=100），gene2=1（高表达，Points=12），age=64（Points≈10），stage=0（Points=0），Total Points=122；在Total Points轴上找到此数值并向Risk of Death概率轴作垂线，则可知该患者死亡风险概率为0.5～0.6。
为了使诺莫图更加清晰，也可以将置信区间去掉。此时将构建的诺莫图的绘图函数和绘图语句可以修改如下：
nom<-nomogram(model,
lp.at=seq(-2,4,by=0.5),
fun=function(x)1/(1+exp(-x)),
funlabel="Risk of Death",
fun.at=c(0.05,seq(0.1,0.9,by=0.1),0.95))
plot(nom,lplabel="linear Predictor",
fun.side=c(3,1,1,1,3,1,3,1,1,1,3),
label.every=2,
col.grid=gray(c(0.8,0.95)))
修改后的诺莫图如图3.20所示，此时置信区间去掉，图片显得更为清晰。
图3.20 Logistic回归模型的诺莫图（去除置信区间）
3.7 多个Logistic回归模型评价的决策曲线分析法
通常情况下，评价一种诊断方法是否好用一般是做ROC（Receiver Operating Characteristic Curve）曲线并计算AUC（Area Under ROC Curve）。但是，ROC只是从该方法的特异性和敏感性考虑，追求的是准确率。当通过某个生物标志物预测患者是否患了某种疾病，无论选取哪个值为临界值，都会遇到假阳性和假阴性的可能，有时避免假阳性时受益更大，而有时则希望能够避免假阴性。既然两种情况都无法避免，那就需要找到一个净收益最大的办法，而这个问题就是临床效用问题。对此，Andrew Vickers博士等研究出另外一种评价方法，也就是决策曲线分析（Decision Curve Analysis，DCA）法。
DCA分析涉及3个基本概念：阈值概率（Threshold Probability）、净收益（Net Benefit）和权重因子（Weighting Factor）。
（1）阈值概率为患者选择治疗的诊断确定性水平。阈值概率考虑患者治疗与否的相对价值，如果治疗效果高且风险低，则选择治疗的阈值概率低；如果治疗效果极低或有很大的风险，则选择治疗的阈值概率将很高。
（2）净收益是指治疗决策所带来的预期收益和预期伤害之和。
（3）权重因子衡量患者不接受治疗（或治疗不足）和过度治疗所带来的风险。
这三者之间的关系是：
例如，发表在2018年Lancet Haematology杂志上的一篇论文就采用了决策曲线分析来评价临床预测模型。该论文用两个队列构建了肿瘤相关的静脉血栓栓塞的预测模型。其中，一个队列的模型评价采用了决策曲线分析，如图3.21所示。
图3.21 决策曲线评价临床预测模型（图片来自Lancet Haematology期刊）
图3.21中的横坐标为阈值概率，纵坐标为净增益。Treat all和Treat none表示两种极端情况。Treat none表示不治疗，所有样本净获益均为0；Treat all表示全面治疗，其净获益是斜率为负值的反斜线。从图中可以看到，相比不治疗或全面治疗，该预测模型有更好的临床应用价值。
下面应用R软件来绘制决策曲线分析图。首先需要安装和加载R软件的rmda软件包。这里采用该软件包自带的数据dcaData，数据结构如图3.22所示。该数据中含有年龄Age（连续变量）、性别Female（二分类变量：1为女性，0为男性）、Smokes（二分类变量：FALSE为不吸烟，TRUE为吸烟）、Marker1和Marker2（两个肿瘤标记物，连续变量）、Cancer（二
图3.22 决策曲线分析的数据结构
分类变量：0为非癌患者，1为癌症患者）。在R窗口中输入语句：
install.packages("rmda")
library(rmda)
data(dcaData)
（安装R软件的rmda软件包）
（加载R软件的rmda软件包）
（加载软件包自带数据dcaData）
下面构建两个Logistic回归模型，其中一个是以癌症（Cancer）为因变量，以年龄（Age）、性别（Female）和吸烟（Smokes）为自变量的模型model1；另一个是以癌症（Cancer）为因变量，以年龄（Age）、性别（Female）、吸烟（Smokes）以及肿瘤标志物1（Marker1）和肿瘤标志物2（Marker2）为自变量的模型model2，在R窗口中输入语句：
model1=Cancer~Age+Female+Smokes
model2=Cancer~Age+Female+Smokes+Marker1+Marker2
（构建模型model1）
（构建模型model2）
下面采用decision_curve()函数分别以公式model1和model2构建基于Logistic回归的DCA模型，在R窗口中输入语句计算model1的决策曲线：
baseline_model<-decision_curve(formula = model1,
data=dcaData, family=binomial(link='logit'),
thresholds=seq(0,1,by=0.01), confidence.intervals=0.95,
study.design = "cohort")
（model1作为基线模型）
（构建Logistic回归模型）
（以0.01为步长取阈值概率）
（计算95%置信区间）
（研究设计为队列研究）
在该语句中family=binomial(link='logit')表示使用Logistic模型进行拟合；threshold设置横坐标阈值概率范围，一般是0～1，但是如果遇到某种具体情况，在临床上一致认为阈值概率达到某个值以上，如40%，则必须采取干预措施，那么0.4以后的研究就没什么意义了，此时可以设阈值概率为0～0.4；study.design为研究类型，根据队列研究和病例对照研究，此值分别可设为“cohort”或“case-control”。
同样，在R窗口中输入语句计算model2的决策曲线：
full_model<-decision_curve(formula=model2,
data=dcaData,
family=binomial(link='logit'),
thresholds=seq(0,1,by=0.01),
confidence.intervals=0.95,
study.design = "cohort")
（model2作为完全模型）
（构建Logistic回归模型）
（以0.01为步长取阈值概率）
（计算95%置信区间）
（研究设计为队列研究）
下面用plot_decision_curve()函数绘制决策曲线，在R窗口中输入语句：
plot_decision_curve(list(baseline_model,full_model),
curve.names=c("Baseline model","Full model"),
col=c("blue","red"),lty=c(1,2),lwd=c(3,2,2,1),
legend.position="topright",
confidence.intervals=FALSE,
cost.benefit.axis=FALSE)
（绘制两个模型的决策曲线）
（设置曲线名称）
（设置不同模型曲线的颜色，样式和宽度）
（设置图注的位置）
（确定是否绘制曲线的置信区间）
（确定是否绘制成本效益坐标轴）
绘制出的决策曲线图如图3.23所示。
图3.23 两个Logistic回归模型的决策曲线分析
图3.23中横坐标为阈值概率，纵坐标为标准净增益。All和None表示两种极端情况：None表示所有样本都净获益为0；All表示所有样本都有净获益，该净获益是斜率为负值的反斜线。从图3.23可以看到，full_model模型的净收益比baseline_model模型高。这说明两个肿瘤标志物对癌症预测起了一定的作用。
26
医学数据挖掘案例与实践（第2版）
27
第3章 Logistic回归分析