UPDATE: 2023-02-16 11:59:37

はじめに

p値とは「特定の統計モデルのもとで(帰無仮説が真であると仮定することも含む)、観察されたデータの統計的要約が観察された値と同じか、それよりも極端である場合の確率」と説明される。p値が利用されるのは、帰無仮説有意検定(NHST)を行うときであり、p値の解釈については誤用が絶えない。これは自分にも当てはまることなので、強い自戒の意味を込めてp値についておさらいする。

検定とp値に関する誤解は絶えず、下記の通り、ASA(American Statistical Association)が声明を出していたりする。京都大学の佐藤先生が日本語版でまとめてくださっているので、詳細は下記を参考に願います。

ここでの目標は、平均値の差の検定を例にp値の計算過程を明らかにし、可視化するところまで。今後は、個人的には最近よく聞くようになった「p値関数」についてまとめたい。

サンプルデータ

下記のような試験のスコアが観測されたとする。ここでは、等分散の正規母集団からのサンプルデータとする。

両グループはある程度、重なっている部分もあるが、赤色のグループAの方が平均値は高そうではある。

set.seed(1989)
A <- floor(rnorm(10, 90, 1))
B <- floor(rnorm(10, 89, 1))
#A <- c(9,7,8,9,8,9,9,10,9,9) 
#B <- c(9,6,7,8,7,9,8,8,8,7)
A_density <- density(A)
B_density <- density(B)

plot(
  range(A_density$x, B_density$x),
  range(A_density$y, B_density$y),
  type = "n",
  xlab = "Score",
  ylab = "Density", 
  xaxt = "n")
axis(1, at = seq(min(A, B)-10, max(A, B)+10, by = 1))

polygon(A_density, col = adjustcolor("#E55A64", alpha.f = 0.5))
polygon(B_density, col = adjustcolor("#7FA6C9", alpha.f = 0.5))
title("Score Density")

summary(A);
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    88.0    89.0    89.5    89.6    90.0    91.0
summary(B)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   86.00   87.25   89.00   88.60   90.00   90.00

仮説検定

天下り的ではあるが、このような2グループの平均値を比較する場合、t検定を行うことになる。つまり、各グループの平均、分散(標準偏差)、サンプルサイズを利用して、t分布(自由度は18)に従うt値を計算することで、2グループ間のスコアの差が0より大きい(小さい)確率を調べる。この過程でp値も計算することができる。

仮説検定では、特定の統計モデルのもとで(帰無仮説が真であると仮定することも含む)、観察されたデータの統計的要約が観察された値と同じか、それよりも極端である場合の確率を計算するが、それがp値である。

ちなみに、ここで利用するt分布(自由度は18)は下記のような形状をしている。

df <- (length(A) + length(B) - 2)
x <- seq(-5, 5, length.out = 100)
y <- dt(x, df = df)
plot(x,
     y,
     type = "l", 
  xaxt = "n")
axis(1, at = seq(min(x), max(x), by = .5))
title("t Distribution")
legend("topright", 
       paste0("t(df=", df, ")"),
       lty = 1)

t検定

計算式に従って、計算してもt値は計算できるが、t.test()でt検定は簡単に実行でき、t値を得られる。

ttest_res <- t.test(A, B, var.equal = TRUE)
ttest_res
## 
##  Two Sample t-test
## 
## data:  A and B
## t = 1.7678, df = 18, p-value = 0.09405
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.188461  2.188461
## sample estimates:
## mean of x mean of y 
##      89.6      88.6

t値は1.767767で、p値は0.094047である。t値はt統計量を計算する計算式から算出され、p値は検定に使用しているt分布から計算できる。

そもそもp値は、「特定の統計モデルのもとで(帰無仮説が真であると仮定することも含む)、観察されたデータの統計的要約が観察された値と同じか、それよりも極端である場合の確率」だった。ここでは、帰無仮説のスコアの差が0であるという仮定は、先ほど紹介した0を中心とするt分布によって表現されている。

つまり、このt分布からp値を計算でき、さきほどt検定によって得られたt値に対応するt分布の両側の赤色部分の確率がp値に対応しているためである。

plot(x,
     y,
     type = "l", 
  xaxt = "n")
axis(1, at = seq(min(x), max(x), by = .5))
title("t Distribution")
abline(v = ttest_res$statistic, lwd = 2, col = "red", lty = 3)
abline(v = -1*ttest_res$statistic, lwd = 2, col = "red", lty = 3)

xx_upr <- seq(ttest_res$statistic, 10, length=100)
xx_lwr <- seq(-1*ttest_res$statistic, -10, length=100)
z <- dt(xx_upr, df = df)
polygon(c(ttest_res$statistic, xx_upr, 8),
        c(0, z, 0),
        col = adjustcolor("#E55A64", alpha.f = 0.5))
polygon(c(-1*ttest_res$statistic, xx_lwr ,8),
        c(0, z, 0),
        col = adjustcolor("#E55A64", alpha.f = 0.5))

下側t値に対応する外側確率はpt()によって得られ、

prob_lwr <- pt(-1 * ttest_res$statistic, df = df)
prob_lwr
##          t 
## 0.04702348

上側t値に対応する外側確率はpt()によって得られた確率から1を引くことで計算できる。

prob_upr <- 1 - pt(ttest_res$statistic, df = df)
prob_upr
##          t 
## 0.04702348

これらの確率を足し合わせると、t検定で得られたp値の確率0.094047が計算できる。

prob_lwr + prob_upr
##          t 
## 0.09404695

帰無仮説のt分布がわかっているので、積分してももちろん同じ結果が得られる。

integral_lwr <- integrate(f = dt, df = df, lower = -Inf, upper = -1 * ttest_res$statistic)$value
integral_upr <- 1 - integrate(f = dt, df = df, lower = -Inf, upper = ttest_res$statistic)$value
sprintf("p-value: %f", (integral_lwr + integral_upr))
## [1] "p-value: 0.094047"

有意水準を慣習的な5％とした仮説検定であれば、対立仮説を棄却できない。この際に、グループA,Bのスコアには差は0であるという帰無仮説を採択するという表現をみかけるが、この帰無仮説を採択するという表現は良くない。仮説検定は「帰無仮説を棄却できるか否か」なので、帰無仮説を「棄却できる」「棄却できない」のいずれかであって、採択することはないためである。

これを可視化すると下記のような形となる。黒線は赤色の部分が合計で5％となるt値である。黒線よりも赤線が内側にある、つまりp値は5%以上あることがわかる。実際、10%近くあるため、観察されたデータが少なくとも同じか、極端な値ではないということになる。

plot(x,
     y,
     type = "l", 
  xaxt = "n")
axis(1, at = seq(min(x), max(x), by = .5))
title("t Distribution")
abline(v = ttest_res$statistic, lwd = 2, col = "red", lty = 3)
abline(v = -1*ttest_res$statistic, lwd = 2, col = "red", lty = 3)

xx_upr <- seq(ttest_res$statistic, 10, length=100)
xx_lwr <- seq(-1*ttest_res$statistic, -10, length=100)
z <- dt(xx_upr, df = df)
polygon(c(ttest_res$statistic, xx_upr, 8),
        c(0, z, 0),
        col = adjustcolor("#E55A64", alpha.f = 0.5))
polygon(c(-1*ttest_res$statistic, xx_lwr ,8),
        c(0, z, 0),
        col = adjustcolor("#E55A64", alpha.f = 0.5))

abline(v = qt(0.025, df = df), lwd = 2, col = "black", lty = 3)
abline(v = qt(0.975, df = df), lwd = 2, col = "black", lty = 3)

まとめ

p値は仮説検定において、「特定の統計モデルのもとで、帰無仮説が真であると仮定した場合に、観察されたデータが少なくとも同じか、もっと極端である場合の確率」である。 仮説検定では、帰無仮説が真であること仮定しているため、帰無仮説に応じた確率分布が必要になる。例えば、t検定の場合はt分布によってスコアの差が0であることを表現する。そこから仮説検定を行って得られた統計量の値を利用し、帰無仮説に利用している確率分布から、観測データから得られた統計量の値よりも極端な確率を計算する。これがp値に対応している。