Likan Zhan

正确理解和使用P值

侃侃迩行 · 2017-03-01

什么是P值?

Informally, a p-value is the probability under a specified statistical model that a statistical summary of the data (e.g., the sample mean di erence between two compared groups) would be equal to or more extreme than its observed value.

P值检验的是研究者采集到的样本与某一个统计假设的相容程度,而不是这个统计假设的正确与否。

1. P值、效力和效应

研究者通常希望确定一个理论或假设正确与否。以独立样本T检验的典型应用情境为例,我们希望确定某种疾病是否影响人的某个生理指标,从而为将来该疾病者的诊断和治疗奠定基础。该研究问题会产生两个假设:

大多数情况下,我们通常无法以一种二元论的方式确定该理论正确还是错误,只能从概率层次上对该问题进行回答:该疾病对某生理指标有影响的概率是 \(\pi_0\),没有影响的概率是 \(1-\pi_0\)。当\(\frac{\pi}{1-\pi}\)的比率足够大时,我们就说该疾病影响该生理指标

实际应用中让这个问题变得更复杂的是:研究者通常没有办法收集这两个被试群体的所有数据,而只能从其中选取一个样本。我们要从收集到的样本数据去推论这两个整体的关系。样本数据也会产生两种结果:

由于整体和样本之间的不同关系,实验会产生四种不同的结果:

第一种情况发生在零假设正确时,即疾病实际上不影响生理指标时。此时如果因为取样误差,样本数据显示疾病影响生理指标,这就产生了所说的虚报或者叫I类错误。虚报发生的概率就是我们通常所说的P值:是当研究总体符合零假设的情况下,我们收集到该样本数据的可能性,\(\alpha = P(A|(1-\pi_0)\)

n <- 10
sd <- 1
p_value <- 0.05
Cohen_D <- 4

Mean0 <- 0
t_value <- qt(p_value / 2, df = n, lower.tail = FALSE)
lw <- Mean0 - t_value * sd
up <- Mean0 + t_value * sd

X0 <- Mean0 + seq(-4, 4, by = 0.01) * sd
Y0 <- dt(X0, df = n)
Mean1 <- Mean0 + Cohen_D * sd
X1 <- X0 + Mean1
Y1 <- Y0

Xmin <- min(c(X0, X1))
Xmax <- max(c(X0, X1))
plot(X0, Y0, type = "n", xlab = "", ylab = "",
     main = NULL, axes = FALSE, xlim = c(Xmin, Xmax))
lines(X0, Y0)
axis(1, at = seq(Xmin, Xmax, sd), pos = 0)
mtext(expression(H[0]), side = 3,
      adj = 4.15 / (Xmax-Xmin), padj = 0.4)
abline(v = t_value, lty = 3)

polygon(c(min(X0[X0 >= up]),
              X0[X0 >= up],
          max(X0[X0 >= up])),
        c(0, Y0[X0 >= up], 0),
        col = rgb(1, 0, 0, 0.4))
polygon(c(min(X0[X0 <= lw]),
              X0[X0 <= lw],
          max(X0[X0 <= lw])),
        c(0, Y0[X0 <= lw], 0),
        col = rgb(1, 0, 0, 0.4))
text(x = t_value + 0.4, y = -0.01,
     expression(alpha/2), pos = 3, col = "blue")

第二种情况也发生在零假设正确,即疾病实际上不影响生理指标时,研究者的样本数据也没发现有影响的证据,即正确拒斥。正确拒斥的计算方法为 \(1-\alpha = P((1-A)|(1-\pi_0))\)

plot(X0, Y0, type = "n",
     xlab = "", ylab = "",
     main = NULL, axes = FALSE,
     xlim = c(Xmin, Xmax))
lines(X0, Y0)
axis(1, at = seq(Xmin, Xmax, sd), pos = 0)
mtext(expression(H[0]), side = 3,
      adj = 4.15 / (Xmax-Xmin), padj = 0.4)
abline(v = t_value, lty = 3)

polygon(c(min(X0[X0 >= lw & X0 <= up]),
              X0[X0 >= lw & X0 <= up],
          max(X0[X0 >= lw & X0 <= up])),
        c(0, Y0[X0 >= lw & X0 <= up], 0),
        col = rgb(0, 0, 1, 0.4))
text(x = Mean0, y = 0.1,
     expression(1-alpha), pos = 3)

第三种情况发生在备择假设正确,即疾病的确影响该生理指标时。此时如果样本数据没能发现疾病对生理指标的影响,我们就说产生了漏报,或者叫II类错误。二类错误产生的概率按如下公式计算:\(\beta=P((1-A)|\pi_0)\)

plot(X1, Y1, type = "n", xlab = "", ylab = "",
     main = NULL, axes = FALSE, xlim = c(Xmin, Xmax))
lines(X0, Y0)
lines(X1, Y1)
mtext(expression(H[0]), side = 3,
      adj = 4.15 / (Xmax-Xmin), padj = 0.4)
abline(v = t_value, lty = 3)
axis(1, at = seq(Xmin, Xmax, sd), pos = 0)
mtext(expression(H[1]), side = 3,
      adj = (4 + Cohen_D) / (Xmax-Xmin), padj = 0.4)
polygon(c(min(X1[X1 <= up]),
              X1[X1 <= up],
          max(X1[X1 <= up])),
        c(0, Y1[X1 <= up], 0),
        col = rgb(0.5, 0, 0.5, 0.4))
text(x = t_value - 0.35, y = 0.00,
        expression(beta), pos = 3, col = "blue")

第四种情况也发生在备择假设正确,即疾病的确影响该生理指标时。此时如果样本数据能正确发现疾病对生理指标的影响,我们就说产生了击中。统计中击中的概率就叫做`统计效力’(statistical power),计算公式如下:\(1-\beta=P(A|\pi_0)\)

plot(X1, Y1, type = "n", xlab = "", ylab = "",
     main = NULL, axes = FALSE, xlim = c(Xmin, Xmax))
lines(X0, Y0)
lines(X1, Y1)
mtext(expression(H[0]), side = 3,
      adj = 4.15 / (Xmax-Xmin), padj = 0.4)
abline(v = t_value, lty = 3)
axis(1, at = seq(Xmin, Xmax, sd), pos = 0)
mtext(expression(H[1]), side = 3,
      adj = (4 + Cohen_D) / (Xmax-Xmin), padj = 0.4)
polygon(c(min(X1[X1 >= up]),
              X1[X1 >= up],
          max(X1[X1 >= up])),
        c(0, Y1[X1 >= up], 0),
        col = rgb(0.5, 0.5, 0.5, 0.4))
text(x = Mean1, y = 0.1,
        expression(1-beta), pos = 3, col = "blue")

上述四种情况可以总结为下面这个表格:

\(H_0\) \(H_1\)
没有拒绝零假设 正确拒斥(\(1-\alpha\) 漏报(II类错误,\(\beta\)值)
拒绝零假设 虚报(I类错误,\(\alpha\)值) 击中(统计效力,\(1-\beta\)

效应值(Effect size)指零假设分布和备择假设分布的标准化距离,如上 图\(H_0\)的平均值和\(H_1\)平均值之间的标准化距离,即\(D = \frac{M_{H0} - M_{H1}}{SD_p}\)\(SD_p\)\(H_0\)\(H_1\) 的合并标准差。根据上图可知,效应值大小影响了零假设分布和备择假设分布之间的分离程度。效应值越大,同样P值情况下,统计的效力就越大。

在效应值不变的情况下,P值越大,统计效力就越大。独立样本t检验中,假设两个样本的样本量均为 \(n=10\),标准差均为 \(\sigma^2=1\),效应值为 \(D=1.32\),那么统计效力和P值的关系如下图:

pvalue <- seq(0, 0.1, by = 0.001)
ptt <- function (pvalue) {
       pt <- power.t.test(
             n = 10,
             delta = 1.32,
             sd = 1,
             sig.level = pvalue,
             type = "two.sample",
             alternative = "two.sided"
             )
      ptt <- pt$power
      return(ptt)
}
power <- sapply(pvalue, FUN = ptt)

plot(pvalue, power,
     ylim = c(0, 1),
     type = "l", lwd = 2,
     font.lab = 2,
     xlab = expression(alpha),
     ylab = "Power")
lines(x = c(0.05, 0.05),
      y = c(-1, ptt(0.05)),
      lty = 3)
lines(x = c(-1, 0.05),
      y = c(ptt(0.05), ptt(0.05)),
      lty = 3)
lines(x = c(0.005, 0.005),
      y = c(-1, ptt(0.005)),
      lty = 3)
lines(x = c(-1, 0.005),
      y = c(ptt(0.005), ptt(0.005)),
      lty = 3)

此时假如两个组别平均值差别为 \(H = 1.2\),两个组别的合并标准差(pooled SD)为 \(S_p = 1.1\),那么此时的P值约为0.025。其计算方式如下P:

n <- 10
HD <- 1.2
Sp <- 1.1
t <- HD / (Sp * sqrt(2/n))
df <- 2 * (n - 1)
pvalue0 <- 2 * pt(t, df, lower.tail = FALSE)
pvalue0
## [1] 0.02529368

所以,P值本身并不能检验两个前提假设的正确与否。但是研究者真正关心的却是前提假设正确与否。在我们的例子中就是:某项疾病是否影响某项生理指标的水平。美国统计学会在2016发文专门阐述P值的文章的附件种,18篇简短评论提供了从P值推论出前提假设正确与否的不同可能性。文章总结出了三种主要的方法。

2. 贝叶斯因子

如前所述P值指的是零假设正确的前提下,我们得到某结果的概率:\(P(A|\pi_0)\);而研究者关心的是在得出某结果的前提下,零假设的正确与否,即\(P(\pi_0|A)\)。它本质上是一个先验概率和后验概率的贝叶斯问题。所以理论上来说可以用贝叶斯方法解决该问题。但是贝叶斯分析必须在零假设和备择假设两种情境下主观的确定一个模型参数的先验分布,这增加了贝叶斯方法的主观性。

为此,在对美国统计学会撰文的评论中, Benjamin 和 Berger提出可以只用P值计算一个贝叶斯因子(Bayes factor)的上限,\(\overline{B}\)。其计算公式如下:

\[ \overline{B}\leq -\frac{1}{e\times P \times ln(P)} \] 其中P是我们计算出来的显著性水平。

这个上限的计算不需要假定一个先验概率。其含义也比较简单。例如 \(\overline{B}=10\) 表示备择假设正确可能性与比零假设正确可能性之比的最大值是10。因为它表示了备择假设和零假设中那一个假设更可能正确,所以贝叶斯成分可以被用来做假设检验。通常认为当贝叶斯因子为20或者更大时,我们就有充分的理由接受备择假设,拒绝零假设。

贝叶斯因子和P值的关系可以用下图表示:

B <- function(P) {
     -1 / ( exp(1) * P * log(P))
  }
plot(pvalue, B(pvalue),
     ylim = c(0, 55),
     type = "l", lwd = 2,
     ylab = expression( bar(B)),
     xlab = "P",
     font.lab = 2
     )
lines(x = c(0.05, 0.05),
      y = c(-100, B(0.05)),
      lty = 3)
lines(x = c(-100, 0.05),
      y = c(B(0.05), B(0.05)),
      lty = 2)
lines(x = c(0.005, 0.005),
      y = c(-100, B(0.005)),
      lty = 3)
lines(x = c(-100, 0.005),
      y = c(B(0.005), B(0.005)),
      lty = 2,
      col = "blue")

按照上述标准,当显著性水平 \(\alpha=0.05\) 时,贝叶斯因子仅为 \(\overline{B} \leq 2.5\)。回到我们之前的例子,当\(P=0.025\)时,备择假设正确与零假设正确的比率最大不会超过3.9,即 \(\overline{B}\leq 3.9\)。如果我们坚持\(\overline{B}\)应该大于20,那么此时的P值应该小于0.0032才行。Johnson认为P值小于0.005就可以了。

贝叶斯因子的计算不涉及测验的统计效力。所以如果我们假设效应值(effect size)仍然为1.32,但在\(\alpha<0.005\)的水平上拒绝零假设,此时的统计效力(power)仅为43%。如下图:

plot(B(pvalue), power,
     ylim = c(0, 1),
     type = "l", lwd = 2,
     xlab = expression( bar(B)),
     font.lab = 2
     )

此时要达到80%的统计效力,我们的样本量应该提高到 \(n=18\)才行。

3. 错误发现率

同样是对美国统计学会撰文的评论中,Altman 认为研究者应该用错误拒绝率(false discovery rate,FDR)来代替P值,来检验前提假设的正确与否。一个测验中,如果零假设的真正概率是\(\pi_0\),那么备择假设的真正概率为\(1-\pi_0\)。假设我们在显著性水平为\(\alpha\)上拒绝零假设,即\(P<\alpha\),此时我们错误拒绝的概率应该是\(\alpha\pi_0\)。假如该检验的统计效力为\(\beta\),那么我们拒绝该测验的概率为\(\beta(1-\pi_0)\)

据此,作者把错误发现率(false discovery rate,FDR)定义为错误拒绝的期望概率与所有拒绝的期望概率的比率。即 \[ eFDR = \frac{\alpha\pi_0}{\alpha\pi_0+\beta(1-\pi_0)} \]

对于先验概率\(\pi_0\)的确定,Altman提出了如下一些拇指原则:在高通量情景中(high-throughput scenario),\(\pi_0\)可以经由数据估计出来;在低通量的实验中( low-throughput testing),\(\pi_0\)的确定应该根据先验证据中零假设的概率来确定:对一个研究设想的首次、第二次验证时,\(\pi_0\)可以等于0.5或0.75;但在事后检验(post hoc test)中,\(\pi_0\)应该设定为0.95。

错误拒绝率和先验零假设概率的关系如下图:

pi0 <- seq(0, 1, by = 0.01)
alpha <- 0.05
power <- ptt(alpha)
eFDR <- (alpha * pi0) / (alpha * pi0 + power * (1 - pi0) )
plot(pi0, eFDR, type = "l",
     xlab = expression(pi[0])
     )

在我们的例子中,当\(P=0.025\)、统计效力为\(80\%\) 时,eFDR在首次、第二次和事后检验中分别为0.03, 0.09和0.38。

该方法的一个缺点是\(\pi_0\)的确定有一定的随意性。回避这种随意性的一个简单的方法是通过让FDR达到某种水平的方式确定\(\pi_0\)。例如在我们的例子中,在统计效力为80%的前提下,要使得\(FD=0.05\)\(\pi_0\)应该不大于0.62才行。

4. 置信区间

在对美国统计学会文章的评论中,大部分研究者提出的建议是:在报告统计结果时,应该既报告P值,也报告置信区间(confidence intervals)。因为置信区间包含了效应值大小(effect size)的信息。

作者最后强调,P值本身是一个随机变量。当数据是连续性数据、而且零假设正确时,P值是在[0, 1]的一个均匀分布(uniform distribution)。这也就是说,P值在不同样本间的变异是非常大的。而且这种变异并不随着实验样本量和统计效力的变化而变化。

当备择假设成立时,P值的变异会随着统计效力的增加而变小。尽管如此,P值仍然是随机分布的。

5. 延伸阅读

Altman, N., & Krzywinski, M. (2017). Points of Significance: Interpreting P values. Nat Meth, 14(3), 213-214. doi:10.1038/nmeth.4210

Nuzzo, B. (2014). Statistical errors. Nature, 506, 150-152. doi:http://doi.org/10.1038/506150a

Wasserstein, R. L., & Lazar, N. A. (2016). The ASA‘s Statement on p-Values: Context, Process, and Purpose. The American Statistician, 70(2), 129-133. doi:10.1080/00031305.2016.1154108