目的

二群の平均値の差（両側検定）を行うときに必要な各群あたりのサンプルサイズを求める。
R には，power.t.test という関数が用意されている。精度的には power.t.test を使う方がよい。

使用法

sample.size(alpha, powd, esize)

引数

alpha	有意水準（α）
powd	望む検出力（1-β）
esize	効果量

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/sample_size.R", encoding="euc-jp")

# 二群の平均値の差（両側検定）を行うときに必要な各群あたりのサンプルサイズを求める
sample.size <- function(alpha,                                       # 有意水準
                        powd,                                   # 検出力
                        esize)                                  # 効果量
{

        gcf <- function(a, x)
        {
                ITMAX <- 100
                EPS <- 3e-7
                FPMIN <- 1e-30
        
                b <- x+1-a
                c <- 1/FPMIN
                d <- 1/b
                h <- d
                for (i in 1:ITMAX) {
                        an <- -i*(i-a)
                        b <- b+ 2
                        d <- an*d+b
                        if (abs(d) < FPMIN) d <- FPMIN
                        c <- b+an/c
                        if (abs(c) < FPMIN) c <- FPMIN
                        d <- 1/d
                        del <- d*c
                        h <- h*del
                        if (abs(del-1) < EPS) {
                                return(exp(-x+a*log(x)-lgamma(a))*h)
                        }
                }
                stop("error")
        }

        gser <- function(a, x)                                       # 不完全ガンマ関数
        {
                ITMAX <- 100
                EPS <- 3e-7
        
                if (x == 0) {
                        return(0)
                }
                else if (x > 0) {
                        ap <- a
                        del <- sum <- 1/a
                        for (n in 1:ITMAX) {
                                ap <- ap+1
                                del <- del*x/ap
                                sum <- sum+del
                                if (abs(del) < abs(sum)*EPS) {
                                        return(sum*exp(-x+a*log(x)-lgamma(a)))
                                }
                        }
                }
                stop("error")
        }

        gammp <- function(a, x)                                      # pgamma(x, a)
        {
                ifelse(x < a+1, gser(a, x), 1-gcf(a, x))
        }

        erff <- function(x)                                  # 1-2*pnorm(x*sqrt(2), lower.tail=FALSE)
        {
                ifelse(x < 0, -gammp(0.5, x*x), gammp(0.5, x*x))
        }

        betacf <- function(a, b, x)
        {
                ITMAX <- 100
                EPS <- 3e-7
                FPMIN <- 1e-30
        
                qab <- a+b
                qap <- a+1
                qam <- a-1
                c <- 1
                d <- 1-qab*x/qap
                if (abs(d) < FPMIN) d <- FPMIN
                d <- 1/d
                h <- d
                for (m in 1:ITMAX) {
                        m2 <- 2*m
                        aa <- m*(b-m)*x/((qam+m2)*(a+m2))
                        d <- 1+aa*d
                        if (abs(d) < FPMIN) d <- FPMIN
                        c <- 1+aa/c
                        if (abs(c) < FPMIN) c <- FPMIN
                        d <- 1/d
                        h <- h*d*c
                        aa <- -(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
                        d <- 1+aa*d
                        if (abs(d) < FPMIN) d <- FPMIN
                        c <- 1+aa/c
                        if (abs(c) < FPMIN) c <- FPMIN
                        d <- 1/d
                        del <- d*c
                        h <- h*del
                        if (abs(del-1) < EPS) return(h)
                }
                stop("error")
        }

        betai <- function(a, b, x)                           # 不完全ベータ関数
        {
                if (x < 0 || x > 1) stop("error")
                bt <- ifelse(x == 0 || x == 1, 0, exp(lbeta(a, b)+a*log(x)+b*log(1-x)))
                ifelse(x < (a+1)/(a+b+2), bt*betacf(a, b, x)/a, 1-bt*betacf(b, a, 1-x)/b)
        }

        sub1 <- function(n, esize, alpha)
        {
                df <- n-2
                t <- qt(alpha/2, df, lower.tail=FALSE)
                dd <- 1-0.25/df+1/(32*df*df)
                0.5*(1+erff((esize*sqrt(n)-sqrt(2)*t*dd) / (2*sqrt(1+t*t*(1-dd*dd)))))
        }

# 関数本体

        n <- 0
        powa <- 0
        INTV <- 200
        EPS <- 0.001
        dir <- -1

        while (powa <= powd) {
                n <- n+100
                powa <- sub1(n, esize, alpha)
        }

        while (abs((powa-powd)/powd) >= EPS) {
                INTV <- INTV*0.5
                n <- n+dir*INTV*0.5
                powa <- sub1(n, esize, alpha)
                dir <- ifelse(powa < powd, 1, -1)
        }
        return(n)
}


使用例

sample.size(0.05, 0.8, 0.5)

出力結果例

> sample.size(0.05, 0.8, 0.5)
[1] 64.84375

同じことを power.t.test でやってみる。

> power.t.test(sig.level=0.05, power=0.8, delta=0.5)

     Two-sample t test power calculation 

              n = 63.76576
          delta = 0.5
             sd = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

 NOTE: n is number in *each* group