【Python】当传统logistic遇上Bayes：逻辑回归的两种面孔

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf
from sklearn.linear_model import BayesianRidge
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import log_loss
from scipy.stats import norm
from statsmodels.stats.contingency_tables import Table2x2
import warnings

warnings.filterwarnings('ignore')

# 设置桌面工作目录和创建结果文件夹
desktop_path = os.path.join(os.path.expanduser("~"), "Desktop")
results_dir = os.path.join(desktop_path, "Result-Bayes")
if not os.path.exists(results_dir):
    os.makedirs(results_dir)

# 设置中文字体显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 显示负号

# 读取示例数据
print("正在读取数据...")
data = pd.read_excel(os.path.join(desktop_path, "示例数据.xlsx"), sheet_name="示例数据")

# 数据预处理
print("正在预处理数据...")
model_data = data.rename(columns={
    '结局': 'outcome',
    '指标1': 'indicator1',
    '指标2': 'indicator2',
    '指标3': 'indicator3',
    '指标4': 'indicator4',
    '指标5': 'indicator5',
    '指标6': 'indicator6',
    '指标7': 'indicator7',
    '指标8': 'indicator8',
    '时间': 'time',
    'BMI': 'bmi',
    '肥胖程度': 'obesity',
    '教育水平': 'education',
    '血型': 'blood_type'
}).copy()

# 转换分类变量
model_data['outcome'] = model_data['outcome'].astype('category')
model_data['indicator8'] = model_data['indicator8'].apply(lambda x: 'Yes'if x == '是'else'No').astype('category')
model_data['obesity'] = model_data['obesity'].astype('category')
model_data['education'] = model_data['education'].astype('category')
model_data['blood_type'] = model_data['blood_type'].astype('category')

# 0. 传统Logistic回归
print("正在进行传统Logistic回归...")
formula = 'outcome ~ indicator1 + indicator2 + indicator3 + indicator4 + indicator5 + indicator6 + indicator7 + indicator8 + time + bmi + obesity + education + blood_type'
logistic_model = smf.glm(formula=formula, data=model_data, family=sm.families.Binomial()).fit()
logistic_summary = logistic_model.summary()

# 提取OR和置信区间
logistic_or = np.exp(logistic_model.params)
logistic_ci = np.exp(logistic_model.conf_int())
logistic_ci.columns = ['Lower_CI', 'Upper_CI']
logistic_pvalues = logistic_model.pvalues

logistic_results = pd.DataFrame({
    'Variable': logistic_model.params.index,
    'OR': logistic_or,
    'Lower_CI': logistic_ci['Lower_CI'],
    'Upper_CI': logistic_ci['Upper_CI'],
    'p_value': logistic_pvalues,
    'Method': 'Traditional'
})

# 1. 使用Scikit-learn的贝叶斯岭回归
print("正在进行Scikit-learn贝叶斯岭回归...")
# 准备设计矩阵
X = pd.get_dummies(model_data.drop('outcome', axis=1), drop_first=True)
y = model_data['outcome'].cat.codes

# 标准化特征
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 拟合贝叶斯岭回归模型 - 使用 max_iter 替代 n_iter
bayesian_ridge = BayesianRidge(max_iter=1000, tol=1e-3, compute_score=True, alpha_1=1e-6, alpha_2=1e-6,
                              lambda_1=1e-6, lambda_2=1e-6)
bayesian_ridge.fit(X_scaled, y)

# 获取系数和不确定性
coef = bayesian_ridge.coef_
intercept = bayesian_ridge.intercept_

# 模拟后验分布
n_samples = 1000
coef_samples = np.random.normal(
    loc=coef,
    scale=np.sqrt(np.diag(bayesian_ridge.sigma_)),  # 修正：使用 np.diag 获取对角线元素
    size=(n_samples, len(coef))
)

# 修正：使用更简单的方法计算截距的标准差
intercept_std = np.sqrt(1.0 / bayesian_ridge.alpha_)
intercept_samples = np.random.normal(
    loc=intercept,
    scale=intercept_std,
    size=n_samples
)

# 计算OR和置信区间
or_samples = np.exp(np.column_stack([intercept_samples.reshape(-1, 1), coef_samples]))
or_mean = np.mean(or_samples, axis=0)
or_lower = np.percentile(or_samples, 2.5, axis=0)
or_upper = np.percentile(or_samples, 97.5, axis=0)

# 创建结果DataFrame
sklearn_results = pd.DataFrame({
    'Variable': ['Intercept'] + list(X.columns),
    'OR': or_mean,
    'Lower_CI': or_lower,
    'Upper_CI': or_upper,
    'p_value': np.nan,
    'Method': 'Sklearn_Bayesian'
})

# 合并结果
all_results = pd.concat([
    logistic_results[logistic_results['Variable'] != 'Intercept'],
    sklearn_results[sklearn_results['Variable'] != 'Intercept']
], ignore_index=True)

# 变量标签
variable_labels = {
    "indicator1": "指标1",
    "indicator2": "指标2",
    "indicator3": "指标3",
    "indicator4": "指标4",
    "indicator5": "指标5",
    "indicator6": "指标6",
    "indicator7": "指标7",
    "indicator8_Yes": "指标8 (是 vs 否)",
    "time": "时间",
    "bmi": "BMI",
    "obesity_超重": "肥胖程度 (超重 vs 正常)",
    "obesity_肥胖": "肥胖程度 (肥胖 vs 正常)",
    "obesity_偏瘦": "肥胖程度 (偏瘦 vs 正常)",
    "education_中学": "教育水平 (中学 vs 小学)",
    "education_大学": "教育水平 (大学 vs 小学)",
    "education_硕士及以上": "教育水平 (硕士及以上 vs 小学)",
    "blood_type_A": "血型 (A vs O)",
    "blood_type_AB": "血型 (AB vs O)",
    "blood_type_B": "血型 (B vs O)"
}

all_results['Variable_Label'] = all_results['Variable'].map(variable_labels)

# 保存结果到Excel
with pd.ExcelWriter(os.path.join(results_dir, "Bayes_Logistic_Comparison_Results.xlsx")) as writer:
    all_results.to_excel(writer, sheet_name="所有方法比较", index=False)
    logistic_results[logistic_results['Variable'] != 'Intercept'].to_excel(
        writer, sheet_name="传统Logistic", index=False)
    sklearn_results[sklearn_results['Variable'] != 'Intercept'].to_excel(
        writer, sheet_name="贝叶斯回归", index=False)

# 1. 两种方法比较图（对数尺度）
plt.figure(figsize=(14, 10))
colors = ['blue', 'red']

for i, method in enumerate(all_results['Method'].unique()):
    method_data = all_results[all_results['Method'] == method]
    plt.scatter(method_data['OR'], range(len(method_data)),
                color=colors[i], label=method, s=50, alpha=0.7)
    for j, row in method_data.iterrows():
        plt.plot([row['Lower_CI'], row['Upper_CI']], [j, j],
                 color=colors[i], linewidth=1.5)

plt.axvline(x=1, linestyle='--', color='gray')
plt.xscale('log')
plt.yticks(range(len(all_results)), all_results['Variable_Label'])
plt.title('贝叶斯 vs Logistic回归OR值比较 (对数尺度)')
plt.xlabel('OR值 (对数尺度)')
plt.ylabel('变量')
plt.legend()
plt.tight_layout()
plt.savefig(os.path.join(results_dir, "两种方法比较图_对数尺度.jpg"), dpi=300)
plt.close()

# 2. 两种方法比较图（原始尺度）
plt.figure(figsize=(14, 10))
for i, method in enumerate(all_results['Method'].unique()):
    method_data = all_results[all_results['Method'] == method]
    plt.scatter(method_data['OR'], range(len(method_data)),
                color=colors[i], label=method, s=50, alpha=0.7)
    for j, row in method_data.iterrows():
        plt.plot([row['Lower_CI'], row['Upper_CI']], [j, j],
                 color=colors[i], linewidth=1.5)

plt.axvline(x=1, linestyle='--', color='gray')
plt.yticks(range(len(all_results)), all_results['Variable_Label'])
plt.title('贝叶斯 vs Logistic回归OR值比较 (原始尺度)')
plt.xlabel('OR值')
plt.ylabel('变量')
plt.legend()
plt.tight_layout()
plt.savefig(os.path.join(results_dir, "两种方法比较图_原始尺度.jpg"), dpi=300)
plt.close()

# 3. 分别绘制两种方法的森林图（对数尺度）
def create_forest_plot_log(data, title, color="steelblue"):
    plt.figure(figsize=(12, 10))
    y_pos = range(len(data))
    plt.errorbar(data['OR'], y_pos,
                 xerr=[data['OR'] - data['Lower_CI'], data['Upper_CI'] - data['OR']],
                 fmt='o', color=color, ecolor=color, elinewidth=1, capsize=5)
    plt.axvline(x=1, linestyle='--', color='red', linewidth=1)
    plt.xscale('log')
    plt.yticks(y_pos, data['Variable_Label'])
    plt.title(f'{title} (对数尺度)')
    plt.xlabel('OR值 (对数尺度)')
    plt.ylabel('变量')
    plt.tight_layout()
    return plt

# 4. 分别绘制两种方法的森林图（原始尺度）
def create_forest_plot_original(data, title, color="steelblue"):
    plt.figure(figsize=(12, 10))
    y_pos = range(len(data))
    plt.errorbar(data['OR'], y_pos,
                 xerr=[data['OR'] - data['Lower_CI'], data['Upper_CI'] - data['OR']],
                 fmt='o', color=color, ecolor=color, elinewidth=1, capsize=5)
    plt.axvline(x=1, linestyle='--', color='red', linewidth=1)
    plt.yticks(y_pos, data['Variable_Label'])
    plt.title(f'{title} (原始尺度)')
    plt.xlabel('OR值')
    plt.ylabel('变量')
    plt.tight_layout()
    return plt

# 创建各个方法的森林图（对数尺度）
logistic_data = logistic_results[logistic_results['Variable'] != 'Intercept'].copy()
logistic_data['Variable_Label'] = logistic_data['Variable'].map(variable_labels)
logistic_plot_log = create_forest_plot_log(logistic_data, "传统Logistic回归结果", "blue")
logistic_plot_log.savefig(os.path.join(results_dir, "传统Logistic回归森林图_对数尺度.jpg"), dpi=300)
logistic_plot_log.close()

sklearn_data = sklearn_results[sklearn_results['Variable'] != 'Intercept'].copy()
sklearn_data['Variable_Label'] = sklearn_data['Variable'].map(variable_labels)
sklearn_plot_log = create_forest_plot_log(sklearn_data, "贝叶斯回归结果", "red")
sklearn_plot_log.savefig(os.path.join(results_dir, "贝叶斯回归森林图_对数尺度.jpg"), dpi=300)
sklearn_plot_log.close()

# 创建各个方法的森林图（原始尺度）
logistic_plot_original = create_forest_plot_original(logistic_data, "传统Logistic回归结果", "blue")
logistic_plot_original.savefig(os.path.join(results_dir, "传统Logistic回归森林图_原始尺度.jpg"), dpi=300)
logistic_plot_original.close()

sklearn_plot_original = create_forest_plot_original(sklearn_data, "贝叶斯回归结果", "red")
sklearn_plot_original.savefig(os.path.join(results_dir, "贝叶斯回归森林图_原始尺度.jpg"), dpi=300)
sklearn_plot_original.close()

# 5. 率比（Risk Ratio）计算函数
def calculate_rr(var_name, data):
    # 连续变量不计算RR
    continuous_vars = ["indicator1", "indicator2", "indicator3", "indicator4",
                       "indicator5", "indicator6", "indicator7", "time", "bmi"]
    if var_name in continuous_vars:
        return pd.DataFrame({'Variable': [var_name], 'RR': [np.nan],
                             'RR_Lower_CI': [np.nan], 'RR_Upper_CI': [np.nan]})

    # 分类变量计算RR
    if var_name == "indicator8_Yes":
        data_temp = data.copy()
        data_temp['temp_var'] = (data_temp['indicator8'] == 'Yes').astype(int)
        tab = pd.crosstab(data_temp['temp_var'], data_temp['outcome'])

    elif'obesity'in var_name:
        level = var_name.split('_')[1]
        data_temp = data.copy()
        data_temp['temp_var'] = (data_temp['obesity'] == level).astype(int)
        tab = pd.crosstab(data_temp['temp_var'], data_temp['outcome'])

    elif'education'in var_name:
        level = var_name.split('_')[1]
        data_temp = data.copy()
        data_temp['temp_var'] = (data_temp['education'] == level).astype(int)
        tab = pd.crosstab(data_temp['temp_var'], data_temp['outcome'])

    elif'blood_type'in var_name:
        level = var_name.split('_')[1]
        data_temp = data.copy()
        data_temp['temp_var'] = (data_temp['blood_type'] == level).astype(int)
        tab = pd.crosstab(data_temp['temp_var'], data_temp['outcome'])

    else:
        return pd.DataFrame({'Variable': [var_name], 'RR': [np.nan],
                             'RR_Lower_CI': [np.nan], 'RR_Upper_CI': [np.nan]})

    if tab.shape == (2, 2) and tab.sum().sum() > 0:
        try:
            table = Table2x2(tab.values)
            rr = table.riskratio
            rr_ci = table.riskratio_confint()
            return pd.DataFrame({'Variable': [var_name], 'RR': [rr],
                                 'RR_Lower_CI': [rr_ci[0]], 'RR_Upper_CI': [rr_ci[1]]})
        except:
            return pd.DataFrame({'Variable': [var_name], 'RR': [np.nan],
                                 'RR_Lower_CI': [np.nan], 'RR_Upper_CI': [np.nan]})
    else:
        return pd.DataFrame({'Variable': [var_name], 'RR': [np.nan],
                             'RR_Lower_CI': [np.nan], 'RR_Upper_CI': [np.nan]})

# 计算所有分类变量的RR
print("正在计算率比 (RR)...")
categorical_vars_for_rr = [
    "indicator8_Yes", "obesity_超重", "obesity_肥胖", "obesity_偏瘦",
    "education_中学", "education_大学", "education_硕士及以上",
    "blood_type_A", "blood_type_AB", "blood_type_B"
]

rr_results = [calculate_rr(var, model_data) for var in categorical_vars_for_rr]
rr_results_df = pd.concat(rr_results, ignore_index=True)

# 保存RR结果
rr_results_df.to_csv(os.path.join(results_dir, "Rate_Ratios.csv"), index=False, encoding='utf-8-sig')

# 合并OR和RR结果
final_results = all_results.merge(rr_results_df, on='Variable', how='left')

# 保存最终结果
final_results.to_csv(os.path.join(results_dir, "Final_Results.csv"), index=False, encoding='utf-8-sig')

# 6. 率比（RR）森林图
print("正在生成率比森林图...")
rr_plot_data = rr_results_df.dropna().copy()
rr_plot_data['Variable_Label'] = rr_plot_data['Variable'].map(variable_labels)

if not rr_plot_data.empty:
    plt.figure(figsize=(10, 8))
    y_pos = range(len(rr_plot_data))
    plt.errorbar(rr_plot_data['RR'], y_pos,
                 xerr=[rr_plot_data['RR'] - rr_plot_data['RR_Lower_CI'],
                       rr_plot_data['RR_Upper_CI'] - rr_plot_data['RR']],
                 fmt='o', color='steelblue', ecolor='steelblue', elinewidth=1, capsize=5)
    plt.axvline(x=1, linestyle='--', color='red')
    plt.yticks(y_pos, rr_plot_data['Variable_Label'])
    plt.title('率比（Risk Ratio）森林图')
    plt.xlabel('率比 (RR)')
    plt.ylabel('变量')
    plt.tight_layout()
    plt.savefig(os.path.join(results_dir, "Rate_Ratio_ForestPlot.jpg"), dpi=300)
    plt.close()
else:
    print("没有可用的率比数据进行可视化")

# 7. 事件率分析图 - 按BMI分组
model_data['bmi_group'] = pd.qcut(model_data['bmi'], q=4, duplicates='drop')

event_rate_data = model_data.groupby('bmi_group')['outcome'].apply(
    lambda x: (x.cat.codes == 1).mean() if x.cat.codes.nunique() > 1 else 0
).reset_index()
event_rate_data.columns = ['bmi_group', 'Event_Rate']
event_rate_data['n'] = model_data.groupby('bmi_group').size().values

plt.figure(figsize=(8, 6))
bars = plt.bar(range(len(event_rate_data)), event_rate_data['Event_Rate'],
               alpha=0.7, color='steelblue')
plt.xticks(range(len(event_rate_data)), event_rate_data['bmi_group'].astype(str), rotation=45)
for i, (event_rate, n) in enumerate(zip(event_rate_data['Event_Rate'], event_rate_data['n'])):
    plt.text(i, event_rate + 0.01, f'{event_rate:.2%}\n(n={n})', ha='center', va='bottom', fontsize=10)
plt.title('按BMI分组的结局事件率')
plt.xlabel('BMI分组')
plt.ylabel('事件率')
plt.tight_layout()
plt.savefig(os.path.join(results_dir, "事件率分析图.jpg"), dpi=300)
plt.close()

# 8. 分类变量的事件率分析
categorical_vars = ['indicator8', 'obesity', 'education', 'blood_type']

for var in categorical_vars:
    event_rate_var = model_data.groupby([var, 'outcome']).size().unstack(fill_value=0)
    event_rate_var_pct = event_rate_var.div(event_rate_var.sum(axis=1), axis=0) * 100

    plt.figure(figsize=(10, 6))
    event_rate_var_pct.plot(kind='bar', alpha=0.7)
    plt.title(f'按 {var} 分组的事件率')
    plt.xlabel(var)
    plt.ylabel('百分比 (%)')
    plt.legend(title='结局')
    plt.xticks(rotation=45)
    plt.tight_layout()
    plt.savefig(os.path.join(results_dir, f"事件率_{var}.jpg"), dpi=300)
    plt.close()

# 计算模型性能指标
# 传统Logistic回归的预测概率
logistic_probs = logistic_model.predict(model_data)
logistic_loss = log_loss(y, logistic_probs)

# Scikit-learn贝叶斯回归的预测概率
sklearn_probs = 1 / (1 + np.exp(-(intercept + X_scaled.dot(coef))))
sklearn_loss = log_loss(y, sklearn_probs)

# 创建性能比较表
performance_comparison = pd.DataFrame({
    'Method': ['Traditional_Logistic', 'Sklearn_Bayesian'],
    'Log_Loss': [logistic_loss, sklearn_loss]
})

performance_comparison.to_csv(os.path.join(results_dir, "Performance_Comparison.csv"), index=False, encoding='utf-8-sig')

# 输出主要结果到控制台
print("\n=== 传统Logistic回归结果 ===")
print(logistic_results.head())

print("\n=== Scikit-learn贝叶斯回归结果 ===")
print(sklearn_results.head())

print("\n=== 率比（RR）结果 ===")
print(rr_results_df.head())

print("\n=== 模型性能比较 ===")
print(performance_comparison)

print("\n分析完成！所有结果已保存到 Result-Bayes 文件夹中。")
print("包含内容：")
print("1. 贝叶斯和Logistic回归比较 (Excel和CSV)")
print("2. 对数尺度和原始尺度比较图")
print("3. 两种方法的单独森林图")
print("4. 率比（RR）森林图")
print("5. 事件率分析图")
print("6. 分类变量的事件率分析图")
print("7. 模型性能比较")