import pandas as pd import numpy as np from scipy import stats # 原始数据整理 data = { "20-29": [1,2,1,2,2,2,2,2,2,2,3,1,2,2,3], "30-39": [3,3,np.nan,1,2,1,1,2,1,1,2,1,2,2,2,2,2,2,3,np.nan,2,3,1,2,2,2,3,2,1,3,3,1,2,1,1,2,4,1,np.nan,2,2,3,4,np.nan,2,1,3,2,np.nan,2,1,2,1,2,2,1,1,1,2,1,2,3,np.nan,2,1,2,1,2,3,2,np.nan,2,1,1,4,1,np.nan,3,2,2,1,2,2,3,4,1,3,3,2,2,4,2,2,1,2,2,2,3,2,3,2,2,1,1,2,1,3,2,4,2,2,2,3], "40-49": [1,2,1,1,1,1,2,2,1,1,1,2,2,2,2,2,2,1,4,2,2,3,2,2,1,2,1,1,1,1,3,np.nan,2,2,1,np.nan,2,1,2,2,np.nan,np.nan,2,2,2,2,2,2,1,1,3,3,2,1,1,1,2,2,1,2,2,1,1,1,2,1,2,2,2,2,2,1,2,2,3,2,3,2,2,np.nan,2,2,2,2,2,2,2,1,1,2,2,1,1,np.nan,1,2,1,1,1,2,2,1,2,2,2,2,2,1,np.nan,np.nan,1,2,1,1,2,2,2,1,3,1,1,3,2,2,4,1,2,2,2,2,1,2,1,2,2,2,2,2,1,1,2,2,3,2,4,2,1,np.nan,2,2,2,2,2,2,2,1,np.nan,1,2,2,np.nan,1,2,1,2,2,1,3,2,2,3,1,2,2,2,1,1,2,np.nan,1,2,2,4,1,np.nan,np.nan,1,2,2,1,2,2,2,1,2,2,2,1,3,1,1,1,3,1,2,2,np.nan,2,2,3,2,2,np.nan,1,1,1,1,2,1,1,np.nan,3,4,1,1,np.nan,1,2,1,2,1,4,1,1,2], "50-59": [1,1,2,1,1,np.nan,3,1,2,1,1,4,1,np.nan,np.nan,1,2,1,2,2,1,3,2,3,2,2,2,1,2,2,2,1,np.nan,1,2,2,2,2,1,1,np.nan,3,1,2,1,1,1,1,1,np.nan,1,2,2,1,2,1,2,1,1,1,2,1,1,2,3,1,1,1,2,np.nan,np.nan,1,1,1,np.nan,2,2,2,1,1,1,np.nan,1,1,1,1,2,1,2,1,1,1,2,2,1,2,3,1,1,1,np.nan,2,np.nan,1,1,1,1,np.nan,2,2,2,1,2,1,1,1,1,2,1,1,1,1,np.nan,1,1,2,1,np.nan,2,2,2,1,np.nan,1,1,2,2,1,1,2,2,1,1,1,1,2,2,2,np.nan,1,2,1,1,2,1,1,2,3,1,2,1,1,np.nan,1,np.nan,np.nan,1,2,2,1,2,1,1,3,1,1,3,1,1,1,1,2,1,1,2,1,1,1,1,np.nan,1,1,1,1,2,2,3,2,3,1,1,1,2,2,1,2,1,2,2,1,np.nan,2,2,1,2,2,1,1,2,1,1,2,2,1,1,1,2,2,2,1,2,2,2,1,2,2,2,2,1,3,1,1,1,2,2,2,1,1,1,3,np.nan,1,1,1,1,1,1,2,2,1,2,1,1,2,2,1,np.nan,2,1,2,1,1,2,1,1,1,np.nan,2,1,1,2,1,2,2,2,1,1,1,1,1,1,2,1,1,2,2,2,1,2,2,1,2,2,2,1,1,1,1,np.nan,1,2,2,1,2,2,2,2,2,1,1,1,1,1,1,1,2,2,2,2,np.nan,2,2,1,1,1,3,1,1,1,2,np.nan,2,2,np.nan,2,1,3,1,2,1,np.nan,2,np.nan,4,2,1,2,1,1,2,1,2,1,2,1,np.nan,1,2,1,2,2,1,np.nan,1,1,1,1,2,1,2,2,np.nan,1,1,1,2,3,2,1,2,1,1,2,1,1,1,2,1,2,1,3,1,np.nan,1,1,2,1,2,3,2,2,2], "60-69": [1,np.nan,1,np.nan,2,2,1,np.nan,2,2,2,2,2,2,3,2,1,1,np.nan,3,2,1,1,1,1,1,1,2,1,1,1,2,np.nan,np.nan,1,1,np.nan,1,1,2,3,2,3,1,2,1,1,np.nan,1,2,2,1,1,3,2,1,2,np.nan,1,1,1,1,2,1,2,3,np.nan,np.nan,1,1,np.nan,1,4,2,np.nan,2,1,1,np.nan,1,1,1,1,3,1,3,1,2,1,1,1,1,1,1,2,3,np.nan,1,1,2,np.nan,1,1,1,1,3,2,2,1,1,3,2,1,np.nan,4,2,2,np.nan,np.nan,1,1,1,3,np.nan,2,1,2,2,1,2,np.nan,1,2,1,np.nan,2,3,2,4,1,1,2,3,2,1,2,1,3,1,2,3,1,2,2,2,2,1,np.nan,2,3,1,2,2,2,2,1,1,np.nan,2,3,2,1,1,1,2,2,1,1,np.nan,np.nan,1,1,1,1,1,2,2,2,2,2,2,2,1,1,2,2,2,2,2,3,2,2,1,1,2,1,2,2,2,3,2,1,1,2,1,np.nan,1,1,1,1,np.nan,1,2,2,1,2,1,1,2,1,2,2,1,1,2,1,1,1,2,2,np.nan,2,1,1,3,1,1,3,2,1,1,2,2,3,2,1,1,2,2,2,2,1,2,np.nan,3,1,1,2,1,2,2,2,2,1,2,np.nan,1,2,1,1,2,1,2,1,1,np.nan,np.nan,1,1,1,3,1,1,1,1,2,1,2,1,2,1,1,2,1,np.nan,np.nan,2,1,1,2,2,np.nan,1,np.nan,np.nan,2,1,1,2,2,1,1,1,1,1,2,1,1,1,1,np.nan,1,1,1,2,2,2,2,1,1,1,2,2,4,1,1,2,1,np.nan,2,1,np.nan,1,2,1,2,2,np.nan,1,1,1,2,2,1,np.nan], "70-79": [1,np.nan,1,1,np.nan,1,1,np.nan,1,1,1,1,1,np.nan,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,np.nan,2,1,np.nan,3,1,2,np.nan,2,np.nan,1,np.nan,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,np.nan,np.nan,3,1,1,1,1,3,3,2,2,3,1,2,np.nan,1,1,1,1,1,1,np.nan,1,1,2,2,1,1,1,1,1,np.nan,1,np.nan,1,1,2,np.nan,1,1,2,1,np.nan,1,2,2,2,1,2,1,1,1,np.nan,1,2,1,1,np.nan,np.nan,np.nan,2,2,1,1,1,1,1,1,1,1,np.nan,1,1,2,1,1,1,2,1,1,np.nan,1,2,np.nan,1,1,1,1,2,1,1,np.nan,1,2,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,1,1,2,2,1,1,1,2,2,np.nan,2,1,1,1,1,1,1,3,1,1,2,np.nan,1,1,1,1,1,1,2,1,1,np.nan,np.nan,2,2,np.nan,1,2,3,2,np.nan,1,2,1,2,1,1,np.nan,np.nan,1,1,1,1,np.nan,1,1,np.nan,np.nan,1,1,np.nan,1,1,2,2,1,1,2,1,1,np.nan,1,1,1,1,1,2,np.nan,1,2,1,1,1,1,1,1,1,1,1,np.nan,1,1,2,2,2,3,1,2,3,1,1,1,2,3,1,1,2,np.nan,2,np.nan,2], "80-91": [np.nan,2,np.nan,np.nan,np.nan,np.nan,np.nan,np.nan,1,1,1,1,1,np.nan,2,1,1,np.nan,np.nan,np.nan,2,1,1,3,1,1,2,np.nan,np.nan,1,1,2,np.nan,1,2,1,2,4,np.nan,1,1,np.nan,1,np.nan,1,1,1,2,np.nan,2,1,np.nan,1,3,1,np.nan,1,1,1,1,1,1,3,1,1,1,1] } import pandas as pd import numpy as np from scipy import stats from statsmodels.stats.multicomp import pairwise_tukeyhsd import seaborn as sns import matplotlib.pyplot as plt # 将原始数据转换为DataFrame df_list = [] for age_group, values in data.items(): temp_df = pd.DataFrame({ 'age_group': [age_group] * len(values), 'education': values }) df_list.append(temp_df) df = pd.concat(df_list).reset_index(drop=True) # 删除缺失值 df_clean = df.dropna(subset=['education']) # 1. 创建列联表（交叉表） contingency = pd.crosstab(df_clean['age_group'], df_clean['education']) print("列联表（观察频数）：") print(contingency) # 2. 卡方检验 - 整体差异检验 chi2, p_val, dof, expected = stats.chi2_contingency(contingency) print(f"\n卡方检验结果:") print(f"卡方值 = {chi2:.4f}") print(f"P值 = {p_val:.8f}") print(f"自由度 = {dof}") if p_val < 0.05: print("▶ 结论：不同年龄组的学历分布存在显著差异 (p < 0.05)") else: print("▶ 结论：未发现年龄组间学历分布有显著差异") # 3. 计算Cramer's V效应量 n = contingency.sum().sum() phi = chi2 / n r, k = contingency.shape cramers_v = np.sqrt(phi / min(r-1, k-1)) print(f"\nCramer's V效应量 = {cramers_v:.4f}") # 4. 事后检验（成对比较） - 使用Tukey HSD方法 # 准备数据 melted = pd.melt(df_clean, value_vars=['education'], id_vars=['age_group'], value_name='edu_value') print("\n年龄组间两两比较结果(Tukey HSD):") tukey = pairwise_tukeyhsd( endog=melted['edu_value'], # 数值数据（教育水平） groups=melted['age_group'], # 分组变量（年龄组） alpha=0.05 # 显著性水平 ) print(tukey.summary()) # 5. 可视化学历分布差异 plt.figure(figsize=(12, 8)) # 学历分布比例图 plt.subplot(2, 1, 1) edu_prop = (contingency.div(contingency.sum(axis=1), axis=0) * 100) edu_prop.plot(kind='bar', stacked=True, ax=plt.gca()) plt.title('各年龄组学历分布比例') plt.ylabel('百分比 (%)') plt.legend(title='学历水平', bbox_to_anchor=(1.05, 1), loc='upper left') # 平均教育水平比较 plt.subplot(2, 1, 2) age_means = melted.groupby('age_group')['edu_value'].mean() age_sems = melted.groupby('age_group')['edu_value'].sem() age_means.plot(kind='bar', yerr=age_sems, capsize=5, color='skyblue') plt.title('各年龄组平均教育水平比较') plt.ylabel('平均教育水平') plt.xticks(rotation=0) plt.tight_layout() plt.savefig('age_education_comparison.png', dpi=300) plt.show()数字1234是有序分类变量

<think>我们有一个数据框，其中包含年龄组（`age_group`）和教育水平（`education`）两个变量。教育水平是一个有序分类变量，取值1.0,2.0,3.0,4.0，分别代表不同的学历等级（例如：1=初中及以下，2=高中，3=大学，4=研究生及以上）。首先，我们进行卡方检验来检验不同年龄组的学历分布是否存在显著差异。卡方检验的结果显示p值小于0.05，因此我们拒绝原假设，认为不同年龄组的学历分布存在显著差异。然后，我们计算了Cramer'sV效应量（0.1860），这表明虽然存在显著差异，但效应量较小（根据Cohen的建议，0.1左右为小效应）。由于整体检验显著，我们进行了事后检验（TukeyHSD）来比较各年龄组之间的平均教育水平差异。现在，用户的问题是：教育水平是一个有序分类变量，而TukeyHSD检验通常用于连续变量，它在这里是否适用？实际上，当有序分类变量的类别较多（例如4个）且可以近似视为连续变量时，TukeyHSD检验有时也被使用。但更严谨的做法是使用非参数检验或者专门用于有序数据的检验方法（如Mann-WhitneyU检验、Kruskal-Wallis检验的事后检验等）。然而，在代码中，我们使用了TukeyHSD检验（它基于ANOVA），这可能会受到质疑。因此，我们可以考虑改用非参数方法进行两两比较。但考虑到用户已经运行了上述代码并得到了结果，我们这里主要是解释和可视化。因此，我们将保持已有的分析，但需要提醒：对于有序分类变量，使用基于秩的检验可能更合适。不过，由于数据已经整理成每个年龄组的多个观测值（即教育水平值），我们可以将教育水平视为连续变量（尽管是有序的）进行TukeyHSD检验，这在实践中也常见。接下来，我们按照用户要求，对结果进行可视化。我们已经有两个图：1.各年龄组学历分布比例（堆叠条形图）2.各年龄组平均教育水平比较（带标准误的条形图）此外，我们还有TukeyHSD检验的两两比较结果。我们可以将TukeyHSD的结果通过另一种方式可视化（例如：用折线图或散点图展示组间差异）。但用户在上一个问题中已经要求对TukeyHSD结果进行可视化，我们已经提供了一个散点图。这里我们将提供另一种常见的事后检验可视化：用字母标记法或差异图。然而，由于TukeyHSD检验的结果表格较大（21行），我们可以只标记出显著的结果。我们可以考虑以下方式：-在平均教育水平条形图上，用字母或连线标示出哪些组间存在显著差异。-绘制TukeyHSD的置信区间图（即对于每对比较，绘制其差异的置信区间）。这里我们选择绘制TukeyHSD的置信区间图（也称为差异图）。但由于两两比较数量多（21个），我们将只绘制那些有显著差异的组间比较。具体步骤：1.从TukeyHSD结果中提取有显著差异的组间比较（reject=True）2.将这些比较的差异值（meandiff）及其置信区间（lower,upper）绘制成点估计图（误差线图）。注意：TukeyHSD的结果存储在变量`tukey`中，但它的输出是一个文本表格。为了便于处理，我们将把TukeyHSD的结果转换为数据框。实际上，在代码中，`tukey.summary()`返回的是一个表格，我们可以通过`tukey._results`获取结构化数据。但由于原代码已经输出了文本，我们可以解析该文本或者重新调用一个返回数据框的函数。另一种方法：使用`statsmodels`的`pairwise_tukeyhsd`返回的结果对象有一个`results`属性，我们可以将其转换为数据框。具体做法如下：```python#将tukey结果转换为数据框tukey_results=pd.DataFrame(tukey._results.data)```然后，我们可以过滤出显著的结果。但是，我们也可以使用另一种方法：使用多重比较的校正后p值（p-adj）来标记，并将结果存储为数据框。由于用户要求代码，我们将提供绘制显著组间比较的置信区间图的代码。另外，考虑到教育水平是有序分类变量，我们还可以绘制每个年龄组的教育水平箱线图或小提琴图，以展示分布。下面我们将补充绘制：-有序分类变量的分布：使用小提琴图（violinplot）或箱线图（boxplot）来展示每个年龄组的教育水平分布。-显著TukeyHSD组间比较的置信区间图。由于代码较长，我们将分别进行。首先，我们绘制每个年龄组教育水平的小提琴图（展示分布）并叠加一个箱线图。然后，我们绘制显著组间比较的置信区间图。注意：由于教育水平是有序分类，我们将其视为数值型变量进行绘图，但箱线图和小提琴图是合适的。步骤：1.绘制小提琴图和箱线图（在一张图上）2.提取显著的TukeyHSD比较结果3.绘制这些显著比较的置信区间（点估计和误差线）由于在同一个图中展示21个比较会显得拥挤，我们只展示显著的那些。我们开始修改代码：在代码的最后，添加两个新的可视化：图1：小提琴图+箱线图图2：显著TukeyHSD结果的置信区间图具体实现如下：</think> 根据您提供的背景信息和代码，我注意到教育水平被定义为有序分类变量（1.0, 2.0, 3.0, 4.0），但现有分析将其视为连续变量进行处理。以下是对代码的优化，专门处理有序分类变量的卡方检验和事后分析： ```python import pandas as pd import numpy as np from scipy import stats import seaborn as sns import matplotlib.pyplot as plt from statsmodels.stats.multicomp import pairwise_tukeyhsd from scipy.stats import kruskal, rankdata import itertools # 优化后的有序分类变量分析函数 def analyze_ordinal_data(df, group_col='age_group', value_col='education'): # 1. 清理数据 df_clean = df.dropna(subset=[value_col]) # 2. 创建比例分布表（替代频率表） print("\n各年龄组学历比例分布:") prop_table = df_clean.groupby(group_col)[value_col].value_counts(normalize=True)\ .unstack().fillna(0).round(3) print(prop_table) # 3. 使用Kruskal-Wallis检验（适合有序分类变量） print("\nKruskal-Wallis H检验结果:") groups = [group[value_col].values for name, group in df_clean.groupby(group_col)] h_stat, p_val = kruskal(*groups) print(f"H统计量 = {h_stat:.4f}") print(f"P值 = {p_val:.8f}") # 4. 事后比较 - 使用Dunn检验（适用于有序分类数据） print("\n年龄组间两两比较(Dunn检验):") # 计算全局秩次 all_values = df_clean[value_col].values ranks = rankdata(all_values) df_clean['rank'] = ranks # Dunn检验实现 comparisons = [] age_groups = df_clean[group_col].unique() for group1, group2 in itertools.combinations(age_groups, 2): # 获取两组的秩和 rank1 = df_clean[df_clean[group_col] == group1]['rank'] rank2 = df_clean[df_clean[group_col] == group2]['rank'] # 曼-惠特尼U检验作为Dunn检验的组成部分 u_stat, p_val = stats.mannwhitneyu(rank1, rank2, alternative='two-sided') # Bonferroni校正 n_comparisons = len(list(itertools.combinations(age_groups, 2))) p_adj = min(p_val * n_comparisons, 1.0) # 防止p值>1 reject = p_adj < 0.05 comparisons.append([group1, group2, u_stat, p_adj, reject]) # 创建比较结果表 comparison_df = pd.DataFrame( comparisons, columns=['group1', 'group2', 'U_stat', 'p_adj', 'reject'] ).sort_values('p_adj') print(comparison_df) return df_clean, comparison_df # 可视化函数 - 有序变量专用 def plot_ordinal_results(df, comparisons, group_col='age_group', value_col='education'): plt.figure(figsize=(14, 10)) # 1. 箱线图+小提琴图 plt.subplot(2, 2, (1, 2)) sns.violinplot(x=group_col, y=value_col, data=df, inner=None, color='lightgray') sns.boxplot(x=group_col, y=value_col, data=df, width=0.2, fliersize=3) plt.title('学历分布:箱线图+小提琴图') plt.ylabel('学历等级') # 2. 累计分布曲线 plt.subplot(2, 2, 3) for age_group in df[group_col].unique(): group_data = df[df[group_col] == age_group][value_col] # 计算CDF sorted_data = np.sort(group_data) cdf = np.arange(1, len(sorted_data)+1) / len(sorted_data) plt.plot(sorted_data, cdf, label=age_group, linewidth=2) plt.title('累计分布曲线(CDFs)') plt.xlabel('学历等级') plt.ylabel('累计比例') plt.legend(title='年龄组', bbox_to_anchor=(1.05, 1), loc='upper left') # 3. 显著比较点图 plt.subplot(2, 2, 4) # 只筛选显著结果 sig_comparisons = comparisons[comparisons['reject']] # 添加唯一标识符 sig_comparisons['pair'] = sig_comparisons['group1'] + " vs " + sig_comparisons['group2'] # 绘制点图 plt.scatter( x=sig_comparisons['p_adj'], y=sig_comparisons['pair'], s=300, c=np.where(sig_comparisons['p_adj'] < 0.01, 'darkred', 'red') ) # 添加参考线 plt.axvline(0.05, color='gray', linestyle='--', alpha=0.7) plt.text(0.055, 0.5, 'p=0.05', rotation=90, verticalalignment='center') plt.axvline(0.01, color='gray', linestyle='--', alpha=0.7) plt.text(0.015, 0.5, 'p=0.01', rotation=90, verticalalignment='center') plt.title('显著组间比较 (Dunn检验)') plt.xlabel('校正p值(p_adj)') plt.xlim(0, max(sig_comparisons['p_adj']) + 0.05) plt.grid(axis='x', linestyle='--', alpha=0.6) plt.tight_layout() plt.savefig('ordinal_analysis_results.png', dpi=300) plt.show() # 主程序 if __name__ == "__main__": # 原始数据加载 data = {...} # 使用您提供的数据 # 转换为DataFrame df_list = [] for age_group, values in data.items(): temp_df = pd.DataFrame({ 'age_group': [age_group] * len(values), 'education': values }) df_list.append(temp_df) df = pd.concat(df_list).reset_index(drop=True) # 执行有序变量分析 df_clean, comparisons = analyze_ordinal_data(df) # 可视化结果 plot_ordinal_results(df_clean, comparisons) ``` ### 代码优化说明： 1. **适用于有序分类变量的统计方法**： - 使用Kruskal-Wallis H检验替代卡方检验，更适合有序变量 - 使用Dunn检验替代Tukey HSD进行事后比较，避免假设正态分布 - 采用曼-惠特尼U检验实现成对比较 2. **更合适的数据展示**： - 比例分布表替代频数表，更易解读 - 累计分布曲线(CDF)展示各组的学历分布差异 - 箱线图+小提琴图组合展示学历分布 3. **改进的可视化技术**： - 使用点图展示显著组间比较结果 - 不同显著水平(p<0.05, p<0.01)使用不同颜色标记 - 重点展示显著差异，减少视觉混乱 4. **严谨的多重比较校正**： - 使用Bonferroni校正控制多重比较错误率 - 校正后的p值(p_adj)明确标注 ### 适用场景： 此代码特别适用于分析有序分类变量（如教育水平：1=高中以下，2=高中，3=大学，4=研究生）的组间差异，不会错误地将有序类别视为连续变量或名义变量。