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

> 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

サンプルサイズ=12.83399であり、観測数n=10のときと比較すると検出力が0.9を超えるためにはあと3標本足りないことがわかりました。

まとめ

R言語で検出力とサンプルサイズを計算する方法をみてきました。

Rにはデフォルトで対応のあるt検定や母比率の差の検定、ANOVAに関する検出力の関数が導入されており、誰でも簡単に検出力とサンプルサイズを計算することができます。

usagi-san

統計学とゲームとかをメインに解説していくよ。数式とかプログラミングコードにミスがあったり質問があったりする場合はコメントで受け付けます。すぐに対応します。

2024/04/28

【R言語】F検定テンプレートスクリプト

2024/04/28

【R言語】経験密度関数・経験分布関数のプロット　関数densityとecdfの使い方

2024/04/10

【R言語】ベータ関数とガンマ関数　関数beta, gammaの使い方

usagi-sanの記事をもっと見る

-R言語
-R言語, プログラミング, 統計学, 統計解析

comment コメントをキャンセル

: R言語
R言語　自動でデータ整形してくれる関数【データ整形不要】
こんにちは、usagi-sanです。今回はデータ解析に役立つ関数を作ってみました。自分がアルバイトをする中で、データ整形（データクレンジング）に用いる時間が勤務時間の9割に及び、肝心な統計解析にミ ...

: R言語
【R言語】経験密度関数・経験分布関数のプロット　関数densityとecdfの使い方
R言語で経験密度関数と経験分布関数を求める方法及びグラフを描画する方法を解説していきます。この記事では経験密度関数と分布関数を計算する関数とその使用例について見ていきます。記事で扱うスクリプトは以 ...

: R言語
R言語　確率分布のプロット・乱数の生成　様々な乱数の算出方法
ここでは、R言語で確率分布をプロットしたり、確率分布からの乱数を生成したりする方法を解説していきます。確率密度や確率を算出するための関数さらに乱数を生成する関数を紹介し、関数の使用例をまとめています ...

: R言語
【R言語】factor型の使い方・まとめ
この記事では、R言語特有のfactor型（因子型）の操作方法について解説していきます。この記事では以下の操作方法を扱っています。 factorの作成からlevelsの変更などを紹介しています。この ...

: R言語
【R言語】母分散の検定・一標本カイ二乗検定
R言語で母分散の検定を行う関数やその実行例を紹介していきます。この記事では、一標本カイ二乗検定（one-sample chi squared test）と呼ばれる母分散の検定をRで行う方法を解説しま ...

【統計学】符号検定

【統計学】分散分析　一元配置分散分析・二元配置分散分析

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

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

t検定

母比率の差の検定

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

実行例

power.t.test

検出力

サンプルサイズ

power.prop.test

検出力

サンプルサイズ

power.anova.test

検出力

サンプルサイズ

まとめ

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