Rで確率･統計:サンプルサイズの決定:t検定の場合(2群、対応なし、異なるサンプルサイズ)

Rでデータサイエンス

サンプルサイズの決定:t検定の場合(2群、対応なし、異なるサンプルサイズ)

lapply(X = c("dplyr", "pwr", "ggplot2", "scales"), require, character.only = T)

t検定(2群、対応なし、異なるサンプルサイズ)の例

購入時期の違う同仕様の部品Z製造装置が2つ(AとB)がある。Aで製造した\(n_1\)個のZの長さの平均値(\(\mu_A\))とBで製造した\(n_2\)個のZの長さの平均値(\(\mu_B\))について、

  • \(\textrm{H}_0\):\(\mu_A\)と\(\mu_B\)に差はない。
  • \(\textrm{H}_1\):\(\mu_A\)と\(\mu_B\)に差がある(\(\mu_A\)の方が\(\mu_B\)より長い又はその逆の両側検定)。

効果量、有意水準、サンプルサイズが与えられた場合の検定力算出

\[d=\dfrac{\mu_A-\mu_B}{\sigma}\]

# 効果量の目安
effect_size <- sapply(c("small", "medium", "large"), function(x) {
    cohen.ES(test = "t", size = x)$effect.size
})
effect_size
 small medium  large 
   0.2    0.5    0.8 
# 検定力算出コードは関数から引用
# 効果量(Effect size)
# 標準偏差に対する平均値の差の寄与が 0.5
d <- 0.5
# 有意水準
sig.level <- 0.05
# サンプルサイズ
n1 <- 40
n2 <- 60
# 検定タイプ
alternative <- "two.sided"
tsample <- 2
ttside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
tside <- switch(alternative, less = 1, two.sided = 2, greater = 1)
# 自由度
nu <- n1 + n2 - 2
nu
[1] 98
# 自由度 nu のt分布における、境界点(両側) sig.level/2 に対応するt値
qu <- qt(p = sig.level/tside, df = nu, lower = FALSE)
qu
[1] 1.984467
# 自由度 nu の非心t分布における検定統計量(qu)に対する上側確率
p_upper <- pt(qu, nu, ncp = d * (1/sqrt(1/n1 + 1/n2)), lower = FALSE)
p_upper
[1] 0.6792001
# 自由度 nu の非心t分布における検定統計量(-qu)に対する下側確率
p_lower <- pt(-qu, nu, ncp = d * (1/sqrt(1/n1 + 1/n2)), lower = TRUE)
p_lower
[1] 0.000005778559
# 両者の合計
power <- p_upper + p_lower
power
[1] 0.6792059

よって、対立仮説が正しく、かつ帰無仮説が棄却される確率は67.921`%

# 改めて関数で検定力を確認
result <- pwr.t2n.test(n1 = n1, n2 = n2, d = d, sig.level = sig.level, alternative = alternative)
result
result %>%
    plot()

     t test power calculation 

             n1 = 40
             n2 = 60
              d = 0.5
      sig.level = 0.05
          power = 0.6792059
    alternative = two.sided
Figure 1

なお、pwr.t2n.test 関数では効果量、有意水準および検定力が与えられた場合、サンプルサイズ n1 および n2 を 2 + 1e-10 から 1e+09 の範囲として求めた検定力と与えられた検定力との差がゼロになる n1 および n2 をニュートン法により求めている。

# n1 <- uniroot(function(n1) eval(p.body) - power, c(2 + 1e-10, 1e+09))$root n2 <- uniroot(function(n2) eval(p.body) - power, c(2 + 1e-10, 1e+09))$root

参考引用資料

最終更新

Sys.time()
[1] "2024-04-20 08:20:52 JST"

R、Quarto、Package

R.Version()$version.string
[1] "R version 4.3.3 (2024-02-29 ucrt)"
quarto::quarto_version()
[1] '1.4.553'
packageVersion(pkg = "tidyverse")
[1] '2.0.0'

著者