竞争风险模型下的Landmark生存分析是一种评估特定时间点之后的生存概率，并考虑到其他事件的竞争效应的重要统计方法。下面我们将介绍如何通过Python编写代码来进行这种分析，以及绘制分段前后的图形。 ### 使用`lifelines`库进行Landmark生存分析 #### 安装依赖包 首先需要安装一些必要的Python包： ```bash pip install lifelines matplotlib pandas numpy scikit-survival ``` #### 导入所需模块并准备数据集 假设我们已经有了一个包含患者信息的数据框df，它应该有类似如下的结构：`(time, event, competing_risk)`，其中`event`指示目标事件是否发生而`competing_risk`则表示是否有竞争性风险发生。 接下来导入所需的库并将数据准备好: ```python import pandas as pd from sksurv.datasets import load_veterans_lung_cancer from sksurv.preprocessing import OneHotEncoder from sksurv.linear_model import CoxPHSurvivalAnalysis from sklearn.model_selection import train_test_split import matplotlib.pyplot as plt from lifelines import KaplanMeierFitter from lifelines.plotting import add_at_risk_counts # 加载示例数据集 (替换为您自己的数据) data_x, data_y = load_veterans_lung_cancer() data_x.loc[:, "Celltype"] = OneHotEncoder().fit_transform(data_x[["Celltype"]]) # 将DataFrame转换成适合sklearn使用的格式 def get_X_y(df): X = df.drop(["Status", "Survival_in_days"], axis=1).values.astype(float) y = [(True if s == 'dead' else False, t) for s, t in zip(df['Status'], df['Survival_in_days'])] return X, np.array(y, dtype=[('status', '?'), ('survival_time', '<f8')]) X, y = get_X_y(pd.DataFrame(np.hstack([data_x.values, data_y.reshape(-1, 1)]), columns=list(data_x.columns)+['Status'])) # 分割训练测试集 train_x, test_x, train_y, test_y = train_test_split(X, y, test_size=.25) kmf = KaplanMeierFitter() times_of_interest = [90] # landmark time points for i, tlmrk in enumerate(times_of_interest): mask_train = train_y["survival_time"] >= tlmrk kmf.fit(train_y[mask_train]["survival_time"] - tlmrk, event_observed=train_y[mask_train]['status'], label=f'Train Data > {tlmrk} days') plt.figure(figsize=(7,4)) ax = plt.gca() _ = kmf.survival_function_.plot(ax=ax) mask_test_pre_landmark = test_y["survival_time"] < max(times_of_interest) prelandmark_data = test_y[mask_test_pre_landmark] post_landmark_times = [] events_after_lm = [] for pt in prelandmark_data: idx = next((i for i,t in enumerate(times_of_interest+[max(prelandmark_data['survival_time'])]) \ if pt['survival_time']<t), None) post_event = ((pt['survival_time'] + times_of_interest[idx-1])-times_of_interest[idx-1], pt['status']) events_after_lm.append(post_event) post_landmark_times.extend([times_of_interest[idx-1]]*len(events_after_lm)) # 创建新的KM估计用于绘图 kmf_post_landmarks = KaplanMeierFitter() kmf_post_landmarks.fit(durations=[x[0]-y for x,y in zip(events_after_lm,max(post_landmark_times))], event_observed=[e for _, e in events_after_lm],label='Post Landmark') # 绘制两组曲线对比图 fig, ax = plt.subplots(1, figsize=(8,6)) _ = kmf.survival_function_[kmf.label].plot(drawstyle="steps-post", linestyle="--", color="#A6CEE3", ax=ax,label='Before Landmark' ) _ = kmf_post_landmarks.survival_function_.plot(color='#B2DF8A', drawstyle="steps-post", ax=ax,label='After Landmark' ) add_at_risk_counts(kmf, kmf_post_landmarks, ax=ax) plt.title("Comparing Survival Curves Before and After a Specified Time Point") plt.xlabel("Time from Start of Study or Post-Landmark Event Time ") plt.ylabel("Estimated Probability of Event-Free Survival") plt.legend(loc='upper right') plt.show() ``` 此代码片段展示了怎样创建两个Kaplan-Meier图表来比较在指定的时间节点之前和之后的存活率差异。请注意实际应用时应调整变量名、路径以及其他细节以适应您的具体情况。