JuliaStats / HypothesisTests.jl

Hypothesis tests for Julia
MIT License
295 stars 87 forks source link

Wrong calculation for p-value for FDist #255

Closed KronosTheLate closed 2 years ago

KronosTheLate commented 2 years ago

from f.jl:

function pvalue(x::VarianceFTest; tail=:both)
    dist = FDist(x.df_x, x.df_y)
    if tail == :both
        return 1 - 2*abs(cdf(dist, x.F) - 0.5)
    elseif tail == :right
        return 1 - cdf(dist, x.F)
    elseif tail == :left
        return cdf(dist, x.F)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

How does this implementation account for the fact that a FDist is not symmetic, and that the cdf up to a point of a given probability density is not the same as the cdf after the next point of equal probability density?

Example: I want to test if two distributions have the same variance. I get very different results from the implemented pvalue, and the one I calculated in the way thought in my class:

julia> using Distributions, Roots, HypothesisTests

julia> using Random; Random.seed!(1234);

julia>

julia> obs1 = rand(Normal(0, 1), 10);

julia> obs2 = rand(Normal(0, 1.3), 10);

julia>

julia> x₀² = var(obs1) / var(obs2)     #* Test statistic
1.1337692792442289

julia> d = FDist(length(obs1) - 1, length(obs2) - 1)
FDist{Float64}(ν1=9.0, ν2=9.0)

julia>

julia> f(x) = pdf(d, x) - pdf(d, x₀²)   # Function whose zero's have the ProbDensity of my test statistic
f (generic function with 1 method)

julia> points = find_zeros(f, mean(d) / 1000, mean(d) * 5)#, no_pts = 20)
2-element Vector{Float64}:
 0.3481301077674779
 1.1337692792442273

julia> pdf(d, (points[1])), pdf(d, (points[2]))   # Showing that the points have the same ProbDensity
(0.5043607626219371, 0.5043607626219371)

julia> likelyhood_until_first_point = cdf(d, points[1])
0.06592645343140012

julia> likelyhood_after_second_point = 1 - cdf(d, points[2])
0.4273541291215681

julia> my_pvalue = likelyhood_until_first_point + likelyhood_after_second_point
0.49328058255296825

julia> VarianceFTest(obs1, obs2) |> pvalue
0.8547082582431333

There is a pretty huge difference in some cases, as seen here. Have I misunderstood something, or is the current implementation of pvalue for FDist wrong? If it is wrong, is it different for other non-symmetric distributions?

Code without outputs for easier copying:

using Distributions, Roots, HypothesisTests
using Random; Random.seed!(1234);

obs1 = rand(Normal(0, 1), 10);
obs2 = rand(Normal(0, 1.3), 10);

x₀² = var(obs1) / var(obs2)     #* Test statistic
d = FDist(length(obs1) - 1, length(obs2) - 1)

f(x) = pdf(d, x) - pdf(d, x₀²)   # Function whose zero's have the ProbDensity of my test statistic
points = find_zeros(f, mean(d) / 1000, mean(d) * 5)#, no_pts = 20)
pdf(d, (points[1])), pdf(d, (points[2]))   # Showing that the points have the same ProbDensity
likelyhood_until_first_point = cdf(d, points[1])
likelyhood_after_second_point = 1 - cdf(d, points[2])
my_pvalue = likelyhood_until_first_point + likelyhood_after_second_point

VarianceFTest(obs1, obs2) |> pvalue
KronosTheLate commented 2 years ago

I just check the wiki article: https://en.wikipedia.org/wiki/P-value There, it is stated that image

I can see that this definition is in agreement with the result provided in this package, and thereby close this issue.

This is however a bit if a sticky point, as this stackexchange discussion elaborates on.