JuliaStats / HypothesisTests.jl

Hypothesis tests for Julia
Other
292 stars 87 forks source link

`confint(FisherExactTest(75, 285, 1, 1140))` fails but works if you change 1 to 0 #307

Open skirsch opened 11 months ago

skirsch commented 11 months ago
julia> using HypothesisTests

julia> print(FisherExactTest(75, 285, 1, 1140))
Fisher's exact test
-------------------
Population details:
    parameter of interest:   Odds ratio
    value under h_0:         1.0
    point estimate:          298.973
ERROR: ArgumentError: The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.

Stacktrace:
  [1] assert_bracket
    @ C:\Users\stk\.julia\packages\Roots\uaf5K\src\Bracketing\bracketing.jl:52 [inlined]
  [2] init_state(::Roots.Bisection, F::Roots.Callable_Function{Val{1}, Val{false}, HypothesisTests.var"#obj#14"{Float64, Symbol, FisherExactTest, HypothesisTests.var"#dist#13"{FisherExactTest}}, Nothing}, x₀::Float64, x₁::Float64, fx₀::Float64, fx₁::Float64; m::Float64, fm::Float64)
    @ Roots C:\Users\stk\.julia\packages\Roots\uaf5K\src\Bracketing\bisection.jl:50
  [3] init_state
    @ C:\Users\stk\.julia\packages\Roots\uaf5K\src\Bracketing\bisection.jl:44 [inlined]
  [4] init_state(M::Roots.Bisection, F::Roots.Callable_Function{Val{1}, Val{false}, HypothesisTests.var"#obj#14"{Float64, Symbol, FisherExactTest, HypothesisTests.var"#dist#13"{FisherExactTest}}, Nothing}, x::Tuple{Float64, Float64})
    @ Roots C:\Users\stk\.julia\packages\Roots\uaf5K\src\Bracketing\bracketing.jl:6
  [5] #init#42
    @ C:\Users\stk\.julia\packages\Roots\uaf5K\src\find_zero.jl:299 [inlined]
  [6] solve(𝑭𝑿  ::Roots.ZeroProblem{HypothesisTests.var"#obj#14"{Float64, Symbol, FisherExactTest, HypothesisTests.var"#dist#13"{FisherExactTest}}, Tuple{Float64, Float64}}, M::Roots.Bisection, p::Nothing; verbose::Bool, kwargs::Base.Iterators.Pairs{Symbol, Roots.NullTracks, Tuple{Symbol}, NamedTuple{(:tracks,), Tuple{Roots.NullTracks}}})
    @ Roots C:\Users\stk\.julia\packages\Roots\uaf5K\src\find_zero.jl:491
  [7] find_zero(f::Function, x0::Tuple{Float64, Float64}, M::Roots.Bisection, p′::Nothing; p::Nothing, verbose::Bool, tracks::Roots.NullTracks, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Roots C:\Users\stk\.julia\packages\Roots\uaf5K\src\find_zero.jl:220
  [8] find_zero (repeats 2 times)
    @ C:\Users\stk\.julia\packages\Roots\uaf5K\src\find_zero.jl:220 [inlined]
  [9] find_zero(f::Function, x0::Tuple{Float64, Float64}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Roots C:\Users\stk\.julia\packages\Roots\uaf5K\src\find_zero.jl:243
 [10] find_zero
    @ C:\Users\stk\.julia\packages\Roots\uaf5K\src\find_zero.jl:243 [inlined]
 [11] confint(x::FisherExactTest; level::Float64, tail::Symbol, method::Symbol)
    @ HypothesisTests C:\Users\stk\.julia\packages\HypothesisTests\r322N\src\fisher.jl:195
 [12] confint(x::FisherExactTest; level::Float64, tail::Symbol, method::Symbol)
    @ HypothesisTests C:\Users\stk\.julia\packages\HypothesisTests\r322N\src\fisher.jl:206
 [13] confint
    @ C:\Users\stk\.julia\packages\HypothesisTests\r322N\src\fisher.jl:183 [inlined]
 [14] show(_io::Base.TTY, test::FisherExactTest)
    @ HypothesisTests C:\Users\stk\.julia\packages\HypothesisTests\r322N\src\HypothesisTests.jl:73
 [15] print(io::Base.TTY, x::FisherExactTest)
    @ Base .\strings\io.jl:35
 [16] print(xs::FisherExactTest)
    @ Base .\coreio.jl:3
 [17] top-level scope
    @ REPL[3]:1

julia>
skirsch commented 11 months ago

If you replace 1 with a 0, there is no error:

julia> print(FisherExactTest(75, 285, 0, 1140))
Fisher's exact test
-------------------
Population details:
    parameter of interest:   Odds ratio
    value under h_0:         1.0
    point estimate:          Inf
    95% confidence interval: (78.46, Inf)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           <1e-48

Details:
    contingency table:
        75   285
         0  1140

julia>
ararslan commented 11 months ago

Huh, this is a fun edge case. The issue here is coming from confint. Mimicking what's happening inside confint, we start with the following setup:

using HypothesisTests, Roots, Distributions
using HypothesisTests: find_brackets

test = FisherExactTest(75, 285, 1, 1140)
distribution((; a, b, c, d), ω) =
    FisherNoncentralHypergeometric(a + b, c + d, a + c, ω)
objective(test, ω, tail) = pvalue(distribution(test, ω), test.a; tail) - 0.05

The default tail is :both, which actually calls confint(...; tail=:left) and confint(...; tail=:right) to construct the limits. Those in turn find a bracketing interval for the zero of objective (the odds ratio such that a p-value from the specified tail of Fisher's noncentral hypergeometric distribution is 0.05) and pass them to the root finding function. This then begins with

find_brackets(ω -> objective(test, ω, :left))

which produces an interval of (5.0e-324, 1.0). However, that... isn't a bracketing interval for the optimum: the objective function doesn't actually have any roots in that interval. Now, if you try again with test = FisherExactTest(75, 285, 0, 1140), you actually get the same interval. Then why does one work and the other does not? It's because we explicitly handle the case where the first cell in the 2x2 table is equal to either extremum of the distribution for an odds ratio of 1. But when the cell is 1 instead of 0, we're 1 away from an extremum, so we don't hit the special case.

I suspect this is a problem with find_brackets.