如何获取线性判别分析(LDA)、二次判别分析(QDA)模型的p值？

嘿，我完全懂你找LDA和QDA模型p值的痛苦——statsmodels的summary()函数确实不会直接给出这个结果，我之前也踩过这个坑！下面给你两个实用的解决思路：

一、整体模型的显著性p值：似然比检验 #

LDA和QDA都是基于似然估计的模型，我们可以通过似然比检验来判断包含所有特征的模型是否比仅含截距的空模型更显著。这里可以结合statsmodels的判别分析模块来实现：

import statsmodels.api as sm
from statsmodels.discriminant_analysis import LinearDiscriminantAnalysis
from scipy.stats import chi2

# 加载示例数据（替换成你的数据集）
X, y = sm.datasets.get_rdataset("iris", "datasets").data.iloc[:, :-1], sm.datasets.get_rdataset("iris", "datasets").data.iloc[:, -1]
X = sm.add_constant(X)  # 添加截距项

# 拟合完整LDA模型
lda_full = LinearDiscriminantAnalysis()
lda_full.fit(X, y)
log_likelihood_full = lda_full.loglike(X, y)

# 拟合仅含截距的空模型
lda_null = LinearDiscriminantAnalysis()
lda_null.fit(X[:, [0]], y)  # 只保留截距列
log_likelihood_null = lda_null.loglike(X[:, [0]], y)

# 计算似然比统计量和自由度
lr_stat = 2 * (log_likelihood_full - log_likelihood_null)
# 自由度 = 完整模型参数数 - 空模型参数数
df = (lda_full.coef_.size + lda_full.intercept_.size) - (lda_null.coef_.size + lda_null.intercept_.size)

# 计算p值
p_value = chi2.sf(lr_stat, df)
print(f"整体模型的似然比检验p值: {p_value:.4f}")

如果是QDA，只需要把LinearDiscriminantAnalysis换成QuadraticDiscriminantAnalysis即可（statsmodels同样支持）。

二、单个特征的显著性p值：置换检验 #

如果需要评估单个特征对模型的贡献，可以用置换检验——这是一种非参数方法，不需要依赖分布假设，步骤是打乱目标特征的取值，对比原模型和置换后模型的系数/性能差异：

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.datasets import load_iris
import numpy as np

# 加载数据（替换成你的数据集）
X, y = load_iris(return_X_y=True)
target_feature_idx = 0  # 要检验的特征索引

# 拟合原模型，记录目标特征的系数
lda_original = LinearDiscriminantAnalysis()
lda_original.fit(X, y)
original_coef = lda_original.coef_[0, target_feature_idx]  # 取第一个类别的特征系数

# 执行置换检验
n_permutations = 1000
permuted_coefs = []
for _ in range(n_permutations):
    X_permuted = X.copy()
    np.random.shuffle(X_permuted[:, target_feature_idx])  # 打乱目标特征的取值
    lda_perm = LinearDiscriminantAnalysis()
    lda_perm.fit(X_permuted, y)
    permuted_coefs.append(lda_perm.coef_[0, target_feature_idx])

# 计算双尾检验的p值
p_value = (np.sum(np.abs(permuted_coefs) >= np.abs(original_coef)) + 1) / (n_permutations + 1)
print(f"特征{target_feature_idx}的置换检验p值: {p_value:.4f}")

这个方法同样适用于QDA，只需要替换模型类即可。

补充一句：之所以statsmodels的summary不给LDA/QDA的p值，是因为这类判别模型的显著性检验没有像线性/逻辑回归那样形成标准化的输出，所以需要我们手动用上述方法计算。

内容的提问来源于stack exchange，提问作者nimi1234