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

11.1 导读

生存分析研究带有删失的事件发生时间。它既能回答“事件是否发生”,也能回答“何时发生”,适用于客户流失、贷款违约、企业退市和产品寿命等问题。本章从风险集与 Kaplan–Meier 估计出发,再讨论 Cox 比例风险模型、假设诊断和结论边界。

11.2 学习目标

完成本章后,学生应能:

  1. 从事件定义、时间原点、观察终点构造 \((Y,\delta)\),并逐行核对右删失编码,错误率为零。
  2. 对 10 个对象的小表写出每个风险集并手算 Kaplan–Meier 曲线,概率误差不超过 \(10^{-3}\)
  3. 解释 Cox 偏似然为何消去基线风险,并区分风险比与绝对风险;完整 score 与信息量推导列为拓展。
  4. 用 Schoenfeld 残差诊断比例风险;从数据收集设计、失访模式与敏感性分析评估独立删失的可信度,并为两类违背匹配不同处理方法。
  5. 在时间锁定的生存任务中报告删失率、事件数、C-index 与失败条件;将未惩罚与 Lasso-Cox 的稳定性比较列为拓展目标。
先修自检

某公司在研究开始后第 3 年退市,另一公司到第 5 年末仍上市。两者的 \((Y,\delta)\) 分别为 \((3,1)\)\((5,0)\)。若把第二家公司写成事件,可回到 式 11.1式 11.2 复习右删失。

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

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

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

11.3 生存时间与删失时间

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

  • 生存时间(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.3.1 删失机制的重要假设

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

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

违反独立删失假设的例子

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

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

比例风险可以借助残差检验寻找反例,但条件独立删失通常不能仅凭已观测的事件时间与删失指标得到证明。实践中应结合随访制度、退出原因、可观测预测变量和失访模式论证其可信度,并报告对未观测删失机制的敏感性分析。时间交互项与分层 Cox 针对的是非比例风险;若删失可能携带结局信息,则需根据问题考虑逆概率删失加权(IPCW)、事件—删失联合模型或模式混合敏感性分析,不能把两类问题混用。

11.3.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.4 Kaplan-Meier 生存曲线

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

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

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

11.4.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.4.2 中国实际案例:A股上市公司退市生存分析

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

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

下面用 stock_basic_data 演示任何原因退市的 Kaplan–Meier 描述。运行前只需把 BOOK_DATA_DIR 指向本地数据根;Cox 案例另读取同一目录下带 info_date 的财务报表。曲线反映当前样本、截止日和混合退市原因下的存续分布,不能单独解释行业因果机制、财务困境或投资价值。

表 11.1: A股退市生存数据摘要
import pandas as pd  # 以表格结构构造事件时间、删失指标与公司协变量
import numpy as np  # 为风险比变换与数值边界提供数组运算
import matplotlib.pyplot as plt  # 绘制生存曲线及其时间尺度标记
plt.rcParams['font.sans-serif'] = ['Source Han Serif SC']  # 使用系统已安装的思源黑体显示中文,避免字体缺失告警
plt.rcParams['axes.unicode_minus'] = False  # 修正负号显示为方块的问题
from lifelines import KaplanMeierFitter  # 估计右删失样本的乘积极限曲线
import os  # 将在线教材的固定数据根同步给本章后续独立代码块

from pathlib import Path  # 使用跨平台路径对象解析显式数据根
BOOK_DATA_DIR = Path('/home/ubuntu/r2_data_mount/data').resolve()  # 明文定义在线教材的BOOK_DATA_DIR绝对路径
DATA_DIR = BOOK_DATA_DIR  # 保留本章后续代码使用的数据根名称
os.environ['BOOK_DATA_DIR'] = str(BOOK_DATA_DIR)  # 为本章后续恢复期与板块案例登记同一路径
if not DATA_DIR.is_dir():  # 在读取前验证数据根
    raise FileNotFoundError(f'BOOK_DATA_DIR 不存在或不是目录: {DATA_DIR}')  # 失败即停止,不回退固定路径
path_basic = DATA_DIR / 'stock' / 'stock_basic_data.h5'  # 拼接股票基本信息文件路径
assert path_basic.is_file(), f'缺少公司基本信息文件: {path_basic}'  # 在HDF读取前核对目标文件
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_1466247/870818981.py:20: 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在独立会话独立运行
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(f'事件数: {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%
事件数: 237

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

kaplan_meier_fitter = KaplanMeierFitter()  # 实例化KM生存函数估计器
kaplan_meier_fitter.fit(  # 拟合整体生存曲线
    company_survival_data['duration'],  # 存续时间
    company_survival_data['event'],  # 事件标记
    label='A股整体生存率 (未退市)'  # 图例标签
)  # 完成总体KM拟合并保留删失信息
fig, ax = plt.subplots(figsize=(10, 6))  # 为总体退市存续曲线保留单一坐标轴
kaplan_meier_fitter.plot(ax=ax)  # 展示风险集递减形成的阶梯生存估计
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 的当次运行输出读取曲线水平与置信区间。只能将它解释为本数据覆盖、观察截止日与“任何原因退市”定义下的描述性存续分布。曲线不能自动识别制度执行、“壳价值”或退市门槛的因果作用;尾部置信区间变宽则提醒风险集中样本减少。

整体生存曲线描述当前数据口径下的上市存续分布。接下来按行业分解,比较样本量前五行业的退市风险差异。

fig, ax = plt.subplots(figsize=(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  # 使用行业名称作为图例
        )  # 完成当前行业的KM拟合
        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.4.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_second = top_industry_categories[1]  # 固定样本量第二的行业作为第二比较组

print(f"比较行业: {industry_first} vs {industry_second}")  # 明示由当前快照规模规则选出的比较族

mask_industry_first = company_survival_data['industry'] == industry_first  # 标记第一行业的公司
mask_industry_second = company_survival_data['industry'] == 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']  # 输入第二行业事件指标
)  # 完成两组同口径log-rank比较

print(f'对数秩检验统计量: {rank_test_results.test_statistic:.4f}')  # 报告预设两组曲线差异统计量
print(f'p值: {rank_test_results.p_value:.4f}')  # 报告同一检验的未调整p值

fig, ax = plt.subplots(figsize=(10, 6))  # 在同一坐标轴比较两个行业的KM曲线
# 用第一行业持续时间与事件拟合KM曲线
kaplan_meier_fitter.fit(company_survival_data.loc[mask_industry_first, 'duration'], company_survival_data.loc[mask_industry_first, 'event'], label=industry_first)  # 保持与log-rank第一组一致
kaplan_meier_fitter.plot(ax=ax)  # 展示第一行业的阶梯生存估计
# 用第二行业持续时间与事件拟合KM曲线
kaplan_meier_fitter.fit(company_survival_data.loc[mask_industry_second, 'duration'], company_survival_data.loc[mask_industry_second, 'event'], label=industry_second)  # 保持与log-rank第二组一致
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.5 风险函数的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)\) 相消,因此偏似然可以在不指定基线风险函数形式的情况下估计 \(\beta\)

然后将所有发生独立事件的时刻 \(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.5.1 准备数据:point-in-time landmark 设计

Cox 模型允许加入连续变量,但报告期末不是信息公开日。这里直接从 BOOK_DATA_DIR 读取财务报表,以 info_date 作为信息可得日,并把公司上市后首份可得财报的日期设为 landmark。缺少可得日或有限财务值的记录不进入模型,也不用“报告期末 + 固定天数”替代披露日。

表 11.3: Cox landmark 数据约定
数据项 本章约定
输入 BOOK_DATA_DIR/stock/financial_statement.h5
必需字段 order_book_id, quarter, info_date, total_assets, total_liabilities
时点规则 特征记录的 info_date 不晚于其 landmark
表 11.4: Cox模型数据准备
financial_statement_path = DATA_DIR / 'stock' / 'financial_statement.h5'  # 定位财务报表文件
assert financial_statement_path.is_file(), f'缺少财务报表: {financial_statement_path}'  # 检查文件存在
financial_statement_raw = pd.read_hdf(financial_statement_path).copy()  # 读取财务报表
required_financial_columns = {'order_book_id', 'quarter', 'info_date', 'total_assets', 'total_liabilities'}  # 定义最小字段集
assert required_financial_columns.issubset(financial_statement_raw.columns), '财务报表缺少 Cox 所需字段'  # 检查字段要求
financial_statement_raw['available_date'] = pd.to_datetime(financial_statement_raw['info_date'], errors='coerce')  # 以披露或修订日定义可得日
assert financial_statement_raw['available_date'].notna().all(), 'info_date 含缺失或无效日期'  # 检查可得日
finite_financial_values = np.isfinite(financial_statement_raw[['total_assets', 'total_liabilities']]).all(axis=1)  # 检查有限财务值
financial_statement_raw = financial_statement_raw.loc[finite_financial_values].copy()  # 只保留有限财务值
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
            })  # 完成当前公司的landmark杠杆记录

下面把逐公司记录转换为特征表,并检查每条记录的信息可得日。

列表 11.3: 构造 landmark 特征表并执行逐行可得日断言
initial_leverage_df = pd.DataFrame(initial_leverage_records, columns=['order_book_id', 'Initial_Leverage', 'Availability_Date', 'Landmark_Date'])  # 保留版本可得日与landmark
assert initial_leverage_df['Landmark_Date'].notna().all(), 'landmark 日期不能为空'  # 检查landmark日期
assert initial_leverage_df['Availability_Date'].le(initial_leverage_df['Landmark_Date']).all(), '特征可得日晚于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()}")  # 杠杆率分布概览
Cox 模型数据准备完成:
  Industry_Group   Region  Initial_Leverage
0          Other  Coastal          0.975951
1          Other  Coastal          0.631962
2            计算机  Coastal          0.368738
3          Other  Coastal          0.515535
4          Other  Coastal          0.712095

杠杆率统计:
count    5263.000000
mean        0.472239
std         0.683149
min         0.009280
25%         0.309223
50%         0.447294
75%         0.579983
max        41.403966
Name: Initial_Leverage, dtype: float64

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

11.5.2 拟合 Cox 比例风险模型

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

from lifelines import CoxPHFitter  # 估计 landmark 协变量与退市风险率的条件关联

# 准备回归数据
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  # 缺少有效杠杆时使用预先声明的中性分界
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模型
cox_proportional_hazards_model.fit(cox_model_data, duration_col='duration', event_col='event')  # 估计landmark协变量的条件风险关联

# 显示结果
cox_proportional_hazards_model.print_summary()  # 打印模型完整回归结果
cox_proportional_hazards_model.check_assumptions(cox_model_data, p_value_threshold=0.05, show_plots=False)  # 输出比例风险假设诊断
print(f'Cox 样本数: {len(cox_model_data)}; 事件数: {int(cox_model_data["event"].sum())}')  # 输出拟合样本概况
模型变量: ['duration', 'event', 'Initial_Leverage', 'Industry_Group_医药', 'Industry_Group_基础化工', 'Industry_Group_机械', 'Industry_Group_电子', 'Industry_Group_计算机', 'Region_Inland']
有效样本量: 5263
The ``p_value_threshold`` is set at 0.05. Even under the null hypothesis of no violations, some
covariates will be below the threshold by chance. This is compounded when there are many covariates.
Similarly, when there are lots of observations, even minor deviances from the proportional hazard
assumption will be flagged.

With that in mind, it's best to use a combination of statistical tests and visual tests to determine
the most serious violations. Produce visual plots using ``check_assumptions(..., show_plots=True)``
and looking for non-constant lines. See link [A] below for a full example.


1. Variable 'Initial_Leverage' failed the non-proportional test: p-value is 0.0011.

   Advice 1: the functional form of the variable 'Initial_Leverage' might be incorrect. That is,
there may be non-linear terms missing. The proportional hazard test used is very sensitive to
incorrect functional forms. See documentation in link [D] below on how to specify a functional form.

   Advice 2: try binning the variable 'Initial_Leverage' using pd.cut, and then specify it in
`strata=['Initial_Leverage', ...]` in the call in `.fit`. See documentation in link [B] below.

   Advice 3: try adding an interaction term with your time variable. See documentation in link [C]
below.


---
[A]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html
[B]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Bin-variable-and-stratify-on-it
[C]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Introduce-time-varying-covariates
[D]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Modify-the-functional-form
[E]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Stratification

Cox 样本数: 5263; 事件数: 208
表 11.6: Cox 比例风险模型回归结果
model lifelines.CoxPHFitter
duration col 'duration'
event col 'event'
penalizer 0.1
l1 ratio 0.0
baseline estimation breslow
number of observations 5263
number of events observed 208
partial log-likelihood -1591.73
time fit was run 2026-09-01 00:55:05 UTC
coef exp(coef) se(coef) coef lower 95% coef upper 95% exp(coef) lower 95% exp(coef) upper 95% cmp to z p -log2(p)
Initial_Leverage 0.15 1.17 0.03 0.09 0.22 1.09 1.24 0.00 4.74 <0.005 18.82
Industry_Group_医药 -0.30 0.74 0.14 -0.57 -0.04 0.57 0.96 0.00 -2.24 0.02 5.33
Industry_Group_基础化工 -0.29 0.75 0.14 -0.56 -0.03 0.57 0.97 0.00 -2.15 0.03 5.00
Industry_Group_机械 -0.22 0.80 0.12 -0.46 0.02 0.63 1.02 0.00 -1.77 0.08 3.69
Industry_Group_电子 -0.20 0.82 0.14 -0.49 0.08 0.61 1.08 0.00 -1.42 0.16 2.68
Industry_Group_计算机 -0.21 0.81 0.16 -0.53 0.11 0.59 1.12 0.00 -1.28 0.20 2.33
Region_Inland 0.17 1.19 0.08 0.02 0.32 1.02 1.38 0.00 2.18 0.03 5.10

Concordance 0.72
Partial AIC 3197.46
log-likelihood ratio test 33.61 on 7 df
-log2(p) of ll-ratio test 15.58
null_distribution chi squared
degrees_of_freedom 1
model <lifelines.CoxPHFitter: fitted with 5263 total...
test_name proportional_hazard_test
test_statistic p -log2(p)
Industry_Group_医药 km 0.00 0.97 0.04
rank 0.00 0.97 0.05
Industry_Group_基础化工 km 0.03 0.87 0.20
rank 0.03 0.86 0.21
Industry_Group_机械 km 0.08 0.78 0.36
rank 0.09 0.76 0.40
Industry_Group_电子 km 0.11 0.74 0.43
rank 0.12 0.73 0.45
Industry_Group_计算机 km 0.03 0.85 0.23
rank 0.04 0.85 0.24
Initial_Leverage km 6.83 0.01 6.80
rank 10.69 <0.005 9.86
Region_Inland km 0.47 0.49 1.02
rank 0.81 0.37 1.45

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

11.5.3 结果解释

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

fig, ax = plt.subplots(figsize=(10, 6))  # 为 Cox 系数点估计与区间共用横向尺度
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()  # 渲染Cox系数点估计与区间

# 提取模型摘要并输出显著变量的风险比解释
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值解释成确认性证据
横向森林图展示各协变量的Cox系数点估计及区间。
图 11.4: 任何原因退市的Cox条件关联森林图。

描述性风险比(非因果效应):
 - Initial_Leverage: HR=1.166
 - Industry_Group_医药: HR=0.738
 - Industry_Group_基础化工: HR=0.747
 - Industry_Group_机械: HR=0.803
 - Industry_Group_电子: HR=0.815
 - Industry_Group_计算机: HR=0.811
 - Region_Inland: HR=1.186

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

11.6 风险函数与累计风险(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.6.1 中国案例:杠杆分组的累计风险

生存函数描述存续超过某一时点的概率,风险函数(Hazard Function)则描述已存续至该时点条件下的瞬时事件率。 下面使用 NelsonAalenFitter 估计累计风险 \(H(t)\)。该估计由 Nelson 的 hazard plotting 与 Aalen 的计数过程框架发展而来 (Nelson 1972年; Aalen 1978年),这里按 landmark 杠杆率分组作描述。累计风险的跳跃不是瞬时风险率 \(h(t)\);若要估计 \(h(t)\),需要额外平滑与带宽选择。 Nelson–Aalen 曲线是累计风险的阶梯估计,跳跃表示累计事件强度增加,并不是某一时点的平滑瞬时风险。曲线本身不能识别制度改革、估值或公司基本面机制;这些解释需要预先规定的协变量、时间模型和额外证据。

from lifelines import NelsonAalenFitter  # 估计各杠杆组的累计事件强度

# 创建杠杆率分组(上限使用 inf 以涵盖杠杆率超过100%的公司)
# 按预设阈值把landmark杠杆率分成三个描述组
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%)'])  # 为三个区间附上可读标签

fig, ax = plt.subplots(figsize=(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'])  # 输入同一组的退市事件指标

    # 直接绘制 Nelson–Aalen 累计风险估计
    cumulative_hazard = nelson_aalen_fitter.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()  # 渲染三组累计风险阶梯曲线
按 landmark 杠杆率分组的多条累计风险阶梯曲线随生存年数上升,图例标明各组。
图 11.5: 不同 landmark 杠杆率分组的累计风险估计

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

11.7 Cox 模型的正则化

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

当 Cox 模型同时包含大量相关财务指标时,系数估计可能不稳定并出现过拟合。第 6 章的 Lasso 惩罚可以扩展到 Cox 偏似然,用于构造稀疏候选;变量是否稳定入选以及预测是否改善仍需训练内选择和未来期评价。 下面的机制演示向上市公司财务数据加入 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))  # 为全部系数使用同一 Lambda 横轴

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()  # 渲染惩罚强度变化下的Cox系数路径

# 只描述当次运行中的收缩路径,不预设噪声变量归零顺序
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=5.9417e-19): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1679: RuntimeWarning: overflow encountered in exp
  scores = weights * exp(dot(X, beta))
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=0): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=3.76969e-17): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1679: RuntimeWarning: overflow encountered in exp
  scores = weights * exp(dot(X, beta))
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=0): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=1.84773e-17): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1679: RuntimeWarning: overflow encountered in exp
  scores = weights * exp(dot(X, beta))
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=0): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1679: RuntimeWarning: overflow encountered in exp
  scores = weights * exp(dot(X, beta))
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=0): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1679: RuntimeWarning: overflow encountered in exp
  scores = weights * exp(dot(X, beta))
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=0): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1679: RuntimeWarning: overflow encountered in exp
  scores = weights * exp(dot(X, beta))
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=0): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1679: RuntimeWarning: overflow encountered in exp
  scores = weights * exp(dot(X, beta))
/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/lifelines/fitters/coxph_fitter.py:1530: LinAlgWarning: Ill-conditioned matrix (rcond=0): result may not be accurate.
  inv_h_dot_g_T = spsolve(-h, g, assume_a="pos", check_finite=False)
多条协变量系数曲线随横轴惩罚参数增大向零收缩,用于比较进入与退出路径。
图 11.6: Cox 模型的 Lasso 正则化路径

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

11.8 生存树模型

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

生存树可以用分割表示非线性与交互,但仍需要规定深度、叶节点样本量并做样本外验证。下面用 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))  # 为排序后的重要性提供统一横轴
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()  # 渲染生存树置换重要性排序

生存树特征重要性:
==================================================
                feature  importance
1     Industry_Group_医药    0.045303
0      Initial_Leverage    0.034602
2   Industry_Group_基础化工    0.033400
7               Noise_0    0.015349
11              Noise_4    0.009062
3     Industry_Group_机械    0.000000
5    Industry_Group_计算机    0.000000
4     Industry_Group_电子    0.000000
6         Region_Inland    0.000000
8               Noise_1    0.000000
9               Noise_2    0.000000
10              Noise_3    0.000000
水平条形图按置换重要性从高到低排列生存树输入特征,横轴为重要性数值。
图 11.7: 生存树可视化

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

11.9 模型评估:C-index

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

C-index(一致性指数,Concordance Index)比较可比样本对的风险排序;0.5 附近表示排序能力接近随机。下面按 landmark 日期划分较早与较晚队列。由于两组结局都随访到同一个 2024 年研究截止日,这只是回顾性队列留出评价,不能解释为在历史切分日真实部署的时间外预测。该指标也不代表概率校准、财务困境预测或因果效应。

表 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'])  # 对齐Cox样本的公司观察起点
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)]  # 以起点时间的80%分位锁定队列边界
training_mask = mature_landmark_dates < cohort_boundary  # 早期公司用于重新估计Cox模型
earlier_cohort_data = mature_cox_data.loc[training_mask]  # 构造较早回顾性训练队列
later_cohort_data = mature_cox_data.loc[~training_mask]  # 构造较晚回顾性评价队列

evaluation_cox_model = CoxPHFitter(penalizer=0.1)  # 使用固定惩罚强度稳定回顾性估计
evaluation_cox_model.fit(earlier_cohort_data, duration_col='duration', event_col='event')  # 只用较早队列拟合
predicted_risk_scores = evaluation_cox_model.predict_partial_hazard(later_cohort_data)  # 对较晚队列生成风险排序分数
print(f'队列切分日期: {cohort_boundary.date()}')  # 报告时间切分以便复核
print(f'较晚回顾性队列样本数: {len(later_cohort_data)}; 事件数: {int(later_cohort_data.event.sum())}')  # 同时报告评价信息量
try:  # 捕获没有可比事件对时指标不可定义的情形
    cox_concordance_index = concordance_index(  # 比较较晚队列中的可比风险排序
        later_cohort_data['duration'], -predicted_risk_scores, later_cohort_data['event'])  # lifelines以较大预测值对应较长生存,故风险分数取负
    print(f'较晚回顾性队列 C-index: {cox_concordance_index:.4f}')  # 输出定义良好的评价值
except ZeroDivisionError:  # 零事件或无可比对时不伪造数值
    cox_concordance_index = np.nan  # 用缺失值明确标记不可评价
    print('较晚回顾性队列 C-index: NA(没有可比较事件对)')  # 解释缺失值的统计原因
队列切分日期: 2017-10-31
较晚回顾性队列样本数: 823; 事件数: 2
较晚回顾性队列 C-index: 0.7500

请以代码实际输出为准,并同时报告切分日期与较晚队列规模。真正的历史部署评价应在每个预测起点把训练结局截断在当时可见范围内,并预先固定预测窗口;存在未完成随访时还需采用适当删失处理。单次回顾性切分仍可能不稳定,应再用多个后续队列或外部样本复核。

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

在投资风险管理中,回撤深度与回撤恢复时间(Time to Recovery)回答不同问题。本节把收盘价首次低于当时历史高点识别为真实水下回撤的开始,时间原点统一设为该水下区间紧邻的前一个峰值交易时段。

  • 目标集合:本地数据快照中,海康威视后复权价格路径从首个可观测交易时段到末个可观测交易时段之间,按下述规则能够完整识别起点的全部正持续时间水下区间。这是一个有限路径内集合,不代表投资者持仓、其他证券或未来回撤总体。
  • 入组:某交易时段首次低于截至当时的路径内累计峰值时,该水下区间入组;计时从紧邻的前一峰值交易时段开始。数据左边界之前已开始的水下状态无法识别,因而不属于目标集合。
  • 事件:水下区间后首次恢复到原峰值或更高收盘价。
  • 时间:从水下区间紧邻的前一峰值交易时段到恢复事件的交易时段数。
  • 删失:样本结束时仍未恢复的末次水下回撤在共同日历终点行政右删失;启动越晚,可获得的最大随访越短。
  • 依赖单位:所有区间来自同一证券价格路径,虽不重叠,却共享波动状态、制度环境和相邻峰值,不能当作独立抽样对象。

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

下面使用海康威视日线数据识别每个连续的水下区间。连续创新高但从未低于历史峰值的相邻时段不是回撤;一旦进入水下区间,就从紧邻的前一峰值时段起计时,首次回到该峰值为恢复事件,样本结束时仍水下的末次区间右删失。代码纳入所有正持续时间的实际水下回撤,不按最终恢复时长事后筛选。

这里不能证明独立删失:晚启动区间受共同样本终点限制,恢复时间又可能随启动日历时期和市场波动状态改变。同一路径区间也存在序列依赖。因此下方 Kaplan–Meier 算法只作为该有限路径内集合的描述性乘积极限汇总;纵轴不赋予外部总体概率含义,不报告基于独立对象近似的常规置信区间。20、60、180 个交易时段节点及中位数仅用于比较同一规则下的路径内恢复经验。

# 在真实 HDF 存储层直接选择海康威视,避免为单股恢复期示例载入全市场行情
drawdown_data_root_value = os.environ.get('BOOK_DATA_DIR')  # 独立读取恢复期案例的数据根配置
assert drawdown_data_root_value, '请先设置 BOOK_DATA_DIR,使其指向包含 stock/ 子目录的数据根'  # 缺失时给出修复方向
drawdown_data_root = Path(drawdown_data_root_value).expanduser().resolve()  # 解析恢复期案例的数据根
assert drawdown_data_root.is_dir(), f'BOOK_DATA_DIR 不存在或不是目录: {drawdown_data_root}'  # 核对数据目录
drawdown_price_path = drawdown_data_root / 'stock' / 'stock_price_post_adjusted.h5'  # 独立定位后复权股价文件
assert drawdown_price_path.is_file(), f'缺少后复权行情文件: {drawdown_price_path}'  # 在HDF读取前核对文件
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_underwater = closing_prices.lt(running_maximum_prices)  # 用严格不等号排除连续持平或创新高时段
is_episode_start = is_underwater & ~is_underwater.shift(fill_value=False)  # 标记每个连续水下区间的首日
episode_start_dates = closing_prices.index[is_episode_start]  # 仅保留真实进入水下的回撤起点

下面的确定性夹具锁定时间单位:周五到下一个交易日周一相隔一个交易时段,而不是三个日历日。

recovery_fixture_dates = pd.to_datetime(['2024-01-05', '2024-01-08'])  # 构造周五与下周一两个相邻交易时段
recovery_fixture_positions = pd.Series(range(2), index=recovery_fixture_dates)  # 用有序样本位置定义交易时段编号
assert recovery_fixture_positions.iloc[1] - recovery_fixture_positions.iloc[0] == 1  # 锁定周末不增加交易时段计数
pd.DataFrame({'date': recovery_fixture_dates, 'session_position': recovery_fixture_positions.to_numpy()})  # 输出可人工核对的夹具
表 11.8: 恢复时间的交易时段计数夹具
date session_position
0 2024-01-05 0
1 2024-01-08 1

基于上述水下指标,下面纳入每个实际回撤区间,并以紧邻的前一峰值时段作为统一时间原点:

# 以有序行情行号度量峰值到恢复或删失的交易时段数
trading_session_positions = pd.Series(np.arange(len(closing_prices)), index=closing_prices.index)  # 将每个观测日映射到连续交易时段位置
recovery_durations = []  # 存储每次回撤的持续交易时段数
recovery_events = []  # 存储事件标记(1=已恢复, 0=删失/未恢复)
recovery_origin_dates = []  # 存储每次回撤紧邻的前一峰值日期以审计入组时期

for episode_start_date in episode_start_dates:  # 逐个处理不重叠的连续水下区间
    start_position = int(trading_session_positions.loc[episode_start_date])  # 定位首个水下交易时段
    origin_position = start_position - 1  # 以紧邻的前一峰值时段作为时间原点
    episode_tail = is_underwater.iloc[start_position:]  # 在当前区间及其后寻找首次恢复
    recovered_dates = episode_tail.index[~episode_tail]  # 收集恢复到原峰值或更高的候选时段
    if len(recovered_dates) > 0:  # 样本内首次离开水下状态即为恢复事件
        end_position = int(trading_session_positions.loc[recovered_dates[0]])  # 锁定首个恢复时段
        recovery_events.append(1)  # 标记完整观测到恢复
    else:  # 末次水下区间可能延伸到样本终点
        end_position = len(closing_prices) - 1  # 以最后可观测交易时段作为删失时点
        recovery_events.append(0)  # 标记右删失而不伪造恢复事件
    recovery_durations.append(end_position - origin_position)  # 保留所有正持续时间的实际回撤
    recovery_origin_dates.append(closing_prices.index[origin_position])  # 登记决定最大潜在随访的日历入组时点

drawdown_recovery_data = pd.DataFrame({'origin_date': recovery_origin_dates, 'duration': recovery_durations, 'event': recovery_events})  # 构建含入组日期的回撤恢复数据框
assert not drawdown_recovery_data.empty, '当前价格样本未识别到实际水下回撤'  # 避免对空风险集拟合KM曲线
assert drawdown_recovery_data['duration'].gt(0).all(), '回撤恢复时间必须为正交易时段数'  # 锁定时间原点与事件顺序
print(f'识别到 {len(drawdown_recovery_data)} 次实际水下回撤')  # 输出无时长筛选的周期总数
recovered_mean_sessions = drawdown_recovery_data.loc[drawdown_recovery_data['event']==1, 'duration'].mean()  # 计算仅限已完成恢复周期的平均交易时段数
print(f'平均恢复时间(仅已恢复周期): {recovered_mean_sessions:.1f} 个交易时段')  # 明示条件均值的样本范围与时间单位
识别到 57 次实际水下回撤
平均恢复时间(仅已恢复周期): 45.0 个交易时段

回撤周期数量与恢复时间以当前数据快照的代码输出为准。这里的周期由“峰值—连续水下区间—恢复或删失”规则构造,不按事后持续时长删除短回撤。周期来自同一价格路径,彼此并非天然独立;复权方式和样本终点也会改变结果,因此它是描述性练习,不是对未来解套时间的承诺。

我们已识别出所有实际水下回撤周期的生存数据。下面使用 Kaplan–Meier 乘积公式汇总无时长筛选的回撤恢复曲线,并标注路径内中位恢复时间:

# 估计仍未恢复概率关于交易时段数的阶梯函数
kaplan_meier_fitter = KaplanMeierFitter()  # 初始化同一删失约定下的乘积极限估计器
kaplan_meier_fitter.fit(drawdown_recovery_data['duration'], drawdown_recovery_data['event'], label='回撤尚未恢复')  # 用交易时段持续时间与恢复指标拟合

fig, ax = plt.subplots(figsize=(10, 6))  # 为恢复曲线和中位数参照线共用坐标轴
kaplan_meier_fitter.plot_survival_function(ax=ax, ci_show=False)  # 路径依赖下撤除常规独立样本置信区间

# 添加中位恢复时间参考线
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')  # 声明横轴采用有序行情时段而非日历日
ax.set_ylabel('描述性乘积极限值', fontsize=12, fontproperties='Source Han Serif SC')  # 避免赋予未经识别的总体概率含义
ax.grid(True, alpha=0.3)  # 辅助读取指定交易时段的未恢复概率
plt.legend(prop={'family': 'Source Han Serif SC'})  # 使用已安装中文字体设置图例
plt.tight_layout()  # 自动调整布局
plt.show()  # 渲染单股回撤恢复的KM曲线

# 在课堂任务约定的交易时段节点读取未恢复概率
recovery_horizon_probabilities = {horizon: kaplan_meier_fitter.predict(horizon) for horizon in [20, 60, 180]}  # 统一计算三个预设交易时段节点
for horizon, probability in recovery_horizon_probabilities.items():  # 逐节点报告同一KM对象的估计
    print(f'全部已识别回撤在{horizon}个交易时段的描述性乘积极限值: {probability:.1%}')  # 明示路径内描述量与时间单位
不含常规置信区间的阶梯曲线汇总路径内水下回撤在各交易时段仍未恢复的比例乘积,并标出中位点。
图 11.8: 海康威视路径内实际水下回撤的描述性乘积极限曲线
全部已识别回撤在20个交易时段的描述性乘积极限值: 29.8%
全部已识别回撤在60个交易时段的描述性乘积极限值: 12.3%
全部已识别回撤在180个交易时段的描述性乘积极限值: 8.8%

图 11.8 描述当前样本中所有已识别水下回撤“尚未恢复”的路径内乘积极限汇总,不以最终持续时间筛选周期。20、60、180 个交易时段及中位恢复时间应直接从本次拟合对象读取;右删失、同一路径依赖与单一标的限制使其不能直接转化为止损规则、系统性风险归因或基本面判断。

为检查日历支持域与恢复状态是否混在一条曲线中,下面按回撤入组日期的中位数把路径分成较早、较晚两段,在每段内用完全相同的事件、时间和删失规则重算三个预设节点。这个切分只是一项透明的路径敏感性检查,不把两个时期视为独立随机样本,也不进行普通显著性检验。

drawdown_period_boundary = drawdown_recovery_data['origin_date'].median()  # 用路径内入组日期中位数形成透明的等规模时期切分
drawdown_recovery_data['start_period'] = np.where(drawdown_recovery_data['origin_date'].le(drawdown_period_boundary), '较早启动', '较晚启动')  # 登记每个周期的启动时期
drawdown_sensitivity_records = []  # 收集两个时期的样本、删失与节点汇总
for start_period, period_data in drawdown_recovery_data.groupby('start_period'):  # 在同一规则下分别重算两个时期
    period_fitter = KaplanMeierFitter().fit(period_data['duration'], period_data['event'])  # 拟合仅用于路径敏感性的乘积极限曲线
    drawdown_sensitivity_records.append({  # 保存可复核的时期内描述量
        'start_period': start_period, 'episodes': len(period_data), 'censored': int(period_data['event'].eq(0).sum()),  # 报告时期样本与行政删失
        'pl_20': period_fitter.predict(20), 'pl_60': period_fitter.predict(60), 'pl_180': period_fitter.predict(180)  # 读取三个预设节点
    })  # 完成当前启动时期的敏感性记录
drawdown_period_sensitivity = pd.DataFrame(drawdown_sensitivity_records).set_index('start_period')  # 汇总为时期对照表
print(f'启动时期切分日期: {drawdown_period_boundary.date()}')  # 报告数据决定的透明切分边界
drawdown_period_sensitivity.round(3)  # 输出各时期节点与删失数量
启动时期切分日期: 2015-01-13
表 11.9: 回撤启动时期的描述性乘积极限敏感性
episodes censored pl_20 pl_60 pl_180
start_period
较早启动 29 0 0.310 0.103 0.069
较晚启动 28 1 0.286 0.143 0.107

表 11.9 的较早、较晚节点差异明显,或删失集中在较晚启动组,应明确写成“总体汇总对日历时期和可随访长度敏感”,而不能选择更有利的一段作为未来恢复概率。更深入的分析可预先规定市场状态、以证券为聚类单位做多证券区块重抽样,或采用滚动起点外推;本单一路径练习不提供足以支持这些总体推断的独立单位。

11.11 回撤恢复时间课堂活动

本活动沿用正文的真实海康威视价格数据,只把连续水下区间识别为回撤。时间原点是水下首日紧邻的前一峰值交易时段,恢复事件是首次回到该峰值;样本结束前未恢复的末次回撤右删失。学生应纳入所有实际回撤,不按最终时长筛选。

11.11.1 数据与解释

用不含常规置信区间的 Kaplan–Meier 计算报告路径内中位恢复时间以及 20、60 个交易时段的描述性乘积极限值,并逐项列出风险集、事件数和删失数。再按 表 11.9 比较较早与较晚启动区间。结论只描述当前标的、当前样本窗口和当前回撤定义;同一价格路径上周期的依赖、晚入组的较短潜在随访、单一证券与市场状态变化都限制外推。

11.12 本章小结

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

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

生存分析的应用建议

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

关键假设

  • 独立删失的可信度主要依靠数据收集设计、退出原因与敏感性分析评估,不能由已观测结局单独证明。
  • 比例风险可用 Schoenfeld 残差寻找违背证据;时间交互或分层 Cox 只处理这类违背。

软件工具

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

11.13 理论来源与前沿

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

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

  1. 高维协变量与正则化:当特征数远大于样本量时,需要用 Lasso/Elastic Net 等正则化构建稀疏Cox 模型,并配套稳定的变量选择与不确定性量化。
  2. 时间变化效应与非比例风险:在真实商业场景(例如客户流失)中,某些协变量的影响可能随时间变化,此时可考虑时间交互项、分层 Cox 或 AFT 等替代。Schoenfeld 残差提供比例风险诊断 (Schoenfeld 1982年),加权残差检验把诊断推广为协变量级与全局检验 (Grambsch 和 Therneau 1994年);有限事件数会降低检验功效,同时检验多个协变量时还应报告检验族和多重性处理,不能把单个未拒绝结果解释为比例风险已被证明。
  3. 机器学习与因果推断的结合:将生存模型与树模型、Boosting、深度学习以及因果推断框架结合,用于处理复杂非线性、异质性处理效应与动态干预策略。

11.14 练习

11.14.1 概念题

  1. [核心|难度:1|分值:4|任务:独立] 解释什么是右删失、左删失和区间删失,并各举一个实际例子。

  2. [核心|难度:2|分值:5|任务:独立] Kaplan-Meier 估计量的核心思想是什么?为什么它比简单地计算生存比例更合理?

  3. [核心|难度:2|分值:5|任务:独立] Cox 比例风险模型中的“比例风险”假设是什么意思?如何检验这个假设?

  4. [核心|难度:2|分值:5|任务:独立] 解释风险函数 \(h(t)\) 和生存函数\(S(t)\) 之间的关系。

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

11.14.2 应用题

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

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

  2. [核心|难度:3|分值:15|任务:独立] 真实 A 股 landmark 退市评价(必做):沿用 小节 11.4.2 的任何原因退市定义和 小节 11.9 的回顾性队列设计,提交一个紧凑证据表,至少包含时间原点、事件日期或行政删失日、\((Y,\delta)\) 逐行检查结果、时间切分日期、早期与晚期队列的样本数、事件数、删失率,以及晚期队列 C-index。若晚期队列零事件或没有可比对,必须输出 NA 与停止原因。最后用两句话分别说明独立删失为何不能只由数据表证明,以及该切分为何不是历史时点真实部署评价。

  3. [拓展|难度:3|分值:15|任务:独立] 模拟信贷违约数据: 假设你是一家网贷平台的风控分析师,请模拟一份贷款数据(包含借款人年龄、收入、信用分、借款金额等特征)。

    要求:

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

    要求:

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

11.14.3 理论题

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

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

  3. [拓展|难度:3|分值:12|任务:独立] 证明在单变量二元协变量的情况下,Cox 模型的score test 等价于对数秩检验。

11.15 练习参考解答

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

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

11.15.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 残差并使用加权残差检验 (Schoenfeld 1982年; Grambsch 和 Therneau 1994年);拒绝结果是反对比例风险的证据,未拒绝则可能来自事件数不足。若同时检查多个协变量,应同时报告全局检验、逐项检验数量与相应多重性处理,再依据偏离形状考虑时间交互、分层 Cox 或其他生存模型。

  4. 风险函数与生存函数

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

  5. 独立删失

    若删失机制包含有关 \(T\) 的信息(例如财务状况恶化导致的数据缺失),则不再独立,标准方法可能产生偏差。这个条件通常不能从已观测的 \((Y,\delta)\) 单独检验成立;应核对随访制度、退出原因、可观测预测变量与失访模式,并对未观测机制做敏感性分析。非比例风险可用时间交互或分层 Cox 处理,而疑似信息性删失应考虑 IPCW、事件—删失联合模型或模式混合分析。

11.15.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. 真实 A 股 landmark 退市评价

    唯一规范数据构造是 表 11.5:时间原点为首份可靠财报的 Landmark_Date,事件为研究截止日前任何原因退市,未退市公司在共同截止日行政右删失,持续时间从 landmark 起算。唯一规范切分与评分实现是 表 11.7;解答不得另开随机切分或重读全市场 HDF。下面只把题目特有的逐行时序检查、队列事件/删失证据和现场 C-index 整理为一个可评分表。

表 11.10: 习题7:真实退市事件、删失、时间切分与 C-index 证据
exercise7_row_audit = company_survival_data[['order_book_id', 'Landmark_Date', 'analysis_end', 'duration', 'event']].copy()  # 复用规范landmark事件时间数据
exercise7_row_audit['is_time_order_valid'] = exercise7_row_audit['analysis_end'].gt(exercise7_row_audit['Landmark_Date'])  # 逐行检查结局晚于入组
exercise7_row_audit['is_status_valid'] = exercise7_row_audit['event'].isin([0, 1])  # 逐行检查删失编码属于二元集合
assert exercise7_row_audit[['is_time_order_valid', 'is_status_valid']].all().all(), '事件时间或删失编码检查失败'  # 任一错误即停止评分链
exercise7_split_evidence = pd.DataFrame({  # 汇总早期训练与晚期评价队列的可观察证据
    'cohort': ['earlier_training', 'later_evaluation'],  # 明示两个时间队列职责
    'n': [len(earlier_cohort_data), len(later_cohort_data)],  # 报告各队列样本数
    'events': [int(earlier_cohort_data['event'].sum()), int(later_cohort_data['event'].sum())],  # 报告可用于排序评价的事件数
    'censoring_rate': [earlier_cohort_data['event'].eq(0).mean(), later_cohort_data['event'].eq(0).mean()]  # 报告各队列删失率
})  # 完成队列审计表
exercise7_c_index_display = 'NA(没有可比较事件对)' if pd.isna(cox_concordance_index) else f'{cox_concordance_index:.4f}'  # 保留指标不可定义状态
print({'time_origin': 'Landmark_Date', 'event_or_censor_end': 'analysis_end', 'split_date': str(cohort_boundary.date()), 'later_c_index': exercise7_c_index_display})  # 输出关键合同字段
print(exercise7_split_evidence.round(3))  # 输出样本、事件和删失率
print(exercise7_row_audit.head(10))  # 展示可人工复核的逐行事件时间证据
{'time_origin': 'Landmark_Date', 'event_or_censor_end': 'analysis_end', 'split_date': '2017-10-31', 'later_c_index': '0.7500'}
             cohort     n  events  censoring_rate
0  earlier_training  3260     205           0.937
1  later_evaluation   823       2           0.998
  order_book_id Landmark_Date analysis_end   duration  event  \
0   000001.XSHE    2006-04-26   2024-01-01  17.683778      0   
1   000002.XSHE    2006-04-25   2024-01-01  17.686516      0   
2   000004.XSHE    2006-04-28   2024-01-01  17.678303      0   
3   000005.XSHE    2006-04-28   2024-01-01  17.678303      0   
4   000006.XSHE    2006-04-22   2024-01-01  17.694730      0   
5   000007.XSHE    2006-04-28   2024-01-01  17.678303      0   
6   000008.XSHE    2006-04-21   2024-01-01  17.697467      0   
7   000009.XSHE    2006-04-29   2024-01-01  17.675565      0   
8   000010.XSHE    2006-04-29   2024-01-01  17.675565      0   
9   000011.XSHE    2006-04-20   2024-01-01  17.700205      0   

   is_time_order_valid  is_status_valid  
0                 True             True  
1                 True             True  
2                 True             True  
3                 True             True  
4                 True             True  
5                 True             True  
6                 True             True  
7                 True             True  
8                 True             True  
9                 True             True  
满分答案必须保留两项限制。第一,行政删失是否在给定协变量后独立于潜在退市时间,不能仅由 $(Y,\delta)$ 表证明,还需结合上市队列成熟度、数据覆盖制度与未观测退出机制;若不可信,应做 IPCW 或敏感性分析。第二,较早、较晚公司都随访到同一 2024 年截止日,@tbl-ex7-real-survival-evidence 只是按 landmark 排序的回顾性队列留出,不是切分日在真实历史中锁定结局并部署的前向评价。晚期队列无事件或无可比对时,`NA` 是正确证据,不能替换成 0.5。
  1. 信贷违约模拟与 Lasso-Cox(完整实现)

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

列表 11.4: 习题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.5: 习题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  # 构造板块级上市持续时间、退市事件与删失记录
from lifelines import KaplanMeierFitter  # 估计各板块任何原因退市的生存曲线
from lifelines.statistics import logrank_test  # 比较预设板块对的事件时间分布
import matplotlib.pyplot as plt  # 在相同坐标尺度叠加板块曲线

from pathlib import Path  # 使用跨平台路径对象解析显式数据根
board_data_root_value = os.environ.get('BOOK_DATA_DIR')  # 独立读取板块比较所需的数据根
assert board_data_root_value, '请先设置 BOOK_DATA_DIR,使其指向包含 stock/ 子目录的数据根'  # 缺失时给出配置提示
DATA_DIR = Path(board_data_root_value).expanduser().resolve()  # 从必需环境变量取得数据根
if not DATA_DIR.is_dir():  # 在读取前验证数据根
    raise FileNotFoundError(f'BOOK_DATA_DIR 不存在或不是目录: {DATA_DIR}')  # 失败即停止,不回退固定路径
board_basic_path = DATA_DIR / 'stock' / 'stock_basic_data.h5'  # 定位板块比较使用的公司基本信息
assert board_basic_path.is_file(), f'缺少公司基本信息文件: {board_basic_path}'  # 在HDF读取前核对文件
stock_basic_data_ex = pd.read_hdf(board_basic_path)  # 加载上市公司基本信息

# 将上市和退市日期转换为 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):  # 按代码前缀复现题面指定的教学板块映射
    """根据股票代码判断所属板块"""  # 文档字符串说明映射对象而非交易所历史制度
    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拟合
    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.9: A股不同板块任何原因退市的 Kaplan–Meier 生存曲线

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

表 11.11: 习题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.11 时必须同时报告事件数;静态退市比例不等于生存概率,也不能替代 图 11.9 的风险集校正。

11.15.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.16 章末回顾

学习自检:定义生存时间、事件与删失;从风险集写出 Kaplan–Meier 乘积极限;解释 Cox 偏似然为何消去基线风险;检查比例风险和 landmark 时点;用 C-index 或预先指定的时间点损失比较模型。完成事件/删失定义、时点检查与评价可定义性后,再解释模型结果。

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

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

展开检索反馈与学习决策

1)不产生事件因子,但会改变后续风险集人数。2)条件化后基线风险相消,偏似然估计相对风险参数;绝对生存概率仍需累计基线风险。3)否则特征或入组会使用结果发生后的信息。需要复习时,可在 小节 11.4 改变删失月份并重算风险集,在 小节 11.5 写出两个风险集的分子分母,并在 小节 11.9 重画 landmark、特征截止和测试窗口。

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

Aalen, Odd O. 1978年. 《Nonparametric Inference for a Family of Counting Processes》. The Annals of Statistics 6 (4): 701~26. https://doi.org/10.1214/aos/1176344247.
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.
Grambsch, Patricia M., 和 Terry M. Therneau. 1994年. 《Proportional Hazards Tests and Diagnostics Based on Weighted Residuals》. Biometrika 81 (3): 515~26. https://doi.org/10.1093/biomet/81.3.515.
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.
Nelson, Wayne. 1972年. 《Theory and Applications of Hazard Plotting for Censored Failure Data》. Technometrics 14 (4): 945~66. https://doi.org/10.1080/00401706.1972.10488991.
Schoenfeld, David. 1982年. 《Partial Residuals for the Proportional Hazards Regression Model》. Biometrika 69 (1): 239~41. https://doi.org/10.1093/biomet/69.1.239.