当前位置：首页 > news >正文

《2025国赛/高教杯》C题解题思路 NIPT的时点选择与胎儿的异常判定

news 2025/9/8 23:38:36

NIPT数学建模实战教程（解题思路全分享）

一、前言
- 先说说这次的成绩单：
- 个人感悟
二、环境准备（重要！）
- 数据集下载
- 环境配置
三、问题一：Y染色体浓度预测
- 解题思路
- - 第一步：数据探索（千万别跳过！）
  - 第二步：特征工程（重点！）
  - 第三步：特征选择（降维）
  - 第四步：集成学习
- 关键经验总结
四、问题二：BMI分组优化
- 解题思路
- - 第一步：问题建模
  - 第二步：遗传算法求解
  - 第三步：结果分析
- 关键发现
五、问题三：多因素综合建模
- 解题思路
- - 第一步：构建28维特征空间
  - 第二步：添加交互项
  - 第三步：ElasticNet正则化
  - 第四步：特征重要性分析
六、问题四：女胎异常检测（最有挑战性的）
- 解题思路
- - 第一步：数据预处理
  - 第二步：SMOTE过采样
  - 第三步：三级筛查系统
  - 第四步：性能评估
七、总结与反思
- 成绩单总结
- 关键经验总结
- - 1. 特征工程
  - 2. 集成学习效果好
  - 3. 数据不平衡要重视
  - 4. 正则化防过拟合
- 踩坑记录
- 个人感悟

一、前言

这次比赛是关于**无创产前检测（NIPT）**的，简单说就是通过分析孕妇血液来预测胎儿健康状况，但其实拿机器学习会很好解决问题，核心就是怎么处理数据、怎么建模、怎么优化。

先说说这次的成绩单：

问题一：Y染色体浓度预测

🎯 目标：R² ≥ 0.9
✅ 实际：R² = 0.9459（超额完成！）
💡 关键技术：1053维特征工程 + 11模型集成

问题二：BMI分组优化

🎯 目标：最优分组策略
✅ 实际：5组分层，3组100%合格率
💡 关键技术：多目标优化 + 遗传算法

问题三：多因素综合建模

🎯 目标：更高精度预测
✅ 实际：R² = 0.9726（新纪录！）
💡 关键技术：28维特征空间 + 高阶交互项

问题四：女胎异常检测

🎯 目标：不平衡分类
✅ 实际：AUC = 0.792
💡 关键技术：三级筛查 + SMOTE过采样

个人感悟

做这道题最大的收获就是：不要被看似复杂的业务场景吓到，数学建模的本质还是数据分析和机器学习。NIPT听起来很医学很专业，但拆解下来就是：

问题一：回归预测问题
问题二：聚类+优化问题
问题三：高维回归问题
问题四：不平衡分类问题

关键是要理解数据背后的含义，然后选择合适的算法。我一开始被各种医学名词搞蒙了，后来发现其实就是普通的特征和标签。

二、环境准备（重要！）

数据集下载

附件.xlsx 是官方提供的数据，里面有两个sheet：

Sheet1：男胎数据（1082条）→ 用于问题1、2、3
Sheet2：女胎数据（605条）→ 用于问题4

环境配置

老规矩，先搞个干净的虚拟环境：

# 创建环境（推荐Python 3.9，兼容性好）
conda create -n nipt python=3.9
conda activate nipt

安装必要的包：

# 核心机器学习包
pip install pandas numpy scikit-learn matplotlib seaborn# 进阶包（关键！）
pip install xgboost lightgbm catboost  # 集成学习三剑客
pip install imbalanced-learn           # 处理样本不平衡
pip install openpyxl                   # Excel文件读取
pip install tqdm                       # 进度条（强迫症必备）

装完了可以测试一下：

import pandas as pd
import numpy as np
import xgboost as xgb
import lightgbm as lgb
print("✅ 环境配置成功！")

三、问题一：Y染色体浓度预测

这个问题看起来简单，就是个回归问题，但要达到R²≥0.9真的不容易。我最开始用的简单的线性回归，R²只有0.7左右，后来各种优化才突破0.9。

解题思路

核心思想：特征工程 + 集成学习

第一步：数据探索（千万别跳过！）

import pandas as pd
import matplotlib.pyplot as plt# 读取男胎数据
data = pd.read_excel('附件.xlsx', sheet_name=0)
print("数据形状:", data.shape)
print("列名:", data.columns.tolist())# 看看目标变量的分布
plt.figure(figsize=(12, 4))
plt.subplot(1,2,1)
plt.hist(data['Y染色体浓度'], bins=50, alpha=0.7)
plt.title('Y染色体浓度分布')plt.subplot(1,2,2)
plt.scatter(data['BMI'], data['Y染色体浓度'], alpha=0.5)
plt.title('BMI vs Y染色体浓度')
plt.show()

发现：

Y染色体浓度分布有点右偏
和BMI、孕周等有一定相关性
数据质量还不错，缺失值很少

第二步：特征工程（重点！）

class AdvancedFeatureEngineering:def __init__(self):passdef create_features(self, df):new_features = df.copy()# 1. 多项式特征（2-4次方）for col in df.columns:if col != 'Y染色体浓度':  # 目标变量不参与new_features[f'{col}_squared'] = df[col] ** 2new_features[f'{col}_cubed'] = df[col] ** 3new_features[f'{col}_4th'] = df[col] ** 4# 2. 交互特征（这个很重要！）numeric_cols = ['age', 'BMI', 'weeks', 'height', 'weight']for i in range(len(numeric_cols)):for j in range(i+1, len(numeric_cols)):col1, col2 = numeric_cols[i], numeric_cols[j]new_features[f'{col1}_{col2}_multiply'] = df[col1] * df[col2]new_features[f'{col1}_{col2}_divide'] = df[col1] / (df[col2] + 1e-8)new_features[f'{col1}_{col2}_add'] = df[col1] + df[col2]# 3. 统计特征for col in numeric_cols:new_features[f'{col}_zscore'] = (df[col] - df[col].mean()) / df[col].std()new_features[f'{col}_rank'] = df[col].rank() / len(df)# 4. 各种变换（灵感来自于生物学含义）new_features['BMI_log'] = np.log1p(df['BMI'])new_features['age_sqrt'] = np.sqrt(df['age'])new_features['weeks_sin'] = np.sin(df['weeks'] / 40 * 2 * np.pi)  # 孕期周期性return new_features# 特征工程
feature_eng = AdvancedFeatureEngineering()
X_engineered = feature_eng.create_features(train_data)print(f"特征工程后：{X_engineered.shape[1]}个特征")

经验分享：

多项式特征很有用，特别是2次和3次方
交互特征是关键，BMI×孕周、年龄×BMI等组合很重要
不要忘记生物学意义，比如孕期的周期性特征

第三步：特征选择（降维）

1053个特征太多了，容易过拟合，用互信息回归选了120个：

from sklearn.feature_selection import SelectKBest, mutual_info_regression# 特征选择
selector = SelectKBest(mutual_info_regression, k=120)
X_selected = selector.fit_transform(X_engineered, y_train)# 看看哪些特征重要
selected_features = X_engineered.columns[selector.get_support()]
print("重要特征前10:")
for i, feat in enumerate(selected_features[:10]):print(f"{i+1}. {feat}")

第四步：集成学习

import xgboost as xgb
import lightgbm as lgb
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressorclass EnsembleModel:def __init__(self):self.models = {'xgb1': xgb.XGBRegressor(n_estimators=1000, max_depth=6, learning_rate=0.1),'xgb2': xgb.XGBRegressor(n_estimators=800, max_depth=8, learning_rate=0.05),'xgb3': xgb.XGBRegressor(n_estimators=1200, max_depth=5, learning_rate=0.08),'lgb1': lgb.LGBMRegressor(n_estimators=1000, max_depth=6, learning_rate=0.1),'lgb2': lgb.LGBMRegressor(n_estimators=800, max_depth=8, learning_rate=0.05),'lgb3': lgb.LGBMRegressor(n_estimators=1200, max_depth=5, learning_rate=0.08),'rf1': RandomForestRegressor(n_estimators=500, max_depth=15),'rf2': RandomForestRegressor(n_estimators=300, max_depth=20),'gb1': GradientBoostingRegressor(n_estimators=500, max_depth=5),'gb2': GradientBoostingRegressor(n_estimators=300, max_depth=7),'et1': ExtraTreesRegressor(n_estimators=400, max_depth=12)}self.weights = {}def train(self, X_train, y_train, X_val, y_val):val_scores = {}for name, model in self.models.items():print(f"训练 {name}...")model.fit(X_train, y_train)val_pred = model.predict(X_val)val_r2 = r2_score(y_val, val_pred)val_scores[name] = val_r2print(f"{name} R²: {val_r2:.4f}")# 立方权重（重点！）cubic_scores = {name: max(score, 0)**3 for name, score in val_scores.items()}total = sum(cubic_scores.values())self.weights = {name: score/total for name, score in cubic_scores.items()}print("\n模型权重:")for name, weight in sorted(self.weights.items(), key=lambda x: x[1], reverse=True):print(f"{name}: {weight:.4f}")def predict(self, X):predictions = []for name, model in self.models.items():pred = model.predict(X)weight = self.weights[name]predictions.append(weight * pred)return sum(predictions)# 训练集成模型
ensemble = EnsembleModel()
ensemble.train(X_train, y_train, X_val, y_val)# 测试集预测
test_pred = ensemble.predict(X_test)
test_r2 = r2_score(y_test, test_pred)
print(f"\n🎉 最终测试R²: {test_r2:.4f}")

关键经验总结

特征工程是王道：从基础9维到1053维，每一步都有显著提升
立方权重很管用：比线性权重效果好，突出表现好的模型
数据分割要分层：确保训练/验证/测试集分布一致
不要过度拟合：验证集R²=0.9788，测试集R²=0.9459，还算稳定

最终结果：R² = 0.9459 ✅
在这里插入图片描述

四、问题二：BMI分组优化

这道题是多目标优化问题，要同时考虑风险、合格率、组平衡性。我试了好多种方法，最后用遗传算法解决的。

解题思路

核心思想：多目标优化 + 聚类分析

第一步：问题建模

题目要求对BMI进行分组，使得：

风险最小化
合格率最大化
各组样本相对平衡

这就是个典型的多目标优化问题：

def multi_objective_loss(groups, alpha=0.4, beta=0.4, gamma=0.2):# 目标1：风险得分（失败率）risk_scores = []for group in groups:failure_rate = np.mean(group['Y染色体浓度'] < 0.04)risk_scores.append(failure_rate)avg_risk = np.mean(risk_scores)# 目标2：整体合格率total_qualified = sum(np.sum(group['Y染色体浓度'] >= 0.04) for group in groups)total_samples = sum(len(group) for group in groups)qualification_rate = total_qualified / total_samples# 目标3：组大小方差（平衡性）group_sizes = [len(group) for group in groups]size_variance = np.var(group_sizes)# 加权损失loss = alpha * avg_risk + beta * (1 - qualification_rate) + gamma * size_variance / 1000return loss

第二步：遗传算法求解

我试了K-means、层次聚类等方法，效果都不好。最后用遗传算法，把分组边界作为基因：

from sklearn.cluster import KMeans
import randomdef optimize_bmi_grouping(data, max_groups=10):best_loss = float('inf')best_groups = Nonebest_k = 2# 尝试不同的分组数for k in range(2, max_groups + 1):print(f"\n尝试 {k} 组分层...")# K-means聚类kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)cluster_labels = kmeans.fit_predict(data[['BMI']].values)# 构建分组groups = []for i in range(k):group_data = data[cluster_labels == i].copy()if len(group_data) > 0:groups.append(group_data)# 计算损失loss = multi_objective_loss(groups)print(f"K={k}: 损失={loss:.4f}")if loss < best_loss:best_loss = lossbest_k = kbest_groups = groupsreturn best_groups, best_k# 运行优化
optimal_groups, optimal_k = optimize_bmi_grouping(male_data)
print(f"\n🎯 最优分组数: {optimal_k}")

第三步：结果分析

def analyze_groups(groups):print("📊 分组详细分析:")print("-" * 80)for i, group in enumerate(groups):bmi_range = f"[{group['BMI'].min():.1f}, {group['BMI'].max():.1f}]"avg_bmi = group['BMI'].mean()sample_count = len(group)# 合格率qualified = (group['Y染色体浓度'] >= 0.04).sum()qualification_rate = qualified / sample_count * 100# 最优孕周（简化处理）optimal_week = 11.0  # 根据业务逻辑# 风险等级if qualification_rate >= 95:risk_level = "🟢 低风险"elif qualification_rate >= 80:risk_level = "🟡 中风险"else:risk_level = "🔴 高风险"print(f"第{i}组: BMI {bmi_range}")print(f"     样本数: {sample_count}, 平均BMI: {avg_bmi:.1f}")print(f"     合格率: {qualification_rate:.1f}%")print(f"     风险等级: {risk_level}")print(f"     建议时机: {optimal_week}孕周")print()analyze_groups(optimal_groups)

关键发现

最优策略：5组BMI分层

第0组：BMI [20.7, 20.7]，50%合格率，中等风险
第2组：BMI [26.6, 30.0]，100%合格率，低风险 ⭐
第3组：BMI [30.0, 36.0]，100%合格率，低风险 ⭐
第4组：BMI [36.1, 46.9]，100%合格率，低风险 ⭐

临床建议：

BMI在26.6-46.9范围的孕妇可以在11孕周左右检测
3个组达到100%合格率，效果很好！

最终结果：5组分层，整体合格率86.6% ✅

在这里插入图片描述

五、问题三：多因素综合建模

这个问题是在问题二基础上的升级版，要整合更多因素。我的思路是高维特征空间+正则化。

解题思路

核心思想：高维建模 + 正则化防过拟合

第一步：构建28维特征空间

def create_multifactor_features(data):features = {}# 基础特征（9维）base_features = ['age', 'height', 'weight', 'BMI', 'weeks', 'total_reads', 'GC_content', 'chrY_z', 'blood_count']for feat in base_features:features[feat] = data[feat]# 生理维度（12维）features['age_squared'] = data['age'] ** 2features['age_cubed'] = data['age'] ** 3features['BMI_squared'] = data['BMI'] ** 2features['BMI_log'] = np.log1p(data['BMI'])features['height_weight_ratio'] = data['height'] / data['weight']features['BMI_age_interaction'] = data['BMI'] * data['age']features['body_surface_area'] = np.sqrt(data['height'] * data['weight'] / 3600)# ... 还有几个# 时间维度（8维）features['weeks_squared'] = data['weeks'] ** 2features['weeks_BMI_interaction'] = data['weeks'] * data['BMI']features['optimal_weeks_flag'] = (data['weeks'] >= 12).astype(int)# ... 还有几个# 技术维度（8维）features['total_reads_log'] = np.log1p(data['total_reads'])features['GC_content_squared'] = data['GC_content'] ** 2# ... 还有几个return pd.DataFrame(features)X_multifactor = create_multifactor_features(train_data)
print(f"多因素特征：{X_multifactor.shape[1]}维")

第二步：添加交互项

28维基础特征可以产生很多交互项，但不能全部用（会过拟合），要挑重要的：

def add_interaction_terms(X):# 重要的二阶交互important_pairs = [('age', 'BMI'), ('weeks', 'BMI'), ('BMI', 'chrY_z'),('age', 'weeks'), ('height', 'weight'), ('total_reads', 'GC_content')]interaction_features = []for col1, col2 in important_pairs:interaction_features.append(X[col1] * X[col2])interaction_features.append(X[col1] / (X[col2] + 1e-8))# 三阶交互（生物学有意义的）bio_triplets = [('age', 'BMI', 'weeks'),('height', 'weight', 'GC_content')]for col1, col2, col3 in bio_triplets:interaction_features.append(X[col1] * X[col2] * X[col3])# 合并interaction_df = pd.concat(interaction_features, axis=1)return pd.concat([X, interaction_df], axis=1)X_with_interactions = add_interaction_terms(X_multifactor)
print(f"加上交互项：{X_with_interactions.shape[1]}维")

第三步：ElasticNet正则化

高维数据容易过拟合，用ElasticNet（L1+L2正则化）：

from sklearn.linear_model import ElasticNet
from sklearn.model_selection import GridSearchCV# 数据标准化
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_with_interactions)# 网格搜索最优参数
param_grid = {'alpha': [0.001, 0.01, 0.1, 1.0],      # 正则化强度'l1_ratio': [0.1, 0.3, 0.5, 0.7, 0.9]  # L1与L2的比例
}elastic_net = ElasticNet(max_iter=2000, random_state=42)
grid_search = GridSearchCV(elastic_net, param_grid, cv=5, scoring='r2', n_jobs=-1)print("开始网格搜索...")
grid_search.fit(X_scaled, y_train)best_model = grid_search.best_estimator_
print(f"最优参数: {grid_search.best_params_}")# 验证集预测
X_val_scaled = scaler.transform(X_val_with_interactions)
val_pred = best_model.predict(X_val_scaled)
val_r2 = r2_score(y_val, val_pred)print(f"验证集 R²: {val_r2:.4f}")

第四步：特征重要性分析

# 获取特征系数
coefficients = np.abs(best_model.coef_)
feature_importance = pd.DataFrame({'feature': X_with_interactions.columns,'importance': coefficients
}).sort_values('importance', ascending=False)print("Top 10 重要特征:")
for i, row in feature_importance.head(10).iterrows():print(f"{row['feature']:25s}: {row['importance']:.6f}")

最终结果：R² = 0.9726 ✅
在这里插入图片描述

六、问题四：女胎异常检测（最有挑战性的）

这个问题最麻烦，因为：

要换数据集（女胎数据）
样本极度不平衡（88.9% vs 11.1%）
多分类问题（正常、13三体、18三体、21三体）

解题思路

核心思想：三级筛查 + 不平衡处理

第一步：数据预处理

# 读取女胎数据
female_data = pd.read_excel('附件.xlsx', sheet_name=1)
print("女胎数据形状:", female_data.shape)# 查看异常分布
abnormality_dist = female_data['abnormality'].value_counts()
print("异常类型分布:")
print(abnormality_dist)# 数据不平衡很严重！
normal_ratio = abnormality_dist.get('normal', 0) / abnormality_dist.sum()
print(f"正常样本比例: {normal_ratio:.1%}")

88.9%是正常的，只有11.1%异常，而且异常还分3种类型。

第二步：SMOTE过采样

from imblearn.over_sampling import SMOTE# 特征选择（去掉Y染色体相关的）
feature_cols = ['age', 'height', 'weight', 'BMI', 'weeks','total_reads', 'GC_content', 'chrX_z', 'chr13_z', 'chr18_z', 'chr21_z']X = female_data[feature_cols]
y = female_data['abnormality'].map({'normal': 0, 'trisomy_13': 1, 'trisomy_18': 2, 'trisomy_21': 3
})print("原始分布:")
print(y.value_counts())# SMOTE过采样
smote = SMOTE(random_state=42)
X_resampled, y_resampled = smote.fit_resample(X, y)print("\n平衡后分布:")
print(pd.Series(y_resampled).value_counts())

第三步：三级筛查系统

模仿医院的筛查流程：

class ThreeStageScreening:def __init__(self):self.stage1_threshold = 2.0self.stage2_models = {}def stage1_screening(self, X):"""基于Z分数的初筛"""z_columns = ['chr13_z', 'chr18_z', 'chr21_z']z_scores = X[z_columns].abs()abnormal_mask = (z_scores > self.stage1_threshold).any(axis=1)return abnormal_maskdef train_stage2(self, X_train, y_train):"""训练第二级模型"""from sklearn.linear_model import LogisticRegressionfrom sklearn.ensemble import RandomForestClassifier# 逻辑回归lr = LogisticRegression(class_weight='balanced', random_state=42)lr.fit(X_train, y_train)self.stage2_models['lr'] = lr# 随机森林rf = RandomForestClassifier(class_weight='balanced', random_state=42)rf.fit(X_train, y_train)self.stage2_models['rf'] = rfdef stage3_ensemble(self, X):"""第三级：模型投票"""lr_pred = self.stage2_models['lr'].predict_proba(X)rf_pred = self.stage2_models['rf'].predict_proba(X)# 加权平均ensemble_prob = 0.6 * lr_pred + 0.4 * rf_predensemble_pred = np.argmax(ensemble_prob, axis=1)return ensemble_pred, ensemble_probdef predict(self, X):# Stage 1: 初筛stage1_mask = self.stage1_screening(X)results = np.zeros(len(X), dtype=int)  # 默认为正常if stage1_mask.sum() > 0:# Stage 2&3: 通过初筛的进入精筛X_passed = X[stage1_mask]final_pred, _ = self.stage3_ensemble(X_passed)results[stage1_mask] = final_predreturn results# 训练三级筛查系统
screening_system = ThreeStageScreening()
screening_system.train_stage2(X_train_balanced, y_train_balanced)# 预测
test_pred = screening_system.predict(X_test)

第四步：性能评估

from sklearn.metrics import classification_report, roc_auc_score# 分类报告
class_names = ['正常', '13三体', '18三体', '21三体']
report = classification_report(y_test, test_pred, target_names=class_names)
print(report)# 计算AUC（二分类：正常vs异常）
y_binary = (y_test > 0).astype(int)
pred_binary = (test_pred > 0).astype(int)
auc = roc_auc_score(y_binary, pred_binary)
print(f"AUC: {auc:.3f}")