1. 项目概述高维数据下的偏均值独立性检验在回归分析的实际应用中我们常常面临一个核心问题在已经知道一部分协变量比如通过历史分析或领域知识确定的基因与响应变量比如某种疾病表型相关后如何科学地判断另一组高维协变量比如成千上万个其他基因探针是否仍然对预测响应变量有显著贡献这个问题在基因表达数据分析、金融风险建模、消费者行为研究等高维场景下尤为突出。传统的F检验或基于线性模型的显著性检验严重依赖于模型形式的正确设定一旦真实关系是非线性的结论就可能产生严重偏差。更一般地这个问题可以抽象为一个统计假设检验在控制了协变量子集Z的非线性效应后检验响应变量Y与另一组协变量W是否条件均值独立。其零假设H0为E(Y|X) E(Y|Z)其中X (Z, W)。如果零假设成立意味着在知道了Z之后W不再提供关于Y均值的额外信息可以从模型中剔除从而实现降维和模型简化。然而在高维p2很大且模型形式自由非参数的设定下直接估计和比较两个条件期望E(Y|X)和E(Y|Z)极具挑战。经典的非参数方法如核平滑或局部多项式回归会遭遇“维度灾难”——随着维度增加估计精度急剧下降所需样本量呈指数级增长导致检验功效近乎为零。近年来机器学习方法如深度神经网络DNN、梯度提升树XGBoost在预测任务上展现了惊人实力。一个很自然的想法是能否借助机器学习强大的预测能力来估计这些高维条件期望进而构建检验统计量这正是本文工作的起点。但直接将机器学习模型“黑箱”式地嵌入传统检验框架会遇到一个根本性难题由此构造的检验统计量其零分布即当H0成立时的分布往往是退化的Degenerate或难以处理的无法用于计算p值。这好比拥有了一把锋利的剑机器学习预测器却不知道挥剑的规则检验的统计分布。本文提出的“偏均值独立性检验”Partial Mean Independence Test, pMIT方法正是为了解决这一矛盾。其核心创新在于通过一种巧妙的假设重述和数据分割策略将问题转化为一个可以构造出具有标准卡方零分布的检验统计量。同时如果检验拒绝了零假设即W是显著的我们还需要一个度量来量化“W在给定Z后对Y的均值贡献有多大”。为此本文进一步提出了“偏广义相关性度量”Partial Generalized Measure of Correlation, pGMC作为传统偏相关系数在非线性、模型自由设定下的自然推广。下面我将结合自己多年在数据科学和统计推断交叉领域的实践经验深入拆解pMIT和pGMC的原理、实现细节、避坑技巧以及实际应用中的考量。2. 核心原理与假设重述化繁为简的关键一步2.1 从比较两个条件期望到一个条件期望与一个无条件期望检验的原假设是 H0: E(Y|X) E(Y|Z)。最直接的想法是构造一个统计量来衡量这两者之间的差异例如估计 E[{E(Y|X) - E(Y|Z)}^2]。然而直接操作会遇到两个“拦路虎”第一两个条件期望都需要估计在高维下误差会累积放大第二更重要的是在H0下这个差异的期望为零导致由此构造的样本统计量即使使用机器学习估计量的极限分布是退化的即方差为零无法用于推断。本文一个关键的洞察在于对原假设进行等价重述。定义一个新的变量 Ȳ Y - E(Y|Z)。那么经过简单的推导可以发现原假设 H0: E(Y|X) E(Y|Z) 等价于新的假设 H0‘: E(Ȳ|X) E(Ȳ)。注意由于E(Ȳ) E[Y - E(Y|Z)] 0所以H0‘实际上是在说 E(Ȳ|X) 0。注意这一步重述是方法的核心。它将比较两个复杂的条件期望E(Y|X) vs E(Y|Z)转化为了比较一个条件期望E(Ȳ|X)和一个已知的常数0。这极大地简化了问题结构。为什么这个转化如此有力因为现在我们可以构造一个更“友好”的统计量。我们不再需要直接估计差异的平方期望而是可以基于样本去评估 E(Ȳ|X) 是否显著地偏离0。具体来说我们可以用机器学习方法基于一部分数据去估计函数 g(X) E(Ȳ|X)然后在另一部分数据上计算这个估计量的平方和并与Ȳ的样本均值理论上应为0的波动进行比较。这种构造方式巧妙地引入了额外的随机性来自Ȳ的样本波动使得最终的检验统计量在H0下收敛到一个标准的卡方分布而非退化分布。2.2 数据分割解耦估计与推断机器学习模型尤其是复杂的深度模型其估计过程本身具有不确定性且估计量通常是有偏的。如果使用同一份数据既做模型训练又做假设检验会导致严重的“过拟合”问题使得检验统计量的分布严重扭曲第一类错误假阳性完全失控。这就好比用同一套试题来训练学生并最终评价他们无法判断其真实水平。因此数据分割是应用机器学习进行统计推断的黄金准则。pMIT方法采用随机分割将总样本量为N的数据集D划分为两个独立的部分一个较大的训练集D1样本量n1 n^γ, γ1和一个较小的测试集D2样本量n2 n。这里γ1意味着训练集比测试集大这是一个关键设计。训练集D1的作用用于训练两个机器学习模型。首先用(Z, Y)数据训练模型估计 h(Z) E(Y|Z)得到估计量 ĥ_D1(Z)。然后计算训练集上的残差 Ȳ Y - ĥ_D1(Z)。接着用(X, Ȳ)数据训练第二个模型估计 g(X) E(Ȳ|X)得到估计量 ĝ_D1(X)。实际上这一步等价于用(X, Y)直接训练一个估计m(X)E(Y|X)的模型然后令 ĝ_D1(X) m̂_D1(X) - ĥ_D1(Z)。但在实操中直接回归残差Ȳ到X上通常更稳定。测试集D2的作用用于构建检验统计量评估假设。我们绝对不使用D2中的任何信息来重新训练或调整模型确保评估的独立性。这种“训练-测试”分离保证了在H0下即使机器学习估计量 ĝ_D1(X) 存在误差只要这个误差的阶数足够小通过使用较大的训练集实现它在测试集上的平方和相对于Ȳ的样本波动就是可以忽略的。从而检验统计量的极限行为由Ȳ的样本波动主导进而推导出标准的分布。2.3 检验统计量的构造与直觉基于训练好的估计量 ĥ_D1 和 ĝ_D1我们在独立的测试集D2上构造如下统计量Tn (1/n) * Σ_{i∈D2} [ ĝ_D1(Xi) - (1/n) Σ_{j∈D2} {Yj - ĥ_D1(Zj)} ]^2这个公式的直觉非常清晰ĝ_D1(Xi)这是我们基于X对“调整后响应变量”Ȳ的最佳预测在训练集上学到的。(1/n) Σ_{j∈D2} {Yj - ĥ_D1(Zj)}这是测试集上Ȳ的样本均值。在H0下E(Ȳ)0所以这个样本均值应该在0附近波动。整个统计量Tn衡量的是我们的预测值 ĝ_D1(X) 与Ȳ的样本均值之间的差异的平方平均。如果H0成立即E(Ȳ|X)0那么ĝ_D1(X)应该只是在估计0因此它与Ȳ的样本均值也在0附近的差异应该很小Tn的值也会很小。如果H1成立那么E(Ȳ|X) ≠ 0我们的预测模型ĝ_D1(X)会捕捉到这种信号导致其值系统地偏离Ȳ的样本均值从而使Tn变大。为了进行标准化我们最终使用的检验统计量是 Vn n * Tn / σ̂²_{Y|Z}其中 σ̂²_{Y|Z} 是测试集上残差 {Y - ĥ_D1(Z)} 的样本方差。理论证明在一定的正则条件下主要关于机器学习估计量的收敛速度当H0成立时Vn 依分布收敛于自由度为1的卡方分布 χ²(1)。这就为我们提供了计算p值的直接依据。3. 实操流程与核心环节实现3.1 算法步骤详解下面我将结合代码片段以Python为例和具体操作详细拆解pMIT的单次数据分割算法流程。假设我们已有数据Y(n_samples,)X(n_samples, p) 其中X可以划分为Z(n_samples, p1) 和W(n_samples, p2)。import numpy as np from sklearn.model_selection import train_test_split from sklearn.neural_network import MLPRegressor from xgboost import XGBRegressor from scipy.stats import chi2 def pMIT_single_split(Y, X, Z, W, test_size0.2, estimator_typeDNN, random_state42): 执行一次数据分割的pMIT检验。 参数: Y: 响应变量数组 X: 全部协变量数组 (包含Z和W) Z: 已知/控制的协变量数组 W: 待检验的协变量数组 test_size: 测试集比例 (D2的比例) estimator_type: 使用的机器学习模型DNN 或 XGBoost random_state: 随机种子 返回: p_value: 检验的p值 test_statistic: 检验统计量Vn # 1. 数据分割D1 (训练集) 和 D2 (测试集) # 总样本量 N N len(Y) # 确保训练集更大这里设置训练集比例为 1 - test_size且 test_size 应较小如0.2 idx_train, idx_test train_test_split(np.arange(N), test_sizetest_size, random_staterandom_state) Y_train, Y_test Y[idx_train], Y[idx_test] Z_train, Z_test Z[idx_train], Z[idx_test] X_train, X_test X[idx_train], X[idx_test] # W_train, W_test W[idx_train], W[idx_test] # 在单次检验中W的划分是自动的因为X包含W # 2. 在训练集D1上估计 h(Z) E(Y|Z) if estimator_type DNN: # 使用深度神经网络需要标准化数据 from sklearn.preprocessing import StandardScaler scaler_Z StandardScaler() Z_train_scaled scaler_Z.fit_transform(Z_train) Z_test_scaled scaler_Z.transform(Z_test) model_h MLPRegressor(hidden_layer_sizes(100, 50), activationrelu, solveradam, max_iter1000, random_staterandom_state) model_h.fit(Z_train_scaled, Y_train) h_pred_train model_h.predict(Z_train_scaled) h_pred_test model_h.predict(Z_test_scaled) elif estimator_type XGBoost: model_h XGBRegressor(n_estimators100, max_depth6, learning_rate0.1, random_staterandom_state) model_h.fit(Z_train, Y_train) h_pred_train model_h.predict(Z_train) h_pred_test model_h.predict(Z_test) else: raise ValueError(estimator_type 必须是 DNN 或 XGBoost) # 3. 计算训练集上的残差 \tilde{Y} Y - \hat{h}(Z) Y_tilde_train Y_train - h_pred_train # 4. 在训练集D1上估计 g(X) E(\tilde{Y}|X) if estimator_type DNN: scaler_X StandardScaler() X_train_scaled scaler_X.fit_transform(X_train) X_test_scaled scaler_X.transform(X_test) model_g MLPRegressor(hidden_layer_sizes(100, 50), activationrelu, solveradam, max_iter1000, random_staterandom_state) model_g.fit(X_train_scaled, Y_tilde_train) g_pred_test model_g.predict(X_test_scaled) # 只在测试集上预测 elif estimator_type XGBoost: model_g XGBRegressor(n_estimators100, max_depth6, learning_rate0.1, random_staterandom_state) model_g.fit(X_train, Y_tilde_train) g_pred_test model_g.predict(X_test) # 只在测试集上预测 # 5. 在测试集D2上计算检验统计量 n_test len(Y_test) # 计算 \tilde{Y} 在测试集上的值 Y_tilde_test Y_test - h_pred_test # 计算 \tilde{Y} 在测试集上的样本均值 Y_tilde_test_mean np.mean(Y_tilde_test) # 计算 T_n T_n np.mean((g_pred_test - Y_tilde_test_mean) ** 2) # 计算 \hat{\sigma}^2_{Y|Z} sigma2_hat np.mean((Y_tilde_test) ** 2) # 注意这里用的是中心化的吗原公式是样本方差但理论中E(\tilde{Y})0所以样本二阶矩即可。实践中常用样本方差。 # 更稳健的做法使用样本方差 sigma2_hat np.var(Y_tilde_test, ddof0) # 使用0自由度的方差与理论推导对应 # 计算 V_n V_n n_test * T_n / sigma2_hat # 6. 计算p值 (基于卡方(1)分布) p_value 1 - chi2.cdf(V_n, df1) return p_value, V_n3.2 模型选择与超参数调优机器学习模型的选择和调优对检验的性能至关重要。DNN和XGBoost是两种强大的选择但各有适用场景。深度神经网络适用于特征间可能存在复杂交互和非线性关系的情况。它不依赖于特征工程能自动学习表征。但其训练需要更多数据对超参数层数、神经元数、学习率敏感且训练过程不稳定。实操心得对于高维小样本数据DNN容易过拟合。务必使用早停法、Dropout、权重衰减等正则化技术。可以先在训练集D1上进一步划分一个验证集用于确定最佳网络结构和训练轮数。一个常见的坑是在测试集D2上做任何形式的模型选择或调参这会严重破坏检验的有效性。所有调优必须严格在训练集D1内部完成。XGBoost基于树的集成方法通常对异构数据混合类型特征、缺失值有较好的鲁棒性且训练速度较快超参数相对容易调节。它通过梯度提升来优化能捕捉复杂的非线性模式。实操心得关键超参数包括max_depth控制树复杂度防止过拟合、n_estimators树的数量太多会过拟合、learning_rate学习率小学习率配合多树通常更好和subsample、colsample_bytree行/列采样引入随机性防止过拟合。同样使用D1内部的交叉验证来调优。重要提示无论选择哪种模型目标都是获得对条件期望h(Z)和g(X)的“最佳可能”估计。这里的“最佳”并非指在独立测试集上的最小预测误差而是指估计量的均方误差MSE收敛速度足够快以满足理论条件(C1)E[(ĝ - g)^2] O(n1^{-a})其中a0n1是训练集大小。实践中我们通过交叉验证选择在D1上预测性能最好的模型配置这通常能保证一个不错的收敛速率。3.3 功效增强与算法稳定性策略单次数据分割的结果可能因随机分割而波动。为了提高检验的稳定性和功效论文提出了两种改进策略。1. 功效增强统计量在计算得到Vn后可以构造一个增强版统计量Vn* Vn τ * Σ_{i∈D2} [ĝ_D1(Xi)]^2其中τ 0是一个调节参数。原理在H0下ĝ_D1(X)估计的是0所以其平方和很小Vn* 与 Vn 渐近等价不影响第一类错误。在H1下ĝ_D1(X)捕捉到了信号其平方和是正的因此 Vn* Vn从而提升了检验功效。实操建议论文建议简单取τ1。在实际代码中可以在计算p值前用Vn* 替代 Vn。但需注意这略微改变了零分布理论上仍是卡方但需要更严格的证明在极端追求第一类错误精确控制时需谨慎。对于大多数应用τ1是一个安全有效的选择。2. 多重数据分割与p值聚合为了消除单次随机分割的偶然性可以进行B次例如B100独立的数据分割每次得到一個p值 p_b (b1,...,B)然后聚合这些p值得到一个更稳健的最终p值。常用聚合方法中位数法Q* 2 * median({p_b})。如果Q* α则拒绝H0。这种方法简单对异常值不敏感。最小p值法经调整Q* min({p_b}) * B。这是Bonferroni校正非常保守。基于均匀分布的方法如果p值在原假设下服从均匀分布可以计算 -2 * Σ log(p_b)它服从卡方分布(2B)。但数据分割导致p值间并非独立此方法可能过于激进。实操建议我推荐使用中位数法。它在控制第一类错误和保持功效之间取得了较好的平衡且计算简单。在代码实现上只需将上述pMIT_single_split函数运行B次收集p值然后计算中位数的两倍即可。def pMIT_multiple_splits(Y, X, Z, W, B100, test_size0.2, estimator_typeDNN, random_seed_start0): 基于多重数据分割的pMIT检验。 p_values [] for b in range(B): p_val, _ pMIT_single_split(Y, X, Z, W, test_sizetest_size, estimator_typeestimator_type, random_staterandom_seed_start b) p_values.append(p_val) # 使用中位数法聚合p值 aggregated_p 2 * np.median(p_values) # 注意聚合p值可能大于1需要截断 aggregated_p min(aggregated_p, 1.0) return aggregated_p, p_values3.4 分割比例ξ的自适应选择训练集大小n1或分割比例ξ n1/N的选择是一个偏差-方差权衡问题。n1越大条件期望估计得越准偏差小但用于构建统计量的测试集n2越小统计量的方差越大功效可能下降。反之亦然。论文提出了一种基于置换检验的数据自适应选择方法设定一个候选比例集合Ξ例如Ξ {0.5, 0.6, 0.7, 0.8, 0.9}。对于每个候选ξ在原数据上随机置换打乱待检验的协变量W破坏其与Y的关系从而模拟H0成立的数据。对每个置换后的数据集应用pMIT算法单次或多次分割计算其p值。重复M次如M200。对于每个ξ计算其经验第一类错误率即p值小于显著性水平α如0.05的比例。选择能满足经验第一类错误率 ≤ α 的最小ξ作为最终的分割比例ξ*。这个方法的逻辑是在H0成立的数据上我们希望检验能控制第一类错误。通过选择能满足这一要求的最小训练集比例我们尽可能为测试集保留了更多的样本以期在H1下获得更高的功效。注意事项这个自适应过程计算量较大因为需要对每个候选ξ进行M次置换检验。在实际应用中如果样本量足够大一个经验性的做法是固定一个较大的训练集比例例如ξ0.8或0.9以确保条件期望的估计足够准确。样本量较小时自适应选择更有价值。4. 偏广义相关性度量从检验到度量当pMIT检验拒绝了H0我们自然想知道“W在给定Z后对Y的均值解释力到底有多强”传统的偏相关系数只能度量线性关联而pGMC提供了一个模型自由的、介于0到1之间的度量。4.1 pGMC的定义与估计pGMC的定义清晰且直观 r²(Y, W|Z) E[{m(X) - h(Z)}²] / E[{Y - h(Z)}²] 其中m(X) E(Y|X), h(Z) E(Y|Z)。分子E[{m(X) - h(Z)}²] 衡量的是当我们从只知道Z进步到知道完整的X时对Y的条件均值预测的改进量均方误差意义下。分母E[{Y - h(Z)}²] 衡量的是如果只用Z来预测Y其不可解释的变异误差方差。因此pGMC解释为W所解释的那部分变异占Z未解释的总变异的比例。它是一个决定系数R²在模型自由和偏相关语境下的推广。基于数据分割和机器学习估计量一个自然的估计量是 r̂² [ Σ_{i∈D2} {m̂_D1(Xi) - ĥ_D1(Zi)}² ] / [ Σ_{i∈D2} {Yi - ĥ_D1(Zi)}² ] 这里D2是独立于模型训练集D1的测试集以确保估计的无偏性或渐近无偏性。4.2 pGMC的性质与置信区间pGMC拥有许多优良性质使其成为一个有吸引力的度量工具取值范围0 ≤ r² ≤ 1。零点意义r² 0当且仅当H0成立即E(Y|X)E(Y|Z)。这意味着给定Z后W不提供任何关于Y均值的额外信息。单位点意义r² 1当且仅当Y m(X) 几乎必然成立即给定Z后Y完全是W通过X的确定性函数。与偏相关系数的关系r² ≥ ρ²。pGMC总是大于等于线性偏相关系数的平方它能捕捉非线性依赖。与GMC的关系r² [GMC(Y|X) - GMC(Y|Z)] / [1 - GMC(Y|Z)]。这提供了另一种计算和理解方式。理论证明在适当的条件下pGMC的估计量 r̂² 是渐近正态的并且收敛速度为根号n测试集大小。这意味着我们可以为其构造渐近置信区间。置信区间构造Wald型 令 θ r²(Y, W|Z)。其估计量为 θ̂ r̂²。 我们需要估计 θ̂ 的渐近方差 σ̂²_θ。这可以通过刀切法或自助法在测试集D2上实现。由于D2独立于训练过程自助法是一个简单可靠的选择。import numpy as np from sklearn.utils import resample def estimate_pGMC_with_CI(Y, X, Z, W, test_size0.2, estimator_typeDNN, n_bootstrap1000, alpha0.05, random_state42): 估计pGMC并计算其自助置信区间。 # 1. 首先进行一次数据分割得到训练好的模型和测试集索引 # 这里简化假设已有训练好的模型 h_estimator, m_estimator 和测试集数据 # 在实际中需要先运行类似 pMIT_single_split 的前半部分来训练模型并得到测试集预测 # 假设我们已经有了 # idx_test: 测试集索引 # h_pred_test: 测试集上 h(Z) 的预测 # m_pred_test: 测试集上 m(X) 的预测 # Y_test: 测试集响应变量 # 计算点估计 numerator np.mean((m_pred_test - h_pred_test) ** 2) denominator np.mean((Y_test - h_pred_test) ** 2) theta_hat numerator / denominator # 2. 在测试集上进行自助法抽样估计标准误 n_test len(Y_test) bootstrap_estimates [] for b in range(n_bootstrap): # 在测试集索引范围内有放回抽样 indices resample(np.arange(n_test), random_staterandom_state b) Y_b Y_test[indices] h_pred_b h_pred_test[indices] m_pred_b m_pred_test[indices] num_b np.mean((m_pred_b - h_pred_b) ** 2) denom_b np.mean((Y_b - h_pred_b) ** 2) theta_b num_b / denom_b bootstrap_estimates.append(theta_b) # 计算自助法标准误和置信区间 se_theta np.std(bootstrap_estimates, ddof1) # 正态近似置信区间 z_alpha 1.96 # for 95% CI ci_lower theta_hat - z_alpha * se_theta ci_upper theta_hat z_alpha * se_theta # 确保置信区间在[0,1]内 ci_lower max(0, ci_lower) ci_upper min(1, ci_upper) # 也可以使用百分位数区间 ci_lower_perc np.percentile(bootstrap_estimates, (alpha/2)*100) ci_upper_perc np.percentile(bootstrap_estimates, (1-alpha/2)*100) return theta_hat, (ci_lower, ci_upper), (ci_lower_perc, ci_upper_perc), se_theta注意当真实pGMC值接近0或1时正态近似可能不佳。百分位数自助区间或偏差校正的自助区间BCa通常更稳健。此外确保分母不为零即Y不能被Z完美预测是定义和估计pGMC的前提。5. 常见问题、陷阱与实战技巧5.1 模型过拟合与欠拟合问题训练集D1上模型过拟合导致 ĥ(Z) 和 m̂(X) 在训练集上表现极好但在独立的测试集D2上预测误差很大。这会严重扭曲检验统计量Tn的分布。欠拟合则会导致估计偏差过大降低检验功效。诊断在D1内部进行交叉验证监控训练误差和验证误差的差距。如果训练误差远小于验证误差可能是过拟合。如果两者都很大可能是欠拟合或模型容量不足。解决对于DNN使用早停法、Dropout、L2正则化、减少网络层数或神经元数量。对于XGBoost降低max_depth、增加min_child_weight、使用更小的learning_rate并配合更多的n_estimators、启用行/列采样。通用策略在D1上使用交叉验证或保留一个验证集来调优超参数。绝对禁止使用测试集D2进行任何模型选择或调参。5.2 样本量不足与分割比例困境问题在高维小样本如n100, p1000场景下数据分割后训练集D1和测试集D2的样本量都可能不足。D1太小导致模型估计不准D2太小导致检验统计量方差大、功效低。应对优先保证训练集规模理论要求γ1/a其中a是估计误差的收敛速率。对于复杂的非线性关系a可能较小如0.2-0.5这意味着需要很大的训练集。在实践中如果总样本量N有限应倾向于给D1分配更大的比例如ξ0.8或0.9牺牲一些测试集样本以换取更准确的估计。考虑使用更简单的模型在样本量很小时过于复杂的模型如很深的DNN弊大于利。可以尝试线性模型带L1/L2正则化、浅层神经网络或树深度受限的XGBoost。利用功效增强启用Vn*统计量τ1可以在不增加样本的情况下提升一些功效。谨慎解读结果当样本量很小时即使检验不显著p值大也不能武断地认为W不重要可能是功效不足。应同时报告pGMC的估计值及其置信区间以量化效应大小及其不确定性。5.3 协变量维度极高时的挑战问题当p尤其是p2即W的维度极高时即使使用DNN或XGBoost估计m(X) E(Y|X)也变得异常困难容易陷入维度灾难。策略特征预筛选在正式应用pMIT之前可以对W进行一轮基于边际关联或简单模型的预筛选保留最相关的一部分特征如top 100或top 500。但这会引入选择偏差需要谨慎。使用专门的高维方法考虑在估计m(X)时使用适用于高维数据的机器学习方法例如带L1正则化的线性模型如果相信关系近似线性或可加。稀疏深度学习使用Dropout、稀疏连接或L1正则化鼓励稀疏性。随机森林/XGBoost它们本身对高维数据有一定鲁棒性通过特征重要性可以进行隐式特征选择。聚焦于Z已知的场景pMIT的优势在于当Z包含已知的重要变量时问题转化为估计E(Ȳ|X)其中Ȳ Y - E(Y|Z)。如果h(Z)能被很好地估计那么Ȳ的变异可能主要由W中少数重要变量驱动从而缓解了高维问题。5.4 计算效率与实现优化问题多重数据分割、自助法置信区间、自适应选择ξ等操作计算开销大。优化建议并行化B次数据分割的pMIT计算、M次置换检验、n_bootstrap次自助抽样这些任务都是完全独立的非常适合并行计算如使用Python的joblib或multiprocessing库。缓存模型在自适应选择ξ或进行多重分割时对于相同的训练集比例ξ和相同的随机种子训练出的模型是相同的。可以设计缓存机制避免重复训练。使用高效库对于XGBoost确保使用GPU加速如果可用。对于DNN使用TensorFlow或PyTorch等框架并利用GPU。设定合理的迭代次数对于探索性分析B50或100可能就够了。对于最终报告B200或500更稳妥。自助法次数n_bootstrap1000通常足够。5.5 与替代方法的比较与选择pMIT不是唯一解决该问题的方法。了解其与相关方法的区别有助于正确选择。方法核心思想优点缺点/注意事项pMIT (本文)重述假设为E(Ȳ|X)0数据分割ML估计构造卡方统计量零分布为标准卡方计算快理论性质清晰提供了pGMC度量需要数据分割损失部分样本信息对ML估计精度有要求Dai et al. (2024)比较两个预测误差的差异引入随机扰动处理退化分布可处理更一般的条件矩限制需要选择扰动大小理论局部备择检出速率可能稍慢Williamson et al. (2023)双数据分割分别估计两个非零参数概念相对直观需要至少两个独立分割样本利用可能更低Lundborg et al. (2022) PCM投影协方差度量基于残差正交化在线性模型下功效可能很高在模型自由设定下需要三个以上子样本计算多个条件期望计算量大传统F检验基于线性模型残差平方和比较计算极其简单在线性情况下最优严重依赖线性模型假设误设时失效选择建议如果追求计算简便和解释直观且相信关系近似线性可以先尝试传统F检验或线性模型的偏F检验。如果处理的是明确的高维线性模型Lundborg et al. (2022)的PCM方法可能是一个强有力的竞争者。如果面对的是高维、非线性、模型自由的设定并且希望有一个统一、有理论保障、附带效应大小度量的框架pMITpGMC的组合是一个非常全面和实用的选择。它的标准卡方推断避免了重抽样计算效率高pGMC提供了有意义的解释。6. 总结与扩展思考偏均值独立性检验及其关联的pGMC度量为我们分析高维数据中变量间的条件依赖关系提供了一套强大的工具。它巧妙地将机器学习的预测能力与严谨的统计推断相结合通过假设重述和数据分割解决了推断中的根本难题。在实际应用中我个人的体会是理解问题的本质比机械地套用算法更重要。pMIT检验回答的是“W是否在给定Z后对Y的均值有贡献”这是一个关于增量贡献的问题。因此如何定义Z已知的、需要控制的变量集至关重要。Z应该包含所有你确信与Y相关且可能与W混淆的变量。例如在基因研究中Z可能是已知的致病通路基因W是待探索的其他基因。此外pMIT检验的结论依赖于机器学习模型对条件期望估计的“足够好”。虽然理论只要求一定的收敛速率但在有限样本下一个糟糕的估计会导致检验功效低下。因此投入精力进行谨慎的模型选择和验证是值得的。不要仅仅把它当作一个黑箱流程。最后pGMC作为一个效应大小度量其价值不亚于检验本身。一个显著的p值告诉我们关联存在而pGMC的估计值及其置信区间则告诉我们这个关联的强度这在许多科学领域对于评估实际意义至关重要。报告结果时应同时提供p值和pGMC估计及置信区间以给出更完整的画面。这个方法可以自然地扩展到其他相关问题例如检验条件独立性通过检验一系列条件均值或者处理分类响应变量通过适当修改损失函数和模型。其核心思想——利用机器学习估计复杂函数再通过巧妙的统计构造进行推断——为高维统计推断开辟了一条富有前景的道路。
高维非线性数据下的偏均值独立性检验：原理、实现与应用
1. 项目概述高维数据下的偏均值独立性检验在回归分析的实际应用中我们常常面临一个核心问题在已经知道一部分协变量比如通过历史分析或领域知识确定的基因与响应变量比如某种疾病表型相关后如何科学地判断另一组高维协变量比如成千上万个其他基因探针是否仍然对预测响应变量有显著贡献这个问题在基因表达数据分析、金融风险建模、消费者行为研究等高维场景下尤为突出。传统的F检验或基于线性模型的显著性检验严重依赖于模型形式的正确设定一旦真实关系是非线性的结论就可能产生严重偏差。更一般地这个问题可以抽象为一个统计假设检验在控制了协变量子集Z的非线性效应后检验响应变量Y与另一组协变量W是否条件均值独立。其零假设H0为E(Y|X) E(Y|Z)其中X (Z, W)。如果零假设成立意味着在知道了Z之后W不再提供关于Y均值的额外信息可以从模型中剔除从而实现降维和模型简化。然而在高维p2很大且模型形式自由非参数的设定下直接估计和比较两个条件期望E(Y|X)和E(Y|Z)极具挑战。经典的非参数方法如核平滑或局部多项式回归会遭遇“维度灾难”——随着维度增加估计精度急剧下降所需样本量呈指数级增长导致检验功效近乎为零。近年来机器学习方法如深度神经网络DNN、梯度提升树XGBoost在预测任务上展现了惊人实力。一个很自然的想法是能否借助机器学习强大的预测能力来估计这些高维条件期望进而构建检验统计量这正是本文工作的起点。但直接将机器学习模型“黑箱”式地嵌入传统检验框架会遇到一个根本性难题由此构造的检验统计量其零分布即当H0成立时的分布往往是退化的Degenerate或难以处理的无法用于计算p值。这好比拥有了一把锋利的剑机器学习预测器却不知道挥剑的规则检验的统计分布。本文提出的“偏均值独立性检验”Partial Mean Independence Test, pMIT方法正是为了解决这一矛盾。其核心创新在于通过一种巧妙的假设重述和数据分割策略将问题转化为一个可以构造出具有标准卡方零分布的检验统计量。同时如果检验拒绝了零假设即W是显著的我们还需要一个度量来量化“W在给定Z后对Y的均值贡献有多大”。为此本文进一步提出了“偏广义相关性度量”Partial Generalized Measure of Correlation, pGMC作为传统偏相关系数在非线性、模型自由设定下的自然推广。下面我将结合自己多年在数据科学和统计推断交叉领域的实践经验深入拆解pMIT和pGMC的原理、实现细节、避坑技巧以及实际应用中的考量。2. 核心原理与假设重述化繁为简的关键一步2.1 从比较两个条件期望到一个条件期望与一个无条件期望检验的原假设是 H0: E(Y|X) E(Y|Z)。最直接的想法是构造一个统计量来衡量这两者之间的差异例如估计 E[{E(Y|X) - E(Y|Z)}^2]。然而直接操作会遇到两个“拦路虎”第一两个条件期望都需要估计在高维下误差会累积放大第二更重要的是在H0下这个差异的期望为零导致由此构造的样本统计量即使使用机器学习估计量的极限分布是退化的即方差为零无法用于推断。本文一个关键的洞察在于对原假设进行等价重述。定义一个新的变量 Ȳ Y - E(Y|Z)。那么经过简单的推导可以发现原假设 H0: E(Y|X) E(Y|Z) 等价于新的假设 H0‘: E(Ȳ|X) E(Ȳ)。注意由于E(Ȳ) E[Y - E(Y|Z)] 0所以H0‘实际上是在说 E(Ȳ|X) 0。注意这一步重述是方法的核心。它将比较两个复杂的条件期望E(Y|X) vs E(Y|Z)转化为了比较一个条件期望E(Ȳ|X)和一个已知的常数0。这极大地简化了问题结构。为什么这个转化如此有力因为现在我们可以构造一个更“友好”的统计量。我们不再需要直接估计差异的平方期望而是可以基于样本去评估 E(Ȳ|X) 是否显著地偏离0。具体来说我们可以用机器学习方法基于一部分数据去估计函数 g(X) E(Ȳ|X)然后在另一部分数据上计算这个估计量的平方和并与Ȳ的样本均值理论上应为0的波动进行比较。这种构造方式巧妙地引入了额外的随机性来自Ȳ的样本波动使得最终的检验统计量在H0下收敛到一个标准的卡方分布而非退化分布。2.2 数据分割解耦估计与推断机器学习模型尤其是复杂的深度模型其估计过程本身具有不确定性且估计量通常是有偏的。如果使用同一份数据既做模型训练又做假设检验会导致严重的“过拟合”问题使得检验统计量的分布严重扭曲第一类错误假阳性完全失控。这就好比用同一套试题来训练学生并最终评价他们无法判断其真实水平。因此数据分割是应用机器学习进行统计推断的黄金准则。pMIT方法采用随机分割将总样本量为N的数据集D划分为两个独立的部分一个较大的训练集D1样本量n1 n^γ, γ1和一个较小的测试集D2样本量n2 n。这里γ1意味着训练集比测试集大这是一个关键设计。训练集D1的作用用于训练两个机器学习模型。首先用(Z, Y)数据训练模型估计 h(Z) E(Y|Z)得到估计量 ĥ_D1(Z)。然后计算训练集上的残差 Ȳ Y - ĥ_D1(Z)。接着用(X, Ȳ)数据训练第二个模型估计 g(X) E(Ȳ|X)得到估计量 ĝ_D1(X)。实际上这一步等价于用(X, Y)直接训练一个估计m(X)E(Y|X)的模型然后令 ĝ_D1(X) m̂_D1(X) - ĥ_D1(Z)。但在实操中直接回归残差Ȳ到X上通常更稳定。测试集D2的作用用于构建检验统计量评估假设。我们绝对不使用D2中的任何信息来重新训练或调整模型确保评估的独立性。这种“训练-测试”分离保证了在H0下即使机器学习估计量 ĝ_D1(X) 存在误差只要这个误差的阶数足够小通过使用较大的训练集实现它在测试集上的平方和相对于Ȳ的样本波动就是可以忽略的。从而检验统计量的极限行为由Ȳ的样本波动主导进而推导出标准的分布。2.3 检验统计量的构造与直觉基于训练好的估计量 ĥ_D1 和 ĝ_D1我们在独立的测试集D2上构造如下统计量Tn (1/n) * Σ_{i∈D2} [ ĝ_D1(Xi) - (1/n) Σ_{j∈D2} {Yj - ĥ_D1(Zj)} ]^2这个公式的直觉非常清晰ĝ_D1(Xi)这是我们基于X对“调整后响应变量”Ȳ的最佳预测在训练集上学到的。(1/n) Σ_{j∈D2} {Yj - ĥ_D1(Zj)}这是测试集上Ȳ的样本均值。在H0下E(Ȳ)0所以这个样本均值应该在0附近波动。整个统计量Tn衡量的是我们的预测值 ĝ_D1(X) 与Ȳ的样本均值之间的差异的平方平均。如果H0成立即E(Ȳ|X)0那么ĝ_D1(X)应该只是在估计0因此它与Ȳ的样本均值也在0附近的差异应该很小Tn的值也会很小。如果H1成立那么E(Ȳ|X) ≠ 0我们的预测模型ĝ_D1(X)会捕捉到这种信号导致其值系统地偏离Ȳ的样本均值从而使Tn变大。为了进行标准化我们最终使用的检验统计量是 Vn n * Tn / σ̂²_{Y|Z}其中 σ̂²_{Y|Z} 是测试集上残差 {Y - ĥ_D1(Z)} 的样本方差。理论证明在一定的正则条件下主要关于机器学习估计量的收敛速度当H0成立时Vn 依分布收敛于自由度为1的卡方分布 χ²(1)。这就为我们提供了计算p值的直接依据。3. 实操流程与核心环节实现3.1 算法步骤详解下面我将结合代码片段以Python为例和具体操作详细拆解pMIT的单次数据分割算法流程。假设我们已有数据Y(n_samples,)X(n_samples, p) 其中X可以划分为Z(n_samples, p1) 和W(n_samples, p2)。import numpy as np from sklearn.model_selection import train_test_split from sklearn.neural_network import MLPRegressor from xgboost import XGBRegressor from scipy.stats import chi2 def pMIT_single_split(Y, X, Z, W, test_size0.2, estimator_typeDNN, random_state42): 执行一次数据分割的pMIT检验。 参数: Y: 响应变量数组 X: 全部协变量数组 (包含Z和W) Z: 已知/控制的协变量数组 W: 待检验的协变量数组 test_size: 测试集比例 (D2的比例) estimator_type: 使用的机器学习模型DNN 或 XGBoost random_state: 随机种子 返回: p_value: 检验的p值 test_statistic: 检验统计量Vn # 1. 数据分割D1 (训练集) 和 D2 (测试集) # 总样本量 N N len(Y) # 确保训练集更大这里设置训练集比例为 1 - test_size且 test_size 应较小如0.2 idx_train, idx_test train_test_split(np.arange(N), test_sizetest_size, random_staterandom_state) Y_train, Y_test Y[idx_train], Y[idx_test] Z_train, Z_test Z[idx_train], Z[idx_test] X_train, X_test X[idx_train], X[idx_test] # W_train, W_test W[idx_train], W[idx_test] # 在单次检验中W的划分是自动的因为X包含W # 2. 在训练集D1上估计 h(Z) E(Y|Z) if estimator_type DNN: # 使用深度神经网络需要标准化数据 from sklearn.preprocessing import StandardScaler scaler_Z StandardScaler() Z_train_scaled scaler_Z.fit_transform(Z_train) Z_test_scaled scaler_Z.transform(Z_test) model_h MLPRegressor(hidden_layer_sizes(100, 50), activationrelu, solveradam, max_iter1000, random_staterandom_state) model_h.fit(Z_train_scaled, Y_train) h_pred_train model_h.predict(Z_train_scaled) h_pred_test model_h.predict(Z_test_scaled) elif estimator_type XGBoost: model_h XGBRegressor(n_estimators100, max_depth6, learning_rate0.1, random_staterandom_state) model_h.fit(Z_train, Y_train) h_pred_train model_h.predict(Z_train) h_pred_test model_h.predict(Z_test) else: raise ValueError(estimator_type 必须是 DNN 或 XGBoost) # 3. 计算训练集上的残差 \tilde{Y} Y - \hat{h}(Z) Y_tilde_train Y_train - h_pred_train # 4. 在训练集D1上估计 g(X) E(\tilde{Y}|X) if estimator_type DNN: scaler_X StandardScaler() X_train_scaled scaler_X.fit_transform(X_train) X_test_scaled scaler_X.transform(X_test) model_g MLPRegressor(hidden_layer_sizes(100, 50), activationrelu, solveradam, max_iter1000, random_staterandom_state) model_g.fit(X_train_scaled, Y_tilde_train) g_pred_test model_g.predict(X_test_scaled) # 只在测试集上预测 elif estimator_type XGBoost: model_g XGBRegressor(n_estimators100, max_depth6, learning_rate0.1, random_staterandom_state) model_g.fit(X_train, Y_tilde_train) g_pred_test model_g.predict(X_test) # 只在测试集上预测 # 5. 在测试集D2上计算检验统计量 n_test len(Y_test) # 计算 \tilde{Y} 在测试集上的值 Y_tilde_test Y_test - h_pred_test # 计算 \tilde{Y} 在测试集上的样本均值 Y_tilde_test_mean np.mean(Y_tilde_test) # 计算 T_n T_n np.mean((g_pred_test - Y_tilde_test_mean) ** 2) # 计算 \hat{\sigma}^2_{Y|Z} sigma2_hat np.mean((Y_tilde_test) ** 2) # 注意这里用的是中心化的吗原公式是样本方差但理论中E(\tilde{Y})0所以样本二阶矩即可。实践中常用样本方差。 # 更稳健的做法使用样本方差 sigma2_hat np.var(Y_tilde_test, ddof0) # 使用0自由度的方差与理论推导对应 # 计算 V_n V_n n_test * T_n / sigma2_hat # 6. 计算p值 (基于卡方(1)分布) p_value 1 - chi2.cdf(V_n, df1) return p_value, V_n3.2 模型选择与超参数调优机器学习模型的选择和调优对检验的性能至关重要。DNN和XGBoost是两种强大的选择但各有适用场景。深度神经网络适用于特征间可能存在复杂交互和非线性关系的情况。它不依赖于特征工程能自动学习表征。但其训练需要更多数据对超参数层数、神经元数、学习率敏感且训练过程不稳定。实操心得对于高维小样本数据DNN容易过拟合。务必使用早停法、Dropout、权重衰减等正则化技术。可以先在训练集D1上进一步划分一个验证集用于确定最佳网络结构和训练轮数。一个常见的坑是在测试集D2上做任何形式的模型选择或调参这会严重破坏检验的有效性。所有调优必须严格在训练集D1内部完成。XGBoost基于树的集成方法通常对异构数据混合类型特征、缺失值有较好的鲁棒性且训练速度较快超参数相对容易调节。它通过梯度提升来优化能捕捉复杂的非线性模式。实操心得关键超参数包括max_depth控制树复杂度防止过拟合、n_estimators树的数量太多会过拟合、learning_rate学习率小学习率配合多树通常更好和subsample、colsample_bytree行/列采样引入随机性防止过拟合。同样使用D1内部的交叉验证来调优。重要提示无论选择哪种模型目标都是获得对条件期望h(Z)和g(X)的“最佳可能”估计。这里的“最佳”并非指在独立测试集上的最小预测误差而是指估计量的均方误差MSE收敛速度足够快以满足理论条件(C1)E[(ĝ - g)^2] O(n1^{-a})其中a0n1是训练集大小。实践中我们通过交叉验证选择在D1上预测性能最好的模型配置这通常能保证一个不错的收敛速率。3.3 功效增强与算法稳定性策略单次数据分割的结果可能因随机分割而波动。为了提高检验的稳定性和功效论文提出了两种改进策略。1. 功效增强统计量在计算得到Vn后可以构造一个增强版统计量Vn* Vn τ * Σ_{i∈D2} [ĝ_D1(Xi)]^2其中τ 0是一个调节参数。原理在H0下ĝ_D1(X)估计的是0所以其平方和很小Vn* 与 Vn 渐近等价不影响第一类错误。在H1下ĝ_D1(X)捕捉到了信号其平方和是正的因此 Vn* Vn从而提升了检验功效。实操建议论文建议简单取τ1。在实际代码中可以在计算p值前用Vn* 替代 Vn。但需注意这略微改变了零分布理论上仍是卡方但需要更严格的证明在极端追求第一类错误精确控制时需谨慎。对于大多数应用τ1是一个安全有效的选择。2. 多重数据分割与p值聚合为了消除单次随机分割的偶然性可以进行B次例如B100独立的数据分割每次得到一個p值 p_b (b1,...,B)然后聚合这些p值得到一个更稳健的最终p值。常用聚合方法中位数法Q* 2 * median({p_b})。如果Q* α则拒绝H0。这种方法简单对异常值不敏感。最小p值法经调整Q* min({p_b}) * B。这是Bonferroni校正非常保守。基于均匀分布的方法如果p值在原假设下服从均匀分布可以计算 -2 * Σ log(p_b)它服从卡方分布(2B)。但数据分割导致p值间并非独立此方法可能过于激进。实操建议我推荐使用中位数法。它在控制第一类错误和保持功效之间取得了较好的平衡且计算简单。在代码实现上只需将上述pMIT_single_split函数运行B次收集p值然后计算中位数的两倍即可。def pMIT_multiple_splits(Y, X, Z, W, B100, test_size0.2, estimator_typeDNN, random_seed_start0): 基于多重数据分割的pMIT检验。 p_values [] for b in range(B): p_val, _ pMIT_single_split(Y, X, Z, W, test_sizetest_size, estimator_typeestimator_type, random_staterandom_seed_start b) p_values.append(p_val) # 使用中位数法聚合p值 aggregated_p 2 * np.median(p_values) # 注意聚合p值可能大于1需要截断 aggregated_p min(aggregated_p, 1.0) return aggregated_p, p_values3.4 分割比例ξ的自适应选择训练集大小n1或分割比例ξ n1/N的选择是一个偏差-方差权衡问题。n1越大条件期望估计得越准偏差小但用于构建统计量的测试集n2越小统计量的方差越大功效可能下降。反之亦然。论文提出了一种基于置换检验的数据自适应选择方法设定一个候选比例集合Ξ例如Ξ {0.5, 0.6, 0.7, 0.8, 0.9}。对于每个候选ξ在原数据上随机置换打乱待检验的协变量W破坏其与Y的关系从而模拟H0成立的数据。对每个置换后的数据集应用pMIT算法单次或多次分割计算其p值。重复M次如M200。对于每个ξ计算其经验第一类错误率即p值小于显著性水平α如0.05的比例。选择能满足经验第一类错误率 ≤ α 的最小ξ作为最终的分割比例ξ*。这个方法的逻辑是在H0成立的数据上我们希望检验能控制第一类错误。通过选择能满足这一要求的最小训练集比例我们尽可能为测试集保留了更多的样本以期在H1下获得更高的功效。注意事项这个自适应过程计算量较大因为需要对每个候选ξ进行M次置换检验。在实际应用中如果样本量足够大一个经验性的做法是固定一个较大的训练集比例例如ξ0.8或0.9以确保条件期望的估计足够准确。样本量较小时自适应选择更有价值。4. 偏广义相关性度量从检验到度量当pMIT检验拒绝了H0我们自然想知道“W在给定Z后对Y的均值解释力到底有多强”传统的偏相关系数只能度量线性关联而pGMC提供了一个模型自由的、介于0到1之间的度量。4.1 pGMC的定义与估计pGMC的定义清晰且直观 r²(Y, W|Z) E[{m(X) - h(Z)}²] / E[{Y - h(Z)}²] 其中m(X) E(Y|X), h(Z) E(Y|Z)。分子E[{m(X) - h(Z)}²] 衡量的是当我们从只知道Z进步到知道完整的X时对Y的条件均值预测的改进量均方误差意义下。分母E[{Y - h(Z)}²] 衡量的是如果只用Z来预测Y其不可解释的变异误差方差。因此pGMC解释为W所解释的那部分变异占Z未解释的总变异的比例。它是一个决定系数R²在模型自由和偏相关语境下的推广。基于数据分割和机器学习估计量一个自然的估计量是 r̂² [ Σ_{i∈D2} {m̂_D1(Xi) - ĥ_D1(Zi)}² ] / [ Σ_{i∈D2} {Yi - ĥ_D1(Zi)}² ] 这里D2是独立于模型训练集D1的测试集以确保估计的无偏性或渐近无偏性。4.2 pGMC的性质与置信区间pGMC拥有许多优良性质使其成为一个有吸引力的度量工具取值范围0 ≤ r² ≤ 1。零点意义r² 0当且仅当H0成立即E(Y|X)E(Y|Z)。这意味着给定Z后W不提供任何关于Y均值的额外信息。单位点意义r² 1当且仅当Y m(X) 几乎必然成立即给定Z后Y完全是W通过X的确定性函数。与偏相关系数的关系r² ≥ ρ²。pGMC总是大于等于线性偏相关系数的平方它能捕捉非线性依赖。与GMC的关系r² [GMC(Y|X) - GMC(Y|Z)] / [1 - GMC(Y|Z)]。这提供了另一种计算和理解方式。理论证明在适当的条件下pGMC的估计量 r̂² 是渐近正态的并且收敛速度为根号n测试集大小。这意味着我们可以为其构造渐近置信区间。置信区间构造Wald型 令 θ r²(Y, W|Z)。其估计量为 θ̂ r̂²。 我们需要估计 θ̂ 的渐近方差 σ̂²_θ。这可以通过刀切法或自助法在测试集D2上实现。由于D2独立于训练过程自助法是一个简单可靠的选择。import numpy as np from sklearn.utils import resample def estimate_pGMC_with_CI(Y, X, Z, W, test_size0.2, estimator_typeDNN, n_bootstrap1000, alpha0.05, random_state42): 估计pGMC并计算其自助置信区间。 # 1. 首先进行一次数据分割得到训练好的模型和测试集索引 # 这里简化假设已有训练好的模型 h_estimator, m_estimator 和测试集数据 # 在实际中需要先运行类似 pMIT_single_split 的前半部分来训练模型并得到测试集预测 # 假设我们已经有了 # idx_test: 测试集索引 # h_pred_test: 测试集上 h(Z) 的预测 # m_pred_test: 测试集上 m(X) 的预测 # Y_test: 测试集响应变量 # 计算点估计 numerator np.mean((m_pred_test - h_pred_test) ** 2) denominator np.mean((Y_test - h_pred_test) ** 2) theta_hat numerator / denominator # 2. 在测试集上进行自助法抽样估计标准误 n_test len(Y_test) bootstrap_estimates [] for b in range(n_bootstrap): # 在测试集索引范围内有放回抽样 indices resample(np.arange(n_test), random_staterandom_state b) Y_b Y_test[indices] h_pred_b h_pred_test[indices] m_pred_b m_pred_test[indices] num_b np.mean((m_pred_b - h_pred_b) ** 2) denom_b np.mean((Y_b - h_pred_b) ** 2) theta_b num_b / denom_b bootstrap_estimates.append(theta_b) # 计算自助法标准误和置信区间 se_theta np.std(bootstrap_estimates, ddof1) # 正态近似置信区间 z_alpha 1.96 # for 95% CI ci_lower theta_hat - z_alpha * se_theta ci_upper theta_hat z_alpha * se_theta # 确保置信区间在[0,1]内 ci_lower max(0, ci_lower) ci_upper min(1, ci_upper) # 也可以使用百分位数区间 ci_lower_perc np.percentile(bootstrap_estimates, (alpha/2)*100) ci_upper_perc np.percentile(bootstrap_estimates, (1-alpha/2)*100) return theta_hat, (ci_lower, ci_upper), (ci_lower_perc, ci_upper_perc), se_theta注意当真实pGMC值接近0或1时正态近似可能不佳。百分位数自助区间或偏差校正的自助区间BCa通常更稳健。此外确保分母不为零即Y不能被Z完美预测是定义和估计pGMC的前提。5. 常见问题、陷阱与实战技巧5.1 模型过拟合与欠拟合问题训练集D1上模型过拟合导致 ĥ(Z) 和 m̂(X) 在训练集上表现极好但在独立的测试集D2上预测误差很大。这会严重扭曲检验统计量Tn的分布。欠拟合则会导致估计偏差过大降低检验功效。诊断在D1内部进行交叉验证监控训练误差和验证误差的差距。如果训练误差远小于验证误差可能是过拟合。如果两者都很大可能是欠拟合或模型容量不足。解决对于DNN使用早停法、Dropout、L2正则化、减少网络层数或神经元数量。对于XGBoost降低max_depth、增加min_child_weight、使用更小的learning_rate并配合更多的n_estimators、启用行/列采样。通用策略在D1上使用交叉验证或保留一个验证集来调优超参数。绝对禁止使用测试集D2进行任何模型选择或调参。5.2 样本量不足与分割比例困境问题在高维小样本如n100, p1000场景下数据分割后训练集D1和测试集D2的样本量都可能不足。D1太小导致模型估计不准D2太小导致检验统计量方差大、功效低。应对优先保证训练集规模理论要求γ1/a其中a是估计误差的收敛速率。对于复杂的非线性关系a可能较小如0.2-0.5这意味着需要很大的训练集。在实践中如果总样本量N有限应倾向于给D1分配更大的比例如ξ0.8或0.9牺牲一些测试集样本以换取更准确的估计。考虑使用更简单的模型在样本量很小时过于复杂的模型如很深的DNN弊大于利。可以尝试线性模型带L1/L2正则化、浅层神经网络或树深度受限的XGBoost。利用功效增强启用Vn*统计量τ1可以在不增加样本的情况下提升一些功效。谨慎解读结果当样本量很小时即使检验不显著p值大也不能武断地认为W不重要可能是功效不足。应同时报告pGMC的估计值及其置信区间以量化效应大小及其不确定性。5.3 协变量维度极高时的挑战问题当p尤其是p2即W的维度极高时即使使用DNN或XGBoost估计m(X) E(Y|X)也变得异常困难容易陷入维度灾难。策略特征预筛选在正式应用pMIT之前可以对W进行一轮基于边际关联或简单模型的预筛选保留最相关的一部分特征如top 100或top 500。但这会引入选择偏差需要谨慎。使用专门的高维方法考虑在估计m(X)时使用适用于高维数据的机器学习方法例如带L1正则化的线性模型如果相信关系近似线性或可加。稀疏深度学习使用Dropout、稀疏连接或L1正则化鼓励稀疏性。随机森林/XGBoost它们本身对高维数据有一定鲁棒性通过特征重要性可以进行隐式特征选择。聚焦于Z已知的场景pMIT的优势在于当Z包含已知的重要变量时问题转化为估计E(Ȳ|X)其中Ȳ Y - E(Y|Z)。如果h(Z)能被很好地估计那么Ȳ的变异可能主要由W中少数重要变量驱动从而缓解了高维问题。5.4 计算效率与实现优化问题多重数据分割、自助法置信区间、自适应选择ξ等操作计算开销大。优化建议并行化B次数据分割的pMIT计算、M次置换检验、n_bootstrap次自助抽样这些任务都是完全独立的非常适合并行计算如使用Python的joblib或multiprocessing库。缓存模型在自适应选择ξ或进行多重分割时对于相同的训练集比例ξ和相同的随机种子训练出的模型是相同的。可以设计缓存机制避免重复训练。使用高效库对于XGBoost确保使用GPU加速如果可用。对于DNN使用TensorFlow或PyTorch等框架并利用GPU。设定合理的迭代次数对于探索性分析B50或100可能就够了。对于最终报告B200或500更稳妥。自助法次数n_bootstrap1000通常足够。5.5 与替代方法的比较与选择pMIT不是唯一解决该问题的方法。了解其与相关方法的区别有助于正确选择。方法核心思想优点缺点/注意事项pMIT (本文)重述假设为E(Ȳ|X)0数据分割ML估计构造卡方统计量零分布为标准卡方计算快理论性质清晰提供了pGMC度量需要数据分割损失部分样本信息对ML估计精度有要求Dai et al. (2024)比较两个预测误差的差异引入随机扰动处理退化分布可处理更一般的条件矩限制需要选择扰动大小理论局部备择检出速率可能稍慢Williamson et al. (2023)双数据分割分别估计两个非零参数概念相对直观需要至少两个独立分割样本利用可能更低Lundborg et al. (2022) PCM投影协方差度量基于残差正交化在线性模型下功效可能很高在模型自由设定下需要三个以上子样本计算多个条件期望计算量大传统F检验基于线性模型残差平方和比较计算极其简单在线性情况下最优严重依赖线性模型假设误设时失效选择建议如果追求计算简便和解释直观且相信关系近似线性可以先尝试传统F检验或线性模型的偏F检验。如果处理的是明确的高维线性模型Lundborg et al. (2022)的PCM方法可能是一个强有力的竞争者。如果面对的是高维、非线性、模型自由的设定并且希望有一个统一、有理论保障、附带效应大小度量的框架pMITpGMC的组合是一个非常全面和实用的选择。它的标准卡方推断避免了重抽样计算效率高pGMC提供了有意义的解释。6. 总结与扩展思考偏均值独立性检验及其关联的pGMC度量为我们分析高维数据中变量间的条件依赖关系提供了一套强大的工具。它巧妙地将机器学习的预测能力与严谨的统计推断相结合通过假设重述和数据分割解决了推断中的根本难题。在实际应用中我个人的体会是理解问题的本质比机械地套用算法更重要。pMIT检验回答的是“W是否在给定Z后对Y的均值有贡献”这是一个关于增量贡献的问题。因此如何定义Z已知的、需要控制的变量集至关重要。Z应该包含所有你确信与Y相关且可能与W混淆的变量。例如在基因研究中Z可能是已知的致病通路基因W是待探索的其他基因。此外pMIT检验的结论依赖于机器学习模型对条件期望估计的“足够好”。虽然理论只要求一定的收敛速率但在有限样本下一个糟糕的估计会导致检验功效低下。因此投入精力进行谨慎的模型选择和验证是值得的。不要仅仅把它当作一个黑箱流程。最后pGMC作为一个效应大小度量其价值不亚于检验本身。一个显著的p值告诉我们关联存在而pGMC的估计值及其置信区间则告诉我们这个关联的强度这在许多科学领域对于评估实际意义至关重要。报告结果时应同时提供p值和pGMC估计及置信区间以给出更完整的画面。这个方法可以自然地扩展到其他相关问题例如检验条件独立性通过检验一系列条件均值或者处理分类响应变量通过适当修改损失函数和模型。其核心思想——利用机器学习估计复杂函数再通过巧妙的统计构造进行推断——为高维统计推断开辟了一条富有前景的道路。
