Play with R 第8期：逻辑回归
目录
8.1 What will this chapter tell me?
8.2 Background to logistic regression
8.3 What are the principles behind logistic regression?
8.4 假设和可能出错的事情
8.5 本章使用的R包
8.6二元逻辑回归分析
8.7 如何报告Logistic回归
8.8 假设检验：另一个例子
8.9多分类逻辑回归（Multinomiallogistic regression）
8.1 What will this chapter tell me?
前面的章节介绍了回归分析，适用于结果变量为连续变量的情况，当结果变量为分类变量时，采用本章所介绍的逻辑回归。
8.2 Background to logistic regression
逻辑回归是预测变量为连续或分类、结果变量为分类变量的多重回归。结果变量为两类的为两分类逻辑回归，多于两类则为多分类逻辑回归。
8.3 What are the principles behind logistic regression?
逻辑回归的方程与一般的回归方程存在很多相似性，当仅有一个预测变量X1时，公式如下：
当有多个预测变量时，公式如下:
对于线性回归，变量关系也是线性的，但是对于逻辑回归需进行对数转换以达到线性表达。
8.3.1 Assessing the model: the log-likelihood statistic
通过观察数据对模型进行估计的原理：
8.3.2 Assessing the model: the deviance statistic
偏差（deviance，也称为-2LL）与对数似然密切相关：
8.3.3 Assessing the model: R and R2
在线性回归中，多重相关系数R和相应的R2是测量数据拟合情况的指标。
以下为不同研究者提出的R2补充:
8.3.4 Assessing the model: information criteria
对于线性模型，可用AIC和BIC作为拟合评价指标。这两个指标的存在用于解决R2的问题：每当模型增加一个变量，R2增加。
8.3.5 Assessing the contribution of predictors: the Z-statistic
在线性模型中，我们不仅想知道模型的整体拟合情况，还有每个预测变量的贡献。线性回归中，计算的是回归系数(b)及其标准为来计算t值。逻辑回归中也有一个类似的统计量——服从正态分布的z值，与t值相似，z说明了预测变量的系数b显著不等于0。
8.3.6 The odds ratio
逻辑回归中的一个重要值是OR值，为B的指数，是预测变量每变化一个单位导致概率变化的指标。
8.3.7 Methods of logistic regression
8.3.7.1 强行进入法
默认的回归方法是将回归模型的预测变量都放入block中，估计每个预测变量的参数。
8.3.7.2 逐步法
可选择向前或向后回归法。向前回归中，模型第一步仅包括一个函数，后续逐步加变量。向后回归则相反，第一步包含全部变量，后续逐步移除变量。
8.3.7.3 如何选择方法？
回归方法的选择最主要考虑的是你是检验一个理论或仅进行探索性分析。有研究者认为逐步法不适用与理论检验，但是对于不考虑因果关系的情况是合适的。对于一般的回归，由于抑制效应，逐步回归法中的向后比向前更优。
8.4 假设和可能出错的事情
8.4.1假设（Assumptions）
1、逻辑回归的一些假设：
（1）线性（Linearity）：一般的线性回归能对连续值结果进行预测；而在逻辑回归中结果变量是分类的，故采用数据的对数（发生概率除以没有发生概率再取对数），使得因变量和自变量之间呈线性关系。因此逻辑回归的线性假设是：任何连续的预测对象和结果变量的对数之间存在线性关系。该假设可以通过观察预测对象和它的对数变换之间的交互项是否显著来检验
（2）独立性误差（Independence of errors）:各观测对象间相互独立，不要相互关联。
（3）多重共线性（Multicollinearity）：预测对象之间的相关性不应该太高。因为逻辑回归对模型中自变量多重共线性较为敏感，例如两个高度相关自变量同时放入模型，可能导致较弱的一个自变量回归符号不符合预期，符号被扭转。需要利用因子分析或者变量聚类分析等手段来选择代表性的自变量，以减少候选变量之间的相关性。
2、逻辑回归的自身特问题：R给出一个完全不正确的结果。这一现象可以通过大的标准误差来体现。主要有两种情况引发该情况：不完全信息和完全分离收敛失败：有时会在输出中得到一条错误信息（如下），意味着我们要忽略R的任何输出，我们的数据是无效的。
8.4.2 预测的不完全信息
1、收集所有变量组合的数据，以保证预测结果的正确
2、在使用交叉表进行分析时需要：收集所有变量组合的数据、检查表中每一个单元格的期望频率，确保它们大于1，且不超过20%的期望频率小于5（逻辑回归中的拟合度测试做出了这样的假设）
8.4.3 完全分离（Completeseparation）
1、逻辑回归从本质上说属于二分类问题，如图8.3结果变量y分类为：1 =Burglar, 0 = Teenager。逻辑回归使用sigmoid函数将预测值映射为(0, 1)上的概率值，帮助判断结果。
2、完全分离：结果变量可以被一个预测变量或多个预测变量组合很好地预测。同时会产生恒大的标准误差。
3、图8.3中的上下三角形数据点有重叠，用R来预测结果概率时，在低权重下，拟合的概率遵循图的底线，在高权重下，它遵循顶线。在中间值时，它试图跟随变化的概率。
图8.4中的上下三角形数据点没有重叠，实现了完美分离，结果可以通过weight完美地进行预测。但是中间部分会由于数据的缺乏而不能确定曲线的斜率导致产生很大的标准误差。该问题经常出现在当太多的变量去适合太少的情况，唯一的解决方案是收集更多的数据（有时也可以使用更简单的模型找到一个简洁的答案）。
8.5 本章使用的R包
在这一章中需要包car（用于重新编码变量和测试多重共线性）和mlogit （用于多项逻辑回归）。如果没有安装这些包，你将需要安装和加载他们。因此，本节所需安装包：
然后通过执行命令来加载这些包
8.6二元逻辑回归分析
例子：测试在肛门塞入鳗鱼是否可以治愈便秘。
结果（因变量）：Cured（Cured/NotCured）；
预测指标1（自变量）：Intervention（Intervention/NoTreatment）；
预测指标2（自变量）：Duration（患者在接受治疗前出现问题的天数，如1、2、3…）。
1. 准备数据
1.1.载入数据——将工作目录设置在数据文档所在位置，并且执行以下代码：eelData<-read.delim("eel.dat",header = TRUE)
由此创建了一个名为eelData的数据集
1.2. 用head(eelData)查看数据集的前六行，如下所示：
1.3. R按文本首字母顺序为变量Cured和Intervention分别设置了两个水平。也就是说，变量Cured的基线类别是Cured而不是NotCured，变量Intervention的基线类别是Intervention而不是NoTreatment。考虑到我们关注的是治愈效果和干预方案的作用性，基线类别应当分别是NotCured和NoTreatment。为了重新设置水平，需要执行以下代码：
如此，现在的基线类别便是NotCured和NoTreatment了。
2. 基本逻辑回归分析
2.1.为了操作逻辑回归分析，我们需要用到glm()的功能，即广义线性模型。它的一般形式是：
newModel<-glm(outcome~predictor(s),data= dataFrame,family = name of a distribution, na.action = an action)
其中newModel包含了模型信息，我们可以通过执行summary(newModel)获得有关该模型的数据。Outcome即因变量，在本例中即Cured。Predictor(s)指的是那些用来预测因变量的预测指标们，在本例中即Intervention和Duration。DataFrame是因变量和自变量来源的数据集。Family指的是分布名称，比如本例中的二项式分布binomial。Na.action主要用于当数据存在缺失时，在本例中，我们可以忽略它。
在本例中，我们实施层次回归：在模型1中，我们只加入Intervention作为预测指标，然后在模型2中我们加上Duration。因此，针对模型1，执行的代码如下：
该模型中的Cured仅由Intervention预测而来。
针对模型2，我们使用以下代码：
该模型的Cured由Intervention和Duration共同预测而来。
2.2. 解释基本逻辑回归分析的结果
为了查看模型，我们需要执行以下代码：
summary(eelModel.1)
summary(eelModel.2)
模型整体拟合性由偏差估算，越大数值的偏差代表着模型拟合越差。R提供两个偏差数据：thenulldeviance 指的是仅包括截距，不包括其他预测指标的模型的偏差，即-2LL(基准)。Theresidualdeviance 指的是包括截距，也包括其他预测指标的模型的偏差，即-2LL(新)。
2.2.1. 当模型中仅加入Intervention
模型1输出的结果如下所示：
模型整体拟合性：由结果输出可知，thenulldeviance数值是154.08，theresidualdeviance 数值是144.16，也就是说，当Intervention被纳入模型之后，偏差变小了，说明加入了Intervention的模型更好地预测了因变量。至于加入了Intervention的模型有多好地提高了模型的拟合性，可以由themodelchi-square计算得到。具体计算方法：154.08-144.16=9.92。根据查阅chi-square的数值分布可知，该数值在.05的水平上是显著的，因此我们可以说，加入Intervention显著地提高了模型对因变量的预测准确性。当然我们也可以通过R直接计算出themodelchi-square 和它的显著性，操作如下：
最后结果即为c2=9.93,p=.002，与之前人工计算的结果一致。因此得出结论，加入了Intervention的新模型可以显著地预测是否被治愈。
模型参数：我们不仅希望知道模型的整体拟合性，也希望知道Intervention作为一个预测指标对该模型的影响，因此我们使用z-statistic来查看它的参数b是否显著地不同于0。从上面的输出结果可以看到，b=1.23，z=3.07，p<.002，因此Intervention是一个显著的预测治愈效果的因子。
R和R2：我们可以利用模型1的结果输出计算得到R，计算公式如下：
因此 R2 = 0.222 = 0.05
而使用R程序，我们可以执行三种不同的R2的计算：
Hosmer and Lemeshow’s measure
Coxand Snell’s measure
Negelkerke’smeasure
虽然最终算出来的R2各不相同，但是我们可以使用它们测量模型的effectsize。
Odds Ratios (OR)：为了知道预测指标的一个单元变化所引起的odds变化，我们必须要先计算病人没有接受干预方案时候的odds，然后计算他们接受干预方案了的odds，最后计算这两个odds的比例变化，如下：
a. 病人没有接受干预方案的odds
b. 病人接受干预方案的odds
c. 两个odds的比例变化
用R程序运算oddsratio结果如下：
那么，该结果说明接受了干预方案的病人治愈的可能性是不接受干预方案病人的3.42倍高。
我们也可以用R程序计算得到该模型中oddsratios的置信区间。执行以下编码可得到结果：
我们需要注意，该结果中的置信区间上限和下限都是大于1的，这个结果让我们可以相信当预测指标Intervention上升的时候，治愈效果的odds也是上升的。
2.2.2. 当模型中同时加入Intervention和Duration
模型2输出的结果如下所示：
将模型1和模型2作比较，我们可以发现两个模型的偏差(residualdeviance)是一致的，这说明模型2并没有在模型1的基础上对模型有所改进。而且AIC在模型2中比模型1中高，这说明模型1是更好的模型。
模型整体拟合性：themodelchi-square计算的方法和之前模型1中提到的类似，当然要将原本使用的空模型的数据用现在要进行计算的模型数据代替，执行编码如下：
P=0.964,数值显著高于0.05，因此我们得出结论，模型2并没有对模型进行显著的提升。
还有另一种更容易地方法也可以帮我们计算到两个模型间的差异，使用anova()功能，如下所示：
可以看到，通过这种方法，我们也得到了两个模型的差异是0.0019835
2.3.逻辑回归中的个案诊断
我们需要检查建立的模型与实际数据的契合度。在逻辑回归分析中，拟合值指的是个案取预测指标的值时，Y发生的预测概率。在R中预测概率可以通过fitted()功能得到。在具体操作数据分析时，我们可以把个案诊断的变量作为新变量加入数据集中，以实现后面的各种计算。以下是代码：
接下来，我们可以通过head()功能查看数据中一些预测概率的例子。
其中，可以看到当病人不接受任何干预时，有0.429的可能性可以被治愈。而当病人接受干预时，有0.719的可能性可以被治愈。可能性为0说明没有可能被治愈，可能性为1说明绝对会被治愈，因此0.719这个数值说明了病人接受干预方案是确实会提高治愈机会的。
最后，我们需要用残差来检测模型是否是个好模型。其中，我们尤其要看thestandardized residuals, studentized residuals, leverage statistics, DFBeta。执行以下代码可以得到它们：
可知，所有的DFBeta都小于1，leveragestatistics接近于0.018，studentizedresiduals都在-2和+2之间。这些残差数据表明模型和实际数据契合得很好，数据中没有对模型产生影响的异常值。
8.7 如何报告Logistic回归
Logistic回归和线性回归的报告模式相同，最精简的方式是报告beta值，标准误，显著性水平以及模型其它的参数（如R方和拟合优度）。强烈建议报告odds ratio及其置信区间。如果报告了常数项，那么读者可以从中得知完整的回归模型。不显著的变量也可以报告，它们和显著的变量一样珍贵。
8.8 假设检验：另一个例子
考虑罚点球是否能中，存在两个影响因素：第一个因素是球员是否担忧（通过PSWQ问卷测量），第二个因素是球员过去的点球成功率（有记录）。本研究随机选取了75名球员完成PSWQ问卷，同时调取了他们从前的踢球数据记录。据此，本研究主要有四个变量：
1.Scored：结果变量，0=点球未中，1=点球击中。
2.PSWQ：这个变量是第一个预测变量，它告诉我们球员担心的程度。
3.Previous：这个变量是特定球员在职业生涯中罚球的百分比。因此，它代表了以前在点球得分上的成功。
4.Anxious：这个变量是我们的第三个预测变量，它以前没有被用来预测惩罚点球的成功。它是在接受处罚之前的一种焦虑状态的量度。
附：我们在前面的线性回归和本章的逻辑回归中学习了如何进行层次回归。尝试对这些数据进行层次逻辑回归分析。在第一个块中输入Previous和PSWQ，在第二个块中输入Anxious。
8.8.1多重共线性检验
首先，将数据读入一个新的dataframe，我们将其称为penalty-Data，可以通过将工作目录设置为文件的位置（参见3.4.4节）并执行：
在第7.7.2.4节中，我们了解了多重共线性如何影响回归模型的标准误差参数。Logistic回归同样容易产生共线性的偏置效应，因此在进行Logistic回归分析后进行共线性检验是十分必要的。我们在逻辑回归中寻找共线性的方法与在线性回归中寻找共线性的方法完全相同。
首先，我们使用来自自助任务的所有三个预测器重新创建模型，通过执行以下步骤来创建模型：
该命令创建一个模型(penaltyModel.2)，其中变量得分由PSWQ、Anxious和Previous(Scored~ Previous + PSWQ +Anxious)预测。创建了这个模型之后，我们可以获得VIF和容差，将模型名称从car包中输入到VIF()函数中。执行:
第一行给出VIF值，第二行给出容差(即VIF的倒数)。
结果显示在Output 8.6中。这些值表明存在共线性问题：一个大于10的VIF值通常被认为是有问题的。这个分析的结果是非常明确的：状态焦虑和以前的处罚经验之间存在共线性关系，这种依赖性导致模型产生偏差。
对于共线性的解决方案之一是忽略其中一个变量(例如，我们可能坚持使用block 1中忽略状态焦虑的模型)。这样做的问题应该很明显：无法知道应该忽略哪个变量。由此得出的理论结论是没有意义的，因为从统计学上讲，任何共线变量都可以省略。没有任何统计依据可以忽略一个变量而忽略另一个变量。即使移除一个预测因子，另一个同样重要的、不具有如此强的多重共线性的预测因子也可以替代它。当多重共线性中涉及多个预测因子时，另一种可能性是对这些预测因子进行因子分析，并使用得到的因子得分作为预测因子。最安全(尽管不令人满意)的补救方法是承认模型的不可靠性。因此，如果我们要报告哪些因素可以预测点球成功，我们可能会承认，在第一个模型中，以往的经验可以显著预测点球成功，但这种经验可能会通过增加状态焦虑来影响点球。这种说法具有很强的推测性，因为“Anxious”和“Previous”之间的相关性并没有告诉我们因果关系的方向，但它承认了这两个预测因素之间无法解释的联系。
8.8.2 检验logit的线性度
在本例中，我们有三个连续变量，因此我们必须检查每个变量是否与结果变量Scored的log线性相关。我在本章前面提到，为了测试这个假设，我们需要运行Logistic回归，但是要包含预测因子，即每个预测因子与自身log之间的相互作用。我们需要使用log()函数创建每个变量及其log的交互项。首先，我们来看PSWQ变量，我们称PSWQ和它的log，也即logPSWQInt的相互作用，通过执行下述命令创建该变量：
该命令在penaltyData dataframe中创建了一个名为logPSWQInt的新变量，该变量是变量PSWQ (penaltyData$PSWQ)乘以该变量的log(log(penaltyData$PSWQ))。这个Dataframe现在是这样的：
注意，有三个新的变量反映了每个预测与该预测的log之间的相互作用。
为了验证这个假设，我们需要像以前一样重新进行分析，只是我们应该把所有的变量放在一个单独的块中（我们不需要分层处理），我们还需要在每个预测及其log中加入三个新的交互项。我们通过执行以下步骤来创建模型：
这个命令创建了一个模型(penaltytest1)，其中得分的变量由PSWQ、Anxious、Previous和我们创建的这些变量与它们的log(logPSWQInt、logAnxInt和logPrevInt)进行预测。然后我们使用summary()函数来显示模型。
Output 8.7显示了假设检验的输出部分。我们只关心交互项是否重要。任何有意义的相互作用都表明主效应违反了对数的线性假设。三种交互作用的显著性值（列Pr(>|z|)中的值）均大于.05，表明PSWQ、Anxious和Previous均满足logit的线性假设。
8.9 多分类逻辑回归（Multinomiallogistic regression）
之前介绍的Logistic回归中，因变量都是二分的，即非此即彼。其实，当因变量是多分类时（如：因变量为教育水平，有高中低三个水平），也可以利用多分类Logistic回归解决。
本质上，多分类Logistic回归还是在进行和Logistic回归相似的两两比较，只不过是讲多分类间的比较转变成了一系列的两两比较（breaks the outcome variable down into a series of comparisonsbetween two categories)。仍以教育水平为例，如果有高、中、低三个水平，多分类Logistic回归就是做了低vs中、低vs高，两次比较。当然这里的比较也可以是高vs中、低vs中等等，如何进行比较全在于你如何设定基线组（baseline）。在R中要利用的包是mlogit。
我们用一个有趣的例子来具体解释这个过程，数据内容见下图：例如某个研究想要研究什么因素能够预测搭讪成功与否。因变量有三个水平，搭讪失败、要到电话、抱得美人归；而自变量有搭讪对象的性别（男、女）、聊天时的性格好坏（0-10），聊天时的有趣程度（0-10），聊天时的“车速”（是否有性爱相关内容）（0-10）。显然，因变量是一个多分类变量，适合利用多分类Logistic回归来解决。
首先，我们需要进行数据整形。需要注意的是，在做多分类Logistic回归时，和通常一个被试一行的数据形式不同。根据因变量水平，一位被试所占的行数应该等于因变量的水平个数。我们需要将因变量的各个水平列出，并在被试真实所在的因变量类别输入TRUE，其余各类别则为FALSE。例如：被试1最终要到了搭讪对象的电话号码，有趣程度为3分，性格好坏为6分，聊天“车速”为7分，搭讪对象为男性，那么他的数据应该形如下图的前三行。后三行为被试2的数据，可以看到被试2抱得美人归，他的有趣程度为5分，性格好坏为2分，聊天“车速”为7分，搭讪对象为男性。当然，在R语言中使用mlogit.data()函数就可以完成这个操作。
newDataframe<-mlogit.data(oldDataFrame,choice = "outcome variable", shape ="wide"/"long")
在本例中，代码为：
mlChat<- mlogit.data(chatData, choice = "Success", shape ="wide")
然后，使用mlogit()函数建立回归方程（具体的设立规则与之前介绍的lm()和glm()
函数类似），特别注意新增的一个参数为reflevel，即设置那一组为基线组，在本例中设置搭讪失败组（No response/Walk off）是比较合理的，所以我们将它设为第三组（reflevel）= 3。完成代码如下：
newModel<-mlogit(outcome ~ predictor(s),data = dataFrame, na.action = an action, reflevel = a number representing thebaseline category for the outcome)
在本例中：
chatModel <- mlogit(Success ~ 1 | Good_Mate+ Funny + Gender + Sex + Gender:Sex + Funny:Gender, data = mlChat, reflevel =3)
最后我们来一起解读这个例子，用summary(chatModel)即可查看回归结果。如下图，我们可以看到R输出了一系列的系数B（SE）和总体模型的相关参数。可以看到在这个例子中，R2 = 0.138，卡方值为278.52，且达到了显著水平。
便于我们后续报告odds ratio，我们可以利用exp()函数输出：
exp(chatModel$coefficients)
data.frame(exp(chatModel$coefficients))
想要输出回归系数的置信区间也可以利用exp()函数：
exp(confint(chatModel))
以性格好坏对搭讪结果的影响为例。见输出结果的第三行和第四行，可以看到性格的好坏可以显著预测搭讪失败或是要到电话，但却对搭讪失败还是抱得美人归的预测作用就不显著。具体的报告格式如下：
本期分享人：
张骊凡（8.1-8.3）
王丹（8.4-8.5）
吴映秋（8.6）
杨雨濛（8.7-8.8）
许岳培（8.9）
排版: 林善焱
REVIEW
往期内容