【Rで統計】ANOVAでは帰無仮説が正しくても有意差は出現する：p値だけでなく効果量を見る

統計検定を使う全員に知っておいてほしいこと、それは帰無仮説が正しい、あるいは、群間に差がないとき、p値は一様分布に従うということです。これは直感に反するかもしれませんが、非常に重要な性質です。つまり、本当はまったく差がなくても、p値が0.999になることも、0.001になることも同様にあり得るのです。
　この事情は当然、ANOVAでも同じです。したがって、ANOVAの結果を見て、次のような思い込みの判断が生じるかもしれません。

p 値が 0.05 以下になったので、群の平均は本当に違っているに違いない

しかし、これは 統計的には正しくありません。ANOVAでも、群の真の平均が本当にすべて等しい場合（帰無仮説が正しい場合）であっても、p 値が 0.05 以下になることは確率的にあり得ます。これは、解析ミスでも、R のバグでもありません。検定というの本質的にそういうものなのです。

p 値以外の何を見ればよいのか：効果量を見る
効果量とは何か：差の「大きさ」を測る指標
- イータ二乗の定義
数値実験で確認
- ■ p値の分布
- ■ p値と効果量の分布
まとめ：追試はとても大切
【付録】数値実験を行うRスクリプト

p 値以外の何を見ればよいのか：効果量を見る

　ここまで読むと、次のような疑問を持つ人もいるかもしれません。

サンプル数を増やせば、 p値の「信頼度」も上がるのではないか？

たしかに素朴な感覚では、「例数が少ないから偶然 p < 0.05 になっただけで、データをたくさん集めれば、もっと正しい判断ができるはずだ」と考えたくなります。しかし、この期待は帰無仮説が正しい場合に限っては、p 値に関して成り立ちません。サンプル数を増やしても、帰無仮説が正しい限り、p 値の分布の性質は変わらないのです。
　ただし、群間の真の平均値に差が存在する場合には状況は異なります。その場合、サンプル数を増やすことで検定の検出力が高まり、p 値は平均的に小さくなる傾向を示します。本記事で扱っているのは、あくまで「帰無仮説が正しい場合」における p 値の性質です。
　前節で述べたとおり「帰無仮説が正しいとき、p値は 0〜1 の一様分布に従う」という性質は、サンプルサイズに依存しません。各群の例数を、「10 にしても」「100 にしても」「1000 にしても」p値の分布は変わらず、一様分布のままです。
　つまり、データを10倍多く集めても、100倍多く集めても、p値が「より信頼できる指標」になることはありません。

　そのように言われると、

それなら、例数を増やすのは無意味ではないか？

と文句を言いたくなるかもしれません。しかし、それは誤りです。サンプル数を増やせば、違いの評価の精度は確実に向上します。ただし、その「評価」は p値ではなく、効果量によって行われるべきなのです。

効果量とは何か：差の「大きさ」を測る指標

　効果量（effect size）は、「群間差が存在するかどうかではなく、群間差がどれくらい大きいのか」を表す指標です。
　ANOVA では、効果量としていくつかの選択肢がありますが、本記事では最も基本的で直感的な指標であるイータ二乗（η²）を用います。

イータ二乗 $\eta^ 2$ の定義

　イータ二乗は、ANOVA において次のように定義されます。

$\displaystyle{ \eta^2 = \frac{\text{群間平方和（SS}_{\text{between}})}{\text{全平方和（SS}_{\text{total}})} }$

これは、

データ全体のばらつきのうち、群平均の違いで説明できる割合

を表しています。 $\eta^ 2$ の解釈は、

$\eta^ 2 = 0$ → 群の違いはまったく意味を持たない
$\eta^ 2 > 0$ が大きい → 群間差が全体のばらつきに大きく寄与している

という、解釈が非常に明確な指標です。

数値実験で確認

　ここでは、理論的に数式を使って説明するのではなく、Rを使って数値実験をした結果をしめします。数値実験の設計（ANOVA）は以下です。

群数：3 群
各群の分布：完全に同一
平均も分散もすべて同じ

つまり、群間に差はないという帰無仮説は常に正しいとします。
　この条件のもとで、各群の例数 $n _ A$ 、 $n _ B$ 、 $n _ C$ を

$n _ A = n _ B = n _ C =10$
$n _ A = n _ B = n _ C =100$
$n _ A = n _ B = n _ C =1000$

と変化させながら、p値の分布と、 $\eta^ 2$ の分布を数値実験で調べます。

■ p値の分布

　私が行った数値実験の結果を下図に示します。下図では、乱数を用いてサンプルデータを生成し、そのデータに対して ANOVA による検定を行う数値実験を 10,000 回 繰り返しました。理論的に帰無仮説が正しい場合、p 値の分布は一様分布に従うことが知られていますが、本数値実験で得られた 10,000 個の p 値のヒストグラムも、ほぼ一様分布となっていることが確認できます。
　この結果からも分かるように、帰無仮説が正しい場合であっても、p 値が特定の値（たとえば 1 に近い値）に偏って現れることはありません。すなわち、p 値は「差がないこと」を示す指標ではなく、帰無仮説のもとで 0 から 1 の間を等しく取りうる量であることが、数値的にも確かめられます。

帰無仮説が正しい場合のANOVAのp値の分布。A、B、Cの３群それぞれの例数 $N _ {\rm A}$ 、 $N _ {\rm B}$ 、 $N _ {\rm C}$ を、10（左）、100（中）、1000（右）として、10000回検定を繰り返した。

■ p値と効果量 $\eta^ 2$ の分布

　ここでは、p 値と効果量 $\eta^ 2$ の関係を、数値実験の結果も示しながら説明します。ポイントは、帰無仮説が正しい場合でも、有意と判定されるデータは必ず一定数現れるという事実です。
　帰無仮説が正しいということは、「本当は群の間に差がない」という状況です。それでも、データは必ずばらつきを持っているため、偶然によって群平均が少し離れて見えることがあります。そのようなデータでは、ANOVA の結果として p 値が小さくなり、「有意」と判定されます。
　このときに重要なのは、有意と判定されたデータだけを見ると、効果量 $\eta^ 2$ が相対的に大きく見えるという点です。これは、「差があったから有意になった」のではなく、「たまたま差が大きく見えたデータが、有意として拾われている」ためです。
　下図は、各群の例数を $n _ A = n _ B = n _ C = 10$ とした数値実験の結果を示しています。この実験では帰無仮説が正しく、群間の真の平均値に差はありません。下図左の赤色のビン（棒）は、 $p \lt 0.05$ となった例を表しています。これらの例に対応する効果量 $\eta^ 2$ の分布が、下図右の赤色の部分です。
　見かけ上は群間の平均値に差があるように見えますが、この程度の差は偶然によって十分に生じうるものです。

p値分布（左）と効果量 $\eta^ 2$ 分布（右）。A、B、Cの３群それぞれの例数が $N _ {\rm A} = 10$ 、 $N _ {\rm B} = 10$ 、 $N _ {\rm C} = 10$ で、帰無仮説が正しい（群間に差がない）場合の数値実験。

　次に、サンプル数の影響を考えます。下図では、各群の例数を 100、1000 に増やした場合の結果を順に示しています。

p値分布（左）と効果量 $\eta^ 2$ 分布（右）。A、B、Cの３群それぞれの例数が $N _ {\rm A} = 100$ 、 $N _ {\rm B} = 100$ 、 $N _ {\rm C} = 100$ で、帰無仮説が正しい（群間に差がない）場合の数値実験。

p値分布（左）と効果量 $\eta^ 2$ 分布（右）。A、B、Cの３群それぞれの例数が $N _ {\rm A} = 1000$ 、 $N _ {\rm B} = 1000$ 、 $N _ {\rm C} = 1000$ で、帰無仮説が正しい（群間に差がない）場合の数値実験。

　各群の例数を 10 倍、さらに 100 倍と増やしていくと、平均値はより安定して推定されるようになります。上図からも分かるように、 $n_A = n_B = n_C = 100$ 、 $n_A = n_B = n_C = 1000$ と例数が増えるにつれて、偶然による群間のずれは急速に小さくなります。その結果、有意と判定されたデータに対応する効果量 $\eta^ 2$ は、例数の増加とともに著しく小さくなり、ほぼ 0 に集中していきます。各群の例数を10から、10倍（100）、100倍（1000）に増やすと、 $\eta^ 2$ の値は、それぞれ、1/10、1/100になっています。
　これは、「本当は群間に差がない」という事実が、サンプル数を増やすことで、より高い精度で数値として確認できるようになることを意味します。このように、サンプル数を増やすことの本質的な役割は、p 値の振る舞いを変えることではなく、効果量を通して差がどれほど小さいかを明確に示すことにあります。

まとめ：追試はとても大切

　本記事では、ANOVA を例に、p 値と効果量の意味と限界について整理しました。
　重要なポイントは、帰無仮説が正しい場合でも、p 値は 0 から 1 の間で一様に分布し、一定の確率で小さな値を取るという事実です。これは ANOVA に限らず、統計的検定一般に共通する性質です。そのため、単一の実験で $p \lt 0.05$ が得られたとしても、それだけで「群間に本当に差がある」と結論づけることはできません。
　そして、サンプル数を増やすことには明確な意味があります。ただし、その効果は p 値に現れるのではなく、効果量に現れます。本記事で用いた効果量 $\eta^ 2$ は「データ全体のばらつきのうち、群間差で説明できる割合」を表す指標です。帰無仮説が正しい場合、例数を増やすにつれて、偶然による見かけの群間差は抑えられ、 $\eta^ 2$ は 0 に近づいていきます。これは「差がない」という構造そのものが、数値としてより明確に表れることを意味します。
　つまり、

p 値は「偶然にどれくらい極端な結果が出たか」を示す指標であり
効果量は「その差がどれくらい大きいのか（あるいは小さいのか）」を示す指標です

という役割の違いがあります。サンプル数を増やす目的は、p 値を安定させることではなく、効果量を通して差の大きさを高精度に評価することにあります。

　ただし、私は、p 値が無意味とまでは言いません。独立した追試が行われ、同様の実験で繰り返し小さな p 値が得られる場合、その確率は積として計算され、急速に小さくなります。このとき、「偶然による結果」である可能性は確率的に極めて低くなり、結果の再現性が強く裏づけられます。
　したがって、統計解析において重要なのは、

単発の p 値に過度な意味を与えないこと
効果量を用いて、差の大きさを定量的に評価すること
独立した追試による再現性を重視すること

です。
　p 値は「有意か否か」あるいは「差があるかどうか」を、一度の実験結果だけから機械的に判断するための指標ではありません。効果量による差の大きさの評価や、独立した追試による再現性の確認と併せて解釈されて初めて、統計的な情報として意味を持ちます。ANOVA の結果を解釈する際には、この点を常に意識することが重要です。
　今回は「真の差がある」場合について示しませんでしたが、下に一例のみ数値実験結果を示しておきます。近いうちに「真の差がある」場合のp値分布の特徴と、多重検定におけるp値の調整について説明したいと思います。

p値分布（左）と効果量 $\eta^ 2$ 分布（右）。A、B、Cの３群それぞれの例数が $N _ {\rm A} = 10$ 、 $N _ {\rm B} = 10$ 、 $N _ {\rm C} = 10$ で、各群の**真の平均**が $\mu _ {\rm A} = 0$ 、 $\mu _ {\rm B} = 0$ 、 $\mu _ {\rm C} = 0.5$ とCにのみ差がある場合。

【付録】数値実験を行うRスクリプト

■ サンプルデータの生成

　以下のスクリプトをRで実行してから、以下の実行例を参照してください。

# ------------------------------------------------------------
# Generate sample data for a 3-group ANOVA demonstration
#
# Inputs:
#   n_vec     : named integer vector, sample sizes per group (e.g., c(A=10,B=10,C=10))
#   mu_vec    : named numeric vector, true means per group (same names as n_vec)
#   var_total : target total variance (between + within)
#   seed      : random seed (NULL to skip)
#
# Output:
#   data.frame with columns: group (factor), y (numeric)
# ------------------------------------------------------------
generate_3group_data <- function(
  n_vec,
  mu_vec,
  var_total,
  seed = 1
) {
  # ---- 0) checks ----
  if (!is.null(seed)) set.seed(seed)

  if (is.null(names(n_vec)) || any(names(n_vec) == "")) {
    stop("n_vec must be a *named* vector (e.g., c(A=10,B=10,C=10)).")
  }
  grp_levels <- names(n_vec)

  if (is.null(names(mu_vec)) || any(names(mu_vec) == "")) {
    stop("mu_vec must be a *named* vector with the same names as n_vec.")
  }

  # Align order by group names
  n_vec  <- n_vec[grp_levels]
  mu_vec <- mu_vec[grp_levels]

  if (length(n_vec) != length(mu_vec)) {
    stop("n_vec and mu_vec must have the same length.")
  }
  if (any(!is.finite(mu_vec))) stop("mu_vec contains non-finite values.")
  if (!is.finite(var_total) || var_total <= 0) stop("var_total must be positive.")

  N <- sum(n_vec)
  if (N <= 0) stop("Total sample size must be positive.")
  if (any(n_vec <= 0)) stop("All group sample sizes must be positive.")

  # ---- 1) decompose variance: var_total = var_between + var_within ----
  mu_bar <- sum(n_vec * mu_vec) / N
  var_between <- sum(n_vec * (mu_vec - mu_bar)^2) / N

  var_within <- var_total - var_between
  if (var_within <= 0) {
    stop(sprintf(
      "Invalid settings: var_between (%.4g) >= var_total (%.4g).\nReduce mean differences (mu_vec) or increase var_total.",
      var_between, var_total
    ))
  }
  sd_within <- sqrt(var_within)

  # ---- 2) generate data ----
  group <- factor(rep(grp_levels, times = n_vec), levels = grp_levels)
  y <- rnorm(N, mean = rep(mu_vec, times = n_vec), sd = sd_within)

  data.frame(group = group, y = y)
}

■ 実行例（帰無仮説が正しい場合）

　以下は、3 群の真の平均値がすべて等しいという条件、すなわち 帰無仮説が正しい場合の数値実験の例です。各群のデータは同一の分布から生成されており、理論的には群間に差は存在しません。

dat <- generate_3group_data(
  n_vec     = c(A = 10, B = 10, C = 10),
  mu_vec    = c(A = 0,  B = 0,  C = 0),
  var_total = 1.0,
  seed      = 1
)

このコードでは、各群の例数を 10 とし、群間の真の平均値はすべて 0 に設定しています。したがって、生成されるデータにおいて観測される群間差は、すべて 偶然によるばらつきに起因するものです。

■ 1万回数値実験してp値と効果量の分布を描く

　今回示したような図を作成するRスクリプトの例を示しておきます。

############################################################
# Simulation of p-values and effect sizes (eta^2) in ANOVA
#
# Scenario:
#   - One-way ANOVA with 3 groups
#   - Null hypothesis is true (all true means are equal)
#
# What is shown:
#   (1) Distribution of p-values (theoretical: Uniform[0,1])
#   (2) Distribution of effect size eta^2
#       - gray: all simulations
#       - red : simulations with p < 0.05
#
# Design choices for publication:
#   - Minimal, clear comments
#   - Effect-size histograms drawn with rect() so that
#     bar width is slightly shrunk, creating visible gaps
############################################################

#-----------------------------------------------------------
# Prerequisite:
#   generate_3group_data(n_vec, mu_vec, var_total, seed)
#   (defined in the previous section)
#-----------------------------------------------------------

#-----------------------------------------------------------
# 1) Run ANOVA simulations and collect p-values and eta^2
#-----------------------------------------------------------
simulate_anova_p_eta2 <- function(n_grp, n_sim = 10000,
                                  var_total = 1.0, seed = 1) {
  set.seed(seed)

  pvals <- numeric(n_sim)
  eta2  <- numeric(n_sim)

  for (i in seq_len(n_sim)) {

    # Data generation under the null hypothesis
    dat <- generate_3group_data(
      n_vec     = c(A = n_grp, B = n_grp, C = n_grp),
      mu_vec    = c(A = 0,     B = 0,     C = 0),
      var_total = var_total,
      seed      = NULL
    )

    fit <- aov(y ~ group, data = dat)
    tab <- summary(fit)[[1]]

    # p-value
    pvals[i] <- tab["group", "Pr(>F)"]

    # Effect size eta^2
    SS_between <- tab["group",     "Sum Sq"]
    SS_within  <- tab["Residuals", "Sum Sq"]
    eta2[i] <- SS_between / (SS_between + SS_within)
  }

  list(pvals = pvals, eta2 = eta2)
}

#-----------------------------------------------------------
# 2) Histogram drawing with a small horizontal gap
#-----------------------------------------------------------
plot_hist_rect_density <- function(h, col, border = NA, gap = 0.92) {
  # h   : output of hist(..., plot = FALSE)
  # gap : < 1 shrinks bar width to create spacing

  bw   <- diff(h$breaks)
  dens <- h$counts / (sum(h$counts) * bw)

  x_l  <- h$breaks[-length(h$breaks)]
  x_r  <- h$breaks[-1]
  x_m  <- (x_l + x_r) / 2

  half_w <- (x_r - x_l) * gap / 2

  rect(x_m - half_w, 0, x_m + half_w, dens,
       col = col, border = border)

  invisible(dens)
}

#-----------------------------------------------------------
# 3) Example run (change n_grp to 100 or 1000 if desired)
#-----------------------------------------------------------
n_grp <- 10
n_sim <- 10000

res   <- simulate_anova_p_eta2(n_grp = n_grp, n_sim = n_sim, seed = 1)
pvals <- res$pvals
eta2  <- res$eta2

# Subset with p < 0.05
eta2_sig <- eta2[pvals < 0.05]
ratio    <- length(eta2_sig) / length(eta2)  # area ratio

#-----------------------------------------------------------
# 4) Visualization (two panels)
#-----------------------------------------------------------
par(mfrow = c(1, 2), mar = c(5, 5, 3, 1))

# ---- (1) p-value distribution ----
hist(
  pvals,
  breaks = seq(0, 1, by = 0.05),
  probability = TRUE,
  col = "gray80",
  border = "white",
  main = sprintf("Distribution of p-values\n(null is true, n = %d per group)", n_grp),
  xlab = "p-value",
  ylab = "Density"
)
abline(h = 1, lty = 2, lwd = 2, col = "lightpink")  # Uniform[0,1]

# ---- (2) Effect size eta^2 ----
nbins <- max(round(n_sim / 5000) * 10, 20)

h_all <- hist(eta2,     breaks = nbins, plot = FALSE)
h_sig <- hist(eta2_sig, breaks = h_all$breaks, plot = FALSE)

plot(
  0, 0, type = "n",
  xlim = range(h_all$breaks),
  ylim = c(
    0,
    max(h_all$counts / (sum(h_all$counts) *
        diff(h_all$breaks)[1])) * 1.05
  ),
  xlab = expression(eta^2),
  ylab = "Density",
  main = expression("Effect size " * eta^2 *
                    " (gray: all, red: p < 0.05)")
)

# All simulations (gray)
plot_hist_rect_density(h_all, col = "gray80", gap = 0.92)

# p < 0.05 subset (red), scaled by area ratio
bw <- diff(h_sig$breaks)
dens_sig <- ratio * (h_sig$counts / (sum(h_sig$counts) * bw))

x_l <- h_sig$breaks[-length(h_sig$breaks)]
x_r <- h_sig$breaks[-1]
x_m <- (x_l + x_r) / 2
half_w <- (x_r - x_l) * 0.92 / 2

rect(x_m - half_w, 0, x_m + half_w, dens_sig,
     col = 2, border = NA)

legend("topright",
       legend = c("all", "p < 0.05"),
       fill = c("gray80", 2),
       border = NA,
       bty = "n")

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。