diff --git a/src/PowerAnalyses.jl b/src/PowerAnalyses.jl index 08f8f44..d93fad8 100644 --- a/src/PowerAnalyses.jl +++ b/src/PowerAnalyses.jl @@ -6,7 +6,9 @@ using Distributions: NoncentralT, UnivariateDistribution, cdf, - quantile + ccdf, + quantile, + cquantile using Roots: find_zero include("types.jl") diff --git a/src/power.jl b/src/power.jl index 63a1974..7763fc8 100644 --- a/src/power.jl +++ b/src/power.jl @@ -5,10 +5,13 @@ Return `power` for which `quantile(d0, alpha) == quantile(d1, beta)` and `beta = """ function _power(d0::UnivariateDistribution, d1::UnivariateDistribution, alpha::Real, tail::Tail) right_tail = tail == one_tail ? alpha : alpha / 2 - except_right_tail = 1 - right_tail - critical_value = quantile(d0, except_right_tail) - beta = cdf(d1, critical_value) - return 1 - beta + upper_critical_value = cquantile(d0, right_tail) + if tail == one_tail + return ccdf(d1, upper_critical_value) + else + lower_critical_value = quantile(d0, right_tail) + return cdf(d1, lower_critical_value) + ccdf(d1, upper_critical_value) + end end """ @@ -17,11 +20,9 @@ end Return `alpha` for which `quantile(d0, alpha) == quantile(d1, beta)` and `beta = 1 - power`. """ function _alpha(d0::UnivariateDistribution, d1::UnivariateDistribution, power::Real, tail::Tail) - beta = 1 - power - critical_value = quantile(d1, beta) - except_right_tail = cdf(d0, critical_value) - right_tail = 1 - except_right_tail - return tail == one_tail ? right_tail : 2 * right_tail + return find_zero((0, 1)) do alpha + return _power(d0, d1, alpha, tail) - power + end end _distribution_parameters(T::IndependentSamplesTTest, n) = 2 * n - 2 # n1 + n2 - 2 @@ -122,9 +123,10 @@ end Return the minimum effect size for some test `T` with significance level `alpha`, power `power` and sample size `n`. """ function get_es(T::StatisticalTest; alpha::Real, power::Real, n) - f(es) = get_alpha(T; es, power, n) - alpha - initial_value = (0, 10) - return find_zero(f, initial_value) + initial_value = (0, 1000) + return find_zero(initial_value) do es + get_power(T; es, alpha, n) - power + end end """ diff --git a/test/runtests.jl b/test/runtests.jl index 1d49d7a..31d335d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,14 +22,15 @@ using Test: @testset, @test @testset "IndependentSamplesTTest" begin @testset "two_tails" begin - # This matches power.t.test(n=50, d=0.5, sig.level=NULL, power=0.95, type="two.sample", alternative="two.sided") - # but not pwr.t.test(n=50, d=0.5, sig.level=NULL, power=0.95, type="two.sample", alternative="two.sided") - # because the latter considers both sides of the distribution. A fix will probably require - # switching to root finding. The pwr.t.test matches # G*Power and is 0.3928954 - @test get_alpha(IndependentSamplesTTest(two_tails); es, power, n) ≈ 0.3950408 rtol=1e-6 - @test get_power(IndependentSamplesTTest(two_tails); es, alpha, n) ≈ 0.6968888 rtol=1e-6 + # This matches pwr.t.test(n=50, d=0.5, sig.level=NULL, power=0.95, type="two.sample", alternative="two.sided") + # but not power.t.test(n=50, d=0.5, sig.level=NULL, power=0.95, type="two.sample", alternative="two.sided") + # because the latter considers both sides of the distribution. The pwr.t.test matches # G*Power and is 0.3928954 + @test get_alpha(IndependentSamplesTTest(two_tails); es, power, n) ≈ 0.3928954 rtol=1e-6 + @test get_power(IndependentSamplesTTest(two_tails); es, alpha, n) ≈ 0.6968934 rtol=1e-6 @test get_es(IndependentSamplesTTest(two_tails); alpha, power, n) ≈ 0.7281209 rtol=1e-4 - @test get_n(IndependentSamplesTTest(two_tails); es, alpha, power) ≈ 104.928 rtol=1e-6 + @test get_n(IndependentSamplesTTest(two_tails); es, alpha, power) ≈ 104.9279 rtol=1e-6 + + @test get_power(IndependentSamplesTTest(two_tails); es = 0.0, alpha, n) ≈ alpha end @testset "one_tail" begin # R gives a warning that full precision might not have been achieved