11  生存分析与删失数据(Survival Analysis and Censored Data)

本章我们将探讨生存分析和删失生存数据。这些方法源于一种独特的结局变量类型:事件发生时间

11.1 学习闭环

先修与可观察目标

完成本章后,学生应能:

  1. 从事件定义、时间原点、观察终点构造 \((Y,\delta)\),并逐行审计右删失编码,错误率为零。
  2. 对 10 个对象的小表写出每个风险集并手算 Kaplan–Meier 曲线,概率误差不超过 \(10^{-3}\)
  3. 写出 Cox 偏似然、score 与信息量,解释基线风险为何消去,并区分风险比与绝对风险。
  4. 诊断比例风险与独立删失假设;若不成立,提出时间交互、分层或替代模型。
  5. 在时间锁定的生存任务中报告删失率、事件数、C-index 与失败条件;将未惩罚与 Lasso-Cox 的稳定性比较列为拓展目标。

入口检查(5 分钟)

某公司在研究开始后第 3 年退市,另一公司到第 5 年末仍上市。写出两者的 \((Y,\delta)\)

展开入口评分、补救与异形复测

答案分别为 \((3,1)\)\((5,0)\)。若把第二家公司写成事件,先回到 式 11.1式 11.2,完成“提前还款/研究结束/失访”三例分类后重测。

无提示检索、渐隐链与迁移

闭书写出“风险集 \(r_j\)—事件数 \(d_j\)—条件生存率—累乘”链。第一轮对照正文小表;第二轮只给事件/删失时点;第三轮把第 5 月事件改到第 6 月并自行更新风险集。陌生迁移任务是为订阅客户定义“从首次付费到流失”的原点、事件、删失和竞争风险,说明为什么不能把静态状态当成事件日期。

目标—评价映射

目标 评价证据 达标标准
1—2 练习 1、2、6、7、9 风险集和时点逐行可复算
3 练习 10—12 核心题 10 给出偏似然与 score/info 骨架;11—12 为拓展推导
4 练习 3—5 能指出假设破坏及修复
5 练习 10 与时间队列案例 核心出口含删失率、事件数、时间外 C-index 和失败条件;练习 8 的惩罚路径为拓展

本章练习均标为 项目:无;生存分析不属于贯穿项目里程碑。但 小节 11.10 是两条标准路线都必交的独立审计 A11:60 分钟、10 分章内证据、占课程成绩 5%;它不能替代任何 M02—M09、M12 或 M13 里程碑,也不能用 A10 替换。

核心概念:生存分析的应用场景

生存分析虽然名称源于医学研究,但其应用远不止于医学领域:

  1. 医学研究:患者生存时间、疾病复发时间
  2. 商业分析:客户流失(churn)时间、产品使用寿命
  3. 工程技术:设备故障时间、系统可靠性
  4. 金融领域:信用卡违约时间、贷款提前还款时间

11.2 生存时间与删失时间

对于每个个体,存在两个真实时间:

  • 生存时间(Survival Time) \(T\):事件发生的真实时间(也称为失败时间或事件时间)
  • 删失时间(Censoring Time) \(C\):删失发生的时刻

我们实际观测到的是:

\[ Y = \min(T, C) \tag{11.1}\]

同时观测到一个状态指示变量:

\[ \delta = \begin{cases} 1 & \text{如果 } T \leq C \text{(观测到事件)} \\ 0 & \text{如果 } T > C \text{(被删失)} \end{cases} \tag{11.2}\]

直观理解:删失的含义

假设你进行一项为期 5 年的客户流失研究:

  • 如果客户在第3年取消服务,你观测到完整的流失时间
  • 如果客户在第 5 年研究结束时仍然活跃,你只知道其流失时间至少为 5 年,但不知道确切值
  • 这第二种情况就是右删失(right censoring)

11.2.1 删失机制的重要假设

为了分析生存数据,我们需要对删失机制做出关键假设:

独立删失假设:在给定特征条件下,事件时间 \(T\) 与删失时间\(C\) 独立。

违反独立删失假设的例子

  1. 病情导致的失访:如果病情严重的患者更倾向于退出研究,会导致对平均生存时间的高估
  2. 选择性失访:如果男性重症患者比女性重症患者更容易失访,可能导致错误的性别生存时间比较

这些情况都违反了独立删失假设,分析结果会有偏差。

11.2.2 删失的类型

  1. 右删失(Right Censoring):在本章 \(\delta=I(T\le C)\) 的约定下,删失个体满足 \(Y=C\) 且真实事件时间 \(T>Y\)\(T=Y\) 会被编码为事件而不是删失
    • 最常见的删失类型
    • 例如:研究结束时患者仍存活
  2. 左删失(Left Censoring):真实事件时间\(T \leq Y\)
    • 例如:调查时事件已经发生,但不知道确切时间
  3. 区间删失(Interval Censoring):只知道事件发生在某个时间区间内
    • 例如:每周随访一次,只知道事件发生在两次随访之间

本章重点讨论右删失

11.3 Kaplan-Meier 生存曲线

生存函数(Survival Function)定义为:

\[ S(t) = \Pr(T > t) \tag{11.3}\]

这是时间 \(t\) 的递减函数,表示生存超过时间\(t\) 的概率。

11.3.1 Kaplan-Meier 估计量

\(t_1 < t_2 < \cdots < t_K\) 为不同的事件时刻,\(d_k\)\(t_k\) 的事件数,\(q_k\) 为该时刻事件发生后退出观察的删失数,\(r_k\)\(t_k\) 之前仍在风险集中的个体数。

Kaplan–Meier 乘积极限估计以各事件时点的条件生存概率连乘来处理右删失 (Kaplan 和 Meier 1958年)

\[ \widehat{S}(t_k) = \prod_{j=1}^{k} \frac{r_j - d_j}{r_j} \tag{11.4}\]

对于介于 \(t_k\)\(t_{k+1}\) 之间的时间 \(t\),设 \(\widehat{S}(t)=\widehat{S}(t_k)\)。删失不会在当时降低生存概率,但会减少后续风险集。

Kaplan-Meier 的直观解释

考虑 10 个产品。第 2 个月有 2 个失效,随后 1 个删失;第 5 个月又有 1 个失效。第 2 月事件发生前 \(r_1=10,d_1=2\),所以 \(\widehat S(2)=8/10=0.8\)。删失后剩 7 个进入后续风险集,因此第 5 月 \(r_2=7,d_2=1\)\[ \widehat S(5)=\frac{8}{10}\times\frac{6}{7}\approx0.686. \] 这张小时间线同时展示了事件降低生存概率、删失只改变后续分母的区别。

11.3.2 中国实际案例:A股上市公司退市生存分析

让我们使得A 股上市公司的基本数据,来进行一个关于公司退市(Delisting) 的生存分析。随着中国推行全面注册制并常态化退市机制,上市公司的“存活状态”成了市场和监管层关注的重点焦点。

  • 时间起点:公司上市日的(listed_date)。
  • 事件de_listed_date 有记录且不晚于研究截止日期的“任何原因退市”;截止日当日事件按 \(T=C\) 编码为事件。本数据没有可靠区分主动退市、并购吸收、监管退市等竞争事件,因此结果不能直接解释成财务困境或破产风险。
  • 删失:截至数据截止,公司仍在正常交易(右删失)。

下面用 stock_basic_data 演示任何原因退市的 Kaplan–Meier 描述。索引所列核心动作的本地 setup 是 表 11.1:运行前设置 BOOK_DATA_DIR,随后按源码执行至 列表 11.1;Cox consumer 另需 表 11.3 的两个 BOOK_M11_* 环境变量,再执行至 表 11.6。核心报告标签分别为 ch11_survival_sample_auditch11_cox_results;缺任一输入时应停止,不能借用交互环境中的旧对象。曲线反映本样本、当前截止日和混合退市原因下的存续分布;它不能单独解释行业因果机制、财务困境或投资价值。

表 11.1: A股退市生存数据摘要
import pandas as pd  # 导入pandas数据分析库
import numpy as np  # 导入numpy数值计算库
import matplotlib.pyplot as plt  # 导入matplotlib绑图库
plt.rcParams['font.sans-serif'] = ['Source Han Serif SC']  # 使用系统已安装的思源黑体显示中文,避免字体缺失告警
plt.rcParams['axes.unicode_minus'] = False  # 修正负号显示为方块的问题
from lifelines import KaplanMeierFitter  # 导入Kaplan-Meier生存函数估计器
import os  # 导入操作系统模块用于跨平台路径处理

from pathlib import Path  # 使用跨平台路径对象解析显式数据根
DATA_DIR = Path(os.environ['BOOK_DATA_DIR']).expanduser().resolve()  # 从必需环境变量取得数据根
if not DATA_DIR.is_dir():  # 在读取前验证数据根
    raise FileNotFoundError(f'BOOK_DATA_DIR 不存在或不是目录: {DATA_DIR}')  # 失败即停止,不回退固定路径
path_basic = DATA_DIR / 'stock' / 'stock_basic_data.h5'  # 拼接股票基本信息文件路径
stock_basic_data = pd.read_hdf(path_basic)  # 直接读取A股上市公司基本信息
stock_basic_data = stock_basic_data.copy()  # 创建副本避免SettingWithCopyWarning
stock_basic_data['listed_date'] = pd.to_datetime(stock_basic_data['listed_date'], errors='coerce')  # 上市日期转datetime
stock_basic_data['de_listed_date'] = pd.to_datetime(stock_basic_data['de_listed_date'], errors='coerce')  # 退市日期转datetime
stock_basic_data = stock_basic_data.dropna(subset=['listed_date'])  # 移除缺少上市日期的记录
/tmp/ipykernel_1432188/1500803791.py:17: UserWarning: Could not infer format, so each element will be parsed individually, falling back to `dateutil`. To ensure parsing is consistent and as-expected, please specify a format.
  stock_basic_data['de_listed_date'] = pd.to_datetime(stock_basic_data['de_listed_date'], errors='coerce')  # 退市日期转datetime

表 11.2 用三个确定性教学行锁定 式 11.2 的边界约定:\(T=C\) 必须归入事件,只有 \(T>C\) 才是右删失。该小表不代表经验样本。

import pandas as pd  # 允许边界fixture在fresh kernel独立运行
censoring_boundary_fixture = pd.DataFrame({'case': ['T<C', 'T=C', 'T>C'], 'event_time': [2.0, 3.0, 5.0], 'censor_time': [4.0, 3.0, 4.0]})  # 构造三个确定性边界情形
censoring_boundary_fixture['observed_time'] = censoring_boundary_fixture[['event_time', 'censor_time']].min(axis=1)  # 按Y=min(T,C)构造观测时间
censoring_boundary_fixture['event'] = censoring_boundary_fixture['event_time'].le(censoring_boundary_fixture['censor_time']).astype(int)  # 按T<=C编码事件
assert censoring_boundary_fixture['event'].tolist() == [1, 1, 0]  # T=C必须与T<C同属事件
assert censoring_boundary_fixture['observed_time'].tolist() == [2.0, 3.0, 4.0]  # 逐行验证观测时间
censoring_boundary_fixture  # 输出可人工审计的三行证据
表 11.2: T<C、T=C 与 T>C 的事件—删失边界断言
case event_time censor_time observed_time event
0 T<C 2.0 4.0 2.0 1
1 T=C 3.0 3.0 3.0 1
2 T>C 5.0 4.0 4.0 0

数据加载完成后,我们需要遍历每一家上市公司,根据其上市日期与退市日期计算存续时间,构建生存分析所需的事件-时间数据集。

observation_end_date = pd.to_datetime('2024-01-01')  # 设定研究观测截止日期
survival_records = []  # 初始化列表,用于收集每家公司的生存记录

for _, row in stock_basic_data.iterrows():  # 遍历每一家上市公司
    listed_date = row['listed_date']  # 获取上市日期
    delisted_date = row['de_listed_date']  # 获取退市日期
    if pd.isna(delisted_date) or delisted_date > observation_end_date:  # 缺失或截止日后退市才右删失
        duration = (observation_end_date - listed_date).days / 365.25  # 从上市到截止日的存续年数
        event = 0  # 标记为删失(未观测到退市事件)
    else:  # 截止日当日或之前退市均为事件
        duration = (delisted_date - listed_date).days / 365.25  # 从上市到退市的存续年数
        event = 1  # 标记为事件(退市)
    if duration > 0:  # 排除存续时间为负的异常记录
        survival_records.append({  # 将计算结果追加到列表
            'order_book_id': row['order_book_id'],  # 股票代码
            'duration': duration,  # 存续年数
            'event': event,  # 事件标记(0=删失,1=退市)
            'industry': row.get('citics_2019_l1_name', '未知')  # 中信一级行业分类
        })  # 完成构建

company_survival_data = pd.DataFrame(survival_records)  # 将列表转换为DataFrame
top_industry_categories = company_survival_data['industry'].value_counts().nlargest(5).index  # 提取样本量前五的行业
company_survival_data['Industry_Group'] = company_survival_data['industry'].apply(  # 对每行/列应用函数
    lambda x: x if x in top_industry_categories else 'Other'  # 非前五行业统一归类为Other
)  # 完成构建
列表 11.1: 退市生存样本量与事件率审计
print('退市生存数据总览:')  # 打印数据预览标题
print(company_survival_data.head())  # 展示前五行数据
print(f'\n总计样本量: {len(company_survival_data)}')  # 输出总样本量
print(f'退市事件发生率: {company_survival_data["event"].mean():.2%}')  # 输出退市率
print({'report_label': 'ch11_survival_sample_audit', 'n': len(company_survival_data), 'events': int(company_survival_data['event'].sum())})  # 输出索引动作的稳定报告标签
退市生存数据总览:
  order_book_id   duration  event industry Industry_Group
0   000001.XSHE  32.747433      0       银行          Other
1   000002.XSHE  32.922656      0      房地产          Other
2   000003.XSHE  10.948665      1      NaN          Other
3   000004.XSHE  33.084189      0      计算机            计算机
4   000005.XSHE  33.059548      0      NaN          Other

总计样本量: 5334
退市事件发生率: 4.44%
{'report_label': 'ch11_survival_sample_audit', 'n': 5334, 'events': 237}

生存数据构建完毕后,我们首先绘制 A 股市场整体的 Kaplan-Meier 生存曲线,展示所有上市公司随时间推移的”存活概率”变化趋势。

kaplan_meier_fitter = KaplanMeierFitter()  # 实例化KM生存函数估计器
kaplan_meier_fitter.fit(  # 拟合整体生存曲线
    company_survival_data['duration'],  # 存续时间
    company_survival_data['event'],  # 事件标记
    label='A股整体生存率 (未退市)'  # 图例标签
)  # 完成构建
fig, ax = plt.subplots(figsize=(10, 6))  # 创建10×6英寸画布
kaplan_meier_fitter.plot(ax=ax)  # 绘制KM生存曲线
ax.set_title('A股上市公司生存曲线 (Time to Delisting)', fontsize=14, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置标题
ax.set_xlabel('上市年数', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 x 轴标签
ax.set_ylabel('生存概率 (未退市)', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 y 轴标签
ax.grid(True, alpha=0.3)  # 添加半透明网格线
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图形
阶梯曲线从一开始,随任何原因退市事件向下跳跃,删失标记不造成当时下降,横轴为上市年数。
图 11.1: A股上市公司整体 Kaplan-Meier 生存曲线

请从 图 11.1 的当次运行输出读取曲线水平与置信区间。只能将它解释为本数据覆盖、观察截止日与“任何原因退市”定义下的描述性存续分布。曲线不能自动识别制度执行、“壳价值”或退市门槛的因果作用;尾部置信区间变宽则提醒风险集中样本减少。

整体生存曲线展示了 A 股市场的”宏观韧性”。接下来,我们进一步按行业分解,比较前五大行业的退市风险差异。

fig, ax = plt.subplots(figsize=(10, 6))  # 创建10×6英寸画布
for ind in top_industry_categories:  # 遍历前五大行业
    if ind != '未知':  # 跳过行业信息缺失的记录
        mask = company_survival_data['industry'] == ind  # 构建当前行业的布尔掩码
        kaplan_meier_fitter.fit(  # 拟合当前行业的KM生存曲线
            company_survival_data.loc[mask, 'duration'],  # 当前行业的存续时间
            company_survival_data.loc[mask, 'event'],  # 当前行业的事件标记
            label=ind  # 使用行业名称作为图例
        )  # 完成构建
        kaplan_meier_fitter.plot(ax=ax)  # 将曲线绘制到同一画布

ax.set_title('不同行业A股的退市生存曲线比较', fontsize=14, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置标题
ax.set_xlabel('上市年数', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 x 轴标签
ax.grid(True, alpha=0.3)  # 添加网格线
plt.legend(prop={'family': 'Source Han Serif SC'})  # 使用已安装中文字体设置图例
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图形
多条按行业着色的 Kaplan-Meier 阶梯曲线在同一上市年数横轴上比较未退市概率。
图 11.2: 不同行业A股上市公司的退市生存曲线比较

图 11.2 的行业、曲线顺序和区间宽度都由当次样本决定。报告时应点名实际出现的行业、给出风险集人数和事件数,并用下一节的预先规定检验补充目测。即使曲线分离,也只是当前混合退市结局的组间关联,不是行业因果效应。

11.3.3 对数秩检验:行业间比较

我们使用对数秩检验来比较不同行业的“仍维持上市”生存曲线是否存在显著差异。

当我们用肉眼观察上一节中不同行业的Kaplan-Meier 生存曲线时,可能会觉得它们之间似乎存在明显的高低差异。但作为严谨的量化分析师,我们绝不能仅仅依靠视觉来做判断——由于各个行业在不同历史时期的上市节奏(即进入风险集的队列)以及最终被删失的比例大相径庭,那些看似分离的曲线,在统计学上是否真的存在不能被随机波动所解释的本质区别? 对数秩检验在各事件时点比较两组的观察事件数与零假设下期望事件数。它给出当前分组下的统计证据,不预设 p 值大小,也不能把组间差异解释为行业“基因”或因果机制。

from lifelines.statistics import logrank_test  # 导入对数秩检验用于比较生存曲线差异

# 假设 company_survival_data 已经存在 (来自上一节)
# 我们选取两个主要行业进行比较,例如"电子" (Tech) 和"房地产" (Real Estate)
# 注意: 具体的行业名称取决于数据中的 industry_name
# 这里我们先打印一下行业名称看看(假设前两个Top行业)

industry_first = top_industry_categories[0]  # 定义industry_first变量
industry_second = top_industry_categories[1]  # 定义industry_second变量

print(f"比较行业: {industry_first} vs {industry_second}")  # 输出结果到控制台

mask_industry_first = company_survival_data['industry'] == industry_first  # 定义mask_industry_first变量
mask_industry_second = company_survival_data['industry'] == industry_second  # 定义mask_industry_second变量

rank_test_results = logrank_test(  # 提取测试集数据
    company_survival_data.loc[mask_industry_first, 'duration'],  # 执行数据处理操作
    company_survival_data.loc[mask_industry_second, 'duration'],  # 执行数据处理操作
    company_survival_data.loc[mask_industry_first, 'event'],  # 执行数据处理操作
    company_survival_data.loc[mask_industry_second, 'event']  # 执行数据处理操作
)  # 完成构建

print(f'对数秩检验统计量: {rank_test_results.test_statistic:.4f}')  # 输出结果到控制台
print(f'p值: {rank_test_results.p_value:.4f}')  # 输出结果到控制台

fig, ax = plt.subplots(figsize=(10, 6))  # 创建子图布局
# 训练/拟合模型
kaplan_meier_fitter.fit(company_survival_data.loc[mask_industry_first, 'duration'], company_survival_data.loc[mask_industry_first, 'event'], label=industry_first)
kaplan_meier_fitter.plot(ax=ax)  # 绘制图形
# 训练/拟合模型
kaplan_meier_fitter.fit(company_survival_data.loc[mask_industry_second, 'duration'], company_survival_data.loc[mask_industry_second, 'event'], label=industry_second)
kaplan_meier_fitter.plot(ax=ax)  # 绘制图形

# 设置子图标题
ax.set_title(f'行业生存曲线比较: {industry_first} vs {industry_second}', fontsize=14, fontproperties='Source Han Serif SC')
ax.set_xlabel('上市年数', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置子图 X 轴标签
ax.grid(True, alpha=0.3)  # 添加子图网格线
plt.legend(prop={'family': 'Source Han Serif SC'})  # 使用已安装中文字体添加图例
plt.show()  # 显示图形
比较行业: 机械 vs 基础化工
对数秩检验统计量: 0.9421
p值: 0.3317
两个样本量最大的行业以两条阶梯生存曲线叠加显示,并配合正文报告 log-rank 检验。
图 11.3: 不同行业的生存曲线比较与检验

图 11.3 上方输出给出本次动态选择行业后的检验统计量和 p 值。应按预先规定的显著性水平报告“拒绝”或“未拒绝”零假设,并同时说明组别是按样本量选择、事件包含所有退市原因。目测曲线分离不能替代检验,未拒绝也不等于两组完全相同。

11.4 风险函数的Cox 比例风险模型

除了比较曲线,我们更希望量化各个因素(如初始财务状况、所在地域)对遭遇退市风险的影响。

这里先定义 Cox 模型使用的风险函数:

\[ h(t\mid x)=\lim_{\Delta t\to0}\frac{\Pr(t<T\le t+\Delta t\mid T>t,x)}{\Delta t}. \]

它是已存续到 \(t\) 时事件在紧接着的瞬间发生的条件率,不是普通概率。本节使用这一定义建模,后文再详述它与生存函数、累计风险的关系。

数学推导:Cox 比例风险模型与偏似然函数 (Partial Likelihood)

Cox 比例风险模型 (Proportional Hazards Model) 的核心思想是将风险函数 \(h(t|x_i)\) 分解为两部分。 \[ h(t|x_i) = h_0(t) \exp(x_i^T \beta) \] 其中 \(h_0(t)\)基线风险函数(Baseline Hazard Function),表示当所有协变量为0 时的潜在基准风险。\(\exp(x_i^T \beta)\) 是相对基准的风险乘数

Cox 比例风险模型通过偏似然估计回归参数 \(\beta\),不必为基准风险函数 \(h_0(t)\) 预先指定参数形式 (Cox 1972年);这项性质不意味着绝对风险已经由偏似然单独识别。

假设所有样本中的\(K\) 个非删失的独特事件(例如退市)发生时间 \(t_1 < t_2 < \dots < t_K\)。令 \(\mathcal{R}(t_i)\) 为在时间 \(t_i\) 刚好处于风险集(Risk Set)中的所有个体的集合(含义是:那些存活到\(t_i\) 并且尚未发生事件、也未被删失的全部个体)。

给定在时的\(t_i\) 确切地一个人发生了事件(暂不考虑 Tied Times 打平),那么这个人恰巧是受试的\(j \in \mathcal{R}(t_i)\) 的条件概率依据风险律直接写出为: \[ P(\text{个体 } j \text{ 在} t_i \text{ 发生事件} \mid \text{在} t_i \text{ 有一个事件发生}) = \frac{h(t_i|x_j)}{\sum_{k \in \mathcal{R}(t_i)} h(t_i|x_k)} \] 清爽地代入Cox 模型的风险函数乘性公式: \[ \frac{h_0(t_i) \exp(x_j^T \beta)}{\sum_{k \in \mathcal{R}(t_i)} h_0(t_i) \exp(x_k^T \beta)} = \frac{\exp(x_j^T \beta)}{\sum_{k \in \mathcal{R}(t_i)} \exp(x_k^T \beta)} \]

奇妙的数学分离在这里发生了:未知的非参数基准常数的\(h_0(t_i)\) 在多项式分子和分母中被平顺地消去了!

然后将所有发生独立事件的时刻 \(t_1, \ldots, t_K\) 的这些离散的条件概率在时间轴上累积相乘,即可得到全样本的偏似然函数: \[ L(\beta) = \prod_{i=1}^K \frac{\exp(x_{(i)}^T \beta)}{\sum_{k \in \mathcal{R}(t_i)} \exp(x_k^T \beta)} \] 其中 \(x_{(i)}\) 就是命运指针抽中在时的\(t_i\) 发生失败(如退市)的那个“倒霉”候选人的特征向量。

对该函数取对数,得到对数偏似然函数 \(\ell(\beta)=\log L(\beta)\),可以用牛顿—拉弗森等数值方法估计 \(\hat{\beta}\)。半参数设计避免预先指定基线风险形状,但其适用性仍依赖比例风险、删失和风险集定义;是否优于其他生存模型必须由相同数据与评价设计决定。

11.4.1 准备数据:point-in-time landmark 设计

Cox 模型允许加入连续变量,但报告期末不是信息公开日。本节与第 12—13 章统一采用 point-in-time 版本合同:课程运行器必须提供含可靠 info_date 的冻结财务切片、合同版本、内容哈希和冻结日;若缺 info_date、哈希不匹配或任一记录的可得日晚于其 landmark,流程必须结构化停止。不得再用“报告期末 + 固定天数”替代披露日,也不得把后续修订回填到较早 landmark。

表 11.3: Chapter 11—13 财务信息的 point-in-time 共同合同
合同项 要求
contract_id ch11-cox-point-in-time-v1
输入接口 BOOK_M11_FINANCIAL_SLICEBOOK_M11_FINANCIAL_SLICE_SHA256
必需 schema order_book_id:str, quarter:str, info_date:datetime64[ns], total_assets:float64, total_liabilities:float64
版本规则 同公司同报告期保留 info_date <= landmark_date 的最后可得版本
结构化停止 missing_info_dateinvalid_info_dateinput_sha256_mismatchavailability_after_landmark
与第 12—13 章的共同不变量 先验冻结切片与哈希;用 info_date 选 point-in-time 版本;断言 availability <= landmark;失败时不拟合、不填补代理日
表 11.4: Cox模型数据准备
import hashlib  # 计算冻结切片的内容哈希
financial_statement_path = Path(os.environ['BOOK_M11_FINANCIAL_SLICE']).expanduser().resolve()  # 消费运行器提供的小切片
expected_financial_sha256 = os.environ['BOOK_M11_FINANCIAL_SLICE_SHA256'].strip().lower()  # 读取预注册哈希
actual_financial_sha256 = hashlib.sha256(financial_statement_path.read_bytes()).hexdigest()  # 在读取数据前校验字节内容
if actual_financial_sha256 != expected_financial_sha256:  # 拒绝未登记或被改写的切片
    raise RuntimeError({'status': 'stopped', 'reason': 'input_sha256_mismatch'})  # 哈希不符即结构化停止
financial_statement_raw = pd.read_hdf(financial_statement_path)  # 只在哈希通过后读取冻结切片
required_financial_columns = {'order_book_id', 'quarter', 'info_date', 'total_assets', 'total_liabilities'}  # 冻结最小schema
if not required_financial_columns.issubset(financial_statement_raw.columns):  # 检查可得日和建模字段
    raise RuntimeError({'status': 'stopped', 'reason': 'missing_info_date_or_required_column'})  # 缺可靠时点即停止
financial_statement_raw['available_date'] = pd.to_datetime(financial_statement_raw['info_date'], errors='coerce')  # 以披露或修订日定义可得日
if financial_statement_raw['available_date'].isna().any():  # 禁止缺失时点绕过版本守卫
    raise RuntimeError({'status': 'stopped', 'reason': 'invalid_info_date'})  # 无法解析即停止
financial_groups = financial_statement_raw.groupby('order_book_id')  # 按公司建立point-in-time候选版本组
列表 11.2: 逐公司提取首份可得财报的 landmark 杠杆率
# ── 2. 遍历每家公司,提取上市后第一份财报的资产负债率(杠杆率) ──
initial_leverage_records = []  # 用于存储每家公司的初始杠杆率

for _, company_row in company_survival_data.iterrows():  # 遍历迭代
    stock_id = company_row['order_book_id']  # 当前公司股票代码
    if stock_id not in financial_groups.groups:  # 该公司在财报数据中不存在,跳过
        continue  # 跳过本次循环
    # 获取该公司的全部财报记录
    company_financial_records = financial_groups.get_group(stock_id).sort_values('available_date')
    # 获取该公司上市日期
    basic_info_row = stock_basic_data[stock_basic_data['order_book_id'] == stock_id]  # 查找基本信息
    if len(basic_info_row) == 0:  # 未找到基本信息则跳过
        continue  # 跳过本次循环
    listed_date = pd.to_datetime(basic_info_row.iloc[0]['listed_date'])  # 上市日期(字符串转为datetime)

    # 以首次可靠披露日作为该公司的landmark候选
    post_listing_financials = company_financial_records[company_financial_records['available_date'] > listed_date]

    if len(post_listing_financials) > 0:  # 存在上市后的财报
        first_financial_record = post_listing_financials.iloc[0]  # 取第一份
        total_assets_value = first_financial_record.get('total_assets', None)  # 总资产
        total_liab_value = first_financial_record.get('total_liabilities', None)  # 总负债
        if pd.notna(total_assets_value) and total_assets_value > 0 and pd.notna(total_liab_value):  # 数据有效性检查
            leverage_ratio = total_liab_value / total_assets_value  # 计算资产负债率(杠杆率)
            initial_leverage_records.append({  # 将计算结果追加到列表
                'order_book_id': stock_id,  # 定义字典键值对条目
                'Initial_Leverage': leverage_ratio,  # landmark 时点杠杆率
                'Availability_Date': first_financial_record['available_date'],  # 保存特征版本的可靠可得日
                'Landmark_Date': first_financial_record['available_date']  # 以该版本首次可得日作为landmark
            })  # 完成构建

下面把逐公司记录转换为稳定表,并单独执行 point-in-time 断言;lst-cox-landmark-extractionlst-cox-landmark-validation 必须在同一 fresh kernel 中依次运行。

列表 11.3: 构造 landmark 特征表并执行逐行可得日断言
initial_leverage_df = pd.DataFrame(initial_leverage_records, columns=['order_book_id', 'Initial_Leverage', 'Availability_Date', 'Landmark_Date'])  # 保留版本可得日与landmark
if initial_leverage_df['Landmark_Date'].isna().any():  # landmark 必须来自可靠可得日
    raise RuntimeError({'status': 'stopped', 'reason': 'missing_landmark_date'})  # 缺时点则不继续建模
if not initial_leverage_df['Availability_Date'].le(initial_leverage_df['Landmark_Date']).all():  # 每行断言版本在landmark时已可得
    raise RuntimeError({'status': 'stopped', 'reason': 'availability_after_landmark'})  # 越界即结构化停止
表 11.5: landmark 合并、事件重编码与地区协变量
# ── 3. 将初始杠杆率合并到生存数据中 ──
company_survival_data = company_survival_data.merge(initial_leverage_df, on='order_book_id', how='left')  # 左连接合并

delisting_mapping = stock_basic_data[['order_book_id', 'de_listed_date']].drop_duplicates('order_book_id')  # 取得事件日期
company_survival_data = company_survival_data.merge(delisting_mapping, on='order_book_id', how='left')  # 合并退市日期
company_survival_data = company_survival_data.dropna(subset=['Landmark_Date']).copy()  # 只保留进入 landmark 风险集的公司
company_survival_data['event'] = (company_survival_data['de_listed_date'].notna() & company_survival_data['de_listed_date'].le(observation_end_date)).astype(int)  # 截止日当日或之前退市为事件
company_survival_data['analysis_end'] = company_survival_data['de_listed_date'].where(company_survival_data['event'].eq(1), observation_end_date)  # 事件或删失终点
company_survival_data['duration'] = (company_survival_data['analysis_end'] - company_survival_data['Landmark_Date']).dt.days / 365.25  # 从协变量可得时点计时
company_survival_data = company_survival_data[company_survival_data['duration'] > 0].copy()  # 排除无有效随访者

# 对未匹配到财报数据的公司,用中位数填充杠杆率
median_leverage = company_survival_data['Initial_Leverage'].median()  # 计算中位数
company_survival_data['Initial_Leverage'] = company_survival_data['Initial_Leverage'].fillna(median_leverage)  # 填充缺失值

# ── 4. 从 stock_basic_data 获取省份信息,构造地区分类变量 ──
province_mapping = stock_basic_data[['order_book_id', 'province']].drop_duplicates('order_book_id')  # 提取省份映射表
company_survival_data = company_survival_data.merge(province_mapping, on='order_book_id', how='left')  # 合并省份

# 定义沿海/内陆分类规则
COASTAL_PROVINCES = ['广东省', '浙江省', '江苏省', '上海市', '北京市', '福建省', '山东省', '天津市']  # 沿海经济发达省份
company_survival_data['Region'] = company_survival_data['province'].apply(  # 对每行/列应用函数
    lambda x: 'Coastal' if x in COASTAL_PROVINCES else 'Inland'  # 按省份分类为沿海或内陆
)  # 完成构建

print("Cox 模型数据准备完成:")  # 输出结果到控制台
print(company_survival_data[['Industry_Group', 'Region', 'Initial_Leverage']].head())  # 预览关键变量
print(f"\n杠杆率统计:\n{company_survival_data['Initial_Leverage'].describe()}")  # 杠杆率分布概览

上方现场输出给出实际样本量、协变量分布和缺失处理结果。这里 Landmark_Date 就是所选记录的可靠 info_date,代码已经断言 available_date <= landmark(取等号);解释 Cox 系数前还应检查极端杠杆值与修订版本分布。正文不固定记录可能随切片版本变化的样本数。

11.4.2 拟合 Cox 比例风险模型

在完成数据预处理后,我们使用 CoxPHFitter 拟合带 L2 惩罚的 Cox 模型。事件是任何原因退市,不是财务困境;惩罚模型的常规 p 值只作描述,不能按未惩罚模型的精确推断方式解释。行业、地区和杠杆的关系是条件关联,不是因果效应。

表 11.6: Cox 比例风险模型回归结果
from lifelines import CoxPHFitter  # 导入Cox比例风险模型拟合器

# 准备回归数据
cox_model_data = company_survival_data[['duration', 'event', 'Industry_Group', 'Region', 'Initial_Leverage']].copy()  # 选择Cox模型所需变量

# ── 填充缺失值,避免 get_dummies 后产生全 NaN 列 ──
leverage_median = cox_model_data['Initial_Leverage'].median()  # 计算杠杆率中位数
if pd.isna(leverage_median):  # 若中位数本身为 NaN(全部缺失),使用行业均值 0.5
    leverage_median = 0.5  # 定义leverage_median变量
cox_model_data['Initial_Leverage'] = cox_model_data['Initial_Leverage'].fillna(leverage_median)  # 填充杠杆率缺失值
cox_model_data['Industry_Group'] = cox_model_data['Industry_Group'].fillna('Other')  # 填充行业缺失值
cox_model_data['Region'] = cox_model_data['Region'].fillna('Inland')  # 填充地区缺失值

# 编码分类变量
cox_model_data = pd.get_dummies(cox_model_data, drop_first=True, dtype=float)  # 将分类变量转为哑变量(float类型确保兼容性)

# 仅删除关键列(duration/event)含 NaN 的行
cox_model_data = cox_model_data.dropna(subset=['duration', 'event'])  # 确保生存时间和事件标记完整

# ── 移除方差接近零的协变量列(防止 Hessian 矩阵奇异导致不收敛) ──
covariate_columns = cox_model_data.columns.difference(['duration', 'event'])  # 提取所有协变量列名
low_variance_threshold = 1e-6  # 设定极低方差阈值
low_variance_columns = [col for col in covariate_columns if cox_model_data[col].var() < low_variance_threshold]  # 识别低方差列
if low_variance_columns:  # 若存在低方差列则移除
    print(f'移除低方差列: {low_variance_columns}')  # 提示用户已移除的低方差列
    cox_model_data = cox_model_data.drop(columns=low_variance_columns)  # 删除低方差列

print('模型变量:', cox_model_data.columns.tolist())  # 输出最终模型使用的变量列表
print(f'有效样本量: {len(cox_model_data)}')  # 确认有效样本量

# 拟合 Cox 模型(添加 L2 正则化惩罚项,防止高共线性导致的矩阵奇异问题)
cox_proportional_hazards_model = CoxPHFitter(penalizer=0.1)  # penalizer=0.1 为正则化强度
# 训练/拟合模型
cox_proportional_hazards_model.fit(cox_model_data, duration_col='duration', event_col='event')

# 显示结果
cox_proportional_hazards_model.print_summary()  # 打印模型完整回归结果
cox_proportional_hazards_model.check_assumptions(cox_model_data, p_value_threshold=0.05, show_plots=False)  # 输出比例风险假设诊断
print({'report_label': 'ch11_cox_results', 'n': len(cox_model_data), 'events': int(cox_model_data['event'].sum())})  # 输出索引动作的稳定报告标签

Cox 表报告 landmark 样本中的条件关联。带惩罚估计的系数和常规 p 值不应机械套用 5% 阈值;风险比只描述比例风险假设成立时、其他已纳入协变量相同的相对退市率。若诊断拒绝比例风险假设,应考虑分层或时变效应,并停止给出全时段不变的风险比。

11.4.3 结果解释

Cox 模型的回归系数可以通过森林图展示。图中的系数描述协变量与“任何原因退市”风险率的条件关联;惩罚拟合下的区间与 p 值只作描述,不应称为财务困境因果效应。

fig, ax = plt.subplots(figsize=(10, 6))  # 创建10×6英寸画布
cox_proportional_hazards_model.plot(ax=ax)  # 绘制Cox模型系数的森林图
ax.set_title('Cox 模型:条件退市关联(Hazard Ratios)', fontsize=14, fontproperties='Source Han Serif SC')  # 事件与解释边界一致
ax.grid(True, alpha=0.3)  # 添加网格线
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图形

# 提取模型摘要并输出显著变量的风险比解释
print('\n描述性风险比(非因果效应):')  # 输出边界明确的标题
cox_model_summary = cox_proportional_hazards_model.summary  # 提取模型摘要表
for var in cox_model_summary.index:  # 遍历每个协变量
    coefficient = cox_model_summary.loc[var, 'coef']  # 获取回归系数
    hazard_ratio = np.exp(coefficient)  # 计算风险比(指数化系数)
    print(f' - {var}: HR={hazard_ratio:.3f}')  # 不把惩罚模型p值解释成确认性证据
图 11.4

图 11.4 用于描述样本中的条件关联,不能揭示因果“核心因子”。行业、地区和杠杆系数还可能受进入 landmark 风险集的选择、退市原因混合及遗漏变量影响;具体数值以当前执行结果为准。

11.5 风险函数与累计风险(Hazard Function)

风险函数(也称为风险率)定义为:

\[ h(t) = \lim_{\Delta t \to 0} \frac{\Pr(t < T \leq t + \Delta t | T > t)}{\Delta t} \tag{11.5}\]

这是在给定生存到时间 \(t\) 的条件下,在紧接着的瞬间发生事件的瞬时率

风险函数 vs 生存函数

  • 生存函数 \(S(t)\):回答生存到时间\(t\) 的概率是多少。
  • 风险函数 \(h(t)\):回答如果已经生存到时间\(t\),在下一瞬间发生事件的概率强度是多少。

关键关系。 \[ S(t) = \exp\left(-\int_0^t h(u) du\right) \tag{11.6}\]

11.5.1 中国案例:杠杆分组的累计风险

如果说生存函数是在俯瞰全局的“存活概率”,那么风险函数(Hazard Function)就是在凝视每一个当下的“死亡威胁”。 下面使用 NelsonAalenFitter 估计累计风险 \(H(t)\),并按 landmark 杠杆率分组。累计风险的跳跃不是瞬时风险率 \(h(t)\);若要估计 \(h(t)\),需要额外平滑与带宽选择。 Nelson–Aalen 曲线是累计风险的阶梯估计,跳跃表示累计事件强度增加,并不是某一时点的平滑瞬时风险。曲线本身不能识别制度改革、估值或公司基本面机制;这些解释需要预先规定的协变量、时间模型和额外证据。

from lifelines import NelsonAalenFitter  # 导入Nelson-Aalen累积风险估计器

# 创建杠杆率分组(上限使用 inf 以涵盖杠杆率超过100%的公司)
# 执行数据处理操作
company_survival_data['leverage_group'] = pd.cut(company_survival_data['Initial_Leverage'],
                                  bins=[0, 0.3, 0.6, float('inf')],  # 按杠杆率分为低、中、高三组
                                  labels=['低杠杆(0-30%)', '中杠杆(30-60%)', '高杠杆(>60%)'])  # 完成模型/Pipeline构建

fig, ax = plt.subplots(figsize=(10, 6))  # 创建10×6英寸画布

nelson_aalen_fitter = NelsonAalenFitter()  # 实例化Nelson-Aalen累积风险估计器

for leverage_group in company_survival_data['leverage_group'].dropna().unique():  # 仅遍历有数据的分组
    mask = company_survival_data['leverage_group'] == leverage_group  # 选取当前分组
    if mask.sum() < 2:  # 跳过观测不足的分组
        continue  # 跳过本次循环
    nelson_aalen_fitter.fit(company_survival_data.loc[mask, 'duration'],  # 训练/拟合模型
            company_survival_data.loc[mask, 'event'])  # 完成模型/Pipeline构建

    # 直接绘制 Nelson–Aalen 累计风险估计
    cumulative_hazard = nelson_aalen_fitter.cumulative_hazard_  # 定义cumulative_hazard变量
    cumulative_hazard.plot(ax=ax, label=leverage_group)  # 绘制当前分组的累计风险

ax.set_title('不同 landmark 杠杆率分组的累计退市风险', fontsize=14, fontproperties='Source Han Serif SC')  # 标题与实际变量一致
ax.set_xlabel('生存时间(年)', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 x 轴标签
ax.set_ylabel('累计风险 H(t)', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 y 轴标签
ax.grid(True, alpha=0.3)  # 添加网格线
plt.legend(title='杠杆率分组', title_fontproperties='Source Han Serif SC', prop={'family': 'Source Han Serif SC'})  # 图例与实际变量一致
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图形
图 11.5

图 11.5 展示各组从 landmark 起累积的退市事件强度。曲线越高表示截至该时点累计风险越大,但它既不是瞬时风险率,也不能单凭组间差异作因果解释。若曲线差异随时间明显改变,应进一步用残差诊断或时间交互项检验比例风险假设。

11.6 Cox 模型的正则化

当特征数量较多时,我们可以使用正则化 Cox 模型进行变量选择。

当我们赋予Cox 比例风险模型处理数百个财务指标的能力时,多重共线性和过拟合的阴影便会随之降临。幸运的是,我们在第 6 章学过的 Lasso 正则化同样可以完美嫁接到生存分析中。 下面的这段“抗噪压力测试”代码非常有趣:我们故意在原本真实的上市财务数据中,人为地注入了 5 列完全由随机数构成的纯噪音特征(Noise_0Noise_4)。接着,我们利了lifelines 库开启了 L1 惩罚(l1_ratio=1.0),让惩罚力度(Lambda)从极小逐渐呈指数级放大,并记录下每一个特征对应系数的挣扎轨迹。 路径图仅展示惩罚增大时各系数如何收缩。在有限样本、相关特征与一次随机噪声注入下,噪声变量不保证比业务变量更早归零,未归零变量也不等于真实风险因子。该图是惩罚机制演示;若要做筛选,还需在时间外验证中选惩罚强度,并用重抽样入选频率评估稳定性。

from sklearn.preprocessing import StandardScaler  # 导入标准化工具用于特征缩放
from lifelines import CoxPHFitter  # 导入Cox比例风险模型拟合器

# 增加一些噪声特征来测试 Lasso 的筛选能力
np.random.seed(1)  # 设定随机种子确保可复现性
for i in range(5):  # 循环生成5个噪声特征列
    cox_model_data[f'Noise_{i}'] = np.random.normal(0, 1, len(cox_model_data))  # 生成标准正态分布随机噪声

predictor_matrix = cox_model_data.drop(['duration', 'event'], axis=1)  # 提取协变量矩阵(排除响应变量)
predictor_matrix = predictor_matrix.select_dtypes(include=[np.number])  # 仅保留数值型列

# 标准化特征(使各变量均值为0、标准差为1)
scaler = StandardScaler()  # 实例化标准化器
scaled_predictor_matrix = scaler.fit_transform(predictor_matrix)  # 拟合并转换特征矩阵
scaled_predictor_df = pd.DataFrame(scaled_predictor_matrix, columns=predictor_matrix.columns)  # 转回DataFrame保留列名
scaled_predictor_df['duration'] = cox_model_data['duration'].values  # 将存续时间列加回
scaled_predictor_df['event'] = cox_model_data['event'].values  # 将事件标记列加回

以下代码使用带L1惩罚的Cox回归(Lasso正则化),通过逐步增大惩罚参数来观察各协变量系数的衰减路径。

# 使用 lifelines 的带惩罚的 Cox 回归
penalty_parameters = np.logspace(-4, 0, 20)  # 生成20个从10^-4到10^0的惩罚参数值
coefficient_path = []  # 初始化列表,用于记录每个惩罚参数下的系数向量

successful_penalty_parameters = []  # 记录成功收敛的惩罚参数值
for lambda_val in penalty_parameters:  # 遍历每个惩罚参数
    try:  # 极小惩罚值可能导致矩阵奇异,需捕获收敛错误
        penalized_cox_model = CoxPHFitter(penalizer=lambda_val, l1_ratio=1.0)  # L1正则化对应Lasso
        penalized_cox_model.fit(scaled_predictor_df, duration_col='duration', event_col='event')  # 拟合带惩罚的Cox模型
        coefficient_path.append(penalized_cox_model.params_.values)  # 记录当前惩罚下的系数向量
        successful_penalty_parameters.append(lambda_val)  # 记录成功的惩罚参数
    except Exception:  # 若收敛失败,跳过该惩罚值
        continue  # 跳过本次循环

penalty_parameters = np.array(successful_penalty_parameters)  # 更新为仅含成功值的数组
coefficient_path = np.array(coefficient_path)  # 将系数路径转换为二维数组

# 绘制正则化路径图
fig, ax = plt.subplots(figsize=(12, 8))  # 创建12×8英寸画布

for i, col in enumerate(predictor_matrix.columns):  # 遍历每个特征变量
    ax.plot(penalty_parameters, coefficient_path[:, i], label=col)  # 绘制该特征的系数随惩罚参数的变化路径

ax.set_xscale('log')  # 将x轴设为对数刻度
ax.set_xlabel('Lambda (正则化参数)', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 x 轴标签
ax.set_ylabel('系数值', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 y 轴标签
ax.set_title('Cox 模型的 Lasso 正则化路径', fontsize=14, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置标题
ax.grid(True, alpha=0.3)  # 添加网格线
plt.legend(prop={'family': 'Source Han Serif SC'}, loc='best', fontsize=8)  # 使用已安装中文字体设置图例
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图形

# 只描述当次运行中的收缩路径,不预设噪声变量归零顺序
图 11.6

图 11.6 展示一次样本、一次噪声注入下的惩罚路径。随着惩罚增大,系数趋向零;但某个噪声变量较早归零并不能证明选择一致性,也不能把其余变量认证为“真实信号”。稳定筛选需要重复生成或重抽样并报告入选频率。

11.7 生存树模型

生存树可以自动发现非线性关系和交互作用。

生存树可以用分割表示非线性与交互,但仍需要规定深度、叶节点样本量并做样本外验证。下面用 scikit-survival 拟合一棵限深树,再用置换重要性描述当前样本中特征与预测表现的关系。重要性没有方向,相关特征会分享或替代重要性,且本练习的训练样本重用结果不是因果“话语权”或外部预测证据。

# 导入生存树模型和数据结构化工具
from sksurv.tree import SurvivalTree  # 导入生存树模型
from sksurv.util import Surv  # 导入生存数据结构化工具

# 准备生存树所需的结构化目标变量
survival_targets = Surv.from_arrays(  # 构造生存目标数组
    event=cox_model_data['event'].astype(bool),  # 事件标记转为布尔型
    time=cox_model_data['duration']  # 存续时间
)  # 完成构建
tree_predictor_matrix = cox_model_data.drop(['duration', 'event'], axis=1)  # 提取特征矩阵

# 拟合生存树模型(限制树深为3以防止过拟合)
survival_tree_model = SurvivalTree(max_depth=3, min_samples_split=10, min_samples_leaf=5, random_state=42)  # 配置生存树参数
survival_tree_model.fit(tree_predictor_matrix, survival_targets)  # 拟合模型

# 计算并排序特征重要性(使用置换重要性替代不可用的impurity-based方法)
from sklearn.inspection import permutation_importance  # 导入置换重要性工具
perm_importance_result = permutation_importance(  # 计算置换重要性
    survival_tree_model, tree_predictor_matrix, survival_targets,  # 执行数据处理操作
    n_repeats=10, random_state=42  # 重复10次取平均
)  # 完成构建
calculated_feature_importance = pd.DataFrame({  # 构建DataFrame存储特征重要性
    'feature': tree_predictor_matrix.columns,  # 特征名称
    'importance': perm_importance_result.importances_mean  # 置换重要性均值
}).sort_values('importance', ascending=False)  # 按重要性降序排列

在完成生存树模型的拟合和特征重要性计算后,下面输出特征重要性排名并以水平条形图可视化展示,直观比较各财务指标对企业生存时间的预测贡献。

print('\n生存树特征重要性:')  # 输出标题
print('='*50)  # 输出分隔线
print(calculated_feature_importance)  # 打印特征重要性表

# 绘制特征重要性水平条形图
fig, ax = plt.subplots(figsize=(10, 6))  # 创建10×6英寸画布
top_feature_importance = calculated_feature_importance.head(10)  # 取前10个最重要的特征
ax.barh(range(len(top_feature_importance)),  # 绘制水平条形图
        top_feature_importance['importance'])  # 条形高度为重要性分数
ax.set_yticks(range(len(top_feature_importance)))  # 设置y轴刻度位置
ax.set_yticklabels(top_feature_importance['feature'], fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 y 轴刻度标签
ax.set_xlabel('特征重要性', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 x 轴标签
ax.set_title('生存树:特征重要性排名', fontsize=14, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置标题
ax.grid(True, alpha=0.3, axis='x')  # 添加x方向网格线
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图形
图 11.7

图 11.7 的重要性是当前训练样本和模型设定下的描述量。排序应从现场输出读取,并通过未来测试期置换重要性检查;它不能单独证明某变量具有经济主导性、因果作用或稳定预测能力。

11.8 模型评估:C-index

使用 C-index 评估模型对“谁先退市”的风险排序能力。

构建完这些花哨的生存模型后,我们该如何评价它们的好坏呢?传统的准确率由于“删失”特征的存在而彻底失效,于是统计学家祭出了专属的裁判指标:C-index(一致性指数,Concordance Index)。 C-index 比较可比样本对的风险排序;0.5 附近表示排序能力接近随机。为避免训练内评价的乐观偏差,下面按 landmark 日期划分较早训练队列与较晚测试队列,并只在后者报告结果。该指标评价的是当前事件定义下的排序,不代表概率校准、财务困境预测或因果效应。

表 11.7: 模型 C-index 评估
from lifelines.utils import concordance_index  # 导入一致性指数用于模型评估

# 先锁定至少具有 3 年潜在随访的成熟队列,再按 landmark 日期划分早期训练与晚期测试
landmark_dates = pd.to_datetime(company_survival_data.loc[cox_model_data.index, 'Landmark_Date'])
study_end_date = observation_end_date  # 使用本章预先锁定的研究截止日期,不从结果数据反推
maturity_cutoff = study_end_date - pd.DateOffset(years=3)
mature_mask = landmark_dates <= maturity_cutoff
mature_cox_data = cox_model_data.loc[mature_mask]
mature_landmark_dates = landmark_dates.loc[mature_mask]
cohort_boundary = mature_landmark_dates.sort_values().iloc[int(len(mature_landmark_dates) * 0.8)]
training_mask = mature_landmark_dates < cohort_boundary
cox_training_data = mature_cox_data.loc[training_mask]
cox_test_data = mature_cox_data.loc[~training_mask]

evaluation_cox_model = CoxPHFitter(penalizer=0.1)
evaluation_cox_model.fit(cox_training_data, duration_col='duration', event_col='event')
predicted_risk_scores = evaluation_cox_model.predict_partial_hazard(cox_test_data)
print(f'队列切分日期: {cohort_boundary.date()}')
print(f'测试队列样本数: {len(cox_test_data)}; 事件数: {int(cox_test_data.event.sum())}')
try:
    cox_concordance_index = concordance_index(
        cox_test_data['duration'], -predicted_risk_scores, cox_test_data['event'])
    print(f'晚期测试队列 C-index: {cox_concordance_index:.4f}')
except ZeroDivisionError:
    cox_concordance_index = np.nan
    print('晚期测试队列 C-index: NA(没有可比较事件对)')

请以代码实际输出为准,并同时报告切分日期与测试队列规模。单次时间队列切分仍可能不稳定;严谨研究还应预先固定特征、事件口径和评价窗口,并用多个后续队列或外部样本复核。

11.9 案例:股价回撤恢复时间分析

在投资风险管理中的最大回撤(Max Drawdown)是一个关键指标。我们不仅关心回撤的深度,还关心回撤恢复时间(Time to Recovery),即股价从跌破高点到重创新高所需的时间。这可以看作是一个生存分析问题。

  • 事件:股价创新高(恢复)。
  • 时间:从上一次创新高到下一次创新高的持续时间。
  • 删失:截至目前仍未创新高。

我们将分析海康威视(002415.XSHE)的历史回撤恢复情况。

在这个章节的最后,我们要打破思维的桎梏:谁说生存分析只能用来研究“死亡”与“退市”?在充满奇迹的交易桌上,一次灾难性的暴跌也可以被视为一次“死亡”,而股价重新跌回并突破前高,则是一场漫长的“复活”。 下面的代码,将生存分析的矛头直指向了一个令无数股民夜不能寐的金融难题:“买在最高点被套牢,到底需要熬多少天才解套?”(即最大回撤恢复时间分析): 我们以智能安防龙头——海康威视的全历史日线数据为手术刀。代码首先计算了股价的历史阻力位与实时回撤幅度,然后用近乎苛刻的条件找出了每一次超了5 天的显著“套牢解套”周期,这构成了我们的生存(持续套牢)时间,而股价重新突破前高就是那激动人心的最终“事件”。同时,对于当前依然被套牢在山顶无法解脱的区间,代码极其专业地将其标记为了“右删失”。 最后输出的这张生存曲线以及底部的硬核文字报告,不再是一份冰冷的医学死亡率,而是一份实打实的交易心理学防线评估:它不仅告诉你海康威视中位数的解套天数,更无情地向你宣告了如果你在错误的时间冲进去,在接下来的半年(180天)内还要继续承受无期徒刑般煎熬的概率。

# 在真实 HDF 存储层直接选择海康威视,避免为单股恢复期示例载入全市场行情
drawdown_price_path = Path(os.environ['BOOK_DATA_DIR']).expanduser().resolve() / 'stock' / 'stock_price_post_adjusted.h5'  # 独立定位后复权股价文件
haikang_data = pd.read_hdf(  # 选择性读取海康威视的全部真实历史行情
    drawdown_price_path,  # 指定后复权股价 HDF 文件
    where="order_book_id='002415.XSHE'",  # 在存储层限定唯一教学标的
    columns=['date', 'order_book_id', 'close'],  # 只读取恢复时间所需字段
).reset_index()  # 恢复日期与证券索引列并重建唯一行索引
haikang_data = haikang_data.sort_values('date').set_index('date')  # 按日期排序并设为索引
closing_prices = haikang_data['close']  # 提取收盘价序列

# 计算滚动最高价与回撤幅度
running_maximum_prices = closing_prices.cummax()  # 计算历史累计最高价
drawdown_series = (closing_prices - running_maximum_prices) / running_maximum_prices  # 计算回撤百分比

# 找出所有创新高的日期
is_new_high = (closing_prices >= running_maximum_prices)  # 标记当日是否创历史新高
new_high_dates = closing_prices.index[is_new_high]  # 提取所有创新高的日期

基于上述回撤计算结果,我们接下来识别每一次显著的回撤-恢复周期(持续时间超过5个交易日的回撤),并计算其持续天数:

# 遍历所有创新高的相邻日期对,计算回撤恢复时间
recovery_durations = []  # 存储每次回撤的持续天数
recovery_events = []  # 存储事件标记(1=已恢复, 0=删失/未恢复)

if len(new_high_dates) > 1:  # 至少需要两个新高日期才能计算间隔
    for i in range(len(new_high_dates) - 1):  # 遍历每对相邻新高日期
        start_date = new_high_dates[i]  # 本次创新高日期
        end_date = new_high_dates[i+1]  # 下次创新高日期
        duration_days = (end_date - start_date).days  # 计算间隔天数
        if duration_days > 5:  # 仅保留超过5天的显著回撤
            recovery_durations.append(duration_days)  # 记录回撤持续天数
            recovery_events.append(1)  # 标记为已恢复事件
    # 处理最后一段(当前是否仍在回撤中)
    last_high_date = new_high_dates[-1]  # 最后一次创新高日期
    last_recorded_date = closing_prices.index[-1]  # 数据最后一天
    if last_recorded_date > last_high_date:  # 若最后一天不是新高,说明当前仍在回撤中
        duration_days = (last_recorded_date - last_high_date).days  # 计算当前回撤持续天数
        if duration_days > 5:  # 仅保留显著回撤
            recovery_durations.append(duration_days)  # 记录当前回撤天数
            recovery_events.append(0)  # 标记为删失(尚未恢复)

drawdown_recovery_data = pd.DataFrame({'duration': recovery_durations, 'event': recovery_events})  # 构建回撤恢复数据框
print(f'识别到 {len(drawdown_recovery_data)} 次显著回撤周期')  # 输出回撤周期总数
recovered_mean_days = drawdown_recovery_data.loc[drawdown_recovery_data['event']==1, 'duration'].mean()  # 计算已恢复周期的平均天数
print(f'平均恢复时间: {recovered_mean_days:.1f} 天')  # 输出平均恢复天数
识别到 38 次显著回撤周期
平均恢复时间: 100.7 天

回撤周期数量与恢复时间以当前数据快照的代码输出为准。这里的周期由“创新高—再次创新高”规则构造,彼此并非天然独立,且价格路径、复权方式和样本终点都会改变结果;因此它是描述性练习,不是对未来解套时间的承诺。

我们已识别出所有显著回撤周期的生存数据。下面使用 Kaplan-Meier 方法估计回撤恢复时间的生存曲线,并标注中位恢复时间:

# 拟合 Kaplan-Meier 生存曲线
kaplan_meier_fitter = KaplanMeierFitter()  # 实例化KM估计器
kaplan_meier_fitter.fit(drawdown_recovery_data['duration'], drawdown_recovery_data['event'], label='回撤恢复')  # 拟合模型

fig, ax = plt.subplots(figsize=(10, 6))  # 创建10×6英寸画布
kaplan_meier_fitter.plot(ax=ax)  # 绘制KM生存曲线

# 添加中位恢复时间参考线
median_recovery_time = kaplan_meier_fitter.median_survival_time_  # 获取中位恢复时间
plt.axvline(median_recovery_time, color='r', linestyle='--', label=f'中位恢复时间: {median_recovery_time:.0f}天')  # 绘制竖直参考线

ax.set_title('股价回撤恢复时间生存曲线 (Time to Recovery)', fontsize=14, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置图标题
ax.set_xlabel('天数', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 x 轴标签
ax.set_ylabel('未恢复概率', fontsize=12, fontproperties='Source Han Serif SC')  # 使用已安装中文字体设置 y 轴标签
ax.grid(True, alpha=0.3)  # 添加网格线
plt.legend(prop={'family': 'Source Han Serif SC'})  # 使用已安装中文字体设置图例
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图形

# 预测特定时间点的未恢复概率
prob_not_recovered_30d = kaplan_meier_fitter.predict(30)  # 30天内未恢复的概率
prob_not_recovered_180d = kaplan_meier_fitter.predict(180)  # 180天内未恢复的概率
print(f'30天内无法恢复的概率: {prob_not_recovered_30d:.1%}')  # 输出30天未恢复概率
print(f'半年(180天)内无法恢复的概率: {prob_not_recovered_180d:.1%}')  # 输出180天未恢复概率
Kaplan-Meier 阶梯曲线显示回撤尚未恢复的概率随天数下降,并标出中位恢复时间。
图 11.8: 海康威视股价回撤恢复时间的生存分析
30天内无法恢复的概率: 44.7%
半年(180天)内无法恢复的概率: 13.2%

图 11.8 描述当前样本与所选回撤定义下“尚未恢复”的经验分布。30 天、180 天及中位恢复时间应直接从本次拟合对象读取;删失、重叠回撤和单一标的限制使其不能直接转化为止损规则、系统性风险归因或基本面判断。

11.10 回撤恢复支线审计任务

本任务为 60 分钟、10 分的可独立描述性课堂支线。它使用与正文海康威视案例不同的 TEACHING_SECURITY_B、8 个营业日纳入阈值和 20/60 个营业日评价窗口;fixture 是无随机数的人工教学轨迹,不代表任何真实证券、交易所行情或投资结果。

表 11.8: Chapter 11 支线 fixture 的版本、模式与资源合同
合同项 冻结值
contract_id ch11-recovery-audit-v1
性质 人工构造的教学 fixture;不得作真实数据主张
合同冻结日 2026-08-12(Asia/Shanghai)
schema fixture_version:str, security_id:str, date:datetime64[ns], adjusted_close:float64, source_status:str
行数与窗口 120 个营业日;回撤周期至少 8 个营业日;报告 20/60 营业日未恢复概率
SHA-256 e2eb5e9e24af7e0d6327efdcb610759147b532091dc73f4679bc6ad564227187
资源上限 普通 CPU 单核、内存 128 MB、总运行 60 秒;仅允许一条 KM 曲线,不拟合预测模型
列表 11.4: Chapter 11 支线的独立 fixture setup 与内容哈希
import hashlib  # 为人工价格轨迹建立内容寻址证据
import numpy as np  # 构造无随机性的周期轨迹
import pandas as pd  # 建立统一日期和数值schema
ch11_fixture_dates = pd.date_range('2024-01-02', periods=120, freq='B')  # 冻结120个营业日
ch11_fixture_positions = np.arange(120)  # 建立确定性位置索引
ch11_fixture_prices = 100 + 0.08 * ch11_fixture_positions + 2 * np.sin(np.pi * ch11_fixture_positions / 10)  # 产生可恢复的教学回撤
ch11_audit_fixture = pd.DataFrame({'fixture_version': 'ch11-recovery-audit-v1', 'security_id': 'TEACHING_SECURITY_B', 'date': ch11_fixture_dates, 'adjusted_close': ch11_fixture_prices, 'source_status': 'synthetic_teaching_fixture'})  # 按合同schema构造fixture
ch11_fixture_bytes = ch11_audit_fixture.to_csv(index=False, date_format='%Y-%m-%d', float_format='%.6f', lineterminator='\n').encode('utf-8')  # 执行规范序列化
ch11_fixture_sha256 = hashlib.sha256(ch11_fixture_bytes).hexdigest()  # 计算实际内容哈希
assert ch11_fixture_sha256 == 'e2eb5e9e24af7e0d6327efdcb610759147b532091dc73f4679bc6ad564227187'  # 内容变化即结构化停止

11.10.1 consumer 与独立证据

独立运行时先执行 列表 11.4。consumer 只能接收 ch11_audit_fixture,先断言版本、列集合、唯一 (security_id,date) 键、日期严格递增、价格为正和 source_status == 'synthetic_teaching_fixture';随后用 adjusted_close.cummax() 找创新高,以相邻创新高间隔至少 8 个营业日定义已恢复周期,末端未恢复区间编码为右删失。不得读取 closing_pricesdrawdown_recovery_data 或正文的 kaplan_meier_fitter

0—15 分钟提交 ch11_event_contract,登记时间原点、恢复事件、行政删失、8 日阈值和 20/60 日窗口;15—30 分钟提交 ch11_period_audit,逐周期报告稳定 period_id、起止日、duration_business_daysevent 与每个事件时点风险集;30—45 分钟提交 ch11_km_evidence,包括周期数、事件数、删失数、KM 表及可定义性;45—60 分钟提交 ch11_window_summarych11_limitations,报告中位数、20/60 日未恢复概率,以及相邻周期依赖、单一人工轨迹和非交易结论。

评分为:事件/删失合同 3 分,逐周期与风险集 3 分,KM 独立输出 2 分,窗口摘要、限制与非真实数据边界 2 分。缺版本/哈希/schema 守卫或事件合同,结构化失败且最高 4 分;风险集不可复算,第二项为 0 分;把曲线写成真实证券发现、未来承诺或交易建议,结论项为 0 分。此支线的证据必须由 ch11-recovery-audit-v1 生成,正文 图 11.8 的图形和数值不能充当提交物。

11.11 本章小结

本章我们系统学习了生存分析的核心方法:

  1. Kaplan-Meier 估计量:非参数估计生存曲线
  2. 对数秩检验:比较两组或多组生存曲线
  3. Cox 比例风险模型:半参数回归模型,评估协变量对生存的影响
  4. 模型评估:使用C-index 评估预测性能
  5. 正则化:对 Cox 模型应用 Lasso/Ridge 正则化

生存分析的应用建议

  1. 医学研究:患者生存时间、无病生存时间
  2. 商业分析:客户流失时间、产品使用寿命
  3. 工程质量:设备故障时间、系统可靠性
  4. 人力资源:员工离职时间、晋升时间

关键假设

  • 独立删失(数据收集时要仔细设计)
  • 比例风险(Cox 模型需要检验)

软件工具

  • Python: lifelines, scikit-survival
  • R: survival, survminer

11.12 理论来源与前沿

生存分析最早在可靠性工程医学随访研究中系统化发展:前者关注设备故障时间,后者关注患者生存时间。Kaplan-Meier 的乘积极限估计提供了在删失存在时对生存函数\(S(t)\) 的非参数估计;Cox 模型则用偏似然把协变量效应与基线风险函数 \(h_0(t)\) 分离,使得回归分析在不指定\(h_0(t)\) 的情况下仍然可行。

近十年的研究前沿主要集中在三类问题上:

  1. 高维协变量与正则化:当特征数远大于样本量时,需要用 Lasso/Elastic Net 等正则化构建稀疏Cox 模型,并配套稳定的变量选择与不确定性量化。
  2. 时间变化效应与非比例风险:在真实商业场景(例如客户流失)中,某些协变量的影响往往随时间变化,此时需要扩展Cox 模型(时间交互项、分段比例风险、AFT 模型等),并用Schoenfeld 残差等方法诊断。
  3. 机器学习与因果推断的结合:将生存模型与树模型、Boosting、深度学习以及因果推断框架结合,用于处理复杂非线性、异质性处理效应与动态干预策略。

11.13 练习

11.13.1 概念题

  1. [核心|难度:1|时间:8分钟|分值:4|项目:无] 解释什么是右删失、左删失和区间删失,并各举一个实际例子。

  2. [核心|难度:2|时间:10分钟|分值:5|项目:无] Kaplan-Meier 估计量的核心思想是什么?为什么它比简单地计算生存比例更合理?

  3. [核心|难度:2|时间:10分钟|分值:5|项目:无] Cox 比例风险模型中的“比例风险”假设是什么意思?如何检验这个假设?

  4. [核心|难度:2|时间:10分钟|分值:5|项目:无] 解释风险函数 \(h(t)\) 和生存函数\(S(t)\) 之间的关系。

  5. [拓展|难度:2|时间:10分钟|分值:5|项目:无] 在什么情况下独立删失假设会被违反?这会对分析结果造成什么影响?

11.13.2 应用题

  1. [核心|难度:2|时间:15分钟|分值:8|项目:无] 内嵌小表的 Kaplan–Meier 手算(必做):10 个产品在第 2 月发生 2 个事件,随后 1 个删失;第 5 月发生 1 个事件;第 7 月有 2 个删失;其余产品在第 10 月研究结束时删失。列出每个时点的风险集、事件数和删失数,计算 \(\widehat S(2)\)\(\widehat S(5)\),并说明第 7 月为何不产生生存曲线下降。

    lungcolon 数据可作为选做扩展,但不是完成本题的外部依赖。 ### 应用于(金融场景)

  2. [拓展|难度:3|时间:40分钟|分值:15|项目:无] 使用 financial_statement.h5 数据,选取两个行业,构建“从首份可得财报到首次亏损”的生存数据集。

    要求:

    • 绘制两个行业的Kaplan-Meier 生存曲线。
    • 使用对数秩检验判断两条曲线是否有显著差异。
    • 解释结果:哪个行业的“财务安全期”更长?
  3. [拓展|难度:3|时间:40分钟|分值:15|项目:无] 模拟信贷违约数据: 假设你是一家网贷平台的风控分析师,请模拟一份贷款数据(包含借款人年龄、收入、信用分、借款金额等特征)。

    要求:

    • 设定真实的违约机制(例如低信用分、高杠杆导致风险指数增加)。
    • 添加删失(例如贷款尚未到期)。
    • 使用带Lasso 正则化的 Cox 模型筛选关键风险因子。
  4. [拓展|难度:3|时间:40分钟|分值:15|项目:无] 综合案例:A股新股生存分析 利用 stock_basic_data.h5 中的 listed_datede_listed_date(退市日期)。

    要求:

    • 定义事件:公司退市。
    • 定义删失:公司截至目前仍在交易。
    • 比较不同板块(主板vs 创业板vs 科创板)的退市风险(生存曲线)。
    • 若退市事件过少,应报告估计不稳定并合并教学分组;不得改用缺少生效日期的静态 ST 状态替换事件。

11.13.3 理论题

  1. [核心|难度:3|时间:20分钟|分值:10|项目:无] 给定三个无并列事件时点及其风险集,写出 Cox 偏似然的一般乘积、score 与观察信息的作用;再说明为什么基线风险相消、风险比不等于绝对风险,以及时间外 C-index 在零事件或全删失测试集为何不可定义。

  2. [拓展|难度:3|时间:25分钟|分值:12|项目:无] 推导 Cox 模型的偏似然函数,并解释为什么不需要估计基线风险\(h_0(t)\)

  3. [拓展|难度:3|时间:25分钟|分值:12|项目:无] 证明在单变量二元协变量的情况下,Cox 模型的score test 等价于对数秩检验。

11.14 练习参考解答

展开完整解答、评分点与常见失败模式

评分以题面分值为准;概率与风险集手算允许 \(10^{-3}\) 绝对误差。应用题若事件、时间原点、删失或观察终点任一未登记,最高得 50%;若把风险比写成绝对概率、静态状态写成事件日期或忽略零事件组,解释项不得分。

11.14.1 概念题解答

  1. 右删失、左删失与区间删失

    • 右删失 (right censoring):只知道事件时间 \(T\) 大于观测结束时刻 \(C\)。例:研究结束时公司仍未退市。
    • 左删失 (left censoring):只知道事件时间 \(T\) 小于某个观测起点。例:只知道公司在数据收集前已经违约,但不知确切日期。
    • 区间删失 (interval censoring):只知道 \(T\) 落在区间 \((L, R]\)。例:只知道违约发生在两次财报披露之间。
  2. Kaplan-Meier 的核心思想

    将总体生存函数分解为一连串条件生存概率的乘积: \[ S(t) = \prod_{t_j \le t} P(T > t_j \mid T \ge t_j) \] 在每个事件时刻 \(t_j\),用风险集大小 \(n_j\) 与事件数 \(d_j\) 估计条件概率为 \(1 - d_j/n_j\)

  3. 比例风险假设与检验

    Cox 模型假设两组个体的风险比与时间无关。常用检验:Schoenfeld 残差与时间的相关性检验。若存在显著相关性,说明违反比例风险假设,应考虑含时间依赖系数的 Cox 模型。

  4. 风险函数与生存函数

    \[ S(t) = \exp\Big(-\int_0^t h(u)\,du\Big) \]

  5. 独立删失

    若删失机制包含有关 \(T\) 的信息(例如财务状况恶化导致的数据缺失),则不再独立,标准方法失效。

11.14.2 应用题解答

  1. 内嵌小表手算

    第 2 月事件前 \(r_1=10,d_1=2\),所以 \(\widehat S(2)=8/10=0.8\);同月随后 1 个删失,故第 5 月事件前 \(r_2=7,d_2=1\)\(\widehat S(5)=0.8(6/7)\approx0.686\)。第 7 月只有删失,条件生存概率没有事件因子,曲线保持 0.686,但第 7 月后风险集减少 2。第 10 月其余对象行政删失,曲线仍不下降。

  2. 行业首次亏损比较(实现契约)

    在公司内按可靠 info_date 排序,将首次 net_profit < 0 的可得日定义为事件日;从首份可得财报开始计时,研究截止仍未亏损者右删失。下面复用 表 11.3 的冻结切片与哈希,不允许固定日数代理;selected_industries 是事前指定的比较族,不按结果选择。

import os
import hashlib
import pandas as pd
from lifelines import KaplanMeierFitter
from lifelines.statistics import logrank_test
from pathlib import Path  # 使用跨平台路径对象解析显式切片
financial_path = Path(os.environ['BOOK_M11_FINANCIAL_SLICE']).expanduser().resolve()  # 读取同一point-in-time切片接口
expected_hash = os.environ['BOOK_M11_FINANCIAL_SLICE_SHA256'].strip().lower()  # 读取冻结内容哈希
if hashlib.sha256(financial_path.read_bytes()).hexdigest() != expected_hash:  # 读取HDF前校验版本内容
    raise RuntimeError({'status': 'stopped', 'reason': 'input_sha256_mismatch'})  # 哈希不符即停止
financial = pd.read_hdf(financial_path).copy()  # 读取已验证的小切片
DATA_DIR = Path(os.environ['BOOK_DATA_DIR']).expanduser().resolve()  # 基本信息仍由显式本地根提供
stock_basic_data = pd.read_hdf(DATA_DIR / 'stock' / 'stock_basic_data.h5').copy()
if 'info_date' not in financial.columns:  # 禁止缺失披露日时继续
    raise RuntimeError({'status': 'stopped', 'reason': 'missing_info_date'})  # 缺时点即停止
financial['available_date'] = pd.to_datetime(financial['info_date'], errors='coerce')  # 以可靠披露或修订日为可得日
if financial['available_date'].isna().any():  # 禁止不可解析时点进入风险集
    raise RuntimeError({'status': 'stopped', 'reason': 'invalid_info_date'})  # 无法解析即停止
observation_end_date = financial['available_date'].max().normalize()
financial = financial.sort_values(['order_book_id', 'available_date'])
first_available = financial.groupby('order_book_id')['available_date'].min().rename('origin')
first_loss = (financial[financial['net_profit'].lt(0)]
              .groupby('order_book_id')['available_date'].min().rename('event_date'))
first_loss_survival = pd.concat([first_available, first_loss], axis=1).reset_index()
first_loss_survival['event'] = first_loss_survival['event_date'].notna().astype(int)
first_loss_survival['end'] = first_loss_survival['event_date'].fillna(observation_end_date)
first_loss_survival['duration'] = (
    (first_loss_survival['end'] - first_loss_survival['origin']).dt.days / 365.25
)
industry_map = stock_basic_data[['order_book_id', 'citics_2019_l1_name']].drop_duplicates('order_book_id')
first_loss_survival = first_loss_survival.merge(industry_map, on='order_book_id').query('duration >= 0')
selected_industries = ['零售', '建筑']
groups = [first_loss_survival[first_loss_survival['citics_2019_l1_name'].eq(name)]
          for name in selected_industries]
if any(group.empty for group in groups):
    raise ValueError('指定行业在当前快照中无样本;请在运行前登记两个实际行业名称。')
for name, group in zip(selected_industries, groups):
    KaplanMeierFitter().fit(group['duration'], group['event'], label=name).plot_survival_function()
图 11.9

图 11.9 的行业曲线必须从现场输出读取;它依赖当前冻结切片与 info_date 版本,不能证明行业因果差异。检验输出单独见 表 11.9

表 11.9: 习题7:两个预注册行业的 log-rank 检验
comparison = logrank_test(groups[0]['duration'], groups[1]['duration'], event_observed_A=groups[0]['event'], event_observed_B=groups[1]['event'])  # 使用与图形相同的两个预注册组
comparison.summary  # 以独立表格报告统计量和p值

补充审计:旧的任何原因退市行业模板(不构成第 7 或第 9 题答案)

#| output: false
import pandas as pd  # 导入pandas用于数据处理
from lifelines import KaplanMeierFitter  # 导入Kaplan-Meier生存函数估计器
from lifelines.statistics import logrank_test  # 导入对数秩检验用于比较生存曲线差异
import matplotlib.pyplot as plt  # 导入matplotlib用于数据可视化
import os  # 导入操作系统模块用于跨平台路径处理

# ---------- 1. 数据准备(复用正文 Cell 1 已生成的 company_survival_data) ----------
# company_survival_data 包含 'industry'、'duration'、'event' 列
# 'industry' 列的值来自 citics_2019_l1_name(中信一级行业分类)

# 动态获取样本量前两大的行业,而非硬编码行业名称
industry_value_counts = company_survival_data['industry'].value_counts()  # 统计各行业样本数
industry_first_name = industry_value_counts.index[0]  # 样本量最多的行业
industry_second_name = industry_value_counts.index[1]  # 样本量第二多的行业
print(f'选取的两个行业: {industry_first_name}, {industry_second_name}')  # 打印实际行业名称

# ---------- 2. 构建行业筛选掩码 ----------
mask_first_industry = company_survival_data['industry'] == industry_first_name  # 第一个行业的布尔掩码
mask_second_industry = company_survival_data['industry'] == industry_second_name  # 第二个行业的布尔掩码

# ---------- 3. 绘制两个行业的 KM 生存曲线并进行 Log-Rank 检验 ----------
fig_exercise, ax_exercise = plt.subplots(figsize=(10, 6))  # 创建画布

kaplan_meier_fitter_exercise = KaplanMeierFitter()  # 实例化 KM 估计器

# 拟合并绘制第一个行业
kaplan_meier_fitter_exercise.fit(  # 训练/拟合模型
    company_survival_data.loc[mask_first_industry, 'duration'],  # 第一行业的存续时间
    company_survival_data.loc[mask_first_industry, 'event'],  # 第一行业的事件标识
    label=industry_first_name  # 图例标签
)  # 完成构建
kaplan_meier_fitter_exercise.plot(ax=ax_exercise)  # 绘制到画布上

完成第一个行业的 Kaplan-Meier 曲线绘制后,接下来对第二个行业执行相同的拟合与绘制流程,并通过 Log-Rank 检验比较两个行业生存曲线之间的差异是否具有统计学显著性。


# 拟合并绘制第二个行业
kaplan_meier_fitter_exercise.fit(  # 训练/拟合模型
    company_survival_data.loc[mask_second_industry, 'duration'],  # 第二行业的存续时间
    company_survival_data.loc[mask_second_industry, 'event'],  # 第二行业的事件标识
    label=industry_second_name  # 图例标签
)  # 完成构建
kaplan_meier_fitter_exercise.plot(ax=ax_exercise)  # 绘制到画布上

# 进行 Log-Rank 检验,比较两个行业生存曲线的差异是否显著
rank_test_result = logrank_test(  # 提取测试集数据
    company_survival_data.loc[mask_first_industry, 'duration'],  # 第一行业持续时间
    company_survival_data.loc[mask_second_industry, 'duration'],  # 第二行业持续时间
    company_survival_data.loc[mask_first_industry, 'event'],  # 第一行业事件
    company_survival_data.loc[mask_second_industry, 'event']  # 第二行业事件
)  # 完成构建

ax_exercise.set_title(f'行业生存曲线比较: {industry_first_name} vs {industry_second_name}')  # 设置标题
ax_exercise.set_xlabel('存续年数')  # 设置 x 轴标签
ax_exercise.set_ylabel('生存概率')  # 设置 y 轴标签
plt.tight_layout()  # 自动调整布局
plt.show()  # 展示图形

print(f'Log-Rank 检验 p 值: {rank_test_result.p_value:.4f}')  # 输出检验 p 值

上述代码动态选取两个行业并输出 Log-Rank 检验。是否存在统计证据必须按现场 p 值和预先阈值判断;它是任何原因退市的辅助练习,不是“首次亏损”答案。

  1. 信贷违约模拟与 Lasso-Cox(完整实现)

本题明确允许模拟,因为目标是检查已知真机制下的变量筛选。设信用分与收入降低风险,杠杆提高风险,年龄与三个噪声变量的真系数为零;独立行政删失并不依赖事件时间。

列表 11.5: 习题8:生成带已知风险机制与删失的信贷数据
import numpy as np  # 使用可复现随机数生成已知真值样本
import pandas as pd  # 组织 Cox 模型所需的列式数据
simulation_rng = np.random.default_rng(20260812)  # 固定随机种子以便逐次复算
borrower_count = 800  # 设置足以支持开发期与验证期的样本数
credit_score_z = simulation_rng.normal(size=borrower_count)  # 标准化信用分
income_z = simulation_rng.normal(size=borrower_count)  # 标准化收入
leverage_z = simulation_rng.normal(size=borrower_count)  # 标准化杠杆
age_z = simulation_rng.normal(size=borrower_count)  # 生成真系数为零的年龄
noise_features = simulation_rng.normal(size=(borrower_count, 3))  # 增加三个无关特征
linear_risk = -0.8 * credit_score_z - 0.4 * income_z + 0.7 * leverage_z  # 写明真风险机制
event_month = simulation_rng.exponential(scale=1 / (0.03 * np.exp(linear_risk)))  # 生成违约时间
censor_month = simulation_rng.uniform(6, 60, size=borrower_count)  # 生成独立行政删失时间
observed_month = np.minimum(event_month, censor_month)  # 构造实际观测时长
default_event = (event_month <= censor_month).astype(int)  # 标记是否观察到违约
credit_survival = pd.DataFrame({'duration': observed_month, 'event': default_event, 'credit_score': credit_score_z, 'income': income_z, 'leverage': leverage_z, 'age': age_z, 'noise_1': noise_features[:, 0], 'noise_2': noise_features[:, 1], 'noise_3': noise_features[:, 2]})  # 汇总建模表
print({'n': len(credit_survival), 'events': int(default_event.sum()), 'censor_rate': 1 - default_event.mean()})  # 审计样本与删失
{'n': 800, 'events': 458, 'censor_rate': 0.4275}

按样本编号锁定前 70% 为开发集、后 30% 为验证集。这里只在验证集选择惩罚强度;验证分数相同到 \(10^{-4}\) 时选择更大的惩罚以偏向稀疏。

列表 11.6: 习题8:验证期选择 Lasso 惩罚并报告筛选稳定性
from lifelines import CoxPHFitter  # 拟合带 L1 惩罚的 Cox 模型
development_end = int(0.7 * len(credit_survival))  # 冻结开发与验证边界
credit_development = credit_survival.iloc[:development_end].copy()  # 提取开发样本
credit_validation = credit_survival.iloc[development_end:].copy()  # 提取锁定验证样本
penalty_grid = [0.001, 0.01, 0.03, 0.1, 0.3]  # 事前登记 L1 惩罚候选值
validation_scores = {}  # 保存每个候选值的验证 C-index
for penalty_value in penalty_grid:  # 仅遍历预注册候选网格
    candidate_cox = CoxPHFitter(penalizer=penalty_value, l1_ratio=1.0)  # 建立纯 Lasso-Cox
    candidate_cox.fit(credit_development, duration_col='duration', event_col='event')  # 只在开发样本拟合
    validation_scores[penalty_value] = candidate_cox.score(credit_validation, scoring_method='concordance_index')  # 锁定验证评价
best_score = max(validation_scores.values())  # 找出现场最佳验证分数
eligible_penalties = [value for value, score in validation_scores.items() if best_score - score <= 1e-4]  # 应用容差规则
selected_penalty = max(eligible_penalties)  # 容差内选更稀疏模型
selected_cox = CoxPHFitter(penalizer=selected_penalty, l1_ratio=1.0)  # 固定最终 Lasso-Cox
selected_cox.fit(credit_development, duration_col='duration', event_col='event')  # 保持验证集未参与拟合
coefficient_audit = selected_cox.params_.rename('coefficient').to_frame()  # 汇总各变量系数
coefficient_audit['selected'] = coefficient_audit['coefficient'].abs() > 1e-4  # 用预设容差标记非零项
print({'penalty': selected_penalty, 'validation_c_index': validation_scores[selected_penalty]})  # 报告调参与评价
print(coefficient_audit)  # 显示关键因子与噪声因子的筛选结果
{'penalty': 0.001, 'validation_c_index': 0.7346171434378335}
              coefficient  selected
covariate                          
credit_score    -0.726608      True
income          -0.347337      True
leverage         0.726567      True
age             -0.006131      True
noise_1         -0.047161      True
noise_2         -0.071014      True
noise_3          0.014402      True

评分证据包括真系数表、删失率、惩罚网格、验证 C-index 和全部系数。允许估计误差,但应检查三项:信用分/收入方向为负、杠杆方向为正、无关项整体收缩更强。若验证集事件少于 20、只含一个事件类别或模型不收敛,应扩大预注册样本而不是报告筛选成功。

  1. A股板块任何原因退市风险
# ---------- A股不同上市板块的退市风险比较:数据准备 ----------
import os  # 导入操作系统模块用于跨平台路径处理
import pandas as pd  # 导入pandas用于数据处理
from lifelines import KaplanMeierFitter  # 导入KM估计器
from lifelines.statistics import logrank_test  # 导入对数秩检验
import matplotlib.pyplot as plt  # 导入matplotlib用于可视化

from pathlib import Path  # 使用跨平台路径对象解析显式数据根
DATA_DIR = Path(os.environ['BOOK_DATA_DIR']).expanduser().resolve()  # 从必需环境变量取得数据根
if not DATA_DIR.is_dir():  # 在读取前验证数据根
    raise FileNotFoundError(f'BOOK_DATA_DIR 不存在或不是目录: {DATA_DIR}')  # 失败即停止,不回退固定路径
stock_basic_data_ex = pd.read_hdf(DATA_DIR / 'stock' / 'stock_basic_data.h5')  # 加载上市公司基本信息

# 将上市和退市日期转换为 datetime 格式
stock_basic_data_ex['listed_date'] = pd.to_datetime(stock_basic_data_ex['listed_date'], errors='coerce')  # 上市日期
stock_basic_data_ex['de_listed_date'] = pd.to_datetime(stock_basic_data_ex['de_listed_date'], errors='coerce')  # 退市日期
stock_basic_data_ex = stock_basic_data_ex.dropna(subset=['listed_date'])  # 去除缺少上市日期的记录

observation_end = pd.to_datetime('2024-01-01')  # 研究观测截止日期

在加载基本信息数据后,下一步根据股票代码后缀规则将各股票归入对应的上市板块(主板、创业板、科创板),以便后续按板块维度进行生存分析比较。

# 根据股票代码后缀判断上市板块
def classify_board(order_book_id):  # 定义函数classify_board
    """根据股票代码判断所属板块"""  # 执行数据处理操作
    code = str(order_book_id).split('.')[0]  # 提取纯数字代码
    if code.startswith('688'):  # 科创板
        return '科创板'  # 返回结果
    elif code.startswith('300'):  # 创业板
        return '创业板'  # 返回结果
    elif code.startswith('60'):  # 上交所主板
        return '主板(沪)'  # 返回结果
    elif code.startswith('00'):  # 深交所主板
        return '主板(深)'  # 返回结果
    else:  # 默认分支
        return '其他'  # 返回结果

stock_basic_data_ex['board_type'] = stock_basic_data_ex['order_book_id'].apply(classify_board)  # 添加板块分类列
stock_basic_data_ex = stock_basic_data_ex[stock_basic_data_ex['board_type'] != '其他']  # 排除无法分类的股票

基于上述板块分类结果,接下来构建各股票的生存数据并绘制 Kaplan-Meier 生存曲线:

# ---------- 构建生存数据并绘制各板块KM曲线 ----------
board_survival_records = []  # 存储各股票的生存记录
for _, row in stock_basic_data_ex.iterrows():  # 遍历迭代
    listed = row['listed_date']  # 上市日期
    delisted = row['de_listed_date']  # 退市日期
    if pd.isna(delisted) or delisted > observation_end:  # 缺失或截止日后退市才删失
        dur = (observation_end - listed).days / 365.25  # 存续年数
        evt = 0  # 事件标记为删失
    else:  # 截止日当日或之前退市均为事件
        dur = (delisted - listed).days / 365.25  # 存续年数
        evt = 1  # 事件标记为退市
    if dur > 0:  # 排除无效记录
        # 将计算结果追加到列表
        board_survival_records.append({'board': row['board_type'], 'duration': dur, 'event': evt})

board_survival_df = pd.DataFrame(board_survival_records)  # 转换为DataFrame

完成生存数据构建和 KM 估计器初始化后,下面遍历各板块分别拟合 Kaplan-Meier 模型并在同一画布上绘制生存曲线,从而直观比较不同板块上市公司的退市风险差异。

fig_board, ax_board = plt.subplots(figsize=(10, 6))  # 为板块曲线创建独立画布
km_board = KaplanMeierFitter()  # 实例化板块曲线共用的KM估计器
for board_name in ['主板(沪)', '主板(深)', '创业板', '科创板']:  # 遍历各板块
    mask_board = board_survival_df['board'] == board_name  # 筛选该板块数据
    if mask_board.sum() < 2:  # 跳过样本量不足的板块
        continue  # 跳过本次循环
    km_board.fit(  # 拟合 KM 模型
        board_survival_df.loc[mask_board, 'duration'],  # 存续时间
        board_survival_df.loc[mask_board, 'event'],  # 事件标识
        label=f'{board_name} (n={mask_board.sum()})'  # 图例含样本量
    )  # 完成构建
    km_board.plot(ax=ax_board)  # 绘制生存曲线

ax_board.set_title('A股不同板块上市公司生存曲线比较')  # 设置标题
ax_board.set_xlabel('存续年数')  # x 轴标签
ax_board.set_ylabel('生存概率')  # y 轴标签
plt.tight_layout()  # 自动调整布局
plt.show()  # 展示图形
主板、创业板与科创板的阶梯生存曲线在同一存续年数横轴上比较,图例显示各组样本数。
图 11.10: A股不同板块任何原因退市的 Kaplan–Meier 生存曲线

图 11.10 只展示当前快照下任何原因退市的描述性曲线;板块样本量与事件率另由 表 11.10 报告,避免把表格输出混入图形对象。

表 11.10: 习题9:各板块样本数、任何原因退市数与事件比例
board_summary = board_survival_df.groupby('board').agg(  # 按板块汇总
    total_count=('event', 'size'),  # 总样本量
    delisted_count=('event', 'sum'),  # 退市数量
    delist_rate=('event', 'mean')  # 退市率
).round(4)  # 保留四位小数
print('\n各板块退市率汇总:')  # 输出表头
print(board_summary)  # 输出汇总表

各板块退市率汇总:
       total_count  delisted_count  delist_rate
board                                          
主板(沪)         1794             102       0.0569
主板(深)         1611             106       0.0658
创业板            977              26       0.0266
科创板            567               2       0.0035

上述代码根据股票代码将 A 股公司分为主板(沪、深)、创业板和科创板四个板块。解释 表 11.10 时必须同时报告事件数;静态退市比例不等于生存概率,也不能替代 图 11.10 的风险集校正。

11.14.3 理论题解答

  1. Cox 偏似然与时间外评价骨架

无并列事件时,若 \(i(j)\)\(t_j\) 发生事件、风险集为 \(R_j\),则

\[L_p(\beta)=\prod_j\frac{\exp(x_{i(j)}^\top\beta)}{\sum_{k\in R_j}\exp(x_k^\top\beta)}.\]

score 是 \(U(\beta)=\partial\log L_p/\partial\beta\),观察信息是 \(I(\beta)=-\partial^2\log L_p/\partial\beta\partial\beta^\top\);前者给出估计方程,后者刻画局部曲率并用于标准误。条件在“风险集中谁发生事件”后,公共的 \(h_0(t_j)\) 从分子分母相消,因此偏似然识别相对风险参数,但不直接给出绝对生存概率。时间外 C-index 需要可比较对以及结果差异;测试期零事件、全删失或没有可比较对时必须报告不可定义,而不能填成 0.5。

  1. Cox 偏似然的连续推导

\(h_i(t)=h_0(t)\exp(x_i^\top\beta)\),在无并列事件的 \(t_j\),已知风险集 \(R_j\) 中恰有一人发生事件,则事件人为 \(i(j)\) 的条件概率为

\[ \frac{h_{i(j)}(t_j)}{\sum_{k\in R_j}h_k(t_j)} =\frac{h_0(t_j)e^{x_{i(j)}^\top\beta}} {\sum_{k\in R_j}h_0(t_j)e^{x_k^\top\beta}} =\frac{e^{x_{i(j)}^\top\beta}}{\sum_{k\in R_j}e^{x_k^\top\beta}}. \]

基线风险在条件概率分子分母中相消。跨事件时点相乘得偏似然

\[ L_p(\beta)=\prod_{j=1}^{D}\frac{e^{x_{i(j)}^\top\beta}}{\sum_{k\in R_j}e^{x_k^\top\beta}}, \]

其对数为 \(\ell_p(\beta)=\sum_j[x_{i(j)}^\top\beta-\log\sum_{k\in R_j}e^{x_k^\top\beta}]\)。因此可以先估计相对风险参数 \(\beta\);若要绝对生存概率,仍需另估累计基线风险。

  1. 二元协变量下 score test 与 log-rank 的等价

\(x_i\in\{0,1\}\)。对上式求导:

\[ U(\beta)=\sum_j\left[x_{i(j)}- \frac{\sum_{k\in R_j}x_ke^{\beta x_k}}{\sum_{k\in R_j}e^{\beta x_k}}\right]. \]

\(H_0:\beta=0\) 下,若 \(Y_{1j}=\sum_{k\in R_j}x_k\)\(Y_j=|R_j|\),且无并列事件,则

\[ U(0)=\sum_j\left[d_{1j}-\frac{Y_{1j}}{Y_j}\right] =\sum_j(O_{1j}-E_{1j}), \]

正是 log-rank 分子。负二阶导给出信息量;允许每时点 \(d_j\) 个事件时,其超几何方差为

\[ I(0)=\sum_j\frac{Y_{1j}Y_{0j}d_j(Y_j-d_j)}{Y_j^2(Y_j-1)}. \]

所以 Cox score 统计量 \(U(0)^2/I(0)\) 与 log-rank 的卡方统计量相同(相同风险集与并列事件处理下)。若使用不同的 tie 近似、权重或时间变化效应,等价关系不再原样成立。

11.15 章末闭环

逐项目标自检:定义生存时间、事件与删失;从风险集写出 Kaplan–Meier 乘积极限;解释 Cox 偏似然为何消去基线风险;检查比例风险和 landmark 时点;用 C-index 或预注册时间点损失比较模型。事件/删失合同、landmark 先于结果窗口及时间外评价可定义性是 must-pass;三项全过且其余两项至少一项有证据,才继续学习。

禁用情境:事件起点、删失机制或风险集无法审计时,不应报告生存概率;比例风险明显不成立且未建模时间变化效应时,不应把单一 Cox 风险比作为全时期摘要。常见误区是把删失当普通缺失删除,或把风险比直接解释为概率差、因果效应或行动收益。

无提示检索:1)删失时点是否在 Kaplan–Meier 乘积中产生事件因子?2)Cox 偏似然消去了什么、没有估计什么?3)训练和验证的 landmark 为什么必须先于各自结果窗口?

展开检索反馈与学习决策

1)不产生事件因子,但会改变后续风险集人数。2)条件化后基线风险相消,偏似然估计相对风险参数;绝对生存概率仍需累计基线风险。3)否则特征或入组会使用结果发生后的信息。第 1 题错:补修 小节 11.3,把第 7 月删失改到第 4 月并重算两次风险集;第 2 题错:补修 小节 11.4,用两个风险集写出分子分母并指出被约去项;第 3 题错:补修 小节 11.8,为“首次逾期”画 landmark、特征截止和测试窗口,要求严格不交叉。三项异形复测均通过才从 小节 11.11 重入。

陌生迁移:为一个贷款提前还款、供应链恢复或客户流失任务,提交时间起点、事件、删失、landmark、风险集和失败条件。下一章转向没有响应变量时的结构探索。

Cox, D. R. 1972年. 《Regression Models and Life-Tables》. Journal of the Royal Statistical Society: Series B (Methodological) 34 (2): 187~220. https://doi.org/10.1111/j.2517-6161.1972.tb00899.x.
Kaplan, E. L., 和 Paul Meier. 1958年. 《Nonparametric Estimation from Incomplete Observations》. Journal of the American Statistical Association 53 (282): 457~81. https://doi.org/10.1080/01621459.1958.10501452.