diff --git a/src/FSharp.Stats/Testing/MultipleTesting.fs b/src/FSharp.Stats/Testing/MultipleTesting.fs index 0a65bee8..308aefd2 100644 --- a/src/FSharp.Stats/Testing/MultipleTesting.fs +++ b/src/FSharp.Stats/Testing/MultipleTesting.fs @@ -63,6 +63,19 @@ module MultipleTesting = |> Seq.ofList + // adapted from: https://en.wikipedia.org/wiki/%C5%A0id%C3%A1k_correction & https://personal.utdallas.edu/~herve/Abdi-Bonferroni2007-pretty.pdf + // checked according to: https://cran.r-project.org/web/packages/aggregation/index.html + // original derivation: https://www.tandfonline.com/doi/abs/10.1080/01621459.1967.10482935 + /// Computes the Dunn-Sidak correction onto a collection of p-values with a given alpha. + /// Returns a critical value. + /// p-values above the critical value should be rejected. Use this correction to reduce family-wise error rate when performing multiple tests on a dataset. p-values above 1, below 0, inf, and nan are excluded. + let dunnSidak alpha (pValues : seq) = + let checkedValues = Seq.filter (fun v -> v >= 0. && v <= 1.) pValues + let m = float (Seq.length checkedValues) + let criticalValue = 1. - (1. - alpha) ** (1. / m) + criticalValue + + // John D. Storey // http://dldcc-web.brc.bcm.edu/lilab/liguow/CGI/R/library/qvalue/html/qvalue.html // http://qvalue.princeton.edu/ diff --git a/tests/FSharp.Stats.Tests/Testing.fs b/tests/FSharp.Stats.Tests/Testing.fs index a3d475b9..c088a3cc 100644 --- a/tests/FSharp.Stats.Tests/Testing.fs +++ b/tests/FSharp.Stats.Tests/Testing.fs @@ -348,6 +348,23 @@ let benjaminiHochbergTests = ] +[] +let dunnSidakTest = + // taken from: https://web.archive.org/web/20200723122809/https://www.real-statistics.com/hypothesis-testing/familywise-error/bonferroni-and-dunn-sidak-tests/ + let testPvalues = [0.01208; 0.00356; 0.11542; 0.02155; 0.03329; 0.01042] + let testAlpha1 = 0.05 + let testAlpha2 = 0.01 + let expectedValue1 = 0.008512 + let expectedValue2 = 0.001674 + let observedValue1 = MultipleTesting.dunnSidak testAlpha1 testPvalues |> round 6 + let observedValue2 = MultipleTesting.dunnSidak testAlpha2 testPvalues |> round 6 + testList "Testing.MultipleTesting.dunnSidak" [ + testCase "Alpha = 5 %" <| fun () -> + Expect.floatClose Accuracy.veryHigh observedValue1 expectedValue1 "corrected p values are about equal." + testCase "Alpha = 1 %" <| fun () -> + Expect.floatClose Accuracy.veryHigh observedValue2 expectedValue2 "corrected p values are about equal." + ] + [] let qValuesTest =