如何计算两样本Hotelling's T²检验的功效与第二类错误

核心计算逻辑 #

检验功效的本质是：当原假设“两组多元均值完全相等”不成立（两组真实存在均值差异）时，检验能正确拒绝原假设的概率；对应的第二类错误率为1 - 检验功效。
整个实现流程完全复用已有计算实际显著性水平的模拟框架，不需要依赖专门的power工具包。

具体实现步骤 #

• 样本生成阶段做唯一调整：不再给两组设置完全相同的均值向量，人为给两组设置已知、非零的均值偏移，保证原假设从生成样本的阶段就不成立；协方差矩阵、样本量、多元维度等其他参数，要和之前测试实际显著性水平的场景完全保持一致，避免场景错位
• 和之前计算第一类错误的流程一致，重复抽样10000次（重复次数越高结果越稳定），每次抽样后运行HotellingsT2()检验获取p值
• 统计所有重复实验中，p值小于等于检验水准（测试中使用的0.05）的次数占总重复次数的比例，这个比例就是对应场景下的实际检验功效
• 第二类错误率可直接用1 - 实际功效计算，不需要额外流程

可直接运行的参考代码 #

library(MASS)
library(ICSNP)

# 全局参数与之前的第一类错误测试场景保持完全一致
set.seed(123)
n_rep = 10000
alpha = 0.05
dim_p = 10
n_per_group = 50
cov_X = diag(rep(1, dim_p))
cov_Y = diag(rep(3, dim_p))
# 可根据测试需求调整均值差大小，对应不同的效应强度
mean_offset = 0.5
reject_count = 0

for (i in 1:n_rep) {
  # 生成样本：X组均值全为0，Y组每个维度偏移mean_offset
  X <- mvrnorm(n_per_group, rep(0, dim_p), cov_X)
  Y <- mvrnorm(n_per_group, rep(mean_offset, dim_p), cov_Y)
  p_val <- HotellingsT2(X, Y)$p.value
  if (p_val <= alpha) {
    reject_count = reject_count + 1
  }
}

empirical_power = reject_count / n_rep
type2_err = 1 - empirical_power

cat("当前场景实际检验功效：", round(empirical_power*100, 2), "%\n")
cat("当前场景第二类错误率：", round(type2_err*100, 2), "%\n")

实用注意技巧 #

• 测试不同效应下的功效变化时，只需要调整代码中mean_offset的取值即可，固定其他参数不动，就能得到对应偏移下的功效值
• 当测试场景严格满足“多元正态分布、协方差矩阵相等”的检验前提时，可通过非中心F分布的理论值做结果对照：Hotelling's T²检验在备择假设成立时，统计量服从非中心F分布，用基础函数pf()即可计算理论功效值，不需要额外推导公式
• 当前测试的“两组协方差不相等”场景，本身是不满足Hotelling's T²检验前提的，不存在现成的理论功效值，所有专门的power工具包计算的都是前提满足时的理想功效，和这类模拟场景不匹配，手动模拟的结果是该场景下最准确的实际性能参考
• 大维度、大重复次数的模拟场景下，可将循环内逐次抽样的逻辑替换为一次性生成全量模拟样本的向量化操作，能大幅提升运行速度

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