Rで確率分布:非心F分布

非心F分布の定義

独立な確率変数 \(U_1\) と \(U_2\) が、それぞれ

• \(U_1 \sim \chi^2(d_1, \lambda)\) （自由度 \(d_1\)、非心性パラメータ \(\lambda\) の非心カイ2乗分布）
• \(U_2 \sim \chi^2(d_2)\) （自由度 \(d_2\) の中心カイ2乗分布）

に従うとします。このとき、それぞれのカイ2乗変数をその自由度で割ったものの比、

\[F = \dfrac{U_1/d_1}{U_2/d_2}\]

が従う分布が、分子の自由度 \(d_1\)、分母の自由度 \(d_2\)、非心性パラメータ \(\lambda\) の非心F分布です。

確率密度関数 (Probability Density Function, PDF)

非心F分布の確率密度関数の分布は、3つのパラメータによってその形状が決定されます。

• \(d_1\): 分子の自由度 (numerator degrees of freedom)
• \(d_2\): 分母の自由度 (denominator degrees of freedom)
• \(\lambda\) (ラムダ): 非心性パラメータ (non-centrality parameter, ncp)
  • 分子の非心カイ2乗分布から引き継がれるパラメータで、分布の中心からの「ズレ」の大きさを表します。
  • \(\lambda=0\) のとき、分子は中心カイ2乗分布となり、分布全体は中心F分布と完全に一致します。

主な特徴

• 定義域: 分散の比であるため、確率変数は常に非負の値 (\(x \ge 0\)) をとります。
• 形状:
  • \(\lambda=0\) のときは、中心F分布（右に歪んだ山型）。
  • \(\lambda > 0\) のときは、中心F分布と比較して分布全体が右にシフトし、ばらつきが大きくなります（山が低く、裾が広くなる）。
• 代表値:
  • 平均 (Mean): \(E[F] = \dfrac{d_2\left(d_1+\lambda\right)}{d_1\left(d_2-2\right)}\) （ただし、\(d_2 > 2\)）
  • 分散も平均と同様に3つのパラメータすべてに依存します。
• 他の分布との関係:
  • 確率変数 \(T\) が自由度 \(k\)、非心性パラメータ \(\delta\) の非心t分布に従うとき、その2乗 \(T^2\) は、自由度(1, k)、非心性パラメータ \(\lambda=\delta^2\) の非心F分布に従います。\(T^2 \sim F(1, k, \delta^2)\)。

Figure 1 の解説

上記のRコードを実行すると、4つの非心F分布が描画されたチャート Figure 1 が生成されます。

• d1=4, d2=30, λ=0 (黒線): 比較の基準となる中心F分布です。ピークが1の近くにあり、右に歪んだ形状をしています。これは帰無仮説（真の差がない）のもとでのF統計量の分布です。
• d1=4, d2=30, λ=5 (青線): 黒線と同じ自由度ですが、非心性パラメータ \(\lambda=5\) が加わっています。分布全体が右にシフトし、山が低く平たくなっている（ばらつきが大きくなっている）ことがわかります。
• d1=4, d2=30, λ=15 (赤線): 非心性パラメータがさらに大きい \(\lambda=15\) になったケースです。青線よりもさらに大きく右にシフトし、より平坦な分布になっています。
• d1=4, d2=10, λ=5 (緑線): 青線と同じ非心性パラメータですが、分母の自由度 \(d_2\) が30から10に小さくなっています。これにより分布のばらつきがさらに大きくなり、青線よりもさらに低く、平坦な形状になっていることが確認できます。

このシミュレーションから、非心性パラメータ \(\lambda\) が大きいほど（＝対立仮説における効果量が大きいほど）、F統計量の分布が右にシフトしていく様子がわかります。ANOVAで棄却域を設定した場合、\(\lambda\) が大きいほど分布の右側の裾が棄却域に大きくかかることになり、検出力が高まるというメカニズムを視覚的に理解することができます。

対立仮説における効果量が大きいほどとは

ある研究者が、3種類の新しい肥料（A, B, C）がトウモロコシの収穫量に与える影響を調べたいと考えています。

• 帰無仮説 (\(H_0\)): 「3種類の肥料の間に効果の差はない」
  • つまり、どの肥料を使っても、トウモロコシの平均収穫量は同じである。(\(\mu_A = \mu_B = \mu_C\))
• 対立仮説 (\(H_1\)): 「少なくとも1種類の肥料は、他の肥料と効果が異なる」
  • 3つの平均収穫量のうち、少なくとも1つは他と異なる。

ここで、「効果量」とは、この対立仮説が「どのくらい本当か」の度合い、つまりグループ間の真の平均値が、実際にどれだけ異なっているかの大きさを指します。

以下に、「効果量が小さい」対立仮説と「効果量が大きい」対立仮説の具体的な例を示します。

1. シナリオ設定

分散分析（ANOVA）を想定した、具体的なシナリオを設定します。

• 検定の種類: 一元配置分散分析 (One-way ANOVA)
• グループ数: 5群
• 各グループのサンプルサイズ: 7人
  • この設定により、分子の自由度 \(d_1 = 5 - 1 = 4\)
  • 分母の自由度 \(d_2 = 5 \times (7-1) = 30\)
  • これは Figure 1 のシミュレーションのパラメータ（黒、青、赤の線）と一致します。
• 有意水準 \(\alpha\): 0.05
• 非心性パラメータ \(\lambda\): 0, 5, 15 の3つのケースを比較

2. 関数 df() と integrate() を利用した検出力の計算

Step 1: 棄却域の決定

まず、帰無仮説（中心F分布, λ=0）のもとで、F統計量が上位5%に入る境界値（臨界値）を求めます。

# 各パラメータを設定
d1 <- 4
d2 <- 30
alpha <- 0.05

# 臨界値を計算 (P(F > critical_value) = alpha となる値)
critical_value <- qf(1 - alpha, df1 = d1, df2 = d2)
cat(paste("臨界値:", round(critical_value, 4)))

臨界値: 2.6896

棄却域は F > 2.6896 となります。

Step 2: 確率密度関数を積分して検出力を計算

次に、各非心性パラメータ \(\lambda\) のもとでの非心F分布 (df(x, df1, df2, ncp)) を、この棄却域（2.6896 から無限大まで）で積分します。

# ケース1: λ = 0 の場合 (これは第1種の過誤の確率αと一致する)
power_lambda0 <- integrate(
  function(x) df(x, df1 = d1, df2 = d2, ncp = 0),
  lower = critical_value,
  upper = Inf
)$value

# ケース2: λ = 5 の場合
power_lambda5 <- integrate(
  function(x) df(x, df1 = d1, df2 = d2, ncp = 5),
  lower = critical_value,
  upper = Inf
)$value

# ケース3: λ = 15 の場合
power_lambda15 <- integrate(
  function(x) df(x, df1 = d1, df2 = d2, ncp = 15),
  lower = critical_value,
  upper = Inf
)$value

# 結果をデータフレームにまとめる
results_df <- data.frame(
  lambda = c(0, 5, 15),
  power_by_df = c(power_lambda0, power_lambda5, power_lambda15)
)

cat("--- 非心性パラメータ毎の検出力 ---\n\n")
print(results_df)

--- 非心性パラメータ毎の検出力 ---

  lambda power_by_df
1      0  0.04999998
2      5  0.34208849
3     15  0.83103854

この結果から、非心性パラメータ \(\lambda\) が0から5、15と大きくなるにつれて、検出力が 5% → 34.2% → 83.1% と上昇していることが数値で明確に確認できました。

2. パッケージ pwr を利用した結果との比較

次に、パッケージ pwr の結果と比較します。

注意点: pwr.anova.test() 関数では、非心性パラメータ \(\lambda\) の代わりに、効果量 (Cohen’s f) を用います。両者の間には以下の関係があります（\(k\): グループ数, \(n\): 1グループのサンプルサイズ）。

\[\lambda = n \cdot k \cdot f^2\]

この関係を使い、\(\lambda\) を \(f\) に変換してから関数に入力します。\[f = \sqrt{\dfrac{\lambda}{n \cdot k}}\]

library(pwr)

# グループ数と1グループのサンプルサイズ
groups_k <- d1 + 1
n_per_group <- (d2 / groups_k) + 1

# ケース1: λ = 0 (効果量f=0)
f_val0 <- sqrt(0 / (n_per_group * groups_k))
pwr_result0 <- pwr.anova.test(
  k = groups_k, n = n_per_group, f = f_val0,
  sig.level = alpha
)$power

# ケース2: λ = 5
f_val5 <- sqrt(5 / (n_per_group * groups_k))
pwr_result5 <- pwr.anova.test(
  k = groups_k, n = n_per_group, f = f_val5,
  sig.level = alpha
)$power

# ケース3: λ = 15
f_val15 <- sqrt(15 / (n_per_group * groups_k))
pwr_result15 <- pwr.anova.test(
  k = groups_k, n = n_per_group, f = f_val15,
  sig.level = alpha
)$power

# 比較結果を元のデータフレームに追加
results_df$power_by_pwr <- c(pwr_result0, pwr_result5, pwr_result15)

cat("--- 関数 pwr.anova.test との比較 ---\n\n")
print(results_df)

--- 関数 pwr.anova.test との比較 ---

  lambda power_by_df power_by_pwr
1      0  0.04999998    0.0500000
2      5  0.34208849    0.3420885
3     15  0.83103854    0.8310385

df関数を積分して求めた検出力 (power_by_df) と、pwrパッケージを用いて計算した検出力 (power_by_pwr) の結果は一致しました。

以上です。