【原】logistic回归中logitP与自变量的线性考察

Memo_Cleon 2023-03-28 发布于上海

展开全文

Logistic回归模型：

Logistic回归属于广义线性模型的范畴，线性仍然是Logistic回归的适用条件之一。只是这里的线性不是因变量P与自变量X间的线性，而是因变量的logit变换值（logitP）与自变量X间的线性。理解了模型的公式，线性的考察似乎就很简单了。但如果按照这个思路进行考察时，你会发现根本行不通，因为就某个个体而言，实际观察结果要么是0，要么是1，即P要么0%，要么100%，实测值的logitP通过ln(P/(1-P))根本就无法计算，从而无从获得logitP残差。通过模型残差来实现线性考察的思路就行不通了。这个问题笔者在《R笔记：二分类资料的logistic回归》一文就有提出。

对于分类自变量而言，一般不需要考察其与模型因变量（logistic回归中是logitP）的线性，因为分类自变量与模型因变量是一定成线性的。分类变量常以哑变量的形式纳入模型，哑变量其实就是将具有K个水平的分类变量转变成K-1个二分类的自变量，两点之间必成直线。很多logistic回归模型中的解释变量全部为分类变量，这或许也是logistic回归常常不在意线性考察的一个原因。所以我们常常说的线性模型中的“线性”要求，其实主要是指连续性自变量与模型因变量之间的线性。

在临床研究当中，年龄和BMI是最为常见的两个人口学特征，我们常常默认它们对模型因变量的效应成线性。但事实往往并非如此。比如相对于正常的BMI，偏小和偏大可能更容易出现不良结局；同样的，在低年龄阶段，年龄与某种疾病可能就没有相关关系。此次笔记采用Framingham心脏研究数据来考察连续性自变量年龄与logitP的线性。

Logistic回归中连续性自变量与logitP的线性考察可以考虑以下几种思路：

既然前面说了通过“模型残差”考察行不通，我们不妨换一种思路。第一种思路是将连续变量处理成有序多分类，按照样本量的大小使类别尽量多，每个类尽量等距，每一类中都有阳性结局和阴性结局，这样我们就可以获得每个类别实际的logitP啦。不过像线性回归中的方差齐性探察一样，这只是一种不得已的变通方法而已。另外，这种方法适合大体考察一个自变量时，当存在多个自变量时不一定很牢靠，因为没有校正其他因素的影响，这些我们在《线性回归中的线性考察》也有介绍。将等级赋值当做连续变量对logitP作图就可以大体判断是否线性。自变量直接使用等级赋值就需要每个类尽量等距。等级赋值虽然可用于有界值的情况，但这样的情况下数据两端可能就是不等距，不等距可能对结果产生误导。本例最小年龄是32对，最高是70岁，所以第一个年龄段应该是32-35，并不跟其他类别等距，可以考虑去掉第一个点后作图。此外，组中值（midpoint）对logitP作图可能更合适一些。

第二种方法：分组法。将连续变量进行分类处理，通常是按四分位数或者五分位数，然后以分类变量形式建立logistic模型，通过模型参数可估计每个分类的组中值对应的logitP值。组中值与logitP简单作图即可判断，本例线性看起来似乎不错。肉眼判断通常可以保证准确性，但毕竟带有一定的主观性，如要确切估计可考虑第三种思路。

第三种方法：加入高次项等非线性项。一个模型加入连续性自变量的高次项（非线性处理）后，与原来的线性模式的模型相比较，如果模型的χ²有了显著提高，说明该自变量在模型中采用了更为合适的尺度，也就意味着该自变量与logitP之间为非线性关系，如果仍按线性关系来处理则意味着模型省略了自变量与logitP之间的非线性关系，可能会导致有偏的系数估计结果。这有点尝试的味道，但是不难理解。结合方法1和方法2的图形可大体判断曲线的形状，必要可以借助excel。

当然评估指标笔者觉得还可以是其他的一些拟合效果与拟合优度的指标，比如通过基于似然比模型检验统计量卡方值、-2LL（使用相同数据的嵌套模型比较）、伪决定系数、信息准则、预测正确率、ROC曲线等。

模型1嵌套于模型2，卡方值与AIC显示模型2更优；但所有模型以AIC准则来看，模型3最优，当然这种差别可能不大。模型1与模型2可以通过似然比检验：Chi2≈3396.610-3391.563=5.047，df=3-2=1，在excel的一个单元格中输入“=CHISQ.DIST.RT(5.047,1)”可获得P=0.025，表明用二次项的加入对模型的改善是有统计学意义的。从系数上来看，二次项Age²的系数也确实具有统计学意义（Wald chi2=4.917，P=0.027）。

第四种方法：借助Box-Tidwell法估计的线性转换参数来判断。Box-Tidwell估计的这个参数是用来对自变量进行幂次转换，以保证转换后的自变量与模型因变量成线性关系。变量x的线性转换参数为a，则表示变量x经a次幂的变换（x^a）后，x^a与因变量y成线性关系【如果a=0，则对自变量x进行自然对数变换：lnx；如果a≠0，则对自变量x进行a次幂的变换：x^a】。如果转换参数=1（统计学检验P>0.05）,即自变量x不需要任何变换就与模型因变量成线性关系。

Box-Tidwell转换参数a的计算是一种迭代估计：（1）y=β₀+β₁x；（2）y=β₀’+β₁’x+γxlnx；（3）a=γ/β₁+1。

将a作为变换参数对x进行变换。变化后的x作为新的自变量重复1-3步，直至得到较为稳定的a值作为自变量x最终的变换参数。很多情况下第一次的迭代就可以获得令人满意的a的估计值。至于为啥是xlnx，涉及用数学中泰勒级数求导啥的笔者就讲不清楚了，应用时其实也不用太在意这些深层次的理论基础。Box-Tidwell法估计时要求y是连续性因变量，对于二分类logistic回归，需要先计算出logitP作为因变量，R语言中的car包提供了boxTidwell函数。

目前借助Box-Tidwell法来检验线性在网络上普遍存在两种方式。一是根据a的计算过程在logistic回归模型中构建xlnx，然后检验其系数γ是否等于0（γ=0则a=1）。这样做的问题在于每次只能计算一次迭代的结果，而且网络上的教程基本上就是利用了首次的计算结果，还好大多数情况下，首次估计的a值往往就是比较理想的结果。第二种方式是利用logistic回归模型获得每个个案的预测概率P，然后求得每个个案的logitP=ln(p/(1-P))，最后再利用Box-Tidwell法来检验logitP与自变量x间的线性关系。这种思路笔者觉得并不合适，因为在建立logistic回归模型的时候，已经默认了logitP与自变量x间的线性关系，在这样的假定下获得的预测logitP理论上肯定与自变量x成线性的。

SPSS操作如下：

Transform>>Compute Variable…

Analyze>>Regression>>Binary Logistic…

从第一次迭代的结果来看，γ=0的检验结果是P=0.022，即Age与logitP成非线性关系。当然你也可以继续进行新一轮的迭代。

以上示例均是单变量的考察。正如我们在《线性回归中的线性考察》提到的：简单的散点图或单因素的线性回归只能大体地判断自变量与因变量是否成线性关系，因为这种方法并没有校正其他因素的影响。在多重线性回归中，自变量间可能相互影响，我们只需要把其他因素放入模型就可以啦，除了第一种方法，后面三种均可实现，此不赘述。

如果发现连续性自变量与logitP是非线性关系，该如何处理？这就涉及到连续性解释变量的编码问题了，分类处理、多项式、数据转换、样条都是常用的处理手法。