为什么完美的预测器会导致p值为0.93？
深入研究完美的预测器
如果您一直在关注我的先前博客文章，您可能会回想起逻辑回归在尝试拟合完全分离的数据时会遇到问题，导致无限的几率比。在Cox回归中，其中风险替换了几率，您可能会想知道是否存在类似的问题与完美的预测器。它确实发生了，但与逻辑回归不同，它在这里发生的方式要不明显得多，甚至什么构成“完美的预测器”也不明确。如后面将更清楚地表明的那样，完美的预测器被定义为其等级与事件时间的等级完全匹配的预测器x（它们的斯皮尔曼相关性为1）。
以前，在“拆箱Cox”中：
拆箱Cox：Cox回归的直观指南
风险和最大似然估计如何预测事件排名？
towardsdatascience.com
我们解释了最大似然估计，并介绍了一个由五个受试者组成的虚构数据集，其中单个预测器x表示延长寿命药物的剂量。为了使x成为事件时间的完美预测器，在此我们交换了受试者C和D的事件时间：
import numpy as npimport pandas as pdimport plotnine as p9from cox.plots import ( plot_subject_event_times, animate_subject_event_times_and_mark_at_risk, plot_cost_vs_beta,)perfect_df = pd.DataFrame({ 'subject': ['A', 'B', 'C', 'D', 'E'], 'time': [1, 3, 4, 5, 6], 'event': [1, 1, 1, 1, 0], 'x': [-1.7, -0.4, 0.0, 0.9, 1.2],})plot_subject_event_times(perfect_df, color_map='x')
为了了解为什么这些“完美的预测器”可能会有问题，让我们从上次离开的地方继续，并检查负对数似然成本针对β的绘图：
negloglik_sweep_betas_perfect_df = neg_log_likelihood_all_subjects_sweep_betas( perfect_df, betas=np.arange(-5, 5, 0.1))plot_cost_vs_beta(negloglik_sweep_betas_perfect_df, width=0.1)
您可以立即看到，现在不存在β的最小值：如果我们使用非常大的负β值，我们最终得到的对数似然拟合对所有事件都接近完美。
现在，让我们深入了解这背后的数学，并查看事件A的似然函数。我们将挖掘当我们调整β时分子和分母如何改变：
当β很高或为一个很大的正数时，分母中的最后一项（具有最大x的1.2），表示受试者E的风险，主导了整个分母并变得非常大。因此，似然函数变得很小，并趋近于零：
这会导致一个非常大的负对数似然。每个单独的似然也会发生同样的情况，因为E主体的最后一个危险将始终超过分子中的任何危险。结果，对于A到D主体，负对数似然增加。在这种情况下，当β很高时，它会降低所有的可能性，导致所有事件的拟合不良。
现在，当β很低或是一个很大的负数时，分母中的第一项，代表A主体的危险，会占据主导地位，因为它具有最低的x值。由于分子中也出现了A主体的相同危险，通过使β越来越负，可能性L（A）可以任意接近1，从而创建几乎完美的拟合：
所有其他单独的可能性也是如此：负的β现在同时增强了所有事件的可能性。基本上，拥有负的β并没有任何不利因素。同时，某些单独的危险增加（具有负x的A和B主体），有些保持不变（x=0的C主体），而其他的减少（具有正x的D主体）。但要记住，这里真正重要的是危险的比率。我们可以通过绘制单独的危险来验证这个草率的数学：
def plot_likelihoods(df, ylim=[-20, 20]): betas = np.arange(ylim[0], ylim[1], 0.5) subjects = df.query("event == 1")['subject'].tolist() likelihoods_per_subject = [] for subject in subjects: likelihoods = [ np.exp(log_likelihood(df, subject, beta)) for beta in betas ] likelihoods_per_subject.append( pd.DataFrame({ 'beta': betas, 'likelihood': likelihoods, 'subject': [subject] * len(betas), }) ) lik_df = pd.concat(likelihoods_per_subject) return ( p9.ggplot(lik_df, p9.aes('beta', 'likelihood', color='subject')) + p9.geom_line(size=2) + p9.theme_classic() ) plot_likelihoods(perfect_df)
将可能性放在一起的方式是，将一个危险除以所有仍处于风险中的主体的危险之和的比率，这意味着负的β值使得每个事件时间等级大于或等于预测器等级的主体的可能性拟合得非常完美！附带一提，如果x与事件时间有完美的负Spearman相关性，情况将会翻转：任意正的β都会给我们带来任意好的拟合。
预测器和时间排名不匹配
实际上，我们可以通过另一个虚构的例子看到这一点，并展示当事件时间排名和预测器排名不匹配时会发生什么：
sample_df = pd.DataFrame({ 'subject': ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H'], 'x': [-1.7, -0.4, 0.0, 0.5, 0.9, 1.2, 1.3, 1.45], 'time': [1, 2, 4, 3, 5, 7, 6, 8], 'rank_x': [1, 2, 3, 4, 5, 6, 7, 8], 'event': [1, 1, 1, 1, 1, 1, 1, 0],})sample_df
在这个特定的例子中，time列的取值范围从1到8，每个值代表自己的排名。我们还有一个x_rank列，对预测器x进行排名。现在，这里是关键观察：对于D和G主体，它们的x_rank实际上高于它们相应的time排名。因此，当我们有大的负β值时，D和G的可能性不会经历分子和分母之间的抵消效应：
它们的可能性现在在β的某些中间有限值上达到最大。 让我们来看一下单个可能性的图以了解这一点：
plot_likelihoods(sample_df)
这些时间和预测器之间的“不匹配”排名发挥了关键作用：当我们有显着负的β时，它们阻止所有可能性从本质上折叠成一个。
总之，在Cox回归中，为了获得预测器x的有限系数β，我们需要至少有一个实例，其中预测器x的排名低于事件时间的排名。
完美确实是好（p值）的敌人
那么，这些完美的预测器在实际情况下如何表现？ 好吧，为了找出答案，让我们再次转向lifelines库进行一些研究：
from lifelines import CoxPHFitterperfect_cox_model = CoxPHFitter()perfect_cox_model.fit( perfect_df, duration_col='time', event_col='event', formula='x')perfect_cox_model.print_summary()
#> /.../coxph_fitter.py:1586: ConvergenceWarning:#> The log-likelihood is getting suspiciously close to 0 and the delta is still large.#> There may be complete separation in the dataset.#> This may result in incorrect inference of coefficients.#> See https://stats.stackexchange.com/q/11109/11867 for more.#> /.../__init__.py:1165: ConvergenceWarning:#> Column x has high sample correlation with the duration column.#> This may harm convergence.#> This could be a form of 'complete separation'.#> See https://stats.stackexchange.com/questions/11109/how-to-deal-with-perfect-separation-in-logistic-regression#> /.../coxph_fitter.py:1611:#> ConvergenceWarning: Newton-Rhaphson failed to converge sufficiently.#> Please see the following tips in the lifelines documentation:#> https://lifelines.readthedocs.io/en/latest/Examples.html#problems-with-convergence-in-the-cox-proportional-hazard-model
就像在逻辑回归中一样，我们遇到了收敛警告，并且我们的预测器系数β的置信区间非常宽。 因此，我们最终得到了一个0.93的p值！
如果我们仅基于p值过滤模型而不考虑或进行进一步的调查，我们将忽略这些完美的预测器。
为了解决这个收敛问题，lifelines库文档和一些有用的StackOverflow线程提出了一个潜在的解决方案：将正则化项纳入成本函数中。 该项有效地增加了大系数值的成本，您可以通过将惩罚器参数设置为大于零的值来激活L2正则化：
perfect_pen_cox_model = CoxPHFitter(penalizer=0.01, l1_ratio=0)perfect_pen_cox_model.fit(perfect_df, duration_col='time', event_col='event', formula='x')perfect_pen_cox_model.print_summary()
这种方法解决了收敛警告，但并没有真正缩小那个讨厌的p值。 即使使用此正则化技巧，完美预测器的p值仍然在相当大的值0.11左右。
时间是相对的：只有排名才重要
最后，我们将使用之前的例子验证事件时间的绝对值对Cox回归拟合没有影响。为此，我们将引入一个新的列称为time2，其中将包含与time列相同顺序的随机数：
sample_df = pd.DataFrame({ 'subject': ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H'], 'x': [-1.7, -0.4, 0.0, 0.5, 0.9, 1.2, 1.3, 1.45], 'time': [1, 2, 4, 3, 5, 7, 6, 8], 'rank_x': [1, 2, 3, 4, 5, 6, 7, 8], 'event': [1, 1, 1, 1, 1, 1, 1, 0],}).sort_values('time')np.random.seed(42)sample_df['time2'] = sorted(np.random.randint(low=-42, high=888, size=8))sample_df
它们的拟合确实相同：
sample_cox_model = CoxPHFitter()sample_cox_model.fit( sample_df, duration_col='time', event_col='event', formula='x')sample_cox_model.print_summary()
sample_cox_model = CoxPHFitter()sample_cox_model.fit( sample_df, duration_col='time2', event_col='event', formula='x')sample_cox_model.print_summary()
结论
我们从中学到了什么？
- 在生存模型中，完美的预测变量是那些排名完全匹配事件时间的预测变量。
- Cox回归无法用有限的系数β拟合这些完美的预测变量，导致广泛的置信区间和大的p值。
- 事件时间的实际值并不重要，一切都取决于它们的排名。
- 当事件时间和预测变量的排名不对齐时，我们在似然度中不会得到大β值的便利抵消效果。因此，我们需要至少有一个情况，其中排名不匹配，才能获得具有有限系数的拟合。
- 即使我们尝试一些花哨的正则化技术，完美的预测变量仍然可能在现实情况下给我们带来那些令人烦恼的广泛置信区间和高p值。
- 就像在逻辑回归中一样，如果我们不太关心这些p值，使用正则化方法仍然可以为我们提供一个方便的模型拟合，以获得正确的预测！
如果您想自行运行代码，请随意使用我的Github上的IPython笔记本：https://github.com/igor-sb/blog/blob/main/posts/cox_perfect.ipynb
下次再见！👋