import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from sklearn.cluster import KMeans
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
import warnings
warnings.filterwarnings('ignore')
# 读取清洗后的数据
df = pd.read_csv('cleaned_data.csv')
# 确保只处理男胎数据:Y染色体浓度非空
df_male = df[df['Y染色体浓度'].notna()].copy()
print(f"男胎数据量: {len(df_male)}")
# 计算孕周数(如果尚未计算)
def extract_weeks(s):
if isinstance(s, str):
if 'w' in s:
parts = s.split('w')
weeks = float(parts[0])
if '+' in parts[1]:
days = float(parts[1].split('+')[1].replace('d', ''))
weeks += days / 7
return weeks
return np.nan
if '孕周数' not in df_male.columns:
df_male['孕周数'] = df_male['检测孕周'].apply(extract_weeks)
# 移除孕周数或Y染色体浓度为NaN的行
df_male = df_male.dropna(subset=['孕周数', 'Y染色体浓度', '孕妇BMI'])
# 定义BMI分组边界(初始分组基于临床经验)
bmi_bins = [20, 28, 32, 36, 40, np.inf]
bmi_labels = ['20-28', '28-32', '32-36', '36-40', '40+']
df_male['BMI分组'] = pd.cut(df_male['孕妇BMI'], bins=bmi_bins, labels=bmi_labels, right=False)
# 检查每组数据量
print("各BMI分组数据量:")
print(df_male['BMI分组'].value_counts())
# 对于每个BMI分组,找到Y染色体浓度达到4%的孕周数
# 由于每个孕妇可能有多次检测,我们首先对每个孕妇找到浓度达标的最早孕周
# 分组按孕妇代码和BMI分组
df_male['达标'] = df_male['Y染色体浓度'] >= 4.0
# 对于每个孕妇,找到达标的最早孕周
earliest_weeks = df_male[df_male['达标']].groupby('孕妇代码')['孕周数'].min().reset_index()
earliest_weeks.rename(columns={'孕周数': '达标孕周'}, inplace=True)
# 合并回原数据,获取每个孕妇的BMI分组
patient_bmi = df_male.groupby('孕妇代码')['孕妇BMI'].first().reset_index()
patient_bmi['BMI分组'] = pd.cut(patient_bmi['孕妇BMI'], bins=bmi_bins, labels=bmi_labels, right=False)
earliest_weeks = earliest_weeks.merge(patient_bmi, on='孕妇代码', how='left')
# 对于每个BMI分组,计算达标孕周的分布
grouped = earliest_weeks.groupby('BMI分组')['达标孕周']
descriptive_stats = grouped.describe(percentiles=[0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95])
print("各BMI分组达标孕周的描述统计:")
print(descriptive_stats)
# 确定最佳NIPT时点:对于每组,选择第95百分位数的达标孕周,以确保95%的孕妇在该时点已达标
best_timing = grouped.quantile(0.95).reset_index()
best_timing.rename(columns={'达标孕周': '最佳NIPT时点(周)'}, inplace=True)
print("各BMI分组的最佳NIPT时点(第95百分位数):")
print(best_timing)
# 可视化达标孕周分布
plt.figure(figsize=(10, 6))
sns.boxplot(x='BMI分组', y='达标孕周', data=earliest_weeks)
plt.axhline(y=12, color='g', linestyle='--', label='12周(早期风险低)')
plt.axhline(y=27, color='r', linestyle='--', label='27周(中期风险高)')
plt.title('各BMI分组Y染色体浓度达标孕周分布')
plt.ylabel('孕周(周)')
plt.legend()
plt.savefig('各BMI分组达标孕周分布.png')
plt.show()
# 风险分析:比较最佳时点与风险周数
best_timing['风险等级'] = best_timing['最佳NIPT时点(周)'].apply(
lambda x: '低风险' if x <= 12 else '高风险' if x <= 27 else '极高风险'
)
print("最佳NIPT时点的风险等级:")
print(best_timing[['BMI分组', '最佳NIPT时点(周)', '风险等级']])
# 检测误差分析
# 假设Y染色体浓度测量误差为±0.5%(根据临床典型误差),孕周估计误差为±0.5周
# 使用蒙特卡洛模拟评估误差对达标时间估计的影响
np.random.seed(42)
n_simulations = 1000
error_results = []
for bmi_group in bmi_labels:
group_data = earliest_weeks[earliest_weeks['BMI分组'] == bmi_group]
if len(group_data) == 0:
continue
# 模拟误差
simulated_weeks = []
for _, row in group_data.iterrows():
base_week = row['达标孕周']
# 模拟孕周误差和浓度误差(浓度误差可能影响达标时间,但这里我们直接模拟达标孕周误差)
# 由于达标孕周是基于浓度≥4%定义的,误差会同时影响浓度和孕周
# 简化:假设达标孕周估计有误差,误差分布为正态(均值0,标准差0.5周)
simulated_week = np.random.normal(base_week, 0.5, n_simulations)
simulated_weeks.extend(simulated_week)
simulated_weeks = np.array(simulated_weeks)
# 计算模拟后的第95百分位数
sim_95perc = np.percentile(simulated_weeks, 95)
error_results.append({
'BMI分组': bmi_group,
'原始最佳时点': best_timing[best_timing['BMI分组'] == bmi_group]['最佳NIPT时点(周)'].iloc[0],
'模拟最佳时点': sim_95perc,
'误差': sim_95perc - best_timing[best_timing['BMI分组'] == bmi_group]['最佳NIPT时点(周)'].iloc[0]
})
error_df = pd.DataFrame(error_results)
print("检测误差对最佳NIPT时点的影响:")
print(error_df)
# 可视化误差影响
plt.figure(figsize=(10, 6))
x_pos = np.arange(len(bmi_labels))
plt.bar(x_pos, error_df['原始最佳时点'], width=0.4, label='原始最佳时点')
plt.bar(x_pos + 0.4, error_df['模拟最佳时点'], width=0.4, label='模拟最佳时点(含误差)')
plt.xticks(x_pos + 0.2, bmi_labels)
plt.xlabel('BMI分组')
plt.ylabel('最佳NIPT时点(周)')
plt.title('检测误差对最佳NIPT时点的影响')
plt.legend()
plt.savefig('检测误差影响.png')
plt.show()
# 输出最终分组和最佳时点
print("\n最终BMI分组和最佳NIPT时点:")
result_df = best_timing.copy()
result_df['BMI区间'] = result_df['BMI分组'].apply(
lambda x: '[20,28)' if x == '20-28' else
'[28,32)' if x == '28-32' else
'[32,36)' if x == '32-36' else
'[36,40)' if x == '36-40' else
'40+')
result_df = result_df[['BMI分组', 'BMI区间', '最佳NIPT时点(周)', '风险等级']]
print(result_df)
这段代码在运行后出现以下问题
C:\Study\anaconda3\envs\MCM\python.exe C:\Study\MCM\2.1.py
男胎数据量: 1082
各BMI分组数据量:
BMI分组
28-32 539
32-36 412
36-40 93
20-28 19
40+ 18
Name: count, dtype: int64
各BMI分组达标孕周的描述统计:
Empty DataFrame
Columns: [count, mean, std, min, 5%, 10%, 25%, 50%, 75%, 90%, 95%, max]
Index: []
各BMI分组的最佳NIPT时点(第95百分位数):
BMI分组 最佳NIPT时点(周)
0 20-28 NaN
1 28-32 NaN
2 32-36 NaN
3 36-40 NaN
4 40+ NaN
最佳NIPT时点的风险等级:
BMI分组 最佳NIPT时点(周) 风险等级
0 20-28 NaN 极高风险
1 28-32 NaN 极高风险
2 32-36 NaN 极高风险
3 36-40 NaN 极高风险
4 40+ NaN 极高风险
检测误差对最佳NIPT时点的影响:
Empty DataFrame
Columns: []
Index: []
Traceback (most recent call last):
File "C:\Study\MCM\2.1.py", line 126, in <module>
plt.bar(x_pos, error_df['原始最佳时点'], width=0.4, label='原始最佳时点')
File "C:\Study\anaconda3\envs\MCM\lib\site-packages\pandas\core\frame.py", line 4102, in __getitem__
indexer = self.columns.get_loc(key)
File "C:\Study\anaconda3\envs\MCM\lib\site-packages\pandas\core\indexes\range.py", line 417, in get_loc
raise KeyError(key)
KeyError: '原始最佳时点'
进程已结束,退出代码为 1