import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
import statsmodels.api as sm # 用于计算P值
# ----------------------
# 1. 生成模拟数据(贴合转运场景)
# ----------------------
np.random.seed(42) # 固定随机数,保证结果可复现
n = 500 # 样本量
data = pd.DataFrame({
# 自变量:候选影响因素
'GCS评分': np.random.randint(3, 16, size=n), # GCS评分3-15分,越低病情越重
'管道数量': np.random.randint(0, 6, size=n), # 0-5根管道
'转运时长': np.random.randint(5, 60, size=n), # 5-60分钟
'镇静程度': np.random.randint(0, 3, size=n), # 0=未镇静,1=浅镇静,2=深镇静
'生命体征波动': np.random.randint(0, 3, size=n), # 0=稳定,1=轻度,2=显著
# 因变量:病情变化(问题1核心)
# 模拟逻辑:GCS低、管道多、转运久、镇静不足、生命体征波动大→更易病情变化
'病情变化': np.where(
(np.random.rand(n) < 0.1) + # 基础概率
(data['GCS评分'] < 9) * 0.3 + # GCS<9(昏迷)增加30%概率
(data['管道数量'] > 2) * 0.2 + # 管道>2根增加20%概率
(data['转运时长'] > 30) * 0.15 + # 转运>30分钟增加15%概率
(data['镇静程度'] == 0) * 0.25 + # 未镇静增加25%概率
(data['生命体征波动'] == 2) * 0.3 # 显著波动增加30%概率
> 0.5, 1, 0
)
})
# ----------------------
# 2. 逻辑回归筛选变量(计算P值,保留P<0.05)
# ----------------------
# 提取自变量和因变量
X = data[['GCS评分', '管道数量', '转运时长', '镇静程度', '生命体征波动']]
y = data['病情变化']
# 用statsmodels计算P值(更适合筛选变量)
X_sm = sm.add_constant(X) # 添加常数项(截距)
logit_model = sm.Logit(y, X_sm)
result = logit_model.fit()
# 输出变量系数和P值
var_pvalues = result.pvalues.drop('const') # 排除常数项
significant_vars = var_pvalues[var_pvalues < 0.05].index.tolist() # 显著变量列表
print("步骤1:逻辑回归筛选的显著影响因素(P<0.05):")
print(significant_vars)
print("\n各变量P值:")
print(var_pvalues.round(4))
# ----------------------
# 3. 可视化:变量重要性(系数绝对值)
# ----------------------
plt.figure(figsize=(8, 5))
coef = result.params.drop('const').abs() # 系数绝对值(代表影响强度)
sns.barplot(x=coef.index, y=coef.values)
plt.title('逻辑回归筛选的变量重要性(系数绝对值)', fontsize=12)
plt.ylabel('系数绝对值(越大影响越显著)')
plt.xticks(rotation=30)
plt.tight_layout()
plt.show()
from pgmpy.models import BayesianNetwork
from pgmpy.estimators import HillClimbSearch, BicScore
from pgmpy.inference import VariableElimination
import networkx as nx
# ----------------------
# 1. 准备数据(仅保留显著变量+新增不良事件因变量)
# ----------------------
# 新增“不良事件”因变量(模拟:由病情变化、管道数量、镇静程度共同影响)
data['不良事件'] = np.where(
(data['病情变化'] == 1) * 0.4 + # 病情变化增加40%概率
(data['管道数量'] > 2) * 0.3 + # 管道多增加30%概率
(data['镇静程度'] == 0) * 0.25 # 未镇静增加25%概率
> 0.5, 1, 0
)
# 仅保留步骤1筛选的显著变量+结局变量
bn_data = data[significant_vars + ['病情变化', '不良事件']]
# ----------------------
# 2. 构建贝叶斯网络(学习结构和参数)
# ----------------------
# 用Hill-Climb算法学习网络结构(基于数据关联)
hc = HillClimbSearch(bn_data)
model = hc.estimate(scoring_method=BicScore(bn_data)) # 用BIC评分优化结构
# 学习条件概率表(CPT)
from pgmpy.estimators import MaximumLikelihoodEstimator
model.fit(bn_data, estimator=MaximumLikelihoodEstimator)
# ----------------------
# 3. 可视化:网络关联路径
# ----------------------
plt.figure(figsize=(10, 6))
pos = nx.spring_layout(model) # 布局
nx.draw(model, pos, with_labels=True, node_size=3000, node_color='lightblue',
font_size=10, arrowsize=20)
# 标注关键路径(示例:手动标注临床关注的路径)
critical_edges = [('GCS评分', '生命体征波动'), ('生命体征波动', '病情变化'),
('管道数量', '不良事件'), ('病情变化', '不良事件')]
nx.draw_networkx_edges(model, pos, edgelist=critical_edges, edge_color='red', width=2)
plt.title('贝叶斯网络:转运相关因素关联路径(红色为关键路径)', fontsize=12)
plt.tight_layout()
plt.show()
# ----------------------
# 4. 推理:关键路径概率验证
# ----------------------
infer = VariableElimination(model)
# 例:计算“GCS评分低(<9)”时,“不良事件”的发生概率
prob_bad = infer.query(variables=['不良事件'],
evidence={'GCS评分': 6}) # GCS=6(低)
print("\n步骤2:贝叶斯网络推理结果:")
print(f"GCS评分低(6分)时,不良事件发生概率:{prob_bad.values[1]:.2f}")
from scipy.stats import spearmanr
# ----------------------
# 1. 提取贝叶斯网络中的核心变量对
# ----------------------
# 基于贝叶斯网络的关键路径,选择待分析的变量对
core_pairs = [
('GCS评分', '生命体征波动'), # GCS低是否与生命体征波动相关
('生命体征波动', '病情变化'), # 波动是否与病情变化相关
('管道数量', '不良事件') # 管道数量是否与不良事件相关
]
# ----------------------
# 2. 斯皮尔曼相关分析
# ----------------------
corr_results = []
for var1, var2 in core_pairs:
rho, p = spearmanr(data[var1], data[var2])
corr_results.append({
'变量对': f"{var1}-{var2}",
'相关系数ρ': round(rho, 2),
'P值': round(p, 4),
'关联强度': '强' if abs(rho) > 0.5 else '中' if abs(rho) > 0.3 else '弱'
})
# 输出结果
corr_df = pd.DataFrame(corr_results)
print("\n步骤3:斯皮尔曼相关分析结果(核心变量对):")
print(corr_df)
# ----------------------
# 3. 可视化:相关系数热力图
# ----------------------
plt.figure(figsize=(6, 4))
corr_matrix = data[significant_vars + ['病情变化', '不良事件']].corr(method='spearman')
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', vmin=-1, vmax=1)
plt.title('斯皮尔曼相关系数热力图(核心变量)', fontsize=12)
plt.tight_layout()
plt.show()怎么改,不能运行