问题1：Y染色体浓度与孕周数、BMI的关系分析

问题重述

分析胎儿Y染色体浓度与孕妇的孕周数、BMI等指标之间的相关性，建立关系模型，并进行显著性检验。

符号定义

• V ：Y染色体浓度（因变量）

• J ：孕周数（自变量）

• K ：BMI（自变量）

• ε ：随机误差项

模型建立

采用多元线性回归模型：

算法与步骤

1. 数据预处理 ：清洗数据，剔除无效值、异常值。

2. 相关性分析 ：计算Pearson相关系数矩阵。

3. 回归建模 ：使用最小二乘法估计参数

4. 显著性检验 ：

• 使用t检验判断每个自变量的显著性（p-value

• 使用F检验判断整体模型的显著性。

模型诊断 ：残差分析、共线性检验（VIF）。

问题2：男胎孕妇BMI分组与最佳NIPT时点

问题重述

根据BMI对男胎孕妇进行分组，确定每组的BMI区间和最佳NIPT时点，使得潜在风险最小，并分析检测误差的影响。

符号定义

模型建立

1. 分组方法 ：使用K-means聚类或基于分位数的分组。

2. 最佳时点确定 ：每组中

3. 风险模型 ：

• 早期发现（≤12周）：风险低

• 中期发现（13–27周）：风险高

• 晚期发现（≥28周）：风险极高

算法与步骤

1. 提取男胎数据，计算每个样本的

2. 根据BMI进行聚类（如K-means），确定分组。

3. 对每组计算

4. 以中位数作为该组的最佳NIPT时点。

5. 分析检测误差（如测量误差、模型误差）对时点选择的影响（敏感性分析）。

问题3：多因素影响下的分组与时点选择

问题重述

综合考虑BMI、年龄、体重等因素，给出男胎孕妇的合理分组和最佳NIPT时点，使得风险最小，并分析误差影响。

符号定义

模型建立

使用 逻辑回归 或 Cox比例风险模型 预测达标概率与时间：

或使用 生存分析模型 ：

算法与步骤

1. 提取男胎数据，构建特征矩阵（BMI、年龄、体重等）。

2. 使用逻辑回归或Cox模型建模。

3. 根据模型输出进行风险分层（如低、中、高风险组）。

4. 为每组确定最佳NIPT时点（如风险最低对应的孕周）。

5. 进行误差传播分析（如Bootstrap法）。

问题4：女胎异常判定方法

问题重述

基于X染色体、21/18/13号染色体的Z值、GC含量、读段数等指标，构建女胎异常的判定模型。

符号定义

模型建立

使用 逻辑回归 或 支持向量机（SVM） 构建分类模型：

算法与步骤

1. 提取女胎数据，标注异常标签（AB列或AE列）。

2. 特征工程：标准化Z值、GC含量、读段数等。

3. 使用逻辑回归、SVM或随机森林进行建模。

4. 交叉验证评估模型性能（准确率、召回率、AUC）。

5. 输出判定规则（如Z值阈值、GC异常范围等）。

准备工作：导入必要库和数据

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport seaborn as snsfrom sklearn.linear_model import LinearRegression, LogisticRegressionfrom sklearn.cluster import KMeansfrom sklearn.model_selection import train_test_splitfrom sklearn.metrics import accuracy_score, confusion_matrix, classification_reportfrom scipy import statsimport statsmodels.api as smfrom statsmodels.formula.api import olsimport warningswarnings.filterwarnings('ignore')
# 设置中文字体plt.rcParams['font.sans-serif'] = ['SimHei']plt.rcParams['axes.unicode_minus'] = False
# 读取数据male_data = pd.read_excel('附件.xlsx', sheet_name='男胎检测数据')female_data = pd.read_excel('附件.xlsx', sheet_name='女胎检测数据')
# 数据预处理# 将孕周转换为数值格式def convert_gestational_week(week_str): if pd.isna(week_str): return np.nan if 'w' in week_str: parts = week_str.split('w') weeks = int(parts[0]) if '+' in parts[1]: days = int(parts[1].split('+')[1]) return weeks + days/7 return weeks return float(week_str)
male_data['检测孕周数值'] = male_data['检测孕周'].apply(convert_gestational_week)female_data['检测孕周数值'] = female_data['检测孕周'].apply(convert_gestational_week)
# 过滤掉Y染色体浓度为空的记录male_data = male_data.dropna(subset=['Y染色体浓度'])

问题1：Y染色体浓度与孕周数、BMI的关系分析

# 问题1: Y染色体浓度与孕周数、BMI的关系分析print("="*50)print("问题1: Y染色体浓度与孕周数、BMI的关系分析")print("="*50)
# 准备数据X = male_data[['检测孕周数值', '孕妇BMI']].dropnay = male_data.loc[X.index, 'Y染色体浓度']
# 多元线性回归X_with_const = sm.add_constant(X)model = sm.OLS(y, X_with_const).fit
print("回归模型摘要:")print(model.summary)
# 绘制关系图fig, axes = plt.subplots(1, 2, figsize=(15, 5))
# 孕周与Y染色体浓度的关系axes[0].scatter(X['检测孕周数值'], y, alpha=0.6)z = np.polyfit(X['检测孕周数值'], y, 1)p = np.poly1d(z)axes[0].plot(X['检测孕周数值'], p(X['检测孕周数值']), "r--", alpha=0.8)axes[0].set_xlabel('孕周')axes[0].set_ylabel('Y染色体浓度')axes[0].set_title('Y染色体浓度与孕周的关系')
# BMI与Y染色体浓度的关系axes[1].scatter(X['孕妇BMI'], y, alpha=0.6)z = np.polyfit(X['孕妇BMI'], y, 1)p = np.poly1d(z)axes[1].plot(X['孕妇BMI'], p(X['孕妇BMI']), "r--", alpha=0.8)axes[1].set_xlabel('BMI')axes[1].set_ylabel('Y染色体浓度')axes[1].set_title('Y染色体浓度与BMI的关系')
plt.tight_layoutplt.savefig('问题1_关系图.png', dpi=300)plt.show
# 计算相关系数corr_week = np.corrcoef(X['检测孕周数值'], y)[0, 1]corr_bmi = np.corrcoef(X['孕妇BMI'], y)[0, 1]
print(f"\nY染色体浓度与孕周的相关系数: {corr_week:.4f}")print(f"Y染色体浓度与BMI的相关系数: {corr_bmi:.4f}")
# 显著性检验_, p_week = stats.pearsonr(X['检测孕周数值'], y)_, p_bmi = stats.pearsonr(X['孕妇BMI'], y)
print(f"\nY染色体浓度与孕周相关性的p值: {p_week:.4e}")print(f"Y染色体浓度与BMI相关性的p值: {p_bmi:.4e}")
if p_week 0.05: print("Y染色体浓度与孕周的相关性在统计上显著")else: print("Y染色体浓度与孕周的相关性在统计上不显著")
if p_bmi 0.05: print("Y染色体浓度与BMI的相关性在统计上显著")else: print("Y染色体浓度与BMI的相关性在统计上不显著")

问题2：男胎孕妇BMI分组与最佳NIPT时点

# 问题2: 男胎孕妇BMI分组与最佳NIPT时点print("\n" + "="*50)print("问题2: 男胎孕妇BMI分组与最佳NIPT时点")print("="*50)
# 计算每个孕妇Y染色体浓度首次达到4%的时间def find_first_above_4(group): group = group.sort_values('检测孕周数值') above_4 = group[group['Y染色体浓度'] >= 0.04] if len(above_4) > 0: return above_4.iloc[0]['检测孕周数值'] return np.nan
first_above_4 = male_data.groupby('孕妇代码').apply(find_first_above_4).reset_indexfirst_above_4.columns = ['孕妇代码', '首次达标孕周']
# 合并数据male_summary = male_data.groupby('孕妇代码').agg({ '孕妇BMI': 'mean', '孕妇年龄': 'first', '孕妇体重': 'mean'}).reset_index
male_summary = pd.merge(male_summary, first_above_4, on='孕妇代码')male_summary = male_summary.dropna
# 基于BMI进行分组# 使用K-means聚类进行分组bmi_data = male_summary[['孕妇BMI']].valueskmeans = KMeans(n_clusters=3, random_state=42)male_summary['BMI分组'] = kmeans.fit_predict(bmi_data)
# 计算每组的BMI范围和最佳NIPT时点group_results = for group_id in range(3): group_data = male_summary[male_summary['BMI分组'] == group_id] bmi_min = group_data['孕妇BMI'].min bmi_max = group_data['孕妇BMI'].max best_time = group_data['首次达标孕周'].median
 group_results.append({ '分组': group_id + 1, 'BMI范围': f"{bmi_min:.1f}-{bmi_max:.1f}", '最佳NIPT时点(周)': f"{best_time:.1f}" })
 print(f"分组 {group_id+1 }: BMI范围 {bmi_min:.1f}-{bmi_max:.1f}, 最佳NIPT时点 {best_time:.1f}周")
# 转换为DataFrame并保存group_df = pd.DataFrame(group_results)print("\nBMI分组及最佳NIPT时点:")print(group_df)
# 分析检测误差对结果的影响# 使用Bootstrap方法估计最佳时点的置信区间bootstrap_results = for group_id in range(3): group_data = male_summary[male_summary['BMI分组'] == group_id] times = group_data['首次达标孕周'].values
 # 进行Bootstrap抽样 n_bootstraps = 1000 bootstrap_estimates = 
 for _ in range(n_bootstraps): sample = np.random.choice(times, size=len(times), replace=True) bootstrap_estimates.append(np.median(sample))
 # 计算95%置信区间 ci_low = np.percentile(bootstrap_estimates, 2.5) ci_high = np.percentile(bootstrap_estimates, 97.5)
 bootstrap_results.append({ '分组': group_id + 1, '最佳时点估计': np.median(times), '95%置信区间': f"{ci_low:.1f}-{ci_high:.1f}" })
 print(f"分组 {group_id+1 }: 最佳时点 {np.median(times):.1f}周, 95%置信区间 {ci_low:.1f}-{ci_high:.1f}周")
# 绘制分组结果plt.figure(figsize=(10, 6))colors = ['red', 'blue', 'green']for group_id in range(3): group_data = male_summary[male_summary['BMI分组'] == group_id] plt.scatter(group_data['孕妇BMI'], group_data['首次达标孕周'],  alpha=0.6, color=colors[group_id], label=f'分组 {group_id+1 }')
plt.xlabel('BMI')plt.ylabel('首次达标孕周')plt.title('BMI分组与Y染色体浓度首次达标时间')plt.legendplt.grid(True, alpha=0.3)plt.savefig('问题2_BMI分组.png', dpi=300)plt.show

问题3：多因素影响下的分组与时点选择

# 问题3: 多因素影响下的分组与时点选择print("\n" + "="*50)print("问题3: 多因素影响下的分组与时点选择")print("="*50)
# 准备数据features = male_summary[['孕妇BMI', '孕妇年龄', '孕妇体重']]target = male_summary['首次达标孕周']
# 使用多元线性回归分析各因素的影响X_with_const = sm.add_constant(features)model = sm.OLS(target, X_with_const).fit
print("多因素回归模型摘要:")print(model.summary)
# 基于多个特征进行聚类分组multi_features = male_summary[['孕妇BMI', '孕妇年龄', '孕妇体重']].values
# 标准化特征from sklearn.preprocessing import StandardScalerscaler = StandardScalerscaled_features = scaler.fit_transform(multi_features)
# 使用K-means进行聚类kmeans_multi = KMeans(n_clusters=3, random_state=42)male_summary['多因素分组'] = kmeans_multi.fit_predict(scaled_features)
# 计算每组的BMI范围和最佳NIPT时点multi_group_results = for group_id in range(3): group_data = male_summary[male_summary['多因素分组'] == group_id] bmi_min = group_data['孕妇BMI'].min bmi_max = group_data['孕妇BMI'].max age_min = group_data['孕妇年龄'].min age_max = group_data['孕妇年龄'].max weight_min = group_data['孕妇体重'].min weight_max = group_data['孕妇体重'].max best_time = group_data['首次达标孕周'].median
 multi_group_results.append({ '分组': group_id + 1, 'BMI范围': f"{bmi_min:.1f}-{bmi_max:.1f}", '年龄范围': f"{age_min:.0f}-{age_max:.0f}", '体重范围': f"{weight_min:.1f}-{weight_max:.1f}", '最佳NIPT时点(周)': f"{best_time:.1f}" })
 print(f"分组 {group_id+1 }: BMI范围 {bmi_min:.1f}-{bmi_max:.1f}, " f"年龄范围 {age_min:.0f}-{age_max:.0f}, " f"体重范围 {weight_min:.1f}-{weight_max:.1f}, " f"最佳NIPT时点 {best_time:.1f}周")
# 转换为DataFrame并保存multi_group_df = pd.DataFrame(multi_group_results)print("\n多因素分组及最佳NIPT时点:")print(multi_group_df)
# 分析检测误差对结果的影响# 计算达标比例def calculate_success_rate(group, threshold=0.04): total = len(group) success = len(group[group['Y染色体浓度'] >= threshold]) return success / total
# 计算每组的达标比例for group_id in range(3): group_data = male_summary[male_summary['多因素分组'] == group_id] patient_codes = group_data['孕妇代码'].unique group_records = male_data[male_data['孕妇代码'].isin(patient_codes)]
 success_rate = calculate_success_rate(group_records) print(f"分组 {group_id+1 } 的Y染色体浓度达标比例: {success_rate:.2%}")
# 误差传播分析 - 使用蒙特卡洛方法print("\n误差传播分析 (蒙特卡洛方法):")for group_id in range(3): group_data = male_summary[male_summary['多因素分组'] == group_id] times = group_data['首次达标孕周'].values
 # 假设测量误差为正态分布，标准差为1周 n_simulations = 1000 error_estimates = 
 for _ in range(n_simulations): # 添加随机误差 times_with_error = times + np.random.normal(0, 1, len(times)) error_estimates.append(np.median(times_with_error))
 # 计算误差影响 original_estimate = np.median(times) error_mean = np.mean(error_estimates) error_std = np.std(error_estimates)
 print(f"分组 {group_id+1 }: 原始估计 {original_estimate:.1f}周, " f"考虑误差后 {error_mean:.1f}±{error_std:.1f}周")

问题4：女胎异常判定方法

# 问题4: 女胎异常判定方法print("\n" + "="*50)print("问题4: 女胎异常判定方法")print("="*50)
# 准备女胎数据# 创建异常标签 - 使用AB列（染色体的非整倍体）或AE列（胎儿是否健康）female_data['异常标签'] = female_data['染色体的非整倍体'].notna | (female_data['胎儿是否健康'] == '否')female_data['异常标签'] = female_data['异常标签'].astype(int)
# 选择特征features = ['X染色体的Z值', '13号染色体的Z值', '18号染色体的Z值', '21号染色体的Z值', 'X染色体浓度', '13号染色体的GC含量', '18号染色体的GC含量', '21号染色体的GC含量', '孕妇BMI', 'GC含量', '被过滤掉读段数的比例']
X = female_data[features].dropnay = female_data.loc[X.index, '异常标签']
print(f"女胎数据中异常样本比例: {y.mean:.2%}")
# 划分训练集和测试集X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42, stratify=y)
# 训练逻辑回归模型model = LogisticRegression(max_iter=1000, class_weight='balanced')model.fit(X_train, y_train)
# 预测并评估模型y_pred = model.predict(X_test)accuracy = accuracy_score(y_test, y_pred)print(f"\n逻辑回归模型准确率: {accuracy:.2%}")
print("\n分类报告:")print(classification_report(y_test, y_pred))
print("混淆矩阵:")print(confusion_matrix(y_test, y_pred))
# 分析特征重要性feature_importance = pd.DataFrame({ '特征': features, '系数': model.coef_[0]}).sort_values('系数', key=abs, ascending=False)
print("\n特征重要性 (按系数绝对值排序):")print(feature_importance)
# 绘制特征重要性图plt.figure(figsize=(10, 6))plt.barh(feature_importance['特征'], feature_importance['系数'])plt.xlabel('系数值')plt.title('逻辑回归特征重要性')plt.tight_layoutplt.savefig('问题4_特征重要性.png', dpi=300)plt.show
# 建立判定规则# 基于Z值的判定规则print("\n基于Z值的判定规则:")z_threshold = 2.5 # 常用的Z值阈值for chrom in ['13', '18', '21']: z_col = f'{chrom}号染色体的Z值' abnormal_count = len(female_data[(female_data[z_col].abs > z_threshold) & (female_data['异常标签'] == 1)]) total_abnormal = len(female_data[female_data['异常标签'] == 1]) print(f"{chrom}号染色体Z值异常(|Z|>{z_threshold})的异常样本: {abnormal_count}/{total_abnormal} ({abnormal_count/total_abnormal:.1%})")
# 基于GC含量的判定规则print("\n基于GC含量的判定规则:")gc_low, gc_high = 0.4, 0.6 # 正常GC含量范围for chrom in ['13', '18', '21']: gc_col = f'{chrom}号染色体的GC含量' abnormal_count = len(female_data[((female_data[gc_col] gc_high)) & (female_data['异常标签'] == 1)]) total_abnormal = len(female_data[female_data['异常标签'] == 1]) print(f"{chrom}号染色体GC含量异常({gc_low}或>{gc_high})的异常样本: {abnormal_count}/{total_abnormal} ({abnormal_count/total_abnormal:.1%})")
# 综合判定方法def integrated_diagnosis(row, z_threshold=2.5, gc_low=0.4, gc_high=0.6): # 检查Z值异常 z_abnormal = any([ abs(row['13号染色体的Z值']) > z_threshold, abs(row['18号染色体的Z值']) > z_threshold, abs(row['21号染色体的Z值']) > z_threshold ])
 # 检查GC含量异常 gc_abnormal = any([ row['13号染色体的GC含量'] or row['13号染色体的GC含量'] > gc_high, row['18号染色体的GC含量'] or row['18号染色体的GC含量'] > gc_high, row['21号染色体的GC含量'] or row['21号染色体的GC含量'] > gc_high ])
 # 检查X染色体异常 x_abnormal = abs(row['X染色体的Z值']) > z_threshold or row['X染色体浓度'] 0.02
 # 综合判定 return int(z_abnormal or gc_abnormal or x_abnormal)
# 应用综合判定方法female_data['综合判定'] = female_data.apply(integrated_diagnosis, axis=1)integrated_accuracy = accuracy_score(female_data['异常标签'], female_data['综合判定'])print(f"\n综合判定方法的准确率: {integrated_accuracy:.2%}")
print("\n综合判定方法分类报告:")print(classification_report(female_data['异常标签'], female_data['综合判定']))

结果保存和总结

# 保存所有结果到Excel文件with pd.ExcelWriter('C题结果汇总.xlsx') as writer: group_df.to_excel(writer, sheet_name='问题2_BMI分组', index=False) multi_group_df.to_excel(writer, sheet_name='问题3_多因素分组', index=False) feature_importance.to_excel(writer, sheet_name='问题4_特征重要性', index=False)
 # 保存问题1的回归结果 problem1_results = pd.DataFrame({ '变量': ['const', '孕周', 'BMI'], '系数': model.params.values, 'p值': model.pvalues.values }) problem1_results.to_excel(writer, sheet_name='问题1_回归结果', index=False)
print("\n所有分析已完成，结果已保存到 'C题结果汇总.xlsx'")

使用说明

1. 确保已安装所需的Python库：pandas, numpy, matplotlib, seaborn, scikit-learn, statsmodels, scipy

2. 将代码保存为Python文件（如 model_c.py ）

3. 确保数据文件 附件.xlsx 在同一目录下

4. 运行代码： python model_c.py

5. 结果将显示在控制台，并保存为图像文件和Excel文件

注意事项

1. 代码中使用了多种统计和机器学习方法，可能需要根据实际数据特点进行调整

2. 对于问题4的女胎异常判定，可以根据实际需求调整判定规则的阈值

3. 如果数据量较大，可能需要优化代码性能

4. 建议在运行前备份原始数据