ANCOVA（协方差分析）是一种结合了方差分析（ANOVA）和线性回归的统计方法，用于在控制一个或多个连续协变量影响的前提下，比较不同组间因变量的均值差异。在Python中，主要可以通过`statsmodels`和`pingouin`这两个库来实现。以下将详细阐述其原理、核心假设、实现步骤，并提供完整的代码示例。 ### 1. ANCOVA 的核心概念与假设 ANCOVA旨在排除协变量（如年龄、基线测量值等）对因变量的影响，从而更纯粹地评估分组变量（处理因素）的效应。其核心假设包括： * **线性关系**：协变量与因变量之间存在线性关系。 * **斜率同质性**：各组的回归线（因变量对协变量的回归）斜率应相同，即不存在分组变量与协变量的交互作用。 * **独立性**：观测值相互独立。 * **正态性**：残差服从正态分布。 * **方差齐性**：各组的残差方差相等。 ### 2. Python 实现方法对比与代码示例 以下通过一个模拟的药物疗效研究案例进行说明：比较两种药物（`Drug_A` 和 `Drug_B`）对患者血压（`BP_After`）的疗效，同时将患者的年龄（`Age`）作为协变量进行控制。 #### 方法一：使用 `statsmodels` 库 `statsmodels` 提供了强大的统计模型接口，可以通过公式 API 构建 ANCOVA 模型。 ```python import pandas as pd import numpy as np import statsmodels.api as sm import statsmodels.formula.api as smf import matplotlib.pyplot as plt import seaborn as sns # 1. 生成模拟数据 np.random.seed(123) n = 50 data = pd.DataFrame({ 'Drug': np.repeat(['Drug_A', 'Drug_B'], n), 'Age': np.random.normal(50, 10, 2*n), # 年龄作为协变量 'BP_Before': np.random.normal(140, 15, 2*n) # 治疗前血压（可选作协变量） }) # 模拟治疗后血压：受药物、年龄和治疗前血压影响，并加入随机噪声 data['BP_After'] = (130 + (data['Drug'] == 'Drug_B') * -8 # Drug_B效果更好 - 0.3 * data['Age'] + 0.6 * data['BP_Before'] + np.random.normal(0, 5, 2*n)) print("数据前5行：") print(data.head()) # 2. 构建并拟合 ANCOVA 模型 # 公式：因变量 ~ 分组变量 + 协变量1 + 协变量2 + ... ancova_model = smf.ols('BP_After ~ Drug + Age + BP_Before', data=data).fit() # 3. 查看模型摘要 print("\nANCOVA 模型摘要：") print(ancova_model.summary()) ``` **关键输出解读**： * **模型整体**：查看 `R-squared` 和 `F-statistic` 及其 `Prob (F-statistic)`，判断模型整体是否显著。 * **分组效应（Drug）**：在 `Drug[T.Drug_B]` 行，`coef` 为 Drug_B 相对于 Drug_A 的均值差异（已调整协变量影响）。`P>|t|` 若小于 0.05，表明在控制年龄和基线血压后，两种药物的疗效存在统计学显著差异 [ref_1][ref_2]。 * **协变量效应**：`Age` 和 `BP_Before` 的系数和 p 值，反映了它们对治疗后血压的独立影响。 #### 方法二：使用 `pingouin` 库 `pingouin` 提供了更简洁、专为心理学和生物统计学设计的接口，并能直接输出方差分析表。 ```python import pingouin as pg # 使用 pingouin 的 ancova 函数 # dv: 因变量， between: 分组变量， covar: 协变量（可以是列表） ancova_result = pg.ancova(data=data, dv='BP_After', between='Drug', covar=['Age', 'BP_Before']) print("ANCOVA 结果表（方差分析格式）：") print(ancova_result) ``` **输出说明**： `pingouin.ancova` 直接返回一个 DataFrame，包含来源（Source）、平方和（SS）、自由度（df）、F值、p值及偏η²（np2）等。这比 `statsmodels` 的回归摘要更贴近传统方差分析报告格式 [ref_5][ref_6]。 ### 3. 假设检验与模型诊断 在执行 ANCOVA 后，必须检验其关键假设是否得到满足。 ```python # 1. 检验斜率同质性（即无交互作用） # 在模型中添加分组变量与协变量的交互项，检验其显著性 interaction_model = smf.ols('BP_After ~ Drug * Age + BP_Before', data=data).fit() print("交互作用模型摘要（检查Drug*Age交互项）：") print(interaction_model.summary()) # 重点关注 'Drug[T.Drug_B]:Age' 交互项的 p 值。若显著，则斜率同质性假设不成立。 # 2. 残差诊断（正态性和方差齐性） residuals = ancova_model.resid fitted = ancova_model.fittedvalues # 绘制残差诊断图 fig, axes = plt.subplots(2, 2, figsize=(12, 10)) # 残差 vs. 拟合值图（检查方差齐性） axes[0, 0].scatter(fitted, residuals, alpha=0.6) axes[0, 0].axhline(y=0, color='r', linestyle='--') axes[0, 0].set_xlabel('拟合值') axes[0, 0].set_ylabel('残差') axes[0, 0].set_title('残差 vs. 拟合值（检查方差齐性）') # Q-Q 图（检查正态性） sm.qqplot(residuals, line='45', ax=axes[0, 1]) axes[0, 1].set_title('Q-Q 图（检查正态性）') # 残差直方图 axes[1, 0].hist(residuals, bins=15, edgecolor='black', alpha=0.7) axes[1, 0].set_xlabel('残差') axes[1, 0].set_ylabel('频数') axes[1, 0].set_title('残差分布直方图') # 按组分组的残差箱线图（检查方差齐性） sns.boxplot(x='Drug', y=residuals, data=data, ax=axes[1, 1]) axes[1, 1].set_title('按药物分组的残差箱线图') plt.tight_layout() plt.show() # 3. 使用统计检验验证假设 # 方差齐性检验（Levene's test） levene_test = pg.homoscedasticity(data, dv='BP_After', group='Drug') print(f"\n方差齐性检验（Levene's test）: p = {levene_test['pval'][0]:.4f}") # 残差正态性检验（Shapiro-Wilk test） shapiro_stat, shapiro_p = pg.normality(residuals) print(f"残差正态性检验（Shapiro-Wilk）: p = {shapiro_p:.4f}") ``` ### 4. 结果可视化：调整后均值的比较 由于 ANCOVA 比较的是调整了协变量影响后的组均值（即最小二乘均值），直接使用原始均值绘图可能产生误导。建议通过模型预测值进行可视化。 ```python # 计算调整后的均值（最小二乘均值） # 创建一个代表“典型”协变量值（如均值）的数据框 typical_covars = pd.DataFrame({ 'Age': [data['Age'].mean()], 'BP_Before': [data['BP_Before'].mean()] }) # 为每种药物生成预测 drugs = ['Drug_A', 'Drug_B'] adjusted_means = {} for drug in drugs: typical_covars['Drug'] = drug adjusted_means[drug] = ancova_model.predict(typical_covars).iloc[0] # 绘制调整后均值的条形图 plt.figure(figsize=(8, 6)) bars = plt.bar(adjusted_means.keys(), adjusted_means.values(), color=['skyblue', 'lightcoral']) plt.ylabel('调整后的平均血压 (BP_After)') plt.title('控制年龄和基线血压后的药物疗效比较（ANCOVA调整后均值）') # 在条形上添加数值标签 for bar, (drug, mean_val) in zip(bars, adjusted_means.items()): plt.text(bar.get_x() + bar.get_width()/2, mean_val + 0.5, f'{mean_val:.1f}', ha='center', va='bottom', fontweight='bold') plt.tight_layout() plt.show() ``` ### 5. 应用场景与注意事项 * **典型场景**： 1. **临床试验**：比较不同治疗方案对结局指标的影响，同时控制基线值、年龄等混杂因素 [ref_1]。 2. **心理学研究**：比较不同干预组在心理量表得分上的差异，控制前测得分或智力等因素 [ref_3]。 3. **观察性研究**：在无法随机分组的条件下，通过统计方法控制已知的混淆变量。 * **注意事项**： 1. **协变量选择**：协变量应与因变量相关，但与分组变量独立（或至少不是处理的结果）。错误地引入中介变量会导致处理效应被低估。 2. **交互作用**：务必检验斜率同质性。如果存在显著的组别×协变量交互作用，说明协变量的影响在不同组中不同，此时标准的 ANCOVA 不再适用，可能需要考虑分层分析或模型中加入交互项。 3. **数据要求**：ANCOVA 对连续型因变量和协变量、分类型分组变量的数据要求严格，需确保满足线性、正态、方差齐性等假设。 综上所述，在Python中执行ANCOVA是一个系统过程，涉及模型构建、假设检验和结果解释。`statsmodels` 提供了更底层的控制和完整的回归诊断，而 `pingouin` 则提供了更直观的方差分析表输出。研究者应根据具体需求和研究设计选择合适的工具，并严格遵守统计分析的规范流程。