p値は何ではないのか？帰無分布をPythonで作り、シミュレーションで理解する

統計解析の結果として、しばしば次のような値を目にします。

$p=0.04$

このとき、次のように読みたくなるかもしれません。

• 「帰無仮説が正しい確率は4%」
• 「96%の確率で差がある」
• 「偶然だけでこの結果が出た確率は4%」

しかし、いずれもp値の定義ではありません。

この記事では、架空の測定データを使った1標本t検定を題材に、帰無仮説のもとで実験を何度も繰り返すシミュレーションを行います。帰無分布を自分で作り、p値が図のどの面積に対応するのかを確認します。

さらに、次の2点もシミュレーションで確かめます。

• 帰無仮説が正しくても、$p<0.05$ は一定の割合で起こる
• 同じ効果の大きさでも、標本サイズによってp値の分布は大きく変わる

コードはGoogle Colabで上から順に実行できます。

TL;DR

• 今回の両側検定におけるp値は、$H_0$を仮定したときに、観測値と同じかそれ以上に極端な検定統計量が得られる確率である
• p値は $P(H_0\mid\mathrm{data})$、つまり「データを見た後に帰無仮説が正しい確率」ではない
• p値は効果量でも、結果の重要性でもない
• $p\geq0.05$ は「差がないことの証明」ではない
• p値だけでなく、効果量、信頼区間、標本サイズ、研究設計を合わせて読む必要がある

この記事で扱う問い

架空の測定法に、基準となる母平均 $50$ があるとします。新しい条件で独立に20回測定し、母平均が基準値と異なるかを調べます。

今回の仮説は次のとおりです。

$$
H_0:\mu=50
$$

$$
H_1:\mu\neq50
$$

• $H_0$：帰無仮説。母平均は基準値50と等しい
• $H_1$：対立仮説。母平均は基準値50と異なる

「大きいか」だけでなく「異なるか」を問うため、左右両方のずれを考える両側検定を使います。

最初に用語を整理する

用語	この記事での意味
母集団	新しい条件で得られうる、すべての測定値
標本	実際に観測した20個の測定値
母平均 $\mu$	母集団全体の平均。実際の研究では通常未知
標本平均 $\bar{x}$	観測した20個の平均
帰無仮説 $H_0$	今回は「母平均が50である」という仮定
対立仮説 $H_1$	今回は「母平均が50ではない」という仮定
検定統計量	データと $H_0$ のずれを、ばらつきも考慮して数値化したもの
帰無分布	$H_0$ が成り立つと仮定したときの検定統計量の分布
p値	帰無分布上で、観測値と同じかそれ以上に極端な領域の確率
有意水準 $\alpha$	検定前に決める判断基準。例として $\alpha=0.05$ を使う

1標本t検定は何を比べているのか

標本平均と帰無仮説の平均との差だけを見ても、その差が大きいかどうかは判断できません。同じ平均差でも、データのばらつきが小さければ目立つ差になり、ばらつきが大きければ珍しくない差になるためです。

そこで1標本t検定では、標本平均と帰無仮説の平均との差を、標準誤差で割ります。

$$
SE=\frac{s}{\sqrt{n}}
$$

$$
t=\frac{\bar{x}-\mu_0}{SE}
=\frac{\bar{x}-\mu_0}{s/\sqrt{n}}
$$

ここで、

• $\bar{x}$：標本平均
• $\mu_0$：帰無仮説で仮定する母平均
• $s$：標本標準偏差
• $n$：標本サイズ
• $SE$：標本平均の標準誤差

です。

$t$は、「標本平均が帰無仮説の値から、標準誤差何個分ずれているか」を表します。

ダミーデータを作り、1標本t検定を行う

教育用のダミーデータなので、生成時の母平均を54、母標準偏差を10に設定します。ただし、実際の研究では真の母平均は分かりません。解析者が利用できるのは、観測された標本と、検討したい帰無仮説です。

from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 図の保存先
OUTPUT_DIR = Path("pvalue_figures")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

# 検定の設定
MU0 = 50.0       # 帰無仮説で仮定する母平均
TRUE_MU = 54.0   # ダミーデータ生成時の母平均
SIGMA = 10.0     # ダミーデータ生成時の母標準偏差
N = 20           # 標本サイズ
ALPHA = 0.05     # 有意水準

# 乱数シードを固定して再現可能にする
rng_observed = np.random.default_rng(28)

# ダミーデータを生成
observed = rng_observed.normal(
    loc=TRUE_MU,
    scale=SIGMA,
    size=N,
)

# 両側の1標本t検定
result = stats.ttest_1samp(
    observed,
    popmean=MU0,
    alternative="two-sided",
)

# 記述統計量と補助的な指標
sample_mean = observed.mean()
sample_sd = observed.std(ddof=1)
standard_error = sample_sd / np.sqrt(N)
cohens_d = (sample_mean - MU0) / sample_sd

# 母平均の95%信頼区間
t_critical = stats.t.ppf(1 - ALPHA / 2, df=N - 1)
ci_low = sample_mean - t_critical * standard_error
ci_high = sample_mean + t_critical * standard_error

print("=== Observed dummy data ===")
print(f"n                  = {N}")
print(f"sample mean        = {sample_mean:.3f}")
print(f"sample SD          = {sample_sd:.3f}")
print(f"standard error     = {standard_error:.3f}")
print(f"t statistic        = {result.statistic:.3f}")
print(f"p value            = {result.pvalue:.4f}")
print(f"95% CI of the mean = [{ci_low:.3f}, {ci_high:.3f}]")
print(f"Cohen's d          = {cohens_d:.3f}")

実行結果は次のようになります。

=== Observed dummy data ===
n                  = 20
sample mean        = 54.655
sample SD          = 9.433
standard error     = 2.109
t statistic        = 2.207
p value            = 0.0398
95% CI of the mean = [50.240, 59.070]
Cohen's d          = 0.493

標本平均は54.655で、帰無仮説の50より4.655大きくなりました。標準誤差は2.109なので、

$$
t=\frac{54.655-50}{2.109}\approx2.207
$$

です。標本平均は、帰無仮説の値から約2.21標準誤差だけ離れています。

観測データを図にする

fig, ax = plt.subplots(figsize=(9, 4.8))

x = np.arange(1, N + 1)

ax.scatter(x, observed, s=55, label="Observed values")
ax.axhline(
    MU0,
    linestyle="--",
    linewidth=2,
    label=f"Null mean = {MU0:.0f}",
)
ax.axhline(
    sample_mean,
    linewidth=2,
    label=f"Sample mean = {sample_mean:.2f}",
)

ax.set(
    xlabel="Observation number",
    ylabel="Measured value",
    title="Observed dummy data and the null-hypothesis mean",
)
ax.set_xticks(np.arange(1, N + 1, 2))
ax.legend()
ax.grid(alpha=0.25)

fig.tight_layout()
fig.savefig(
    OUTPUT_DIR / "fig1_observed_data.png",
    dpi=180,
    bbox_inches="tight",
)
plt.show()

図1. 観測した20個のダミーデータ、帰無仮説の平均50、標本平均54.65。

個々の値はかなりばらついています。標本平均が50を上回っていることだけではなく、このばらつきに対して平均差がどれくらい大きいかを$t$統計量で評価します。

帰無分布をシミュレーションで作る

p値を理解するために、いったん「$H_0:\mu=50$ が成り立っている」と仮定します。

その世界で、同じ標本サイズ20の実験を10万回繰り返します。各実験について$t$統計量を計算すると、$H_0$のもとで起こりうる$t$の分布、つまり帰無分布が得られます。

N_SIM_NULL = 100_000
rng_null = np.random.default_rng(202604)

# H0: 母平均50のもとで、標本サイズ20の実験を10万回生成
null_data = rng_null.normal(
    loc=MU0,
    scale=SIGMA,
    size=(N_SIM_NULL, N),
)

# 各実験のt統計量を計算
null_means = null_data.mean(axis=1)
null_sds = null_data.std(axis=1, ddof=1)
t_null = (null_means - MU0) / (null_sds / np.sqrt(N))

# 観測された|t|以上に極端な割合
simulated_p = np.mean(
    np.abs(t_null) >= abs(result.statistic)
)

print("=== p value from the null distribution ===")
print(f"SciPy p value      = {result.pvalue:.5f}")
print(f"Simulated p value  = {simulated_p:.5f}")

# 帰無分布を描画
fig, ax = plt.subplots(figsize=(9, 5.2))

bins = np.linspace(-5, 5, 121)
ax.hist(
    t_null,
    bins=bins,
    density=True,
    color="0.80",
    edgecolor="white",
    label="Simulated null distribution",
)

# 理論的なt分布
x_grid = np.linspace(-5, 5, 1000)
t_density = stats.t.pdf(x_grid, df=N - 1)

ax.plot(
    x_grid,
    t_density,
    linewidth=2,
    label=f"Theoretical t distribution (df={N - 1})",
)

# 両側検定なので、左右両方の裾を塗る
threshold = abs(result.statistic)
right_tail = x_grid >= threshold
left_tail = x_grid <= -threshold

ax.fill_between(
    x_grid[right_tail],
    0,
    t_density[right_tail],
    alpha=0.45,
    color="tab:red",
    label="As or more extreme than observed",
)
ax.fill_between(
    x_grid[left_tail],
    0,
    t_density[left_tail],
    alpha=0.45,
    color="tab:red",
)

ax.axvline(
    result.statistic,
    linestyle="--",
    linewidth=2,
    color="tab:red",
)
ax.axvline(
    -result.statistic,
    linestyle="--",
    linewidth=2,
    color="tab:red",
)

ax.set(
    xlabel="t statistic",
    ylabel="Density",
    title=(
        "Null distribution and two-sided p value "
        f"(p ≈ {simulated_p:.4f})"
    ),
)
ax.legend()
ax.grid(alpha=0.25)

fig.tight_layout()
fig.savefig(
    OUTPUT_DIR / "fig2_null_distribution.png",
    dpi=180,
    bbox_inches="tight",
)
plt.show()

実行結果は次のようになります。

=== p value from the null distribution ===
SciPy p value      = 0.03982
Simulated p value  = 0.03904

図2. 帰無仮説のもとで作った$t$統計量の分布。赤い左右の裾が、観測値と同じかそれ以上に極端な領域。

シミュレーションで得たp値は0.03904、SciPyが理論的なt分布から計算したp値は0.03982でした。有限回の乱数シミュレーションなので完全には一致しませんが、よく対応しています。

p値の定義

今回の両側1標本t検定では、p値は次のように表せます。

$$
p=P\left(|T|\geq|t_{\mathrm{obs}}|\mid H_0\right)
$$

言葉にすると、次の意味です。

帰無仮説と検定モデルの仮定が成り立つとしたときに、観測された$t$統計量と同じか、それ以上に帰無仮説から離れた$t$統計量が得られる確率。

今回の $p\approx0.0398$ は、

母平均が50で、今回の正規モデルや独立性などの仮定が成り立つ世界で同様の実験を繰り返すと、$|t|\geq2.207$ となる実験が約4%ある

という意味です。

条件の向きを入れ替えてはいけない

p値が表すのは、

$$
P(\mathrm{data\ as\ or\ more\ extreme}\mid H_0)
$$

です。次の確率ではありません。

$$
P(H_0\mid\mathrm{data})
$$

前者は「$H_0$を仮定したときのデータ側の確率」、後者は「データを見た後の仮説側の確率」です。条件付き確率の向きが逆です。

後者を求めるには、仮説に対する事前情報や、対立仮説側のモデルなどが追加で必要です。p値だけからは計算できません。

シミュレーション1：帰無仮説が正しくてもp<0.05は起こる

先ほどは、$H_0$が成り立つ世界で10万回の実験を生成しました。各実験のp値を計算して、その分布を見てみます。

# t分布の右裾確率を使って、各実験の両側p値を計算
p_values_h0 = 2 * stats.t.sf(
    np.abs(t_null),
    df=N - 1,
)

false_positive_rate = np.mean(p_values_h0 < ALPHA)

print("=== Repeating experiments when H0 is true ===")
print(
    f"Fraction with p < {ALPHA:.2f} = "
    f"{false_positive_rate:.4f}"
)

fig, ax = plt.subplots(figsize=(9, 4.8))

p_bins = np.linspace(0, 1, 21)
weights_h0 = np.ones_like(p_values_h0) / p_values_h0.size

ax.hist(
    p_values_h0,
    bins=p_bins,
    weights=weights_h0,
    alpha=0.70,
    edgecolor="black",
)
ax.axhline(
    0.05,
    linestyle="--",
    linewidth=2,
    label="Expected fraction per 0.05-wide bin = 0.05",
)
ax.axvspan(
    0,
    ALPHA,
    alpha=0.25,
    color="tab:red",
    label=f"p < {ALPHA:.2f}",
)

ax.set(
    xlabel="p value",
    ylabel="Fraction of simulations",
    title=(
        "Distribution of p values when H0 is true\n"
        f"p < {ALPHA:.2f}: "
        f"{false_positive_rate * 100:.2f}%"
    ),
)
ax.set_xlim(0, 1)
ax.legend()
ax.grid(alpha=0.25)

fig.tight_layout()
fig.savefig(
    OUTPUT_DIR / "fig3_pvalues_under_h0.png",
    dpi=180,
    bbox_inches="tight",
)
plt.show()

=== Repeating experiments when H0 is true ===
Fraction with p < 0.05 = 0.0495

図3. $H_0$が正しいときのp値の分布。0から1まで、ほぼ一様に現れる。

今回のように検定の仮定が満たされ、連続的な検定統計量を使う場合、$H_0$のもとでp値は0から1までほぼ一様に分布します。

そのため、$H_0$が正しい実験だけを集めても、約5%では $p<0.05$ になります。今回のシミュレーションでは4.95%でした。

これは、

$$
P(p<0.05\mid H_0)\approx0.05
$$

を表しています。しかし、次を表しているわけではありません。

$$
P(H_0\mid p<0.05)=0.05
$$

この2つは別の条件付き確率です。

有意水準を $\alpha=0.05$ にするとは、$H_0$が正しい状況で同じ検定手順を長期的に繰り返したとき、誤って棄却する割合を約5%に制御する、という設計です。この誤りを第一種過誤と呼びます。

シミュレーション2：p値は効果量ではない

次に、母集団の標準化効果量を同じ $d=0.5$ に固定し、標本サイズだけを変えます。

1標本の場合の母集団レベルの標準化効果量を、ここでは次のように置きます。

$$
d=\frac{\mu-\mu_0}{\sigma}
$$

シミュレーションでは $\mu_0=0$、$\sigma=1$、$\mu=0.5$ とするため、どちらの条件でも真の効果量は $d=0.5$ です。

標本サイズ $n=10$ と $n=100$ で、それぞれ5万回の実験を行います。

def simulate_one_sample_pvalues(effect_size, n, n_sim, seed):
    """
    母平均0に対する両側1標本t検定を繰り返す。

    母標準偏差を1にしているため、生成時の母平均が
    そのまま標準化効果量dになる。
    """
    rng = np.random.default_rng(seed)

    data = rng.normal(
        loc=effect_size,
        scale=1.0,
        size=(n_sim, n),
    )

    means = data.mean(axis=1)
    sds = data.std(axis=1, ddof=1)
    t_stats = means / (sds / np.sqrt(n))

    return 2 * stats.t.sf(
        np.abs(t_stats),
        df=n - 1,
    )


EFFECT_SIZE = 0.5
N_SIM_EFFECT = 50_000

p_small = simulate_one_sample_pvalues(
    effect_size=EFFECT_SIZE,
    n=10,
    n_sim=N_SIM_EFFECT,
    seed=202605,
)

p_large = simulate_one_sample_pvalues(
    effect_size=EFFECT_SIZE,
    n=100,
    n_sim=N_SIM_EFFECT,
    seed=202606,
)

sig_small = np.mean(p_small < ALPHA)
sig_large = np.mean(p_large < ALPHA)

print("=== Same true effect size, different sample sizes ===")
print(
    f"d = {EFFECT_SIZE}, n = 10 : "
    f"p < 0.05 in {sig_small * 100:.2f}%"
)
print(
    f"d = {EFFECT_SIZE}, n = 100: "
    f"p < 0.05 in {sig_large * 100:.2f}%"
)
print(
    f"d = {EFFECT_SIZE}, n = 10 : "
    f"p >= 0.05 in {(1 - sig_small) * 100:.2f}%"
)

fig, axes = plt.subplots(
    1,
    2,
    figsize=(11, 4.6),
    sharey=True,
)

p_bins = np.linspace(0, 1, 21)

for ax, p_values, n_value, sig_rate in [
    (axes[0], p_small, 10, sig_small),
    (axes[1], p_large, 100, sig_large),
]:
    weights = np.ones_like(p_values) / p_values.size

    ax.hist(
        p_values,
        bins=p_bins,
        weights=weights,
        alpha=0.70,
        edgecolor="black",
    )
    ax.axvspan(
        0,
        ALPHA,
        alpha=0.25,
        color="tab:red",
    )
    ax.axvline(
        ALPHA,
        linestyle="--",
        linewidth=2,
        color="tab:red",
    )
    ax.set(
        xlabel="p value",
        title=(
            f"Same true effect d={EFFECT_SIZE}, n={n_value}\n"
            f"p < 0.05: {sig_rate * 100:.1f}%"
        ),
    )
    ax.set_xlim(0, 1)
    ax.grid(alpha=0.25)

axes[0].set_ylabel("Fraction of simulations")
fig.suptitle(
    "The p-value distribution changes with sample size",
    y=1.03,
)

fig.tight_layout()
fig.savefig(
    OUTPUT_DIR / "fig4_same_effect_different_n.png",
    dpi=180,
    bbox_inches="tight",
)
plt.show()

=== Same true effect size, different sample sizes ===
d = 0.5, n = 10 : p < 0.05 in 29.24%
d = 0.5, n = 100: p < 0.05 in 99.88%
d = 0.5, n = 10 : p >= 0.05 in 70.76%

図4. 真の標準化効果量を同じ $d=0.5$ にしても、標本サイズによってp値の分布は大きく変わる。

どちらも真の効果量は同じです。それでも、$n=10$ では $p<0.05$ になった実験は29.24%だけでした。一方、$n=100$ では99.88%が $p<0.05$ になりました。

この結果から、少なくとも次の2点が分かります。

1. p値は効果の大きさそのものではない
2. $p\geq0.05$ は、効果が存在しないことの証明ではない

$n=10$ の条件では、真の効果量が $d=0.5$ であっても、70.76%の実験で $p\geq0.05$ でした。標本が小さいと、実在する効果を検出できないことがあります。

特定の効果が存在するときに $p<\alpha$ となる確率を検出力と呼びます。今回の推定では、$n=10$ の検出力は約29%、$n=100$ の検出力は約99.9%です。

p値は「偶然だけで起きた確率」でもない

観測データには通常、個体差、測定誤差、標本抽出による変動など、確率的なばらつきが含まれます。p値は、観測結果の何%が「偶然」で、何%が「本物」かを分解する指標ではありません。

p値が行っているのは、あくまで次の計算です。

1. $H_0$と統計モデルを仮定する
2. その仮定のもとで検定統計量の分布を考える
3. 観測値と同じか、それ以上に極端な領域の確率を求める

したがって、$p=0.04$ を「偶然で起きた確率が4%」と言い換えるのは適切ではありません。

よくある誤解を整理する

読み方	判定	より正確な理解
$p=0.04$ なら $H_0$ が正しい確率は4%	誤り	p値は $P(H_0\mid\mathrm{data})$ ではない
$p=0.04$ なら96%の確率で差がある	誤り	対立仮説が正しい確率もp値だけでは分からない
$p=0.04$ なら偶然だけで起きた確率は4%	誤り	$H_0$のもとで同等以上に極端な統計量が出る確率
$p<0.05$ なら効果が大きい	誤り	小さな効果でも標本が大きいとp値は小さくなりうる
$p\geq0.05$ なら差はない	誤り	検出力不足や大きな不確かさの可能性がある
$p=0.04$ なら追試も96%成功する	誤り	再現確率は効果量、標本サイズ、ばらつき、設計などに依存する

米国統計学会の声明でも、p値は仮説が正しい確率や効果の大きさを測るものではなく、科学的結論を単一の閾値だけで決めるべきではないと整理されています。[1, 2]

p<0.05のとき、何と言えばよいのか

今回、あらかじめ有意水準を $\alpha=0.05$ と決めていたなら、

$$
p=0.0398<0.05
$$

なので、検定手順上は $H_0$ を棄却します。

ただし、より丁寧には次のように表現できます。

観測データは、母平均が50であるという帰無仮説と比較的整合しにくかった。

「帰無仮説が偽であると証明した」「差がある確率は96%」とまでは言えません。

反対に $p\geq0.05$ だった場合は、通常「帰無仮説を棄却できなかった」と表現します。「帰無仮説を証明した」「差がないと確定した」ではありません。

また、0.049と0.051の間に、科学的に突然大きな断絶が生じるわけでもありません。0.05は自然法則ではなく、あらかじめ定める判断基準です。

p値と一緒に何を報告するか

p値だけでなく、少なくとも次の情報を合わせて示すと、結果を解釈しやすくなります。

• 標本サイズ
• 個々のデータや分布
• 平均差など、元の単位での効果量
• 標準化効果量
• 信頼区間
• 使用した検定と、その仮定
• 除外や欠損、多重比較などの解析条件

今回のダミーデータなら、たとえば次のようにまとめられます。

20個の測定値の平均は54.65、標準偏差は9.43だった。基準値50との差は4.65で、母平均の95%信頼区間は50.24から59.07だった。両側1標本t検定では $t(19)=2.21$、$p=0.0398$、標本Cohen's $d=0.49$ だった。

この書き方なら、読者はp値だけでなく、差の大きさと不確かさも確認できます。

実データに適用するときの注意

今回の例は、仕組みを見やすくするため、独立な正規分布からダミーデータを生成しています。実データでは、少なくとも次を確認する必要があります。

• 観測値は本当に独立か
• 外れ値や強い歪みが結果を支配していないか
• 仮説や片側・両側の選択を、結果を見た後で変更していないか
• 多数の仮説を同時に検定していないか
• 欠損、除外、前処理を含めて解析手順を透明に説明できるか
• 統計的な差が、研究分野で意味のある大きさか

p値の計算が正しくても、研究設計やデータ生成過程に問題があれば、科学的な結論は信頼できません。

試してみる演習

コードの値を変更すると、p値の性質をさらに確認できます。

1. TRUE_MU = 50.0 にして、観測データ自体も $H_0$ から生成する
2. 観測データの乱数シードを変え、同じ生成条件でもp値が変動することを見る
3. N = 10, 50, 100 に変え、標準誤差とp値の変化を見る
4. SIGMA を大きくし、ばらつきが検定結果へ与える影響を見る
5. ALPHA = 0.01 に変え、第一種過誤率がどう変わるか確認する
6. 科学的に事前に方向を限定できる場合だけ、片側検定との違いを調べる

まとめ

今回の両側1標本t検定におけるp値は、

$$
p=P\left(|T|\geq|t_{\mathrm{obs}}|\mid H_0\right)
$$

です。

重要な点をまとめます。

• p値は「帰無仮説が正しい確率」ではない
• p値は「偶然だけで結果が起きた確率」ではない
• p値は効果量や科学的重要性を表さない
• $H_0$が正しくても、$p<0.05$ は約5%起こる
• 同じ効果量でも、標本サイズによってp値の分布は大きく変わる
• $p\geq0.05$ は、効果がないことの証明ではない
• p値は、効果量、信頼区間、データ分布、研究設計と一緒に解釈する

p値は、仮説の正しさを直接確率化する数字ではありません。帰無仮説を仮定した世界から、観測データがどれくらい外れて見えるかを測るための一つの指標です。