import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 读取处理后的数据
df = pd.read_excel('分层去重后数据.xlsx')
# 筛选中孕期数据 (15-20周)
mid_pregnancy_data = df[(df['检测孕周'] >= 15) & (df['检测孕周'] <= 20)].copy()
print(f"中孕期数据记录数: {len(mid_pregnancy_data)}")
# 场景1: 静态关联模型 - 中孕期数据
print("="*50)
print("场景1: 静态关联模型 (中孕期数据)")
print("="*50)
# 检查Y染色体浓度数据
print("Y染色体浓度描述统计:")
print(mid_pregnancy_data['Y染色体浓度'].describe())
# 检查是否有异常值或缺失值
print(f"Y染色体浓度缺失值数量: {mid_pregnancy_data['Y染色体浓度'].isnull().sum()}")
# 移除缺失值
mid_pregnancy_data = mid_pregnancy_data.dropna(subset=['Y染色体浓度', '检测孕周', '孕妇BMI'])
# 绘制散点图矩阵查看变量关系
sns.pairplot(mid_pregnancy_data[['Y染色体浓度', '检测孕周', '孕妇BMI', '年龄']])
plt.suptitle('变量关系散点图矩阵', y=1.02)
plt.show()
# 计算相关系数矩阵
corr_matrix = mid_pregnancy_data[['Y染色体浓度', '检测孕周', '孕妇BMI', '年龄']].corr()
print("相关系数矩阵:")
print(corr_matrix)
# 绘制热图 - 修正颜色映射名称
plt.figure(figsize=(8, 6))
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', center=0) # 修正为 'coolwarm'
plt.title('变量相关系数热图')
plt.show()
# 建立多元线性回归模型
# 准备自变量和因变量
X = mid_pregnancy_data[['检测孕周', '孕妇BMI', '年龄']]
X = sm.add_constant(X) # 添加常数项
y = mid_pregnancy_data['Y染色体浓度']
# 拟合模型
model = sm.OLS(y, X).fit()
# 输出模型摘要
print("多元线性回归模型摘要:")
print(model.summary())
# 绘制残差图
plt.figure(figsize=(10, 6))
plt.scatter(model.fittedvalues, model.resid)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('预测值')
plt.ylabel('残差')
plt.title('残差图')
plt.show()
# 绘制Q-Q图检验残差正态性
plt.figure(figsize=(8, 6))
sm.qqplot(model.resid, line='s')
plt.title('Q-Q图 (残差正态性检验)')
plt.show()
# 场景2: 动态关联模型 - 混合效应模型
print("="*50)
print("场景2: 动态关联模型 (混合效应模型)")
print("="*50)
# 筛选有多个孕周记录的孕妇
pregnancy_counts = df['孕妇代码'].value_counts()
multi_measure_patients = pregnancy_counts[pregnancy_counts > 1].index
longitudinal_data = df[df['孕妇代码'].isin(multi_measure_patients)].copy()
print(f"有多条记录的孕妇数量: {len(multi_measure_patients)}")
print(f"纵向数据集记录数: {len(longitudinal_data)}")
# 检查每个孕妇的记录数分布
print("每个孕妇记录数分布:")
print(longitudinal_data['孕妇代码'].value_counts().value_counts())
# 建立混合效应模型
# 将孕妇代码转换为分类变量
longitudinal_data['孕妇代码'] = longitudinal_data['孕妇代码'].astype('category')
# 拟合混合效应模型
mixed_model = smf.mixedlm("Y染色体浓度 ~ 检测孕周 + 孕妇BMI + 年龄",
longitudinal_data,
groups=longitudinal_data["孕妇代码"]).fit()
# 输出模型摘要
print("混合效应模型摘要:")
print(mixed_model.summary())
# 可视化模型结果
# 1. 绘制每个孕妇的Y染色体浓度随孕周的变化
plt.figure(figsize=(12, 8))
for patient in longitudinal_data['孕妇代码'].unique()[:20]: # 只显示前20个孕妇
patient_data = longitudinal_data[longitudinal_data['孕妇代码'] == patient]
plt.plot(patient_data['检测孕周'], patient_data['Y染色体浓度'], 'o-', label=patient, alpha=0.7)
plt.xlabel('检测孕周')
plt.ylabel('Y染色体浓度')
plt.title('Y染色体浓度随孕周变化趋势 (前20个孕妇)')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()
# 2. 绘制Y染色体浓度与BMI的关系
plt.figure(figsize=(10, 6))
sns.scatterplot(data=longitudinal_data, x='孕妇BMI', y='Y染色体浓度', hue='检测孕周', palette='viridis')
plt.title('Y染色体浓度与BMI的关系')
plt.show()
# 3. 绘制模型预测值与实际值对比
longitudinal_data['预测值'] = mixed_model.predict(longitudinal_data)
plt.figure(figsize=(10, 6))
plt.scatter(longitudinal_data['Y染色体浓度'], longitudinal_data['预测值'])
plt.plot([longitudinal_data['Y染色体浓度'].min(), longitudinal_data['Y染色体浓度'].max()],
[longitudinal_data['Y染色体浓度'].min(), longitudinal_data['Y染色体浓度'].max()], 'r--')
plt.xlabel('实际值')
plt.ylabel('预测值')
plt.title('混合效应模型预测值与实际值对比')
plt.show()
# 模型诊断 - 残差分析
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(mixed_model.fittedvalues, mixed_model.resid)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('预测值')
plt.ylabel('残差')
plt.title('混合效应模型残差图')
plt.subplot(1, 2, 2)
sm.qqplot(mixed_model.resid, line='s', ax=plt.gca())
plt.title('混合效应模型Q-Q图')
plt.tight_layout()
plt.show()
# 输出模型系数和显著性
print("="*50)
print("模型系数和显著性")
print("="*50)
print("多元线性回归模型系数:")
print(model.params)
print("\n多元线性回归模型p值:")
print(model.pvalues)
print("\n混合效应模型固定效应系数:")
print(mixed_model.fe_params)
print("\n混合效应模型固定效应p值:")
print(mixed_model.pvalues)
# 保存模型结果
results_df = pd.DataFrame({
'模型': ['多元线性回归', '混合效应'],
'R²': [model.rsquared, mixed_model.rsquared],
'调整R²': [model.rsquared_adj, None],
'AIC': [model.aic, mixed_model.aic],
'BIC': [model.bic, mixed_model.bic]
})
results_df.to_excel('模型比较结果.xlsx', index=False)
print("模型比较结果已保存到 '模型比较结果.xlsx'")
# 保存预测结果
predictions = longitudinal_data[['孕妇代码', '检测孕周', '孕妇BMI', '年龄', 'Y染色体浓度', '预测值']]
predictions.to_excel('模型预测结果.xlsx', index=False)
print("预测结果已保存到 '模型预测结果.xlsx'")
扫一扫
2513
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
