FR: Dunn-Šidák

src/FSharp.Stats/Testing/MultipleTesting.fs

@@ -63,6 +63,16 @@ 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
+    // Old: /// Computes the Dunn-Šidák correction onto a collection of p-values with a given alpha. Returns a list of booleans which state if the p-value should be rejected or not, and the new critical value for this alpha.
+    /// Computes the Dunn-Šidák correction onto a collection of p-values with a given alpha. Returns a critical value. p-values above the critical value should be rejected.
+    let dunnSidak alpha (pValues : seq<float>) = 
+        let m = float (Seq.length pValues)
+        let criticalValue = 1. - (1. - alpha) ** (1. / m)
+        // List.map (fun p -> if p > criticalValue then true, p else false, p) (List.ofSeq pValues),
+        criticalValue
+
+
     // John D. Storey
     // http://dldcc-web.brc.bcm.edu/lilab/liguow/CGI/R/library/qvalue/html/qvalue.html
     // http://qvalue.princeton.edu/

tests/FSharp.Stats.Tests/Testing.fs

@@ -348,6 +348,23 @@ let benjaminiHochbergTests =
     ]
 
 
+[<Tests>]
+let dunnSidakTest =
+    // taken from: 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."
+    ]
+
 
 [<Tests>]
 let qValuesTest =