统计

理解统计效力

LOSSES

在实际的数据分析场景中,我们总是希望得到一个确定性的结果。先前我们讨论到的 Fisher 先生提出的显著性系统作为一种证据评估框架不能帮助我们做出决策,它只告诉我们证据的强度为何,具体的决策需要依赖研究者的心证。这中语带保留的态度体现出了一名统计学家面对未知的谦逊态度,但它给实际的应用场景带来了很多麻烦。为了应对这个问题, Neyman 和 Pearson 延续了 Fisher 「反证」的基本精神,发展出了一套新的假设检验分析体系。尽管 Fisher 先生对于这个演绎有极大的不满,但这套体系已经近乎成为了社科领域的事实标准。

这套体系最重要的贡献是引入了假阳性和假阴性[1]这两个重要的概念。通常期刊里面只看假阳性,也就是原本没有效应的东西被你推断成了有效应。这相当于向学界引入了一个假知识,自然需要被编辑警惕。但是期刊并不 Care 假阴性,也就是原本有效应的东西被你算出了没效应。毕竟实验花得是你的钱不是人家杂志社的钱。

但对于研究者而言,原本应该找到的效应找不到是一件让人万分惋惜的事。另外一方面,或许很多研究者都没有意识到这个问题,很多研究在开始执行的那一刻就已经注定做不出阳性了。

让我们来做一个简单的统计学实验来验证这件事情。你和我将共同扮演一个全知全能的上帝,我们创建一个由小玩偶构成的世界,这个世界有两万个活灵活现的玩偶一起快乐生活。每个玩偶有自己的身高和种族(玩具小熊和玩具浣熊)。我们假定玩具小熊和玩具考拉的平均身高有差异,一者是 50cm,一者是 52cm。

# SHARE
set.seed(2026)
n <- 20000

race <- sample(c("bear", "raccoon"), size = n, replace = TRUE, prob = c(0.5, 0.5))

height <- ifelse(race == "bear",
                 rnorm(n, mean = 50, sd = 8),
                 rnorm(n, mean = 52, sd = 8))

height <- round(height, 0)

population <- data.frame(
  id = 1:n,
  race = race,
  height = height
)

#NOSHARE
summary(population)

作为上帝,我们明确地规定了,小熊和小浣熊的身高的确是有差异的。接下来,我们来做模拟抽样的动作。既然我们是上帝,那么为什么不把阵仗搞得大一点。我们假设有一千个来自不同星球的外星人团队来到了这个奇妙世界,分别测量了这些小生物的身高。每个团队都整齐划一地抽了 50 个小熊和 50个小浣熊。一百个样本!算是大样本研究了吧!

下面让我们来看看双样本 t 检验是否能够正确地发现小熊和小浣熊的身高「真的有差」。为了控制假阳性的概率,我们规定,如果小熊和小浣熊的身高之间没有差异的话,我们可以接受 20 次抽出一次假阳性的风险。

alpha   <- 0.05
n_sim   <- 1000
n_each  <- 50
df      <- 2 * n_each - 2
t_crit  <- qt(1 - alpha / 2, df = df)   # 双尾临界值

bears_pool    <- population$height[population$race == "bear"]
raccoons_pool <- population$height[population$race == "raccoon"]

run_one <- function() {
  b <- sample(bears_pool,    n_each)
  r <- sample(raccoons_pool, n_each)
  t <- t.test(b, r)$statistic
  c(diff = mean(b) - mean(r), t = t)
}

res      <- as.data.frame(t(replicate(n_sim, run_one())))
res$rej  <- abs(res$t) > t_crit   # 落入拒绝域 = 拒绝 H0

cat("临界值 ±", round(t_crit, 3), "(α =", alpha, ",df =", df, ")\n")
cat("检出差异:", sum(res$rej),  "次,占比", round(mean(res$rej) * 100, 1), "%\n")
cat("伪阴性:  ", sum(!res$rej), "次,占比", round(mean(!res$rej) * 100, 1), "%\n")
cat("均值差均值:", round(mean(res$diff), 2), ",标准差:", round(sd(res$diff), 2), "\n")

hist(res$diff,
     breaks = 40,
     col    = "#aaaaaa",
     border = "white",
     main   = paste0("1000次抽样:两组均值之差(α = ", alpha, ")"),
     xlab   = "bear 均值 − raccoon 均值",
     ylab   = "频次")

hist(res$diff[res$rej],
     breaks = 40,
     col    = "#e05c5c",
     border = "white",
     add    = TRUE)

abline(v = 0, lty = 2, lwd = 1.5)

legend("topleft",
       legend = c(paste0("阳性(", sum(res$rej),  "次)"),
                  paste0("阴性(", sum(!res$rej), "次)")),
       fill   = c("#e05c5c", "#aaaaaa"),
       border = "white",
       bty    = "n")

wow,令人印象深刻,竟然超过七成的检验都没能发现这个真实存在的差异。现在我们非常直观地看到了,我们的统计方案没能发现这重要的差异。对于这种情况,我们将其称作「统计效力不彰」。实际研究中,我们只是那一千个研究团队中的一个,在此等统计效力之下盲目进行实验无异于一场赌博。

在公式层面,有四个因素影响统计效力:样本量(外星人们测量多少个小熊和小浣熊的身高)、效应量(在这里是小熊和小浣熊身高的实际差异)、判断标准(你能接受伪阳性的风险是多少)、噪声水平(在这里是两种小动物身高的标准差)。

效应量和噪声水平是数据的固有性质,你控制不了。判断标准和样本量你的确能控制。但是放松标准真的是一个正确的决定吗?为了演示这点,让我们再来扮演一次上帝,这次我们建立一个新的世界,这个世界里面依有小鸭子和小鸡两种玩偶生物。他们的平均身高一模一样,我们完整地模拟之前的实验,把标准放松成如果小鸡和小鸭的身高之间没有差异的话,我们可以接受 5 次抽出一次假阳性的风险。

set.seed(2026)
n <- 20000

race <- sample(c("duck", "chicken"), size = n, replace = TRUE, prob = c(0.5, 0.5))

height <- ifelse(race == "duck",
                 rnorm(n, mean = 50, sd = 40),
                 rnorm(n, mean = 50, sd = 40))

height <- round(height, 0)

population <- data.frame(
  id = 1:n,
  race = race,
  height = height
)

alpha   <- 0.2
n_sim   <- 1000
n_each  <- 50
df      <- 2 * n_each - 2
t_crit  <- qt(1 - alpha / 2, df = df)   # 双尾临界值

duck_pool    <- population$height[population$race == "duck"]
chicken_pool <- population$height[population$race == "chicken"]

run_one <- function() {
  b <- sample(duck_pool,    n_each)
  r <- sample(chicken_pool, n_each)
  t <- t.test(b, r)$statistic
  c(diff = mean(b) - mean(r), t = t)
}

res      <- as.data.frame(t(replicate(n_sim, run_one())))
res$rej  <- abs(res$t) > t_crit   # 落入拒绝域 = 拒绝 H0

cat("临界值 ±", round(t_crit, 3), "(α =", alpha, ",df =", df, ")\n")
cat("阴性:  ", sum(!res$rej), "次,占比", round(mean(!res$rej) * 100, 1), "%\n")
cat("伪阳性:", sum(res$rej),  "次,占比", round(mean(res$rej)  * 100, 1), "%\n")
cat("均值差均值:", round(mean(res$diff), 2), ",标准差:", round(sd(res$diff), 2), "\n")

hist(res$diff,
     breaks = 40,
     col    = "#aaaaaa",
     border = "white",
     main   = paste0("1000次抽样:两组均值之差(α = ", alpha, ")"),
     xlab   = "duck 均值 − chicken 均值",
     ylab   = "频次")

hist(res$diff[res$rej],
     breaks = 40,
     col    = "#e05c5c",
     border = "white",
     add    = TRUE)

abline(v = 0, lty = 2, lwd = 1.5)

legend("topleft",
       legend = c(paste0("阳性(", sum(res$rej),  "次)"),
                  paste0("阴性(", sum(!res$rej), "次)")),
       fill   = c("#e05c5c", "#aaaaaa"),
       border = "white",
       bty    = "n")

瞧,字面意义上我们得到了 20% 的伪阳性概率。

这时你可能会问,是不是我大力出奇迹,采超多数据就能做出统计效力更高的分析了?当然,但现实生活中收数据很有可能非常贵,做实验是要钱的,出于对钱包的怜悯我们也不能无脑堆高数据量。

这时可能又会问,那我能不能分阶段采数据,先收集一小部分数据(比如先采集十个人的数据),看看有没有差异,如果没有的话,追加一部分数据(比如说追加十个人),再算一次,差异还不够大的话就再追加数据直到算出阳性为止。看起来很聪明,但实际上这个做法偷换了「当效应不存在」时均值差异的分布,让我们实际演示这个行为造成的影响。我们就按照这个方案,在小鸡和小鸭的世界里(身高没有实质差异的世界)模拟一千次实验,看看均值的分布会发生什么变化。

alpha   <- 0.05
n_sim   <- 1000
n_max   <- 50
n_start <- 10
n_step  <- 10

race <- sample(c("duck", "chicken"), size = n, replace = TRUE, prob = c(0.5, 0.5))

population <- data.frame(
  id = 1:n,
  race = race,
  height = height
)

duck_pool    <- population$height[population$race == "duck"]
chicken_pool <- population$height[population$race == "chicken"]


# 固定样本量方案:每次直接抽满 n_max 个
run_fixed <- function() {
  b <- sample(duck_pool,    n_max)
  r <- sample(chicken_pool, n_max)
  p <- t.test(b, r)$p.value
  c(diff = mean(b) - mean(r), rej = as.numeric(p < alpha))
}

# 追加数据方案:看到阳性就停,否则追加数据直到 n_max
run_sequential <- function() {
  n <- n_start
  repeat {
    b <- sample(duck_pool,    n)
    r <- sample(chicken_pool, n)
    p <- t.test(b, r)$p.value
    if (p < alpha || n >= n_max) {
      return(c(diff = mean(b) - mean(r), rej = as.numeric(p < alpha)))
    }
    n <- n + n_step
  }
}

res_fixed <- as.data.frame(t(replicate(n_sim, run_fixed())))
res_seq   <- as.data.frame(t(replicate(n_sim, run_sequential())))

cat("固定采样伪阳性:", sum(res_fixed$rej), "次,占比", round(mean(res_fixed$rej) * 100, 1), "%\n")
cat("追加采样伪阳性:", sum(res_seq$rej),   "次,占比", round(mean(res_seq$rej)   * 100, 1), "%\n")

xlim_range  <- range(c(res_fixed$diff, res_seq$diff))
breaks_shared <- seq(xlim_range[1], xlim_range[2], length.out = 41)

y_max <- max(
  hist(res_fixed$diff, breaks = breaks_shared, plot = FALSE)$counts,
  hist(res_seq$diff,   breaks = breaks_shared, plot = FALSE)$counts
) * 1.15

par(mfrow = c(2, 1), mar = c(4, 4, 3, 2))

# 固定采样
hist(res_fixed$diff,
     breaks = breaks_shared,
     col    = "#aaaaaa", border = "white",
     xlim   = xlim_range, ylim = c(0, y_max),
     main   = paste0("固定采样(伪阳性 ", sum(res_fixed$rej),
                     " 次,", round(mean(res_fixed$rej) * 100, 1), "%)"),
     xlab   = "均值之差", ylab = "频次")

hist(res_fixed$diff[as.logical(res_fixed$rej)],
     breaks = breaks_shared,
     col    = "#e05c5c", border = "white", add = TRUE)

abline(v = 0, lty = 2, lwd = 1.5)
legend("topleft",
       legend = c(paste0("伪阳性(", sum( res_fixed$rej),        " 次)"),
                  paste0("阴性(",   sum(!as.logical(res_fixed$rej)), " 次)")),
       fill   = c("#e05c5c", "#aaaaaa"), border = "white", bty = "n")

# 追加采样
hist(res_seq$diff,
     breaks = breaks_shared,
     col    = "#aaaaaa", border = "white",
     xlim   = xlim_range, ylim = c(0, y_max),
     main   = paste0("追加采样(伪阳性 ", sum(res_seq$rej),
                     " 次,", round(mean(res_seq$rej) * 100, 1), "%)"),
     xlab   = "均值之差", ylab = "频次")

hist(res_seq$diff[as.logical(res_seq$rej)],
     breaks = breaks_shared,
     col    = "#e05c5c", border = "white", add = TRUE)

abline(v = 0, lty = 2, lwd = 1.5)
legend("topleft",
       legend = c(paste0("伪阳性(", sum( res_seq$rej),        " 次)"),
                  paste0("阴性(",   sum(!as.logical(res_seq$rej)), " 次)")),
       fill   = c("#e05c5c", "#aaaaaa"), border = "white", bty = "n")

par(mfrow = c(1, 1))

我们可以清晰地看到伪阳性率膨胀了非常多。我们可以从两个角度解释这件事,首先,对于每一个「外星人研究团队」,他们并不是只做了一次 t 检验。从 10 人开始测验,如果阴性就追加样本直直样本有 50 人才放弃,这意味着一个团队最多做了五次测验,每次测验都要独立面对 5% 的伪阳性风险,那么连续做五次测验的伪阳性风险就会提升到 1 − (1 − 0.05)⁵ ≈ 22.6% [2]。此外,相信你一定记得 t 分布当中有一个和样本量有关的参数:自由度。这个参数是固定的,但是在这种「全新的抽样检测技术」之下,样本量不再是固定的了,它有可能是 10、20、30、40、50。此时实验真正的抽样参数分布变成了上述五种分布的联合分布,这也是为什么你看到最后做出来的图,均值的分布变得又矮又宽。

上述现象进一步解释了,为什么现代社科实验都鼓励研究者做预注册,一旦固定了样本量之后就不允许在不经统计学处理的情况下胡乱添加样本,除非你有能力使用正确的统计工具处理这个复杂情况。因为在这种情况下依然使用朴素的分析技术,在事实层面上膨胀了伪阳性率[3],哪怕你的初衷是为了提高统计效力。

在 Neyman 和 Pearson 的统计体系中,实验设计是数据分析品质控制的重要一环,不可以被剥离出去独立思考。

这时你可能又会问,样本量小会导致统计效力低,样本量大又会让实验变得很贵。那我到底该怎么做才好?学界有一种流行的做法,在正式大规模实验之前,先做小样本预实验(比如先采集三十分数据),计算效应量(标准化组间、相关系数、beta* 等)然后告诉统计软件你可以接受的伪阴性风险是多少,统计软件就会告诉你最为经济实惠的样本量了。

但这里有一个需要注意的地方,通常来讲小样本实验的效力天然很弱,所以在填写效应量的时候,你应当注意填写预设置信区间的下限结果,否则很有可能会得到一个偏小的样本量,导致错失本应能发现的效应。值得注意的是,有研究指出基于这样的方式得出来的统计效应依然会虚高导致正式实验的时候得不出阳性结果。这在学界是一个缺乏统一见解的灰色区域。作为一名同样深陷泥潭的研究者,我只能姑且建议你把标准设置的严格一点,让置信区间比较宽,以得到更为保守的效应量,进一步得到更大的样本总数估计值。

另外,实际上你可以自己决定「最小且有意义的效应量」是多少,如果效应量小于这个数字,那你发现出来了阳性也没啥价值。如果你对这个数字心理有谱,那么直接放到统计效力分析工具里面就好。

如果你在做的分析已经被充分复现过、有领域知识,那么搜搜元分析,看看大家的效应量都是多少,再来决定对自己研究效应量的期待也是一个不错的做法。


  1. 大多数教材包括理论提出者本人将这二者称作一类错误和二类错误。但是本文拒绝这种称呼。一来理论构建者自己也坦然承认过这就是按顺序给的两个数字,这个名字本身没有表义功能,二来,这个随意的命名几乎造成了认知灾难,研究者每次都要在脑袋里面做一次人工转换。从人因的角度来看,人们会倾向于逃避复杂的概念,这种因随意带来的认知负担应当被避免。我写这段内容也只是为了让你在读教材的时候能把这些名词对起来。 ↩︎

  2. 这里只是一个粗略估计,实际情况要比这个公式更复杂一点。 ↩︎

  3. 这种看情况添加样本的方法被称作序贯分析,有专门属于它的数据分析方法,如果你真的遇到了需要使用这种样本添加方法的情况,必须转而使用正确的分析工具,不可以使用最简单的分析工具。 ↩︎

Comments

Loading animation

Loading comments...