临床数据分析与回归建模（PICU临床数据应用）

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from collections import Counter
import statsmodels.api as sm
from sklearn.metrics import roc_auc_score, roc_curve, auc
import matplotlib.pyplot as plt

path = "icu_first24hours.csv"
picu_data = pd.read_csv(path)
picu_data
# 筛选感兴趣的特征集（年龄、结局变量、关键实验室指标）
core_features = [
    'age_month', 'HOSPITAL_EXPIRE_FLAG',
    'lab_5097_max', 'lab_5099_min', 'lab_5237_min',
    'lab_5227_min', 'lab_5225_range', 'lab_5235_max', 'lab_5257_min'
]
new_picu_data = picu_data[core_features].copy()

# 查看各列的非空值数量 (count)
new_picu_data.count()
# 计算均值 (mean)
new_picu_data.mean()
# 计算中位数 (median)
new_picu_data.median()
# 计算方差 (var)
new_picu_data.var()
#new_picu_data.to_excel("data.xlsx", index=False)

# === 直方图绘制代码 ===
# 定义要绘制直方图的列名列表，包含8个临床指标
colname =['age_month','lab_5097_max','lab_5099_min','lab_5237_min','lab_5227_min','lab_5225_range','lab_5235_max','lab_5257_min']

# 创建子图，共4行2列，使用seaborn的histplot函数绘制直方图
fig, axs = plt.subplots(int(len(colname)/2), 2, constrained_layout=True, figsize=(8, 6), dpi=150)

# 使用seaborn的histplot函数绘制直方图
for i in range(len(colname)):
    sns.histplot(x=colname[i], data=new_picu_data, alpha=0.4, kde=True, ax=axs[i//2, i%2])
plt.suptitle("Bar Plot")
# === 箱线图绘制代码） ===
# 创建一个新的图形对象和子图数组用于绘制箱线图
fig, axs = plt.subplots(
    int(len(colname) / 2), 2, 
    constrained_layout=True, 
    figsize=(8, 6), 
    dpi=150
)
# 使用seaborn的boxplot函数绘制分组箱线图
for i in range(len(colname)): 
    sns.boxplot(
        data=new_picu_data, 
        x="HOSPITAL_EXPIRE_FLAG", 
        y=colna

#数据清洗，填补缺失值
picu_data_with_median = new_picu_data.fillna(new_picu_data.median())

# 自变量 X 和 因变量 y
X = picu_data_with_median.drop(columns="HOSPITAL_EXPIRE_FLAG")
y = picu_data_with_median["HOSPITAL_EXPIRE_FLAG"]

# 加截距项
X = sm.add_constant(X)

# Logistic 回归
logit_model = sm.Logit(y, X)
result = logit_model.fit()

# 提取回归系数和置信区间
params = result.params
conf = result.conf_int()
conf.columns = ['2.5%', '97.5%']

# 合并并指数化（得到 OR）
or_table = pd.concat([params, conf], axis=1)
or_table = np.exp(or_table)

# 重命名列
or_table.columns = ['OR', '2.5%', '97.5%']
# 去掉截距（
or_table = or_table.drop(index='const', errors='ignore')
# 保留两位小数
or_table = or_table.round(2)
# 按 OR 排序
or_table = or_table.sort_values(by='OR', ascending=False)
or_table

# 预测概率
y_pred_prob = result.predict(X)

# AUC
auc = roc_auc_score(y, y_pred_prob)
print("AUC =", auc)

# ROC 曲线
fpr, tpr, _ = roc_curve(y, y_pred_prob)

plt.figure(figsize=(5, 5), dpi=150)
plt.plot(fpr, tpr, label=f"AUC = {auc:.3f}")
plt.plot([0, 1], [0, 1], linestyle="--", color="gray")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title("ROC Curve")
plt.legend()
plt.show()

selected_features = [
    'lab_5237_min',
    'lab_5227_min',
    'lab_5225_range',
    'lab_5235_max',
    'lab_5257_min'
]

# 构建新数据集
X = picu_data_with_median[selected_features]
y = picu_data_with_median['HOSPITAL_EXPIRE_FLAG']
X_const = sm.add_constant(X)
# 拟合模型
logit_model_reduced = sm.Logit(y, X_const)
result_reduced = logit_model_reduced.fit()

# 输出结果
y_pred_prob = result_reduced.predict(X_const)

params = result_reduced.params
conf = result_reduced.conf_int()
conf.columns = ['2.5%', '97.5%']

# 转换为 OR
or_df = pd.DataFrame({
    'OR': np.exp(params),
    '2.5%': np.exp(conf['2.5%']),
    '97.5%': np.exp(conf['97.5%'])
})

# 去掉截距
or_df = or_df.drop(index='const', errors='ignore')

# 保留两位小数
or_df = or_df.round(2)

or_df

# y_true: 真实标签
# y_pred_prob: 预测概率（logit_model.predict()）

fpr, tpr, thresholds = roc_curve(y, y_pred_prob)
roc_auc = auc(fpr, tpr)

plt.figure(figsize=(5, 5), dpi=150)

# ROC 曲线 —— 棕色
plt.plot(
    fpr, tpr,
    color='saddlebrown',
    lw=2,
    label=f'ROC curve (AUC = {roc_auc:.3f})'
)

# 对角线 —— 灰色
plt.plot(
    [0, 1], [0, 1],
    color='gray',
    lw=1.5,
    linestyle='--',
    label='Random'
)

plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve')
plt.legend(loc='lower right')
plt.grid(alpha=0.3)

plt.show()