Julia に翻訳--070 サンプルサイズが異なる t 検定の検出力

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

サンプルサイズが異なる t 検定の検出力
http://aoki2.si.gunma-u.ac.jp/R/power_t_test2.html

ファイル名: powerttest2.jl　　関数名: powerttest2

翻訳するときに書いたメモ

noncentral t distribution

R:
pt(q, df, ncp, lower.tail=TRUE) == ntdistcdf(df, ncp, q)
pt(1,2,3, lower.tail=TRUE)  # 0.028897773687030438

Julia:
ntdistcdf(df, ncp, q)
ntdistcdf(2, 3, 1)          # 0.028897773687030427

R:
pt(q, df, ncp, lower.tail=FALSE) == ntdistccdf(df, ncp, q)
pt(1,2,3, lower.tail=FALSE) # 0.97110222631296961

Julia:
ntdistccdf(df, ncp, q)
ntdistccdf(2, 3, 1)         # 0.9711022263129696

==========#

using Distributions, StatsFuns

function powerttest2(n1, n2, delta, siglevel=0.05)
    phi = n1+n2-2
    lambda = sqrt(n1*n2/(n1+n2))*delta
    q = cquantile(TDist(phi), siglevel/2)
    return ntdistcdf(phi, lambda, -q) + ntdistccdf(phi, lambda, q)
end

powerttest2(10, 8, 0.6) # 0.2214264393396599

設問

全体で 100 のサンプルを採るとき，
サンプルサイズをどのように配分したら最も検出力が高くなるか。

using Plots
n = 100
delta = 0.6
power = [powerttest2(n1, n-n1, delta) for n1 = 3:n-3]
pyplot()
plot(3:n-3, power, ylabel="Power", xlabel="n1", label="")
vline!([50], label="n1 = n2 = 50")
savefig("fig1.png")

解答

n1 = n2 = 50 にするともっとも検出力が高くなる。

 

Julia の小ネタ--010 非心 t 分布

非心 t 分布　julia ntdistcdf

検索しても，1件しか見つからない。しかも，関数名のリストだけ。

https://github.com/JuliaStats/StatsFuns.jl

で，やってみて，以下のようであるらしいことがわかる。

noncentral t distribution

R での，非心 t 分布（下側確率）

pt(q, df, ncp, lower.tail=TRUE)
pt(1, 2, 3, lower.tail=TRUE)  # 0.028897773687030438

Julia では

ntdistcdf(df, ncp, q)
ntdistcdf(2, 3, 1)            # 0.028897773687030427

R での，非心 t 分布（上側確率）

pt(q, df, ncp, lower.tail=FALSE)
pt(1, 2, 3, lower.tail=FALSE) # 0.97110222631296961

Julia では

ntdistccdf(df, ncp, q)
ntdistccdf(2, 3, 1)           # 0.9711022263129696

Julia に翻訳--069 サンプルサイズが異なる二群の比率の差の検定の検出力

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

サンプルサイズが異なる二群の比率の差の検定の検出力
http://aoki2.si.gunma-u.ac.jp/R/power_prop_test2.html

ファイル名: powerproptest.jl　　関数名: powerproptest2, powerproptest3

翻訳するときに書いたメモ

淡々と

==========#

using Distributions

function powerproptest2(Nc, Nt, Pc, Pt; siglevel=0.05)
    P = (Nc*Pc+Nt*Pt)/(Nc+Nt)
    Zalpha = cquantile(Normal(), siglevel/2)
    Zbeta = (abs(Pc-Pt)-Zalpha*sqrt(P*(1-P)/Nc+P*(1-P)/Nt)) / sqrt(Pc*(1-Pc)/Nc+Pt*(1-Pt)/Nt)
    return cdf(Normal(), Zbeta)
end

function powerproptest3(Pc, Pt; r=1, siglevel=0.05, power=0.8)
    P = (r*Pc+Pt) / (r+1)
    Zalpha = cquantile(Normal(), siglevel/2)
    Zbeta  = cquantile(Normal(), 1-power)
    Nt = (Zalpha*sqrt((r+1)*P*(1-P))+Zbeta*sqrt(r*Pt*(1-Pt)+Pc*(1-Pc)))^2 / r / (Pt-Pc)^2
    Dict(:Nc => Nt*r, :Nt => Nt)
end

powerproptest2(200, 200, 0.6, 0.5) # 0.5200849014176392

powerproptest2(54, 269, 0.5, 0.7)  # 0.8011748762749364

powerproptest3(0.5, 0.7)
# Dict{Symbol,Float64} with 2 entries:
#   :Nt => 92.9988
#   :Nc => 92.9988

powerproptest3(0.5, 0.7, r=1/5)
# Dict{Symbol,Float64} with 2 entries:
#   :Nt => 268.977
#   :Nc => 53.7955

Julia に翻訳--068 平均値の差の検定に必要なサンプルサイズ

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

平均値の差の検定に必要なサンプルサイズ
http://aoki2.si.gunma-u.ac.jp/R/sample-size2.html

ファイル名: samplesize.jl　　関数名: samplesize

翻訳するときに書いたメモ

SpecialFunctions にある関数だけでは書けないようだ。

==========#

using SpecialFunctions, Distributions

function samplesize(alpha, powd, esize; EPS=1e-10)
    n = 0
    powa = 0
    INTV = 200
    dir = -1
    while powa <= powd
        n += 100
        powa = sub1(n, esize, alpha)
    end
    while abs((powa-powd)/powd) >= EPS
        INTV *= 0.5
        n += dir*INTV*0.5
        powa = sub1(n, esize, alpha)
        dir = powa < powd ? 1 : -1
    end
    return n
end

function gcf(a, x; EPS=1e-10, ITMAX=1000)
    FPMIN = 1e-30
    b = x+1-a
    c = 1/FPMIN
    d = 1/b
    h = d
    for i = 1:ITMAX
        an = -i*(i-a)
        b += 2
        d = an*d+b
        if abs(d) < FPMIN
            d = FPMIN
        end
        c = b+an/c
        if abs(c) < FPMIN
            c = FPMIN
        end
        d = 1/d
        del = d*c
        h *= del
        if abs(del-1) < EPS
            return exp(-x+a*log(x)-loggamma(a))*h
        end
    end
    error("not converged in 'gcf'")
end

function gser(a, x; EPS=1e-10, ITMAX=1000) # incomplete gamma function
    if x == 0
        return 0
    elseif x > 0
        ap = a
        del = sm = 1/a
        for n = 1:ITMAX
            ap += 1
            del *= x/ap
            sm += del
            if abs(del) < abs(sm)*EPS
                return sm*exp(-x+a*log(x)-loggamma(a))
            end
        end
    end
    error("not converged in 'gser'")
end

gammp(a, x) = x < a+1 ? gser(a, x) : 1-gcf(a, x)

erff(x) = x < 0 ? -gammp(0.5, x*x) : gammp(0.5, x*x)

function betacf(a, b, x; EPS=1e-10, ITMAX=1000)
    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
    end
    d = 1/d
    h = d
    for m = 1:ITMAX
        m2 = 2*m
        aa = m*(b-m)*x/((qam+m2)*(a+m2))
        d = 1+aa*d
        if abs(d) < FPMIN
            d = FPMIN
        end
        c = 1+aa/c
        if abs(c) < FPMIN
            c = FPMIN
        end
        d = 1/d
        h *= d*c
        aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
        d = 1+aa*d
        if abs(d) < FPMIN
            d = FPMIN
        end
        c = 1+aa/c
        if abs(c) < FPMIN
            c = FPMIN
        end
        d = 1/d
        del = d*c
        h *= del
        if abs(del-1) < EPS
            return h
        end
    end
    stop("not converged in 'betacf'")
end

function betai(a, b, x)
    0 <= x <= 1 || error("0 <= x <= 1")
    bt = x == 0 || x == 1 ? 0 : exp(lbeta(a, b)+a*log(x)+b*log(1-x))
    x < (a+1)/(a+b+2) ? bt*betacf(a, b, x)/a : 1-bt*betacf(b, a, 1-x)/b
end

function sub1(n, esize, alpha)
    df = n-2
    t = cquantile(TDist(df), alpha/2)
    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)))))
end

samplesize(0.05, 0.8, 0.5)

出力結果例
samplesize(0.05, 0.8, 0.5)
[1] 64.84375

EPSILON を小さくして
64.77003814652562
となった

同じことを 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

Julia に翻訳--067 母相関係数の差の検定に必要なサンプルサイズ・検出力

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

母相関係数の差の検定に必要なサンプルサイズ・検出力
http://aoki2.si.gunma-u.ac.jp/R/power-cor-test.html

ファイル名: powercortest.jl　　関数名: powercortest

翻訳するときに書いたメモ

cor0, cor1 を求めるときは，二通り求まる場合がある。

==========#

using Distributions, Roots

function powercortest(; n=NaN, cor0=NaN, cor1=NaN, siglevel=NaN, power=NaN,
                        alternative="two.sided")
    sum(map(isnan, [n, cor0, cor1, siglevel, power])) == 1 ||
        error("n, cor0, cor1, power, siglevel のどれか一つだけを NaN にする")
    alternative in ["one.sided", "two.sided"] ||
        error("'alternative' は 'one.sided' か 'two.sided'")
    tside = alternative=="two.sided" ? 2 : 1

    if isnan(power)
        power = cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
                cquantile(Normal(), siglevel/tside))
    elseif isnan(n)
        n = find_zero(n -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
            cquantile(Normal(), siglevel/tside))-power, [4, 1e7], Bisection())
    elseif isnan(cor0)
        a = []
        try
            b = find_zero(cor0 -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
                cquantile(Normal(), siglevel/tside))-power, [-1, cor1], Bisection())
            append!(a, b)
        catch
            println("No cor0 in [-1, cor1] can be found to achieve the desired power")
        end
        try
            b = find_zero(cor0 -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
                cquantile(Normal(), siglevel/tside))-power, [cor1, 1], Bisection())
            append!(a, b)
        catch
            println("No cor0 in [cor1, 1] can be found to achieve the desired power")
        end
        cor0 = length(a) == 1 ? a[1] : (a[1], a[2])
    elseif isnan(cor1)
        a = []
        try
            b = find_zero(cor1 -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
                cquantile(Normal(), siglevel/tside))-power, [-1, cor0], Bisection())
            append!(a, b)
        catch
            println("No cor1 in [-1, cor0] can be found to achieve the desired power")
        end
        try
            b = find_zero(cor1 -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
                cquantile(Normal(), siglevel/tside))-power, [cor0, 1], Bisection())
            append!(a, b)
        catch
            println("No cor0 in [cor0, 1] can be found to achieve the desired power")
        end
        cor1 = length(a) == 1 ? a[1] : (a[1], a[2])
    elseif isnan(siglevel)
        siglevel = find_zero(siglevel -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
            cquantile(Normal(), siglevel/tside))-power, [1e-5, 0.99999], Bisection())
    else
        error("internal error")
    end
    println("          n = $n")
    println("       cor0 = $cor0")
    println("       cor1 = $cor1")
    println(" sig. level = $siglevel")
    println("      power = $power")
    println("alternative = $alternative")
    Dict(:n => n, :cor0 => cor0, :cor1 => cor1, :siglevel => siglevel,
         :power => power, :alternative => alternative)
end

powercortest(siglevel=0.05, power=0.8, cor0=0.0, cor1=0.4) # find n
#           n = 46.73160799446543
#        cor0 = 0.0
#        cor1 = 0.4
#  sig. level = 0.05
#       power = 0.8
# alternative = two.sided

powercortest(siglevel=0.05, power=0.8, cor0=0.3, cor1=0.4) # find n
#           n = 605.5778584973955
#        cor0 = 0.3
#        cor1 = 0.4
#  sig. level = 0.05
#       power = 0.8
# alternative = two.sided

powercortest(n=605, siglevel=0.05, cor0=0.3, cor1=0.4) # find power
#           n = 605
#        cor0 = 0.3
#        cor1 = 0.4
#  sig. level = 0.05
#       power = 0.7996236163502786
# alternative = two.sided

powercortest(n=605, power=0.8, cor0=0.3, cor1=0.4) # find siglevel
#           n = 605
#        cor0 = 0.3
#        cor1 = 0.4
#  sig. level = 0.050157266442876956
#       power = 0.8
# alternative = two.sided

powercortest(n=605, power=0.8, siglevel=0.05, cor0=0.3) # find cor1
#           n = 605
#        cor0 = 0.3
#        cor1 = (0.19288845143335967, 0.40004600000298063)
#  sig. level = 0.05
#       power = 0.8
# alternative = two.sided

powercortest(n=605, power=0.8, siglevel=0.05, cor1=0.4) # find cor0
#           n = 605
#        cor0 = (0.2999501647530865, 0.4913458910993233)
#        cor1 = 0.4
#  sig. level = 0.05
#       power = 0.8
# alternative = two.sided

Julia に翻訳--066 母平均値の差の検定に必要なサンプルサイズ・検出力

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

母平均値の差の検定に必要なサンプルサイズ・検出力
http://aoki2.si.gunma-u.ac.jp/R/postt.html

ファイル名: poweronesamplettest.jl　　関数名: poweronesamplettest

翻訳するときに書いたメモ

==========#

using Distributions, Roots

function poweronesamplettest(; n=NaN, delta=NaN, sd=NaN, siglevel=NaN,
                    power=NaN, alternative="two.sided")
    sum(map(isnan, [n, delta, sd, power, siglevel])) == 1 ||
        error("n, delta, sd, power, and siglevel のどれか一つだけが NaN でなければならない")
    alternative in ["one.sided", "two.sided"] ||
        error("'alternative' は 'one.sided' か 'two.sided'")
    tside = alternative=="two.sided" ? 2 : 1
    if isnan(power)
        power = cdf(Normal(), (sqrt(n)*delta/sd-cquantile(Normal(), siglevel/tside)))
    elseif isnan(n)
        n = find_zero(n -> cdf(Normal(), (sqrt(n)*delta/sd-cquantile(Normal(), siglevel/tside)))-power, [1, 1e7], Bisection())
    elseif isnan(delta)
        delta = find_zero(delta -> cdf(Normal(), (sqrt(n)*delta/sd-cquantile(Normal(), siglevel/tside)))-power, [-1e7, 1e7], Bisection())
    elseif isnan(sd)
        sd = find_zero(sd -> cdf(Normal(), (sqrt(n)*delta/sd-cquantile(Normal(), siglevel/tside)))-power, [1e-7, 1e7], Bisection())
    elseif isnan(siglevel)
        siglevel = find_zero(siglevel -> cdf(Normal(), (sqrt(n)*delta/sd-cquantile(Normal(), siglevel/tside)))-power, [1e-5, 0.99999], Bisection())
    else
        error("internal error")
    end
    println("      n = $n")
    println("      delta = $delta")
    println("         sd = $sd")
    println(" sig. level = $siglevel")
    println("      power = $power")
    println("alternative = $alternative")
    Dict(:n => n, :delta => delta, :sd => sd, :siglevel => siglevel, :power => power, :alternative => alternative)
end

poweronesamplettest(delta=0.5, sd=sqrt(3), power=0.8, siglevel=0.05, alternative="two.sided") # find n
#       n = 94.18655681218937
#       delta = 0.5
#          sd = 1.7320508075688772
#  sig. level = 0.05
#       power = 0.8
# alternative = two.sided

poweronesamplettest(delta=0.5, sd=sqrt(3), n=46, siglevel=0.05, alternative="two.sided") # find power
#       n = 46
#       delta = 0.5
#          sd = 1.7320508075688772
#  sig. level = 0.05
#       power = 0.49917260874732433
# alternative = two.sided

poweronesamplettest(sd=sqrt(3), n=46, power=0.8, siglevel=0.05, alternative="two.sided") # find delta
#       n = 46
#       delta = 0.7154603140187524
#          sd = 1.7320508075688772
#  sig. level = 0.05
#       power = 0.8
# alternative = two.sided

poweronesamplettest(delta=0.5, n=46, power=0.8, siglevel=0.05, alternative="two.sided") # find sd
#       n = 46
#       delta = 0.5
#          sd = 1.2104450614737239
#  sig. level = 0.05
#       power = 0.8
# alternative = two.sided

poweronesamplettest(delta=0.5, sd=sqrt(3), n=46, power=0.8, alternative="two.sided") # find siglevel
#       n = 46
#       delta = 0.5
#          sd = 1.7320508075688772
#  sig. level = 0.2643070979829777
#       power = 0.8
# alternative = two.sided