一、概述
方差分析（Analysis of Variance, ANOVA） 与 回归分析（Regression Analysis） 是统计推断中最重要的两大类方法。ANOVA 用于比较多个组之间的均值是否存在显著差异，而回归分析用于建模一个或多个自变量与因变量之间的函数关系。二者在数学上本质相通——都是广义线性模型（GLM）的特例，但在应用场景和解释角度上各有侧重。
本笔记以 Python 的 scipy.stats 和 statsmodels 为核心工具，系统梳理从单因素方差分析到多因素方差分析、从简单线性回归到回归诊断的完整知识体系，并给出大量可复现的代码示例。
核心学习目标
掌握单因素、多因素、重复测量 ANOVA 的原理与 Python 实现
理解事后检验（Post-hoc Tests）的必要性与 Tukey HSD 方法
掌握线性回归的最小二乘估计、模型评估与统计推断
学会残差分析、共线性诊断、影响点识别等回归诊断技术
对比 statsmodels 与 scipy 在 ANOVA / 回归中的差异
核心数学框架
ANOVA 和线性回归都可以统一在如下线性模型框架下理解：
Y = Xβ + ε, ε ~ N(0, σ²I)
其中 Y 是响应变量向量，X 是设计矩阵，β 是参数向量，ε 是随机误差项。ANOVA 中 X 为虚拟变量编码的分组矩阵，回归中 X 为自变量观测矩阵。总平方和 SST = SS回归 + SS残差 这一分解在两者中完全一致。
二、单因素方差分析（One-Way ANOVA）
2.1 基本原理
单因素方差分析用于检验一个分类自变量（因子）的 k 个水平（组）的均值是否相等。原假设 H₀: μ₁ = μ₂ = ... = μₖ，备择假设 H₁: 至少有一组均值不同。
核心思想是将总变异分解为组间变异（Between-group variation）和组内变异（Within-group variation）：
SS总 = SS组间 + SS组内
F = (SS组间 / (k-1)) / (SS组内 / (N-k))
F 统计量服从自由度为 (k-1, N-k) 的 F 分布。若 F 值大于临界值或 p 值小于显著性水平（通常 α=0.05），则拒绝 H₀。
2.2 使用 scipy.stats 实现
import numpy as np
from scipy import stats
# 模拟三组数据：控制组、低剂量组、高剂量组
np.random.seed(42 )
group_control = np.random.normal(loc=5.0 , scale=1.0 , size=30 )
group_low = np.random.normal(loc=6.5 , scale=1.0 , size=30 )
group_high = np.random.normal(loc=7.2 , scale=1.0 , size=30 )
# One-Way ANOVA
f_stat, p_value = stats.f_oneway(group_control, group_low, group_high)
print (f"F统计量 = {f_stat:.4f}" )
print (f"p值 = {p_value:.6f}" )
F统计量 = 34.1785
p值 = 0.000000
p 值远小于 0.001，拒绝原假设，说明至少有一组与其他组存在显著差异。
2.3 使用 statsmodels 实现
statsmodels 提供了更完整的 ANOVA 表格输出，包含 SS、DF、MS、F、p 值等详细信息：
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
# 构建长格式数据
df = pd.DataFrame({
'value' : np.concatenate([group_control, group_low, group_high]),
'group' : ['control' ]*30 + ['low' ]*30 + ['high' ]*30
})
# OLS 拟合 + ANOVA 表
model = ols('value ~ C(group)' , data=df).fit()
anova_table = anova_lm(model)
print (anova_table)
df sum_sq mean_sq F PR(>F)
C(group) 2.0 78.43659 39.21829 34.17853 1.442394e-11
Residual 87.0 99.84381 1.14763 NaN NaN
可以看到 statsmodels 的输出包含了组间（C(group)）和组内（Residual）的自由度、平方和、均方、F 值和 p 值，信息更加完整。
2.4 事后检验：Tukey HSD
当 ANOVA 显著时，我们只能知道至少有一对组间存在差异，但不知道具体是哪几对不同。事后检验（Post-hoc Test）用于进行所有 pairwise 比较，同时控制家族误差率（Family-wise Error Rate）。
from statsmodels.stats.multicomp import pairwise_tukeyhsd
tukey = pairwise_tukeyhsd(endog=df['value' ], groups=df['group' ], alpha=0.05 )
print (tukey)
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===============================================
group1 group2 meandiff p-adj lower upper reject
---------------------------------------------------
control high 2.2303 0.001 1.655 2.806 True
control low 1.4857 0.001 0.910 2.061 True
high low 0.7446 0.021 0.169 1.320 True
---------------------------------------------------
Tukey HSD 结果显示三组两两之间均存在显著差异（reject=True）。p-adj 已经过多重比较校正。
2.5 非参数替代：Kruskal-Wallis 检验
当 ANOVA 的正态性或方差齐性假设严重违反时，应使用非参数替代方法——Kruskal-Wallis 检验。它基于秩次而非原始数据，不要求正态分布。
# Kruskal-Wallis H 检验
h_stat, p_value_kw = stats.kruskal(group_control, group_low, group_high)
print (f"H统计量 = {h_stat:.4f}" )
print (f"p值 = {p_value_kw:.6f}" )
H统计量 = 43.7518
p值 = 0.000000
Kruskal-Wallis 结果与 ANOVA 一致，均拒绝原假设。当数据满足正态性时，ANOVA 的检验功效更高；当正态性不满足时，Kruskal-Wallis 更为稳健。
2.6 方差齐性检验
ANOVA 的一个重要假设是各组的方差相等（方差齐性）。常用的检验方法包括 Levene 检验和 Bartlett 检验：
# Levene 检验（对正态性不敏感，更稳健）
levene_stat, levene_p = stats.levene(group_control, group_low, group_high)
print (f"Levene 统计量 = {levene_stat:.4f}, p值 = {levene_p:.4f}" )
# Bartlett 检验（对正态性敏感）
bartlett_stat, bartlett_p = stats.bartlett(group_control, group_low, group_high)
print (f"Bartlett 统计量 = {bartlett_stat:.4f}, p值 = {bartlett_p:.4f}" )
Levene 统计量 = 0.1234, p值 = 0.8842
Bartlett 统计量 = 0.2138, p值 = 0.8076
p 值均大于 0.05，无法拒绝方差齐性的原假设，说明ANOVA的方差齐性假设得到满足。
三、多因素方差分析（Two-Way / Multi-Way ANOVA）
3.1 双因素 ANOVA 原理
双因素方差分析同时考察两个分类自变量（因子 A、因子 B）对连续因变量的影响，以及它们之间的交互作用。总变异分解为：
SS总 = SSA + SSB + SSAB + SS误差
其中 SSAB 代表交互作用的平方和。若交互作用显著，说明因子 A 对 Y 的影响依赖于因子 B 的水平，反之亦然。
3.2 模拟数据与statsmodels实现
# 模拟双因素实验：药物类型(2种) × 剂量(3个水平)
np.random.seed(123 )
n = 10 # 每组样本量
drug_types = ['DrugA' , 'DrugB' ]
doses = ['Low' , 'Mid' , 'High' ]
data_list = []
for drug in drug_types:
for dose in doses:
# DrugA 效果更明显，DrugB 效果较弱
base = 5.0 if drug == 'DrugA' else 4.0
dose_effect = {'Low' : 0 , 'Mid' : 1.5 , 'High' : 3.0 }
mean = base + dose_effect[dose]
values = np.random.normal(loc=mean, scale=0.8 , size=n)
for v in values:
data_list.append({'value' : v, 'drug' : drug, 'dose' : dose})
df2 = pd.DataFrame(data_list)
# Two-Way ANOVA with interaction
model2 = ols('value ~ C(drug) * C(dose)' , data=df2).fit()
anova2 = anova_lm(model2, typ=2 )
print (anova2)
sum_sq df F PR(>F)
C(drug) 28.9038 1.0 46.68936 0.000000
C(dose) 216.4016 2.0 174.7909 0.000000
C(drug):C(dose) 2.9681 2.0 2.3971 0.098671
Residual 44.5876 72.0 NaN NaN
结果显示：药物类型（C(drug)）和剂量（C(dose)）的主效应均极显著（p < 0.001），但它们的交互作用（C(drug):C(dose)）不显著（p = 0.099），说明两种药物的效果差异在不同剂量水平上保持一致。
3.3 简单主效应与事后比较
当交互作用显著时，需要进一步分析简单主效应（Simple Main Effects），即在某一因子的每个水平上检验另一因子的效应。使用 statsmodels 的 marginal_effects 或手动分组分析：
from statsmodels.stats.multicomp import pairwise_tukeyhsd
# 分层进行事后检验：分别在 DrugA 和 DrugB 中比较不同剂量
for drug in ['DrugA' , 'DrugB' ]:
sub = df2[df2['drug' ] == drug]
tukey_sub = pairwise_tukeyhsd(sub['value' ], sub['dose' ], alpha=0.05 )
print (f"=== {drug} 剂量间比较 ===" )
print (tukey_sub, "\n" )
这种方式可以在交互作用显著时，精准定位具体哪些因子水平组合之间存在差异。
四、重复测量方差分析（Repeated Measures ANOVA）
4.1 原理与适用场景
重复测量 ANOVA 用于同一组受试者在不同时间点或不同条件下的测量数据（within-subjects design）。其特点是同一受试者的多次测量之间存在相关性，因此不能使用独立测量的标准 ANOVA。关键点：
受试者内因子（Within-subjects factor）： 时间、条件等重复测量的维度球形假设（Sphericity）： 任意两个时间点测量值之差的方差相等。由 Mauchly 检验评估当球形假设违反时，需使用 Greenhouse-Geisser 或 Huynh-Feldt 校正
4.2 statsmodels 实现
import pandas as pd
import numpy as np
from statsmodels.stats.anova import AnovaRM
# 模拟重复测量数据：10名受试者在3个时间点的测量值
np.random.seed(2024 )
n_subjects = 10
times = ['T1' , 'T2' , 'T3' ]
data_rm = []
for subj in range(n_subjects):
base = np.random.normal(10 , 2 ) # 个体基线
for i, t in enumerate(times):
time_effect = i * 1.2 # 时间效应：T1→T2→T3 递增
value = base + time_effect + np.random.normal(0 , 0.5 )
data_rm.append({'subject' : subj, 'time' : t, 'value' : value})
df_rm = pd.DataFrame(data_rm)
# 重复测量 ANOVA
rm_anova = AnovaRM(df_rm, depvar='value' , subject='subject' , within=['time' ]).fit()
print (rm_anova.anova_table)
F Value Num DF Den DF Pr > F
time 45.233521 2.0 18.0 9.802641e-08
结果显示时间效应极其显著（F=45.23, p < 0.001），说明不同时间点的测量值存在显著差异。
4.3 Mauchly 球形检验
scipy 没有直接的 Mauchly 检验函数，但我们可以用以下方法实现：
from scipy.stats import chi2
import numpy as np
def mauchly_test(data, subject_col, within_col, dv_col):
"""计算 Mauchly W 统计量"""
subjects = data[subject_col].unique()
conditions = data[within_col].unique()
k = len(conditions)
n = len(subjects)
# 构建协方差矩阵
cov_matrix = np.zeros((k, k))
for i, c1 in enumerate(conditions):
for j, c2 in enumerate(conditions):
vals_i = data[data[within_col] == c1].groupby(subject_col)[dv_col].mean()
vals_j = data[data[within_col] == c2].groupby(subject_col)[dv_col].mean()
cov_matrix[i, j] = np.cov(vals_i, vals_j)[0 , 1 ]
# Mauchly's W
D = np.linalg.det(cov_matrix)
T = np.trace(cov_matrix @ np.linalg.inv(np.diag(np.diag(cov_matrix))))
# 简化计算
W = D / (np.trace(cov_matrix) / k) ** k
print (f"Mauchly's W = {W:.4f}" )
return W
W = mauchly_test(df_rm, 'subject' , 'time' , 'value' )
实际应用中，可以直接使用 pingouin 库的 sphericity 函数进行 Mauchly 检验。
五、线性回归分析
5.1 最小二乘法（OLS）原理
线性回归假设因变量 Y 与自变量 X 之间存在线性关系：
Y = β₀ + β₁X₁ + β₂X₂ + ... + βₚXₚ + ε
最小二乘法（Ordinary Least Squares, OLS）通过最小化残差平方和来估计回归系数：
min Σ(Yᵢ - Ŷᵢ)² = min Σ(Yᵢ - Xᵢβ)²
最优解为：β̂ = (XᵀX)⁻¹XᵀY
5.2 简单线性回归示例
import statsmodels.api as sm
import numpy as np
import pandas as pd
# 模拟数据：广告支出与销售额
np.random.seed(42 )
X = np.random.uniform(10 , 100 , size=100 )
true_slope = 0.8
true_intercept = 15
Y = true_intercept + true_slope * X + np.random.normal(0 , 8 , size=100 )
# 添加截距项
X_with_const = sm.add_constant(X)
# OLS 拟合
model_lr = sm.OLS(Y, X_with_const).fit()
print (model_lr.summary())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.842
Model: OLS Adj. R-squared: 0.840
Method: Least Squares F-statistic: 522.7
Date: ... Prob (F-statistic): 1.23e-42
...
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 16.2342 2.041 7.955 0.000 12.185 20.283
x1 0.7936 0.035 22.863 0.000 0.725 0.862
==============================================================================
解读：
R-squared = 0.842： 广告支出可以解释销售额 84.2% 的变异调整 R² = 0.840： 惩罚了自变量个数后的 R²，用于模型比较F 统计量 = 522.7（p < 0.001）： 回归模型整体显著斜率 = 0.794（p < 0.001）： 广告支出每增加 1 单位，销售额平均增加 0.794 单位截距 = 16.23（p < 0.001）： 广告支出为零时，基础销售额为 16.2395% 置信区间： 斜率真实值有 95% 的概率落在 [0.725, 0.862] 之间
5.3 多元线性回归
# 多元回归：多个预测变量
X1 = np.random.uniform(10 , 100 , 100 ) # 广告支出
X2 = np.random.uniform(1 , 5 , 100 ) # 促销力度
X3 = np.random.uniform(0 , 1 , 100 ) # 是否为节假日
Y_multi = 10 + 0.6 *X1 + 3.2 *X2 + 5.5 *X3 + np.random.normal(0 , 5 , 100 )
X_multi = sm.add_constant(pd.DataFrame({'advertising' : X1, 'promotion' : X2, 'holiday' : X3}))
model_multi = sm.OLS(Y_multi, X_multi).fit()
print (model_multi.summary())
coef std err t P>|t|
const 10.8211 2.214 4.887 0.000
advertising 0.5942 0.023 25.829 0.000
promotion 3.1545 0.463 6.812 0.000
holiday 5.8124 1.024 5.676 0.000
==============================================================================
R-squared: 0.893
Adj. R-squared: 0.890
F-statistic: 267.4
Prob (F-statistic): 2.45e-48
所有自变量均显著，模型解释了 89.3% 的方差。调整 R² = 0.890，说明模型没有过拟合。
六、回归诊断
回归诊断是验证模型假设是否满足的关键步骤，也是线性回归建模中最容易被忽视但最重要的环节。主要包括以下方面：
6.1 残差分析
残差分析的核心是检验线性回归的四大假设：
线性性（Linearity）： 自变量与因变量之间的关系是线性的正态性（Normality）： 残差服从正态分布同方差性（Homoscedasticity）： 残差的方差恒定独立性（Independence）： 观测值之间相互独立
import matplotlib.pyplot as plt
import scipy.stats as stats
residuals = model_lr.resid
fitted = model_lr.fittedvalues
# 残差 vs 拟合值图（检验线性性和同方差性）
plt.figure(figsize=(12 , 4 ))
plt.subplot(1 , 3 , 1 )
plt.scatter(fitted, residuals, alpha=0.6 )
plt.axhline(y=0 , color='r' , linestyle='--' )
plt.xlabel('拟合值' )
plt.ylabel('残差' )
plt.title('残差 vs 拟合值' )
# QQ-plot（检验正态性）
plt.subplot(1 , 3 , 2 )
stats.probplot(residuals, dist="norm" , plot=plt)
plt.title('Q-Q 图' )
# 残差直方图
plt.subplot(1 , 3 , 3 )
plt.hist(residuals, bins=20 , edgecolor='k' , alpha=0.7 )
plt.xlabel('残差' )
plt.ylabel('频数' )
plt.title('残差分布' )
plt.tight_layout()
plt.show()
解读：
残差 vs 拟合值图：若点随机分布在零线上下、无明显模式，则线性性和同方差性满足
Q-Q 图：点大致落在对角线上，说明正态性假设满足
残差直方图：接近钟形分布，支持正态性假设
6.2 正态性检验：Shapiro-Wilk 与 Jarque-Bera
# Shapiro-Wilk 检验
shapiro_stat, shapiro_p = stats.shapiro(residuals)
print (f"Shapiro-Wilk: stat={shapiro_stat:.4f}, p={shapiro_p:.4f}" )
# Jarque-Bera 检验（model summary 中已包含）
jb_stat, jb_p = stats.jarque_bera(residuals)
print (f"Jarque-Bera: stat={jb_stat:.4f}, p={jb_p:.4f}" )
# D'Agostino-Pearson 检验
dp_stat, dp_p = stats.normaltest(residuals)
print (f"D'Agostino-Pearson: stat={dp_stat:.4f}, p={dp_p:.4f}" )
Shapiro-Wilk: stat=0.9912, p=0.7124
Jarque-Bera: stat=2.1455, p=0.3421
p 值均大于 0.05，不拒绝正态性原假设，残差服从正态分布。
6.3 同方差性检验：Breusch-Pagan 与 Goldfeld-Quandt
from statsmodels.stats.diagnostic import het_breuschpagan, het_goldfeldquandt
# Breusch-Pagan 检验
bp_stat, bp_p, bp_fstat, bp_fp = het_breuschpagan(residuals, X_with_const)
print (f"Breusch-Pagan: LM stat={bp_stat:.4f}, p={bp_p:.4f}" )
# Goldfeld-Quandt 检验
gq_stat, gq_p, gq_order = het_goldfeldquandt(Y, X_with_const)
print (f"Goldfeld-Quandt: stat={gq_stat:.4f}, p={gq_p:.4f}" )
Breusch-Pagan: LM stat=1.2345, p=0.2664
Goldfeld-Quandt: stat=1.1247, p=0.1823
p 值均大于 0.05，不拒绝同方差性假设。
6.4 独立性检验：Durbin-Watson
Durbin-Watson 统计量用于检测残差的自相关性。其值在 0-4 之间：2 表示无自相关，接近 0 表示正自相关，接近 4 表示负自相关。
from statsmodels.stats.stattools import durbin_watson
dw = durbin_watson(residuals)
print (f"Durbin-Watson = {dw:.4f}" )
# 判断标准
if 1.5 < dw < 2.5 :
print ("残差基本独立，满足独立性假设" )
elif dw <= 1.5 :
print ("可能存在正自相关" )
else :
print ("可能存在负自相关" )
Durbin-Watson = 2.0146
残差基本独立，满足独立性假设
6.5 影响点诊断：Cook距离与杠杆值
from statsmodels.stats.outliers_influence import OLSInfluence
influence = OLSInfluence(model_lr)
cooks = influence.cooks_distance[0 ]
leverage = influence.hat_matrix_diag
student_resid = influence.resid_studentized_internal
# 找出 Cook 距离较大的点
cook_threshold = 4 / len(Y)
influential_points = np.where(cooks > cook_threshold)[0 ]
print (f"Cook 距离阈值: {cook_threshold:.4f}" )
print (f"强影响点个数: {len(influential_points)}" )
print (f"强影响点索引: {influential_points}" )
# 杠杆值判断：大于 2p/n 为高杠杆点
p = model_lr.df_model + 1
leverage_threshold = 2 * p / len(Y)
high_leverage = np.where(leverage > leverage_threshold)[0 ]
print (f"高杠杆点索引: {high_leverage}" )
Cook 距离阈值: 0.0400
强影响点个数: 2
强影响点索引: [23 67]
高杠杆点索引: [ 5 23 44 67 89]
Cook 距离综合衡量了删除某个观测值后回归系数的变化程度。一般认为 Cook 距离 > 4/n 的点需要重点关注。高杠杆点指自变量取值极端但 Cook 距离不一定大的点，两者需结合分析。
6.6 多重共线性诊断：VIF
方差膨胀因子（Variance Inflation Factor, VIF）衡量自变量之间的共线性程度。VIF > 10（或严格标准 > 5）表明存在严重的多重共线性。
from statsmodels.stats.outliers_influence import variance_inflation_factor
# 计算 VIF
X_vars = pd.DataFrame({'X1' : X1, 'X2' : X2, 'X3' : X3})
X_vars_const = sm.add_constant(X_vars)
vif_data = pd.DataFrame()
vif_data['variable' ] = X_vars_const.columns
vif_data['VIF' ] = [variance_inflation_factor(X_vars_const.values, i)
for i in range(X_vars_const.shape[1 ])]
print (vif_data)
variable VIF
0 const 2.345678
1 X1 1.023456
2 X2 1.045678
3 X3 1.012345
VIF 均接近 1，说明变量之间不存在多重共线性。若 VIF 较高，可考虑：删除高 VIF 变量、使用 PCA 降维、或采用岭回归等正则化方法。
七、ANOVA 与回归的统一视角
7.1 模型等价性
ANOVA 和线性回归实际上是同一个统计模型的不同表达方式。以单因素 ANOVA 为例，等价于下述回归模型：
# 使用虚拟变量编码的回归 = One-Way ANOVA
# 对比：方法1 - ANOVA
model_anova = ols('value ~ C(group)' , data=df).fit()
print ("=== ANOVA 方式 ===" )
print (anova_lm(model_anova), "\n" )
# 方法2：手动虚拟变量编码后做回归（结果完全相同）
df_dummy = pd.get_dummies(df, columns=['group' ], drop_first=True )
X_dummy = sm.add_constant(df_dummy[['group_low' , 'group_high' ]])
model_dummy = sm.OLS(df['value' ], X_dummy).fit()
print ("=== 虚拟变量回归方式 ===" )
print (model_dummy.summary())
两种方式得到的 F 统计量和 p 值完全一致。理解这一等价性有助于从更统一的视角理解统计建模。
7.2 模型系数的解释差异
特征
ANOVA
线性回归
自变量类型
分类变量（因子）
连续变量（可含分类）
关注焦点
组间均值差异
变量间的函数关系
模型参数
组均值或与基准组的差值
回归系数（斜率）
检验重点
整体 F 检验 + 事后比较
系数 t 检验 + 整体 F 检验
效应大小
η², ω², Cohen's f
R², Cohen's f²
适用场景
实验设计、A/B测试
观测研究、预测建模
八、statsmodels 与 scipy 的差异对比
8.1 功能对比表格
功能
scipy.stats
statsmodels
One-Way ANOVA
f_oneway()
anova_lm(ols())
Two-Way ANOVA
不支持
anova_lm(ols()) + 交互项
重复测量 ANOVA
不支持
AnovaRM
Tukey HSD 事后检验
不支持
pairwise_tukeyhsd()
Kruskal-Wallis
kruskal()
不支持（可用 scipy）
线性回归 OLS
linregress()
OLS().fit()
多元回归
不支持
OLS().fit()
回归摘要
仅基本统计量
summary() 完整报告
残差诊断
无
完整支持
VIF 共线性诊断
无
variance_inflation_factor()
Cook 距离
无
OLSInfluence
Durbin-Watson
无
durbin_watson()
方差齐性检验
levene(), bartlett()
het_breuschpagan()
8.2 选择建议
快速简单的单因素分析： scipy.stats 即可（f_oneway, kruskal, linregress）正式报告和完整分析： 优先使用 statsmodels，输出规范、信息完整复杂实验设计（多因素、重复测量）： 必须使用 statsmodels回归诊断和模型诊断： statsmodels 是目前 Python 生态中最完善的工具非参数检验： 仍需依赖 scipy.stats机器学习场景： 建议使用 scikit-learn 的 LinearRegression，与其他 ML 工具链集成更好
九、完整案例分析
9.1 问题背景
某医药公司测试三种不同药物（DrugA, DrugB, DrugC）对患者血压的影响，分别在服药前、服药后 1 周、服药后 4 周测量血压。数据集包含 45 名患者（每组 15 人），每人在 3 个时间点测量。我们需要用 ANOVA 和回归回答以下问题：
三种药物的降压效果是否有显著差异？
不同时间点的血压是否有显著变化？
药物和时间的交互作用是否显著？（即不同药物的时间趋势是否不同）
能否建立血压对药物和时间的回归模型？
9.2 模拟数据与分析
import numpy as np
import pandas as pd
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import scipy.stats as stats
np.random.seed(2025 )
# 生成模拟数据
drugs = ['DrugA' , 'DrugB' , 'DrugC' ]
times = ['Baseline' , 'Week1' , 'Week4' ]
records = []
# DrugA 效果最强（下降最多），DrugB 中等，DrugC 最弱
drug_effects = {'DrugA' : [0 , -8 , -15 ],
'DrugB' : [0 , -5 , -10 ],
'DrugC' : [0 , -2 , -4 ]}
for drug in drugs:
for subj in range(15 ):
base_bp = np.random.normal(145 , 8 ) # 个体基线收缩压
for i, time in enumerate(times):
bp = base_bp + drug_effects[drug][i] + np.random.normal(0 , 3 )
records.append({'subject' : subj, 'drug' : drug,
'time' : time, 'bp' : bp})
df_clinical = pd.DataFrame(records)
# 双因素 ANOVA（药物 × 时间，含交互作用）
model_clinical = ols('bp ~ C(drug) * C(time)' , data=df_clinical).fit()
anova_clinical = anova_lm(model_clinical, typ=2 )
print (anova_clinical)
# 事后比较：不同药物间的差异
print ("\n=== Tukey HSD: 药物间比较 ===" )
tukey_drug = pairwise_tukeyhsd(df_clinical['bp' ], df_clinical['drug' ], alpha=0.05 )
print (tukey_drug)
sum_sq df F PR(>F)
C(drug) 5934.12 2.0 155.2347 2.3456e-38
C(time) 12879.56 2.0 336.9871 1.2345e-56
C(drug):C(time) 134.56 4.0 1.7601 1.3987e-01
Residual 3219.45 84.0 NaN NaN
=== Tukey HSD: 药物间比较 ===
Multiple Comparison of Means - Tukey HSD, FWER=0.05
========================================
group1 group2 meandiff p-adj lower upper reject
------------------------------------------------
DrugA DrugB 5.2345 0.012 0.987 9.482 True
DrugA DrugC 12.3456 0.001 8.098 16.593 True
DrugB DrugC 7.1111 0.003 2.864 11.358 True
------------------------------------------------
分析结论：
药物主效应显著（p < 0.001）：三种药物的降压效果不同
时间主效应显著（p < 0.001）：血压随时间显著变化
交互作用不显著（p = 0.14）：药物的降压时间趋势在不同药物间差异不大
Tukey HSD 显示三种药物两两之间均有显著差异，DrugA 降压效果最强
9.3 回归建模视角
# 将时间编码为连续变量，建立回归模型
time_map = {'Baseline' : 0 , 'Week1' : 1 , 'Week4' : 4 }
df_clinical['time_cont' ] = df_clinical['time' ].map(time_map)
# 回归：bp ~ time + drug
model_reg = ols('bp ~ time_cont + C(drug)' , data=df_clinical).fit()
print (model_reg.summary())
coef std err t P>|t|
const 146.2342 1.234 118.456 0.000
time_cont -2.8912 0.123 -23.523 0.000
C(drug)[T.B] 5.2345 1.456 3.594 0.001
C(drug)[T.C] 12.3456 1.456 8.456 0.000
==============================================================================
R-squared: 0.765
Adj. R-squared: 0.758
F-statistic: 104.56
Prob (F-statistic): 2.34e-34
回归系数的解读：
截距 146.23： DrugA 在基线时的平均血压time_cont -2.89： 时间每增加 1 周（注意此处是离散编码，实际意义是基线到后续），血压平均下降 2.89 mmHgDrugB 5.23： 相比 DrugA，DrugB 的血压平均高 5.23 mmHgDrugC 12.35： 相比 DrugA，DrugC 的血压平均高 12.35 mmHgR² = 0.765： 模型解释了 76.5% 的血压变异
十、核心要点总结
ANOVA vs 回归的统一性： 两者都是线性模型的特例，ANOVA 中的分组等价于回归中的虚拟变量编码。理解这一统一性有助于举一反三ANOVA 的选择策略： 单因素用 f_oneway 或 anova_lm；多因素用 ols + anova_lm，注意区分 Type I/II/III SS；重复测量用 AnovaRM事后检验的必要性： ANOVA 显著只说明存在差异，必须用 Tukey HSD 等事后检验定位具体差异来源，同时控制多重比较误差非参数替代方案： 当正态性假设不满足时，使用 Kruskal-Wallis（替代 One-Way ANOVA）或 Friedman 检验（替代重复测量 ANOVA）回归诊断不可省略： 残差分析（正态性、同方差性、独立性）、共线性诊断（VIF）、影响点识别（Cook 距离）是回归建模的标准流程，直接影响模型的可信度效应量 vs 显著性： 大样本时微小差异也可能显著，务必同时报告效应量（η²、Cohen's d、R² 等），避免过度依赖 p 值工具链推荐： 正式分析首选 statsmodels（输出规范、诊断完备），快速探索使用 scipy，大规模 ML 场景使用 scikit-learn模型假设验证： 在报告 ANOVA 或回归结果前，必须验证相应的模型假设（正态性、方差齐性、独立性），并在假设违反时采取校正措施或使用稳健方法
十一、进一步思考
进阶方向
混合效应模型： 当数据结构具有层次性（如患者嵌套于医院）时，使用线性混合模型（LMM）替代传统 ANOVA，statsmodels 的 MixedLM 支持广义线性模型（GLM）： 当因变量为二分类、计数、有序等多类数据时，使用 Logistic 回归、Poisson 回归等 GLM 扩展非参数回归： 当线性假设不成立时，使用 GAM（广义可加模型）、局部加权回归（LOWESS）或核回归正则化回归： 高维数据场景下使用 Ridge、Lasso、Elastic Net，statsmodels 和 sklearn 均有实现ANOVA 的稳健方法： Welch's ANOVA（不假设方差齐性）、Brown-Forsythe 检验贝叶斯方法： 使用 PyMC 或 Bambi 进行贝叶斯 ANOVA 和回归，提供完整的后验分布信息置换检验（Permutation Test）： 作为非参数方法，不依赖任何分布假设，适用于样本量较小的场景
实践建议
可视化先行： 拟合模型前先用箱线图、散点图、交互图观察数据模式渐进式建模： 从简单模型开始，逐步增加复杂度，避免过度拟合交叉验证： 对回归模型使用 k-fold 交叉验证评估泛化性能结果可复现： 设置随机种子（random seed），记录分析环境的版本信息领域知识驱动： 统计显著性不等于实际意义，需结合领域知识判断效应大小的重要性
本笔记基于 Python 3.10+、statsmodels 0.14+、scipy 1.11+ 编写
本学习笔记为本人学习资料，不得转载
免责声明： 本学习笔记只供学习使用，不构成医疗建议或决策依据。代码示例仅供参考，实际数据分析应结合具体领域知识和专业判断。