import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import os
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.formula.api import ols
import statsmodels.api as sm
from scipy.stats import chi2_contingency, fisher_exact
import warnings
warnings.filterwarnings('ignore')
# 设置matplotlib参数,使图表更美观,支持中文显示
plt.rcParams['font.sans-serif'] = ['SimHei'] # 设置中文字体
plt.rcParams['axes.unicode_minus'] = False # 解决保存图像时负号'-'显示为方块的问题
plt.rcParams['figure.dpi'] = 300 # 设置图像DPI为300
plt.rcParams['savefig.dpi'] = 300 # 设置保存图像的DPI为300
# 创建保存结果的文件夹
if not os.path.exists('问题1结果'):
os.makedirs('问题1结果')
# 读取数据
file_path = '附件1:两个医院临床受试者及抑郁症的基本数据.xlsx'
# 读取两个子表
df1 = pd.read_excel(file_path, sheet_name='一院临床受试者及抑郁症的基本数据', header=1)
df2 = pd.read_excel(file_path, sheet_name='二院临床受试者及抑郁症的基本数据', header=1)
# 重命名列以匹配实际数据
columns = ['序号', '组别', '年龄_岁', '未婚', '已婚', '离异', '丧偶',
'无', '使用过抗抑郁药', '其它', '轻度', '中度', '重度']
df1.columns = columns
df2.columns = columns
# 添加医院标识
df1['医院'] = '医院1'
df2['医院'] = '医院2'
# 合并数据集
combined_df = pd.concat([df1, df2], ignore_index=True)
# 数据预处理
# 将婚姻状况转换为单列分类变量
combined_df['婚姻状况'] = combined_df[['未婚', '已婚', '离异', '丧偶']].idxmax(axis=1)
# 将既往抗抑郁药使用情况转换为单列分类变量
combined_df['既往用药'] = combined_df[['无', '使用过抗抑郁药', '其它']].idxmax(axis=1)
# 处理抑郁程度列
for col in ['轻度', '中度', '重度']:
combined_df[col] = combined_df[col].astype(str)
combined_df[col] = combined_df[col].apply(lambda x: 1 if x == '1' else 0)
combined_df['抑郁程度'] = combined_df[['轻度', '中度', '重度']].idxmax(axis=1)
# 转换组别为字符串类型
combined_df['组别'] = combined_df['组别'].astype(str)
# 添加药物名称对应关系
drug_map = {'1': '药物A', '2': '药物B', '3': '药物C'}
combined_df['药物名称'] = combined_df['组别'].map(drug_map)
# 年龄分析
# 描述性统计
age_stats = combined_df.groupby('组别')['年龄_岁'].describe()
age_stats.to_excel('问题1结果/问题1_年龄描述性统计.xlsx')
# 绘制年龄分布箱线图
plt.figure(figsize=(10, 6))
sns.boxplot(x='药物名称', y='年龄_岁', data=combined_df, palette='Set3')
plt.title('各组受试者年龄分布箱线图')
plt.xlabel('药物名称')
plt.ylabel('年龄(岁)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig('问题1结果/问题1_年龄分布箱线图.png')
plt.show()
# 绘制年龄分布核密度图
plt.figure(figsize=(10, 6))
for group in sorted(combined_df['药物名称'].unique()):
group_data = combined_df[combined_df['药物名称'] == group]['年龄_岁']
sns.kdeplot(group_data, label=group, shade=True, alpha=0.5)
plt.title('各组受试者年龄分布核密度图')
plt.xlabel('年龄(岁)')
plt.ylabel('密度')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig('问题1结果/问题1_年龄分布核密度图.png')
plt.show()
# ANOVA分析检验年龄差异
formula = '年龄_岁 ~ C(组别)'
model = ols(formula, data=combined_df).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print("年龄ANOVA分析结果:")
print(anova_table)
# 计算效应量
eta_squared = anova_table['sum_sq'][0] / (anova_table['sum_sq'][0] + anova_table['sum_sq'][1])
print(f"年龄方差分析效应量 (Eta Squared): {eta_squared:.4f}")
# 事后检验
tukey_results = pairwise_tukeyhsd(combined_df['年龄_岁'], combined_df['组别'])
print("年龄Tukey HSD事后检验结果:")
print(tukey_results)
# 将Tukey结果转换为DataFrame并保存
tukey_df = pd.DataFrame(data=tukey_results._results_table.data[1:],
columns=tukey_results._results_table.data[0])
tukey_df.to_excel('问题1结果/问题1_年龄Tukey事后检验.xlsx', index=False)
# 婚姻状况分析
marriage_counts = pd.crosstab(combined_df['药物名称'], combined_df['婚姻状况'])
marriage_percentages = pd.crosstab(combined_df['药物名称'], combined_df['婚姻状况'], normalize='index') * 100
marriage_counts.to_excel('问题1结果/问题1_婚姻状况分布频数.xlsx')
marriage_percentages.to_excel('问题1结果/问题1_婚姻状况分布百分比.xlsx')
# 绘制婚姻状况分布堆积柱状图
plt.figure(figsize=(10, 6))
marriage_counts.plot(kind='bar', stacked=True, colormap='Paired')
plt.title('各药物组受试者婚姻状况分布')
plt.xlabel('药物名称')
plt.ylabel('人数')
plt.legend(title='婚姻状况')
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig('问题1结果/问题1_婚姻状况分布堆积柱状图.png')
plt.show()
# 绘制婚姻状况分布热图
plt.figure(figsize=(10, 6))
sns.heatmap(marriage_percentages, annot=True, fmt='.1f', cmap='YlGnBu')
plt.title('各药物组受试者婚姻状况分布百分比热图')
plt.xlabel('婚姻状况')
plt.ylabel('药物名称')
plt.tight_layout()
plt.savefig('问题1结果/问题1_婚姻状况分布热图.png')
plt.show()
# 绘制婚姻状况分布堆积柱状图(带比例标注)
fig, ax = plt.subplots(figsize=(10, 6))
marriage_counts.plot(kind='bar', stacked=True, colormap='Paired', ax=ax)
for idx, bar in enumerate(ax.patches):
if idx % 4 == 0:
height = bar.get_height()
width = bar.get_width()
x = bar.get_x() + width / 2
y = height + 0.1
ax.text(x, y, f'{height/sum(marriage_counts.iloc[int(idx/4)])/4:.1%}', ha='center', fontsize=9)
plt.title('各药物组受试者婚姻状况分布(带比例标注)')
plt.xlabel('药物名称')
plt.ylabel('人数')
plt.legend(title='婚姻状况')
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig('问题1结果/问题1_婚姻状况分布堆积柱状图优化.png')
plt.show()
# 卡方检验婚姻状况差异
chi2, p, dof, expected = chi2_contingency(marriage_counts)
print("婚姻状况卡方检验结果:")
print(f"卡方值: {chi2:.4f}, p值: {p:.4f}, 自由度: {dof}")
# 计算Cramer's V效应量
n = marriage_counts.sum().sum()
cramer_v = np.sqrt(chi2 / (n * (min(marriage_counts.shape) - 1)))
print(f"Cramer's V效应量: {cramer_v:.4f}")
# 保存婚姻状况卡方检验结果
chi2_results = pd.DataFrame({
'统计量': ['卡方值', 'p值', '自由度', 'Cramer\'s V'],
'数值': [chi2, p, dof, cramer_v]
})
chi2_results.to_excel('问题1结果/问题1_婚姻状况卡方检验结果.xlsx', index=False)
# 既往抗抑郁药使用情况分析
med_history_counts = pd.crosstab(combined_df['药物名称'], combined_df['既往用药'])
med_history_percentages = pd.crosstab(combined_df['药物名称'], combined_df['既往用药'], normalize='index') * 100
med_history_counts.to_excel('问题1结果/问题1_既往用药分布频数.xlsx')
med_history_percentages.to_excel('问题1结果/问题1_既往用药分布百分比.xlsx')
# 绘制既往用药分布堆积柱状图
plt.figure(figsize=(10, 6))
med_history_counts.plot(kind='bar', stacked=True, colormap='Set2')
plt.title('各药物组受试者既往抗抑郁药使用情况分布')
plt.xlabel('药物名称')
plt.ylabel('人数')
plt.legend(title='既往用药情况')
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig('问题1结果/问题1_既往用药分布堆积柱状图.png')
plt.show()
# 绘制既往用药分布饼图组合
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
groups = sorted(combined_df['药物名称'].unique())
for i, group in enumerate(groups):
group_data = med_history_counts.loc[group]
axes[i].pie(group_data, labels=group_data.index, autopct='%1.1f%%',
startangle=90, colors=sns.color_palette('Set2', n_colors=len(group_data)))
axes[i].set_title(f'{group} 既往用药情况')
plt.tight_layout()
plt.savefig('问题1结果/问题1_既往用药分布饼图组合.png')
plt.show()
# 绘制既往用药分布堆积柱状图(带比例标注)
fig, ax = plt.subplots(figsize=(10, 6))
med_history_counts.plot(kind='bar', stacked=True, colormap='Set2', ax=ax)
for idx, bar in enumerate(ax.patches):
if idx % 3 == 0:
height = bar.get_height()
width = bar.get_width()
x = bar.get_x() + width / 2
y = height + 0.1
ax.text(x, y, f'{height/sum(med_history_counts.iloc[int(idx/3)])/3:.1%}', ha='center', fontsize=9)
plt.title('各药物组受试者既往抗抑郁药使用情况分布(带比例标注)')
plt.xlabel('药物名称')
plt.ylabel('人数')
plt.legend(title='既往用药情况')
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig('问题1结果/问题1_既往用药分布堆积柱状图优化.png')
plt.show()
# 卡方检验既往用药差异
chi2, p, dof, expected = chi2_contingency(med_history_counts)
print("既往抗抑郁药使用情况卡方检验结果:")
print(f"卡方值: {chi2:.4f}, p值: {p:.4f}, 自由度: {dof}")
# 计算Cramer's V效应量
n = med_history_counts.sum().sum()
cramer_v = np.sqrt(chi2 / (n * (min(med_history_counts.shape) - 1)))
print(f"Cramer's V效应量: {cramer_v:.4f}")
# 保存既往用药卡方检验结果
chi2_results = pd.DataFrame({
'统计量': ['卡方值', 'p值', '自由度', 'Cramer\'s V'],
'数值': [chi2, p, dof, cramer_v]
})
chi2_results.to_excel('问题1结果/问题1_既往用药卡方检验结果.xlsx', index=False)
# 初始抑郁程度分析
depression_counts = pd.crosstab(combined_df['药物名称'], combined_df['抑郁程度'])
depression_percentages = pd.crosstab(combined_df['药物名称'], combined_df['抑郁程度'], normalize='index') * 100
depression_counts.to_excel('问题1结果/问题1_抑郁程度分布频数.xlsx')
depression_percentages.to_excel('问题1结果/问题1_抑郁程度分布百分比.xlsx')
# 绘制抑郁程度分布堆积柱状图
plt.figure(figsize=(10, 6))
depression_counts.plot(kind='bar', stacked=True, colormap='viridis')
plt.title('各药物组受试者初始抑郁程度分布')
plt.xlabel('药物名称')
plt.ylabel('人数')
plt.legend(title='抑郁程度')
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig('问题1结果/问题1_抑郁程度分布堆积柱状图.png')
plt.show()
# 绘制抑郁程度分布马赛克图
plt.figure(figsize=(10, 6))
from statsmodels.graphics.mosaicplot import mosaic
depression_counts_df = pd.DataFrame(depression_counts)
depression_index = depression_counts_df.index.tolist()
depression_columns = depression_counts_df.columns.tolist()
mosaic_data = {}
for idx in depression_index:
for col in depression_columns:
mosaic_data[(idx, col)] = depression_counts.loc[idx, col]
mosaic(mosaic_data, gap=0.02, title='各药物组受试者初始抑郁程度分布马赛克图',
properties=lambda key: {'color': plt.cm.viridis(depression_columns.index(key[1]) / len(depression_columns))})
plt.tight_layout()
plt.savefig('问题1结果/问题1_抑郁程度分布马赛克图.png')
plt.show()
# 绘制抑郁程度分布堆积柱状图(带比例标注)
fig, ax = plt.subplots(figsize=(10, 6))
depression_counts.plot(kind='bar', stacked=True, colormap='viridis', ax=ax)
for idx, bar in enumerate(ax.patches):
if idx % 3 == 0:
height = bar.get_height()
width = bar.get_width()
x = bar.get_x() + width / 2
y = height + 0.1
ax.text(x, y, f'{height/sum(depression_counts.iloc[int(idx/3)])/3:.1%}', ha='center', fontsize=9)
plt.title('各药物组受试者初始抑郁程度分布(带比例标注)')
plt.xlabel('药物名称')
plt.ylabel('人数')
plt.legend(title='抑郁程度')
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig('问题1结果/问题1_抑郁程度分布堆积柱状图优化.png')
plt.show()
# 卡方检验抑郁程度差异
chi2, p, dof, expected = chi2_contingency(depression_counts)
print("初始抑郁程度卡方检验结果:")
print(f"卡方值: {chi2:.4f}, p值: {p:.4f}, 自由度: {dof}")
# 计算Cramer's V效应量
n = depression_counts.sum().sum()
cramer_v = np.sqrt(chi2 / (n * (min(depression_counts.shape) - 1)))
print(f"Cramer's V效应量: {cramer_v:.4f}")
# 保存抑郁程度卡方检验结果
chi2_results = pd.DataFrame({
'统计量': ['卡方值', 'p值', '自由度', 'Cramer\'s V'],
'数值': [chi2, p, dof, cramer_v]
})
chi2_results.to_excel('问题1结果/问题1_抑郁程度卡方检验结果.xlsx', index=False)
# 综合基线特征雷达图
categories = ['平均年龄', '未婚比例(%)', '无既往用药比例(%)', '重度抑郁比例(%)']
drug_labels = sorted(combined_df['药物名称'].unique())
# 计算各组的值
baseline_values = {}
for drug in drug_labels:
group_data = combined_df[combined_df['药物名称'] == drug]
age_mean = group_data['年龄_岁'].mean()
unmarried_pct = marriage_percentages.loc[drug, '未婚'] if drug in marriage_percentages.index else 0
no_med_history_pct = med_history_percentages.loc[drug, '无'] if drug in med_history_percentages.index else 0
severe_depression_pct = depression_percentages.loc[drug, '重度'] if drug in depression_percentages.index else 0
baseline_values[drug] = [age_mean, unmarried_pct, no_med_history_pct, severe_depression_pct]
# 绘制雷达图
angles = np.linspace(0, 2 * np.pi, len(categories), endpoint=False).tolist()
angles += angles[:1] # 闭合图形
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(polar=True))
for drug, values in baseline_values.items():
values += values[:1] # 闭合图形
ax.plot(angles, values, linewidth=2, label=f'{drug}')
ax.fill(angles, values, alpha=0.1) # 修正变量名
ax.set_thetagrids(np.degrees(angles[:-1]), categories)
ax.set_ylim(0, 100)
ax.set_title('各药物组基线特征雷达图比较', size=15)
ax.legend(loc='upper right', bbox_to_anchor=(0.1, 0.1))
plt.tight_layout()
plt.savefig('问题1结果/问题1_基线特征雷达图.png')
plt.show()
# 基线特征总结表格
age_summary = combined_df.groupby('组别')['年龄_岁'].agg(['mean', 'std']).round(2)
age_summary.columns = ['平均年龄', '年龄标准差']
marriage_summary = marriage_percentages.round(1)
med_history_summary = med_history_percentages.round(1)
depression_summary = depression_percentages.round(1)
summary_columns = pd.MultiIndex.from_tuples([
('年龄', '平均值'), ('年龄', '标准差'),
('婚姻状况', '未婚%'), ('婚姻状况', '已婚%'), ('婚姻状况', '离异%'), ('婚姻状况', '丧偶%'),
('既往用药', '无%'), ('既往用药', '使用过%'), ('既往用药', '其它%'),
('抑郁程度', '轻度%'), ('抑郁程度', '中度%'), ('抑郁程度', '重度%')
])
baseline_summary = pd.DataFrame(index=sorted(combined_df['组别'].unique()), columns=summary_columns)
for group in sorted(combined_df['组别'].unique()):
baseline_summary.loc[group, ('年龄', '平均值')] = age_summary.loc[group, '平均年龄']
baseline_summary.loc[group, ('年龄', '标准差')] = age_summary.loc[group, '年龄标准差']
for status in marriage_percentages.columns:
baseline_summary.loc[group, ('婚姻状况', f'{status}%')] = marriage_percentages.loc[group, status] if group in marriage_percentages.index else 0
for history in med_history_percentages.columns:
baseline_summary.loc[group, ('既往用药', f'{history}%')] = med_history_percentages.loc[group, history] if group in med_history_percentages.index else 0
for level in depression_percentages.columns:
baseline_summary.loc[group, ('抑郁程度', f'{level}%')] = depression_percentages.loc[group, level] if group in depression_percentages.index else 0
print("基线特征总结表:")
print(baseline_summary)
baseline_summary.to_excel('问题1结果/问题1_基线特征总结表.xlsx')
# 统计检验结果总结表
age_anova = ols('年龄_岁 ~ C(组别)', data=combined_df).fit()
age_anova_table = sm.stats.anova_lm(age_anova, typ=2)
marriage_chi2, marriage_p, marriage_dof, _ = chi2_contingency(marriage_counts)
med_history_chi2, med_history_p, med_history_dof, _ = chi2_contingency(med_history_counts)
depression_chi2, depression_p, depression_dof, _ = chi2_contingency(depression_counts)
stat_tests = pd.DataFrame({
'基线特征': ['年龄', '婚姻状况', '既往用药', '抑郁程度'],
'统计方法': ['ANOVA', '卡方检验', '卡方检验', '卡方检验'],
'检验统计量': [age_anova_table['F'][0], marriage_chi2, med_history_chi2, depression_chi2],
'p值': [age_anova_table['PR(>F)'][0], marriage_p, med_history_p, depression_p],
'效应量': [eta_squared, cramer_v, cramer_v, cramer_v],
'效应量类型': ['Eta Squared', 'Cramer\'s V', 'Cramer\'s V', 'Cramer\'s V'],
'组间差异显著性': ['显著' if p < 0.05 else '不显著' for p in
[age_anova_table['PR(>F)'][0], marriage_p, med_history_p, depression_p]]
})
stat_tests['检验统计量'] = stat_tests['检验统计量'].map(lambda x: f"{x:.4f}")
stat_tests['p值'] = stat_tests['p值'].map(lambda x: f"{x:.4f}")
stat_tests['效应量'] = stat_tests['效应量'].map(lambda x: f"{x:.4f}")
print("统计检验结果总结:")
print(stat_tests)
stat_tests.to_excel('问题1结果/问题1_统计检验结果总结.xlsx', index=False)
将他转为matlab代码,。同时保证所有功能不变