R言語　検出力・サンプルサイズ　t検定や母比率の検定の検出力の計算

R言語で様々な検定の検出力やサンプルサイズの計算方法を紹介します。

この記事では、二標本t検定や母比率の差の検定、ANOVAの検出力とサンプルサイズの計算方法をみていきます。

観測数やパラメータ、第一の過誤が与えられたときの検出力や、反対にあるパラメータの下で検出力がある値を超えるために必要なサンプルサイズを計算していきます。

この記事で紹介するソースコードは以下からダウンロードできます。

R言語　検出力・サンプルサイズ

1 ファイル 1.43 KB

ダウンロード

Rによるやさしい統計学

統計学初学者向けの教材です

¥2,970 （2022/06/01 17:41時点 | Amazon調べ）

Amazon

楽天市場

Yahooショッピング

ポチップ

検出力・サンプルサイズの計算

t検定

t検定の検出力及びサンプルサイズの計算方法をみていきます。

対応のあるt検定における検出力及びサンプルサイズは以下のようにして計算されます。この後紹介するR言語の関数では対応のある二標本検定だけでなく一標本検定や二標本検定の検出力とサンプルサイズを計算できます。

t検定の検出力・サンプルサイズ

$x_{11}, \ldots, x_{1n}$と$x_{21}, \ldots, x_{2n}$をそれぞれ$N(\mu_1, \sigma^2)$、$N(\mu_2, \sigma^2)$からの大きさ$n$の無作為標本とする。次の「2群間の平均は等しい」という仮説を考える。

\begin{align}&H_0:\ \delta =0\\&H_1:\ \delta \neq = 0\end{align}

ここに、$\delta = \mu_1 - \mu_2$。第一の過誤を$\alpha$、第二の過誤を$\beta$とすると、この仮説検定の検出力とサンプルサイズは以下の式より決定される。

\begin{align}n = (t_{n - 1, \alpha / 2} - t_{n - 1, \beta} )^2 / \delta^2,\end{align}

ここに、$t_{n- 1, \varepsilon}$は自由度$n-1$のt分布の上側$\varepsilon$点である。上式より、第一の過誤$\alpha$、2つの母集団の分散$\sigma^2$、観測数$n_1$、$n_2$、母平均の差$\delta$が与えられているときのサンプルサイズ$n$及び検出力$1-\beta$を求めることができる。

R言語では次の関数を用いることでt検定の検出力とサンプルサイズを計算することができます。

power.t.test(n = NULL, delta = NULL, sd = 1, sig.level = 0.05, power = NULL, type = c("two.sample", "one.sample", "paired"), alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25)

関数power.t.testの引数
n	観測数、標本数。
delta	母平均の差。
sd	標準偏差。
sig.level	第一の過誤。検定の有意水準。
type	検定の種類。一標本（"one.sample"）、2標本（"two.sample"）、対応のある（"paired"）検定をcharacter型で指定できる。
alternative	両側検定（"two.sided"）か片側検定（"one.sided"）かをcharacter型で指定できる。
strict	TRUEのとき両側検定の際に効果（true effect）の符号とは反対の棄却確率も検出力に含める。厳密な計算を行う。
tol	求根アルゴリズムにおける有効桁数。

関数power.t.testの使用例はこの後の実行例で解説しています。

母比率の差の検定

次に母比率の差の検定の検出力及びサンプルサイズの計算方法をみていきます。

母比率の差の検定における検出力及びサンプルサイズは以下のようにして計算されます。

母比率の差の検定の検出力・サンプルサイズ

$x_{11}, \ldots, x_{1n}$と$x_{21}, \ldots, x_{2n}$をそれぞれ$Beroulli(p_1)$、$Bernoulli(p_2)$からの大きさ$n$の無作為標本とする。次の「2群の母比率は等しい」という仮説を考える。

\begin{align}&H_0:\ p_1 - p_2 =0\\&H_1:\ p_1 - p_2 \neq = 0\end{align}

第一の過誤を$\alpha$、第二の過誤を$\beta$とすると、この仮説検定の検出力とサンプルサイズは以下の式より決定される。

\begin{align}n = (Z_{\alpha/2} + Z_{\beta})^2\bigl\{ p_1(1-p_1) + p_2(1-p_2)\bigr\}/(p_1 - p_2)^2, \end{align}

ここに、$Z_{\varepsilon}$は標準正規分布の上側$\varepsilon$点である。上式より、第一の過誤$\alpha$、2つの母集団の分散$\sigma^2$、観測数$n_1$、$n_2$、母平均の差$\delta$が与えられているときのサンプルサイズ$n$及び検出力$1-\beta$を求めることができる。

R言語では次の関数を用いることで母比率の差の検定の検出力とサンプルサイズを計算することができます。

power.prop.test(n = NULL, p1 = NULL, p2 = NULL, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25)

関数power.prop.testの引数
n	観測数、標本数。
p1	第一群の母比率
p2	第二群の母比率。
sig.level	第一の過誤。検定の有意水準。
power	検出力。
alternative	両側検定（"two.sided"）か片側検定（"one.sided"）かをcharacter型で指定できる。
strict	TRUEのとき両側検定の際に効果（true effect）の符号とは反対の棄却確率も検出力に含める。厳密な計算を行う。
tol	求根アルゴリズムにおける有効桁数。

関数power.prop.testの使用例はこの後の実行例で解説しています。

ANOVA（多群間の母平均の差の検定）

最後にANOVAの検出力及びサンプルサイズの計算方法をみていきます。

ANOVAの検出力及びサンプルサイズは以下のようにして計算されます。

ANOVAの検出力・サンプルサイズ

$x_{i1}, \ldots, x_{1n_i}$を$N(\mu_g, \sigma^2)$からの大きさ$n$の第$i = 1, \ldots, g$群の無作為標本とする。次の「第$1$から第$g$群の母比率は等しい」という仮説を考える。

\begin{align}&H_0:\ \mu1 = \mu_2 = \cdots = \mu_k =0\\&H_1:\ \mathrm{not}\ H_0\end{align}

非心パラメータを次で与える。

\begin{align}\lambda = (k-1)n \sigma_B^2 / \sigma_W^2,\end{align}

ここに、$\sigma_B^2$と$\sigma_W^2$はそれぞれ群間、群内分散である。第一の過誤を$\alpha$、第二の過誤を$\beta$とすると、この仮説検定の検出力とサンプルサイズは以下の式より決定される。

\begin{align}1 - \beta = 1 - F(f_{n-1, k(n-1), \alpha} ; k- 1, k(n-1) , \lambda), \end{align}

ここに、$f_{n-1, k(n-1), \alpha}$は自由度$n-1$、$k(n-1)$のF分布の上側$\alpha$点であり、$F(f_{\alpha}; n-1, k(n-1), \lambda)$は非心パラメータ$\lambda$、自由度$n-1$、$k(n-1)$のF分布の分布関数である。

R言語では次の関数を用いることでANOVAの検出力とサンプルサイズを計算することができます。

power.anova.test(groups = NULL, n = NULL, between.var = NULL, within.var = NULL, sig.level = 0.05, power = NULL)

関数power.anova.testの引数
goups	群数、グループ数。
n	観測数、標本数。
beetween.var	群間の分散。
within.var	群内の分散。
sig.level	第一の過誤。検定の有意水準。
power	検出力。

関数power.anova.testの使用例をこの後みていきます。

実行例

power.t.test

関数power.t.testを用いて検出力とサンプルサイズの計算を行っていきます。

検出力

次を実行することで、観測数$n=10$、母平均の差$\delta = 0.1$、標準偏差$\sigma = 0.05$、第一の過誤$\alpha = 0.05$のときのt検定の検出力を計算することができます。

n <- 10
alpha <- 0.05
delta <- 0.1
sd <- 0.05

power <- power.t.test(n = n, delta = delta, sd = sd, sig.level = alpha)

n <- 10

alpha <- 0.05

delta <- 0.1

sd <- 0.05

power <- power.t.test(n = n, delta = delta, sd = sd, sig.level = alpha)

powerを参照すると検出力=0.988179であることが分かります。

> power

     Two-sample t test power calculation 

              n = 10
          delta = 0.1
             sd = 0.05
      sig.level = 0.05
          power = 0.988179
    alternative = two.sided

NOTE: n is number in *each* group

> power

Two-sample t test power calculation

n = 10

delta = 0.1

sd = 0.05

sig.level = 0.05

power = 0.988179

alternative = two.sided

NOTE: n is number in *each* group

power.t.testの返り値はリスト型となっており、観測数$n=20$や検出力$1-\beta$の値がリストの要素として格納されています。

次のようにして検出力やそのほかのパラメータを取得することが可能です。

power$n         #観測数:n
power$delta     #仮説検定の2群間の差:δ
power$sd        #標準偏差
power$sig.level #有意水準:α
power$power     #検出力:1-β

power$n #観測数:n

power$delta #仮説検定の2群間の差:δ

power$sd #標準偏差

power$sig.level #有意水準:α

power$power #検出力:1-β

次のように、power.t.testによって計算された値をデータフレームに格納すると、csvなどに表として出力するとき便利です。

> powerResult <- data.frame(n = power$n, delta = power$delta, sd = power$sd,
+                           alpha = power$sig.level, power = power$power, row.names = NULL)
> powerResult
   n delta   sd alpha    power
1 10   0.1 0.05  0.05 0.988179

> powerResult <- data.frame(n = power$n, delta = power$delta, sd = power$sd,

+ alpha = power$sig.level, power = power$power, row.names = NULL)

> powerResult

n delta sd alpha power

1 10 0.1 0.05 0.05 0.988179

write.csvを使ってpowerResultを出力すると次の画像のcsvファイルが保存されているのが確認できます。

write.csv("検出力_t検定", powerResult, row.names = FALSE)

1	write.csv("検出力_t検定", powerResult, row.names = FALSE)

サンプルサイズ

次にサンプルサイズの計算方法についてみていきます。

検出力が与えられたときのサンプルサイズはpower.t.testのn以外のパラメータを与えることで計算できます。

次を実行することで、母平均の差$\delta = 0.1$、標準偏差$\sigma = 0.05$、第一の過誤$\alpha = 0.05$のとき検出力が$0.9$を超えるようなサンプルサイズを計算することができます。

beta <- 0.1

power <- power.t.test(delta = delta, sd = sd, sig.level = alpha, power = 1 - beta)
power$n #サンプルサイズ

beta <- 0.1

power <- power.t.test(delta = delta, sd = sd, sig.level = alpha, power = 1 - beta)

power$n #サンプルサイズ

powerを参照するとサンプルサイズは6.386756であることが分かりました。

> power

     Two-sample t test power calculation 

              n = 6.386756
          delta = 0.1
             sd = 0.05
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

NOTE: n is number in *each* group

> power

Two-sample t test power calculation

n = 6.386756

delta = 0.1

sd = 0.05

sig.level = 0.05

power = 0.9

alternative = two.sided

NOTE: n is number in *each* group

余談ですが、実際に検出力と観測数についてのグラフを描くと上で計算したサンプルサイズが正しいことが確認できます。

powers <- sapply(seq_len(50), function(x) {
                  power <- power.t.test(n = x, delta = delta, sd = sd, sig.level = alpha)
                  power$power
                 })

plot(powers, type = "l", xlab = "n", ylab= "Power")

powers <- sapply(seq_len(50), function(x) {

power <- power.t.test(n = x, delta = delta, sd = sd, sig.level = alpha)

power$power

})

plot(powers, type = "l", xlab = "n", ylab= "Power")

power.prop.test

次にpower.prop.testの使い方についてみていきます。

検出力

次を実行することで、観測数$n=10$、第一群の母比率$p_1= 0.5$、第二群の母比率$p_2 = 0.1$、第一の過誤$\alpha = 0.05$のときの母比率の差の検定の検出力を計算することができます。

n <- 10
alpha <- 0.05
p1 <- 0.5
p2 <- 0.1

power <- power.prop.test(n = n, p1 = p1, p2 = p2, sig.level = alpha)

n <- 10

alpha <- 0.05

p1 <- 0.5

p2 <- 0.1

power <- power.prop.test(n = n, p1 = p1, p2 = p2, sig.level = alpha)

power.prop.testの使い方はpower.t.testとほとんど同じであり、引数を母比率の検定の者に換えてあげる必要があります。

powerを参照すると検出力=0.4963802であることが分かります。

> power

     Two-sample comparison of proportions power calculation 

              n = 10
             p1 = 0.5
             p2 = 0.1
      sig.level = 0.05
          power = 0.4963802
    alternative = two.sided

NOTE: n is number in *each* group

> power

Two-sample comparison of proportions power calculation

n = 10

p1 = 0.5

p2 = 0.1

sig.level = 0.05

power = 0.4963802

alternative = two.sided

NOTE: n is number in *each* group

power.t.testと同様に観測数や検出力の値を以下のようにして取得することができます。

power$n         #観測数:n
power$p1        #1つ目のグループの母比率（確率）
power$p2        #2つ目のグループの母比率（確率）
power$sig.level #有意水準:α
power$power     #検出力:1-β

power$n #観測数:n

power$p1 #1つ目のグループの母比率（確率）

power$p2 #2つ目のグループの母比率（確率）

power$sig.level #有意水準:α

power$power #検出力:1-β

> powerResult <- data.frame(n = power$n, p1 = power$p1, p2 = power$p2,
+                           alpha = power$sig.level, power = power$power, row.names = NULL)
> powerResult
   n  p1  p2 alpha     power
1 10 0.5 0.1  0.05 0.4963802

> powerResult <- data.frame(n = power$n, p1 = power$p1, p2 = power$p2,

+ alpha = power$sig.level, power = power$power, row.names = NULL)

> powerResult

n p1 p2 alpha power

1 10 0.5 0.1 0.05 0.4963802

サンプルサイズ

次を実行することで、第一群の母比率$p_1= 0,5$、第二群の母比率$p_2 = 0.1$。第一の過誤$\alpha = 0.05$のとき検出力が$0.9$を超えるようなサンプルサイズを計算することができます。

power <- power.prop.test(p1 = p1, p2 = p2, sig.level = alpha, power = 1- beta)
power$n #サンプルサイズ

1 2	power <- power.prop.test(p1 = p1, p2 = p2, sig.level = alpha, power = 1- beta) power$n #サンプルサイズ

powerを参照するとサンプルサイズは25.43862であることが分かります。

> power

     Two-sample comparison of proportions power calculation 

              n = 25.43862
             p1 = 0.5
             p2 = 0.1
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

NOTE: n is number in *each* group

> power

Two-sample comparison of proportions power calculation

n = 25.43862

p1 = 0.5

p2 = 0.1

sig.level = 0.05

power = 0.9

alternative = two.sided

NOTE: n is number in *each* group

検出力の計算のときは観測数10でやっており検出力は0.4963802でした。検出力が0.9を超えるような検定を実施したい場合、あと最低でも16標本必要だといえます。

power.anova.test

最後にpower.anova.testの使い方を紹介します。

検出力

power.anova.testも他の関数と同じで、引数に検定で用いるパラメータや第一の過誤を設定することで実行できます。

次を実行すると、群の数$k = 4$、観測数$n=10$、群間の分散$\sigma_B^2= 0.01$、群内の分散$\sigma_W^2= 0.025$、第一の過誤$\alpha = 0.05$のときの母比率の差の検定の検出力を計算することができます。

n <- 10
alpha <- 0.05
numberOfGroups <- 4
sigmaB <- 0.01 #群間の分散
sigmaW <- 0.025 #群内の分散

power <- power.anova.test(groups = numberOfGroups, n = n, between.var = sigmaB,
                          within.var = sigmaW, sig.level = alpha)

n <- 10

alpha <- 0.05

numberOfGroups <- 4

sigmaB <- 0.01 #群間の分散

sigmaW <- 0.025 #群内の分散

power <- power.anova.test(groups = numberOfGroups, n = n, between.var = sigmaB,

within.var = sigmaW, sig.level = alpha)

powerを参照すると検出力=0.7942482であることが確認できます。

> power

     Balanced one-way analysis of variance power calculation 

         groups = 4
              n = 10
    between.var = 0.01
     within.var = 0.025
      sig.level = 0.05
          power = 0.7942482

NOTE: n is number in each group

> power

Balanced one-way analysis of variance power calculation

groups = 4

n = 10

between.var = 0.01

within.var = 0.025

sig.level = 0.05

power = 0.7942482

NOTE: n is number in each group

他の関数と同様に検出力などの値はリストの要素として格納されている多め、powerのうしろに$を付けることで各値を取得することが可能です。

power$groups      #群の数
power$n           #観測数
power$between.var #群間の分散
power$within.var  #群内の分散
power$sig.level   #有意水準:α
power$power       #検出力:1-β

power$groups #群の数

power$n #観測数

power$between.var #群間の分散

power$within.var #群内の分散

power$sig.level #有意水準:α

power$power #検出力:1-β

> powerResult <- data.frame(No.groups = power$groups, n = power$n, sigmaB = power$between.var,
+                           sigmaW = power$within.var, alpha = power$sig.level, power = power$power, row.names = NULL)
> powerResult
  No.groups  n sigmaB sigmaW alpha     power
1         4 10   0.01  0.025  0.05 0.7942482

> powerResult <- data.frame(No.groups = power$groups, n = power$n, sigmaB = power$between.var,

+ sigmaW = power$within.var, alpha = power$sig.level, power = power$power, row.names = NULL)

> powerResult

No.groups n sigmaB sigmaW alpha power

1 4 10 0.01 0.025 0.05 0.7942482

サンプルサイズ

サンプルサイズも他の関数と同様に観測数以外の引数を与えてあげることで計算することができます。

以下、群の数$k = 4$、群間の分散$\sigma_B^2= 0.01$、群内の分散$\sigma_W^2= 0.025$、第一の過誤$\alpha = 0.05$のとき検出力が0.9を超えるようなサンプルサイズを計算する例です。

power <- power.anova.test(groups = numberOfGroups, between.var = sigmaB,
                          within.var = sigmaW, sig.level = alpha, power = 1 - beta)
power$n #サンプルサイズ

power <- power.anova.test(groups = numberOfGroups, between.var = sigmaB,

within.var = sigmaW, sig.level = alpha, power = 1 - beta)

power$n #サンプルサイズ

> power

     Balanced one-way analysis of variance power calculation 

         groups = 4
              n = 12.83399
    between.var = 0.01
     within.var = 0.025
      sig.level = 0.05
          power = 0.9

NOTE: n is number in each group