1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

はじめに

2群の測定値を比較したとき、解析結果を次のようにまとめることがあります。

$p < 0.05$ だったので、2群には有意差があった。

これは間違いとは限りません。しかし、p値だけでは「差がどれくらい大きいのか」も「推定値がどれくらい不確かなのか」も分かりません

そこで本記事では、Pythonでダミーデータを作り、次の3つを役割の異なる情報として整理します。

  • p値:帰無仮説のもとで、観測結果がどれくらい極端か
  • 効果量:2群の差がどれくらい大きいか
  • Bootstrap信頼区間:推定値がどれくらい不確かか

Google Colabでコードを順番に実行すると、次の4種類の図が得られます。

  1. 生データと各群の平均値
  2. p値に対応するt分布の裾の面積
  3. 平均差とCohen's $d$ のBootstrap分布
  4. 各群のサンプルサイズによって「$p < 0.05$ になる割合」が変わる様子

本記事の数値は、統計の考え方を説明するために作ったダミーデータです。実在する実験、疾患、処置などを表すものではありません。

この記事の結論

先に要点をまとめると、次のようになります。

  • $p < 0.05$ は、差が大きいことや実用上重要であることを意味しない
  • $p \ge 0.05$ は、差が存在しないことの証明ではない
  • 平均差は、元の単位で差の大きさを示す
  • Cohen's $d$ は、差を標準偏差単位に直した標準化効果量である
  • Bootstrap信頼区間を併記すると、推定値の幅を確認できる
  • p値、効果量、信頼区間、生データ、研究デザインを合わせて判断する必要がある

想定する状況

本記事では、母集団から得た一まとまりのデータを「標本」と呼び、その標本に含まれる観測値の個数 $n$ を「サンプルサイズ(標本サイズ)」と呼びます。Group AとGroup Bはそれぞれ1つの標本であり、今回の各群のサンプルサイズは $n=30$ です。

独立した2群、Group AとGroup Bがあり、各群で連続量を30個ずつ測定したとします。

今回の設定は次のとおりです。

  • Group Aを平均50、標準偏差10の正規分布から生成
  • Group Bを平均55、標準偏差10の正規分布から生成
  • 各群のサンプルサイズは30($n=30$)
  • 乱数シードを固定し、誰が実行しても同じ結果になるようにする

ここで指定する50、55、10は、データを生成する母集団側の設定です。実際に得られる標本平均や標本標準偏差は、乱数によって少しずれます。

1. Google Colabの準備

Google Colabの新しいノートブックを開き、最初のセルで次のコードを実行してください。

from pathlib import Path

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

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

# 図の基本設定
plt.rcParams.update({
    "figure.dpi": 120,
    "savefig.dpi": 180,
    "font.size": 11,
    "axes.grid": True,
    "grid.alpha": 0.25,
})

Google Colabには、通常、NumPy、pandas、Matplotlib、SciPyがあらかじめ入っています。そのため、追加のインストールは不要です。

2. ダミーデータを作る

DATA_SEED = 2024
rng = np.random.default_rng(DATA_SEED)

n_per_group = 30

group_a = rng.normal(
    loc=50,      # 母平均
    scale=10,    # 母標準偏差
    size=n_per_group,
)

group_b = rng.normal(
    loc=55,
    scale=10,
    size=n_per_group,
)

df = pd.DataFrame({
    "group": ["A"] * n_per_group + ["B"] * n_per_group,
    "value": np.concatenate([group_a, group_b]),
})

summary = (
    df.groupby("group")["value"]
      .agg(n="count", mean="mean", sd="std")
      .round(3)
)

summary

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

group n mean sd
A 30 50.295 9.556
B 30 55.995 9.596

標本平均の差は、およそ $55.995-50.295=5.70$ です。

ただし、この表だけでは個々の値の分布や重なり方が見えません。まずは生データを図にします。

生データを可視化する

PLOT_SEED = 100
plot_rng = np.random.default_rng(PLOT_SEED)

fig, ax = plt.subplots(figsize=(7, 5))

for x, (label, values) in enumerate([
    ("A", group_a),
    ("B", group_b),
]):
    color = f"C{x}"

    # 点が重ならないように、横方向へ小さくずらす
    jitter = plot_rng.normal(
        loc=0,
        scale=0.055,
        size=len(values),
    )

    ax.scatter(
        np.full(len(values), x) + jitter,
        values,
        alpha=0.75,
        s=40,
        color=color,
    )

    mean_value = np.mean(values)

    # 平均値を横線で表示
    ax.hlines(
        mean_value,
        x - 0.22,
        x + 0.22,
        linewidth=3,
        color=color,
    )

    ax.text(
        x + 0.26,
        mean_value,
        f"mean = {mean_value:.2f}",
        va="center",
        color=color,
    )

ax.set_xticks([0, 1], ["Group A", "Group B"])
ax.set_ylabel("Observed value")
ax.set_title("Raw data and group means")
ax.set_xlim(-0.5, 1.65)

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

fig1_raw_data.png

Group Bの平均値はGroup Aより高くなっています。一方で、個々の測定値にはかなりの重なりがあります。

このように、平均値だけでなく、生データの分布も確認することが大切です。

3. p値とは何か

帰無仮説と対立仮説

2群の母平均を、それぞれ $\mu_A$、$\mu_B$ とします。

今回の両側検定では、帰無仮説を次のように置きます。

$$
H_0:\mu_A=\mu_B
$$

対立仮説は次のとおりです。

$$
H_1:\mu_A\ne\mu_B
$$

帰無仮説は、今回の文脈では「母平均に差がない」という仮説です。

Welchのt検定

独立した2群の平均値を比較するため、ここではWelchのt検定を使います。

Welchのt統計量は、概念的には次の形です。

$$
t=
\frac{\bar{x}_B-\bar{x}_A}
{\sqrt{\frac{s_A^2}{n_A}+\frac{s_B^2}{n_B}}}
$$

ここで、

  • $\bar{x}_A,\bar{x}_B$:各群の標本平均
  • $s_A,s_B$:各群の標本標準偏差
  • $n_A,n_B$:各群のサンプルサイズ

です。

分子は観測された平均差、分母は平均差の標準誤差です。したがって、t統計量は「観測された差が、その標準誤差の何倍に相当するか」を表します。

SciPyでは、equal_var=Falseを指定するとWelchのt検定になります。

def welch_degrees_of_freedom(a, b):
    """Welch-Satterthwaiteの近似自由度を計算する。"""
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)

    n_a, n_b = len(a), len(b)
    v_a = np.var(a, ddof=1) / n_a
    v_b = np.var(b, ddof=1) / n_b

    numerator = (v_a + v_b) ** 2
    denominator = (
        v_a**2 / (n_a - 1)
        + v_b**2 / (n_b - 1)
    )

    return numerator / denominator


test = stats.ttest_ind(
    group_b,
    group_a,
    equal_var=False,
)

t_stat = float(test.statistic)
p_value = float(test.pvalue)
df_welch = float(
    welch_degrees_of_freedom(group_a, group_b)
)

print(f"Welch's t = {t_stat:.3f}")
print(f"df = {df_welch:.3f}")
print(f"p = {p_value:.6f}")

実行結果は次のとおりです。

Welch's t = 2.305
df = 57.999
p = 0.024754

p値の定義

両側検定のp値は、概念的には次のように書けます。

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

これは、

帰無仮説と統計モデルの仮定が成り立つとしたとき、今回観測されたt統計量と同じか、それ以上に絶対値が大きいt統計量が得られる確率

です。

今回の $p=0.024754$ は、帰無仮説のもとで $|t|\ge 2.305$ となる確率が約2.48%であることを表します。

p値を図で見る

x = np.linspace(-4.5, 4.5, 2000)
density = stats.t.pdf(x, df=df_welch)

tail_mask = np.abs(x) >= abs(t_stat)

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

ax.plot(
    x,
    density,
    linewidth=2,
    label=f"t distribution under H0 (df = {df_welch:.1f})",
)

ax.fill_between(
    x,
    0,
    density,
    where=tail_mask,
    alpha=0.35,
    label=f"Two-sided p-value = {p_value:.4f}",
)

ax.axvline(
    t_stat,
    linestyle="--",
    linewidth=2,
    label=f"Observed t = {t_stat:.2f}",
)
ax.axvline(
    -t_stat,
    linestyle="--",
    linewidth=2,
)

ax.set_xlabel("t statistic")
ax.set_ylabel("Probability density")
ax.set_title(
    "The p-value is the tail area under the null hypothesis"
)
ax.legend()

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

fig2_pvalue_tail.png

色が付いた左右の裾の面積を合計したものが、両側検定のp値です。

p値が直接は教えてくれないこと

$p=0.024754$ から、次のように読むことはできません。

  • 帰無仮説が正しい確率は2.48%である
  • Group BがGroup Aより高い確率は97.52%である
  • 2群の差は大きい
  • 2群の差は実用上重要である
  • 同じ実験を繰り返せば97.52%の確率で再現する

p値は、効果の大きさそのものを表す指標ではありません

4. 効果量とは何か

効果量は、差や関連の大きさを数値として表したものです。

2群の平均値を比較するとき、最も直接的な効果量は平均差です。

$$
\Delta=\bar{x}_B-\bar{x}_A
$$

平均差は元の測定単位を保つため、実験や応用上の意味を考えやすいという長所があります。

一方、測定単位の異なる研究間でも比較しやすいように、差を標準偏差で割った標準化効果量も使われます。代表例がCohen's $d$ です。

Cohen's d

本記事では、2群のプールした標準偏差 $s_p$ を使います。

$$
s_p=
\sqrt{
\frac{
(n_A-1)s_A^2+(n_B-1)s_B^2
}{
n_A+n_B-2
}
}
$$

Cohen's $d$ は次のように計算します。

$$
d=
\frac{
\bar{x}_B-\bar{x}_A
}{s_p}
$$

今回の符号はGroup BからGroup Aを引く向きに統一します。したがって、$d>0$ はGroup Bの平均値が高いことを表します。

def pooled_sd(a, b):
    """2群のプールした標準偏差を計算する。"""
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)

    n_a, n_b = len(a), len(b)
    s_a2 = np.var(a, ddof=1)
    s_b2 = np.var(b, ddof=1)

    return np.sqrt(
        (
            (n_a - 1) * s_a2
            + (n_b - 1) * s_b2
        )
        / (n_a + n_b - 2)
    )


def cohens_d(a, b):
    """Group B - Group AのCohen's dを計算する。"""
    return (
        np.mean(b) - np.mean(a)
    ) / pooled_sd(a, b)


mean_diff = float(
    np.mean(group_b) - np.mean(group_a)
)
d_value = float(
    cohens_d(group_a, group_b)
)

print(f"Mean difference (B - A) = {mean_diff:.3f}")
print(f"Cohen's d = {d_value:.3f}")

実行結果は次のとおりです。

Mean difference (B - A) = 5.699
Cohen's d = 0.595

今回の標本では、

  • 元の単位での平均差:$5.699$
  • 標準化平均差:$d=0.595$

となりました。

Cohen's $d$ については、$0.2$、$0.5$、$0.8$ をそれぞれ小、中、大の目安として紹介することがあります。しかし、これは固定的な判定基準ではありません。

$d=0.2$ でも重要な現象はあり得ますし、$d=0.8$ でも測定コスト、安全性、再現性などを考えると重要でない場合があります。効果量は、分野固有の基準や元の測定単位と合わせて解釈する必要があります

5. Bootstrapとは何か

観測データから再標本化する

Bootstrapは、観測データから復元抽出を繰り返し、統計量のばらつきを調べる方法です。

復元抽出では、一度選ばれた値を元へ戻してから次の値を選びます。そのため、1回のBootstrap標本の中に、同じ観測値が複数回含まれることがあります。

今回の手順は次のとおりです。

  1. Group Aから30個を復元抽出する
  2. Group Bから30個を復元抽出する
  3. 平均差とCohen's $d$ を計算する
  4. 1〜3を20,000回繰り返す
  5. 得られた値の2.5パーセンタイルと97.5パーセンタイルを95%区間とする

ここでは、最も直感的なパーセンタイル法を使います。

BOOTSTRAP_SEED = 2025
boot_rng = np.random.default_rng(BOOTSTRAP_SEED)

n_boot = 20_000

boot_diff = np.empty(n_boot)
boot_d = np.empty(n_boot)

for i in range(n_boot):
    sample_a = boot_rng.choice(
        group_a,
        size=len(group_a),
        replace=True,
    )

    sample_b = boot_rng.choice(
        group_b,
        size=len(group_b),
        replace=True,
    )

    boot_diff[i] = (
        np.mean(sample_b) - np.mean(sample_a)
    )

    boot_d[i] = cohens_d(
        sample_a,
        sample_b,
    )

# 単純なパーセンタイルBootstrap区間
diff_ci = np.percentile(
    boot_diff,
    [2.5, 97.5],
)

d_ci = np.percentile(
    boot_d,
    [2.5, 97.5],
)

print(f"Mean difference = {mean_diff:.3f}")
print(
    "95% bootstrap CI for mean difference = "
    f"[{diff_ci[0]:.3f}, {diff_ci[1]:.3f}]"
)

print(f"Cohen's d = {d_value:.3f}")
print(
    "95% bootstrap CI for Cohen's d = "
    f"[{d_ci[0]:.3f}, {d_ci[1]:.3f}]"
)

実行結果は次のとおりです。

Mean difference = 5.699
95% bootstrap CI for mean difference = [0.900, 10.488]
Cohen's d = 0.595
95% bootstrap CI for Cohen's d = [0.096, 1.172]

Bootstrap分布を可視化する

fig, axes = plt.subplots(
    1,
    2,
    figsize=(12, 4.5),
)

# 平均差のBootstrap分布
axes[0].hist(
    boot_diff,
    bins=45,
    alpha=0.8,
)
axes[0].axvline(
    0,
    linestyle="--",
    linewidth=2,
    label="No difference",
)
axes[0].axvline(
    mean_diff,
    linewidth=2,
    label="Observed",
)
axes[0].axvspan(
    diff_ci[0],
    diff_ci[1],
    alpha=0.18,
    label="95% bootstrap CI",
)
axes[0].set_xlabel("Mean difference (B - A)")
axes[0].set_ylabel("Bootstrap count")
axes[0].set_title(
    "Bootstrap distribution: mean difference"
)
axes[0].legend()

# Cohen's dのBootstrap分布
axes[1].hist(
    boot_d,
    bins=45,
    alpha=0.8,
)
axes[1].axvline(
    0,
    linestyle="--",
    linewidth=2,
    label="No standardized difference",
)
axes[1].axvline(
    d_value,
    linewidth=2,
    label="Observed",
)
axes[1].axvspan(
    d_ci[0],
    d_ci[1],
    alpha=0.18,
    label="95% bootstrap CI",
)
axes[1].set_xlabel("Cohen's d (B - A)")
axes[1].set_ylabel("Bootstrap count")
axes[1].set_title(
    "Bootstrap distribution: Cohen's d"
)
axes[1].legend()

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

fig3_bootstrap_distributions.png

左図は平均差、右図はCohen's $d$ のBootstrap分布です。

  • 破線:差が0の位置
  • 実線:元の標本から計算した点推定値
  • 淡い帯:95%パーセンタイルBootstrap区間

今回の平均差の点推定値は $5.699$ ですが、95%区間は $[0.900,10.488]$ です。

Cohen's $d$ の点推定値は $0.595$ ですが、95%区間は $[0.096,1.172]$ とかなり幅があります。

ここが重要です。

$p<0.05$ という二値的な結果だけを見ると、推定値に残っている幅広い不確かさが見えにくくなります。

なお、頻度論的な95%信頼区間は、計算後の特定の区間について「真値が95%の確率で入っている」と読むものではありません。同じ手続きを繰り返したとき、作られた区間がおおむね95%の割合で真の値を含む、という考え方です。

6. 同じ効果でも、サンプルサイズでp値は変わる

p値は、効果の大きさだけで決まりません。サンプルサイズとデータのばらつきにも影響されます。

そこで、母集団を次のように固定します。

$$
A\sim N(0,1)
$$

$$
B\sim N(0.4,1)
$$

母平均の差は0.4、母標準偏差は1なので、母集団の標準化平均差は $d=0.4$ です。

この真の効果を変えずに、各群のサンプルサイズだけを10、20、30、50、100、200と変えます。各サンプルサイズについて、独立した2群のデータを新たに生成するシミュレーションを3,000回反復し、次を調べます。

  • $p<0.05$ になった実験の割合
  • 推定されたCohen's $d$ の分布

コードでは計算を効率化するため、3,000回分のデータを配列としてまとめて生成しています。配列の各行が、1回分のシミュレーションに対応します。

SIMULATION_SEED = 2026
sim_rng = np.random.default_rng(SIMULATION_SEED)

sample_sizes = np.array([
    10,
    20,
    30,
    50,
    100,
    200,
])

n_sim = 3000
true_d = 0.4
alpha = 0.05

rows = []

for n in sample_sizes:
    # 各行を1回のシミュレーションとして、3,000回分の2群データをまとめて生成する
    a = sim_rng.normal(
        loc=0.0,
        scale=1.0,
        size=(n_sim, n),
    )

    b = sim_rng.normal(
        loc=true_d,
        scale=1.0,
        size=(n_sim, n),
    )

    # 各行を1回のシミュレーションとしてWelchのt検定を行う
    test_result = stats.ttest_ind(
        b,
        a,
        axis=1,
        equal_var=False,
    )
    p_values = test_result.pvalue

    # 各シミュレーションでCohen's dを計算する
    mean_a = a.mean(axis=1)
    mean_b = b.mean(axis=1)
    var_a = a.var(axis=1, ddof=1)
    var_b = b.var(axis=1, ddof=1)

    pooled = np.sqrt(
        (
            (n - 1) * var_a
            + (n - 1) * var_b
        )
        / (2 * n - 2)
    )

    d_values = (
        mean_b - mean_a
    ) / pooled

    rows.append({
        "n_per_group": n,
        "p_lt_0.05_rate": np.mean(
            p_values < alpha
        ),
        "median_d": np.median(d_values),
        "d_2.5%": np.percentile(
            d_values,
            2.5,
        ),
        "d_97.5%": np.percentile(
            d_values,
            97.5,
        ),
    })

simulation_summary = pd.DataFrame(rows)

print(
    simulation_summary
    .round(3)
    .to_string(index=False)
)

実行結果は次のとおりです。

 n_per_group  p_lt_0.05_rate  median_d  d_2.5%  d_97.5%
          10           0.132     0.406  -0.532    1.414
          20           0.237     0.409  -0.205    1.051
          30           0.323     0.403  -0.111    0.935
          50           0.510     0.402   0.017    0.806
         100           0.810     0.399   0.117    0.683
         200           0.980     0.397   0.205    0.602

p_lt_0.05_rateは、このシミュレーション条件における検出力の推定値に相当します。

結果を図にする

fig, axes = plt.subplots(
    1,
    2,
    figsize=(12, 4.5),
)

# p < 0.05になった割合
axes[0].plot(
    simulation_summary["n_per_group"],
    simulation_summary["p_lt_0.05_rate"],
    marker="o",
)
axes[0].set_xlabel("Sample size per group")
axes[0].set_ylabel("Proportion with p < 0.05")
axes[0].set_ylim(0, 1.02)
axes[0].set_title(
    "Same true effect, different sample sizes"
)

# Cohen's dの中央値と、反復実験間の2.5〜97.5パーセンタイル
lower_error = (
    simulation_summary["median_d"]
    - simulation_summary["d_2.5%"]
)
upper_error = (
    simulation_summary["d_97.5%"]
    - simulation_summary["median_d"]
)

axes[1].errorbar(
    simulation_summary["n_per_group"],
    simulation_summary["median_d"],
    yerr=[lower_error, upper_error],
    marker="o",
    capsize=4,
)
axes[1].axhline(
    true_d,
    linestyle="--",
    linewidth=1.5,
    label="True d = 0.4",
)
axes[1].set_xlabel("Sample size per group")
axes[1].set_ylabel("Estimated Cohen's d")
axes[1].set_title(
    "Effect-size estimates become more precise"
)
axes[1].legend()

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

fig4_sample_size_simulation.png

左図を見ると、真の効果を $d=0.4$ に固定しているにもかかわらず、$p<0.05$ になる割合は各群のサンプルサイズによって大きく変わります。

  • 各群10例:約13.2%
  • 各群50例:約51.0%
  • 各群100例:約81.0%
  • 各群200例:約98.0%

右図の点は、3,000回のシミュレーションから得られたCohen's $d$ の中央値です。縦線は、3,000個のCohen's $d$ の推定値について求めた2.5〜97.5パーセンタイルです。これは、1回のデータから計算した信頼区間ではなく、シミュレーションを繰り返したときに推定値がどの程度ばらつくかを表しています。

中央値はどのサンプルサイズでも真の値 $d=0.4$ 付近にあります。一方、サンプルサイズが小さいと推定値は大きくばらつき、サンプルサイズが増えると真の値の周囲へ狭く集まります。

この結果から、次の2点が分かります。

  1. 同じ大きさの効果でも、サンプルサイズが大きいほど $p<0.05$ になりやすい
  2. サンプルサイズが大きいほど、効果量を精密に推定しやすい

したがって、

サンプルサイズが小さいデータで $p\ge0.05$ だったから差がない

とも、

サンプルサイズが大きいデータで $p<0.05$ だったからといって、重要な差であるとは限らない

とも、p値だけからは結論できません。

7. 結果をどう報告するか

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

Group Aは平均 $50.30$、標準偏差 $9.56$、Group Bは平均 $55.99$、標準偏差 $9.60$ であった。平均差(B−A)は $5.70$、95%パーセンタイルBootstrap信頼区間は $[0.90,10.49]$ であった。標準化平均差はCohen's $d=0.60$、95%Bootstrap信頼区間は $[0.10,1.17]$ であった。Welchのt検定では $t(58.0)=2.31$、$p=0.0248$ であった。

この書き方には、少なくとも次の情報が含まれています。

  • 各群の中心とばらつき
  • 元の単位での平均差
  • 標準化効果量
  • それぞれの推定の不確かさ
  • 仮説検定の結果

ただし、実際の研究では、研究デザイン、欠測、外れ値、多重比較、事前に定めた解析計画、測定の信頼性なども確認する必要があります。

8. よくある誤解

誤解 より適切な考え方
$p<0.05$ なら重要な差である p値は効果の大きさや実用的重要性を直接表さない
$p\ge0.05$ なら差は存在しない データが不十分で差を検出できなかった可能性もある
p値は帰無仮説が正しい確率である p値は帰無仮説のもとでのデータ側の極端さに関する確率である
95%信頼区間には95%の確率で真値が入る 頻度論では、区間を作る手続きの長期的な被覆率として解釈する
Cohen's $d=0.5$ はどの分野でも「中程度」である 目安は文脈依存であり、元の単位や分野固有の基準も必要である
Bootstrapなら仮定は不要である 独立性や標本の代表性など、再標本化の前提は残る

9. 実データへ適用するときの注意

独立性

本記事では、各観測値が互いに独立であると仮定しています。

同一対象を繰り返し測定したデータ、同じ個体から得た複数細胞、同一施設に属する複数症例などを、すべて独立な値として扱うことはできません。

対応のあるデータでは対応差を再標本化し、階層構造のあるデータでは階層Bootstrapや混合効果モデルなどを検討します。

時系列データ

時系列には自己相関があるため、各時点を独立に再標本化すると時間構造が失われます。必要に応じてブロックBootstrapなどを使います。

外れ値と分布

平均値、t検定、Cohen's $d$ は外れ値の影響を受けます。生データの図を確認し、必要なら中央値、ロバストな効果量、置換検定なども検討します。

サンプルサイズが小さい場合

サンプルサイズが小さい場合(小標本)には、効果量の推定値が大きくばらつきやすくなります。また、Cohen's $d$ には小標本バイアスがあるため、それを補正したHedges' $g$ が使われることもあります。

Bootstrap区間の方法

本記事では理解しやすさを優先し、単純なパーセンタイル法を使いました。実務では、統計量やサンプルサイズに応じてBCa法などの利用も検討します。

「差がない」を示したい場合

通常の有意差検定で $p\ge0.05$ になっただけでは、2群が十分に近いとは結論できません。

「実質的に無視できる差の範囲」を事前に定めたうえで、同等性検定などを検討する必要があります。

10. 発展課題

コードの値を変更すると、統計量の振る舞いをさらに確認できます。

課題1:サンプルサイズを変える

n_per_group = 10

または、

n_per_group = 200

として、p値、Cohen's $d$、Bootstrap区間がどう変わるか確認してください。

課題2:平均差を小さくする

group_b = rng.normal(
    loc=51,
    scale=10,
    size=n_per_group,
)

サンプルサイズが大きいと、小さな平均差でも $p<0.05$ になることがあります。

課題3:ばらつきを大きくする

group_b = rng.normal(
    loc=55,
    scale=20,
    size=n_per_group,
)

同じ平均差でも、ばらつきが大きいほど推定が難しくなることを確認できます。

課題4:外れ値を加える

group_b_with_outlier = group_b.copy()
group_b_with_outlier[0] = 120

平均値、p値、Cohen's $d$、Bootstrap分布がどのように変わるか比較してみてください。

まとめ

本記事では、ダミーデータを使ってp値、効果量、Bootstrap信頼区間を比較しました。

重要な点は、3つの指標が別々の問いに答えていることです。

指標 主に答える問い
p値 帰無仮説のもとで、今回の結果はどれくらい極端か
平均差・Cohen's $d$ 差はどれくらい大きいか
Bootstrap信頼区間 効果量の推定にはどれくらい幅があるか

統計解析では、結果を $p<0.05$ と $p\ge0.05$ の2種類だけに分けるのではなく、

  1. 生データを見る
  2. 元の単位で差を見る
  3. 標準化効果量を見る
  4. 信頼区間で不確かさを見る
  5. p値とモデルの仮定を確認する
  6. 分野上の重要性を考える

という順序で考えると、より多くの情報を保った解釈ができます。

1
0
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?