import pandas as pd
import numpy as np
from lifelines import KaplanMeierFitter
from lifelines.statistics import logrank_test
from sklearn.ensemble import RandomForestRegressor
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
import seaborn as sns
import re
import warnings
warnings.filterwarnings('ignore')
# 设置中文字体和样式
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False
sns.set_style("whitegrid")
def load_data():
"""加载数据"""
excel_path = r"C:\Users\de'l'l\AppData\Local\Temp\wpscompress\{adbe2bb2-d895-47c2-b31a-c60db3c874a0}\{200e1ff9-59ec-4064-b2c6-d8c338897282}\附件.xlsx"
try:
df = pd.read_excel(excel_path, sheet_name="男胎检测数据")
print("成功从Excel文件加载数据")
return df
except Exception as e:
print(f"加载Excel文件失败: {e}")
print("使用模拟数据...")
return generate_simulated_data()
def generate_simulated_data():
"""生成模拟数据"""
np.random.seed(42)
n_patients = 300
data = []
for pid in range(n_patients):
bmi = np.random.normal(22, 3)
age = np.random.randint(20, 40)
height = np.random.normal(162, 5)
weight = bmi * (height/100)**2
# 模拟3次检测
for week in [12, 16, 20]:
y_concentration = np.clip(
(week-10)/20 * 0.1 * (1 + 0.1*(bmi-22)) * (1 + np.random.normal(0, 0.1)),
0.01, 0.1
)
data.append({
'孕妇代码': f'P{pid:03d}',
'检测孕周': f"{week}w+{np.random.randint(0,6)}d",
'Y染色体浓度': y_concentration,
'孕妇BMI': bmi,
'年龄': age,
'身高': height,
'体重': weight
})
return pd.DataFrame(data)
def convert_week_to_numeric(week_str):
"""将孕周字符串转换为数值"""
if pd.isna(week_str):
return np.nan
if isinstance(week_str, (int, float)):
return week_str
week_str = str(week_str).lower().replace('周', '').replace('w', '')
match = re.search(r'(\d+)(?:\+(\d+))?', week_str)
if match:
weeks = int(match.group(1))
days = int(match.group(2)) if match.group(2) else 0
return weeks + days / 7
return np.nan
def prepare_survival_data(df):
"""准备生存分析数据"""
df['检测孕周数值'] = df['检测孕周'].apply(convert_week_to_numeric)
df['达标'] = df['Y染色体浓度'] >= 0.04
# 计算每个孕妇的最早达标时间
earliest_pass = df[df['达标']].groupby('孕妇代码')['检测孕周数值'].min().reset_index()
earliest_pass.rename(columns={'检测孕周数值': '达标孕周'}, inplace=True)
# 计算每个孕妇的最后检测时间
last_measure = df.groupby('孕妇代码')['检测孕周数值'].max().reset_index()
last_measure.rename(columns={'检测孕周数值': '最后检测孕周'}, inplace=True)
# 合并数据
survival_df = df.drop_duplicates('孕妇代码').merge(
earliest_pass, on='孕妇代码', how='left'
).merge(
last_measure, on='孕妇代码', how='left'
)
# 设置生存时间和事件状态
survival_df['T'] = np.where(
survival_df['达标孕周'].notna(),
survival_df['达标孕周'],
survival_df['最后检测孕周']
)
survival_df['E'] = survival_df['达标孕周'].notna().astype(int)
return survival_df
def analyze_with_lifelines():
"""使用lifelines进行生存分析"""
# 1. 加载和准备数据
df = load_data()
survival_data = prepare_survival_data(df)
# 2. 特征重要性分析
features = ['孕妇BMI', '年龄', '身高', '体重']
X = survival_data[features].fillna(0)
y = survival_data['T'].fillna(0)
rf = RandomForestRegressor(random_state=42)
rf.fit(X, y)
importance = pd.DataFrame({
'特征': features,
'重要性': rf.feature_importances_
}).sort_values('重要性', ascending=False)
# 3. 基于BMI聚类分组
kmeans = KMeans(n_clusters=3, random_state=42)
survival_data['分组'] = kmeans.fit_predict(X)
# 4. 分组生存分析
results = []
for group in sorted(survival_data['分组'].unique()):
group_data = survival_data[survival_data['分组'] == group]
# Kaplan-Meier估计
kmf = KaplanMeierFitter()
kmf.fit(
group_data['T'],
event_observed=group_data['E'],
label=f'分组{group+1}'
)
# 找到达标比例达到95%的孕周
target_survival = 0.05 # 1-95%
if any(kmf.survival_function_.values <= target_survival):
rec_week = kmf.survival_function_.index[
np.where(kmf.survival_function_.values <= target_survival)[0][0]
]
else:
rec_week = kmf.survival_function_.index[-1]
# 计算该时点的实际达标率
actual_rate = 1 - kmf.predict(rec_week)
results.append({
'分组': group + 1,
'推荐孕周': rec_week,
'预期达标率': actual_rate,
'人数': len(group_data),
'BMI范围': f"{group_data['孕妇BMI'].min():.1f}-{group_data['孕妇BMI'].max():.1f}"
})
# 绘制生存曲线
plt.figure(figsize=(8, 5))
kmf.plot_survival_function()
plt.axvline(x=rec_week, color='r', linestyle='--', label=f'推荐时点: {rec_week:.1f}周')
plt.axhline(y=target_survival, color='g', linestyle=':', label='5%未达标')
plt.title(f'分组{group+1}的达标曲线')
plt.xlabel('孕周')
plt.ylabel('未达标概率')
plt.legend()
plt.tight_layout()
plt.savefig(f'survival_group_{group+1}.png', dpi=300)
plt.close()
return pd.DataFrame(results), importance
if __name__ == "__main__":
# 运行分析
recommendations, importance = analyze_with_lifelines()
# 输出结果
print("\n=== 特征重要性分析结果 ===")
print(importance.to_string(index=False))
print("\n=== 分组推荐检测时点 ===")
print(recommendations.to_string(index=False))
# 保存结果
recommendations.to_excel('推荐结果_lifelines.xlsx', index=False)
print("\n分析完成，结果已保存到以下文件:")
print("- survival_group_*.png")
print("- 推荐结果_lifelines.xlsx")
多因素影响下的分组与优化（多因素）
· 传统方法：多元回归预测达标时间，再分组。
· 创新性算法：主成分分析（PCA） + 混合整数规划（MIP） 或 神经网络（DNN） + 遗传算法（GA）
· 分析：因素增多（身高、体重、年龄、BMI）导致维度灾难和多重共线性。需要降维并同时优化分组和时点，问题更复杂。
· 算法融合应用：
1. 降维与特征构建：首先使用PCA将身高、体重、BMI等高度相关的指标融合成几个独立的“体格综合指标”主成分，消除共线性，简化问题。
2. 优化模型构建：将问题直接建模为一个混合整数规划（MIP） 模型。
· 决策变量：① 二进制变量，表示每个孕妇属于哪个组；② 连续变量，表示每组的最佳NIPT时点。
· 目标函数：最小化所有孕妇的潜在风险总和（如：预测达标时间与组内时点差值的绝对值的加权和）。
· 约束：每组人数上下限、达标比例约束（如95%的孕妇在该时点已达标）。
3. 使用遗传算法（GA） 或模拟退火（SA） 来高效求解这个复杂的MIP模型，它们特别擅长处理这种组合优化问题。
· 误差分析：同样采用蒙特卡罗模拟，在PCA降维后的综合得分和达标时间预测中引入误差，多次运行优化模型，观察分组方案和推荐时点的鲁棒性。
"F:\附件男胎.xlsx"（数据文件地址）按照问题要求对代码进行整体优化
本文深入探讨了Objective-C与Swift在iOS开发领域的应用与对比，旨在为开发者提供选择合适的编程语言的指南。
扫一扫
1637
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
