add Holm and Hochberg FWER correction

Author: ZimmerD

src/FSharp.Stats/Testing/MultipleTesting.fs

@@ -1,11 +1,83 @@
 namespace FSharp.Stats.Testing
 
+open System
 open FSharp.Stats
 
 /// This module contains functions to adjust for multiple testing errors in statistical tests.
 module MultipleTesting = 
 
 
+    /// Holm–Bonferroni (step-down) FWER adjustment (NaN-safe)
+    /// NaN p-values are ignored in the computation and preserved at their original indices.
+    let inline holmFWER (p : float[]) : float[] =
+        let m = p.Length
+        let indexed = p |> Array.mapi (fun i v -> i, v)
+
+        // keep only finite p-values for the adjustment
+        let valid  = indexed |> Array.filter (fun (_, v) -> not (Double.IsNaN v))
+        let nValid = valid.Length
+
+        // start with an all-NaN result and fill only valid positions
+        let result = Array.create m Double.NaN
+
+        if nValid = 0 then result else
+            let sorted = valid |> Array.sortBy snd
+
+            // (nValid - i)·p_(i) for i = 0…nValid-1
+            let raw =
+                sorted
+                |> Array.mapi (fun i (_, pv) -> float (nValid - i) * pv)
+
+            // running max from the left, capped at 1.0
+            let adjAsc =
+                raw
+                |> Array.scan (fun runningMax r -> max runningMax r) 0.0
+                |> Array.tail
+                |> Array.map (min 1.0)
+
+            // write adjusted values back into their original indices
+            Array.zip (sorted |> Array.map fst) adjAsc
+            |> Array.iter (fun (i, adj) -> result.[i] <- adj)
+
+            result
+
+    /// Hochberg (step-up) FWER adjustment (NaN-safe)
+    /// NaN p-values are ignored in the computation and preserved at their original indices.
+    let inline hochbergFWER (p : float[]) : float[] =
+        let m = p.Length
+        let indexed = p |> Array.mapi (fun i v -> i, v)
+
+        // keep only finite p-values for the adjustment
+        let valid  = indexed |> Array.filter (fun (_, v) -> not (Double.IsNaN v))
+        let nValid = valid.Length
+
+        // start with an all-NaN result and fill only valid positions
+        let result = Array.create m Double.NaN
+
+        if nValid = 0 then result else
+            let sorted = valid |> Array.sortBy snd
+
+            // (nValid - i)·p_(i), same raw values pattern as Holm
+            let raw =
+                sorted
+                |> Array.mapi (fun i (_, pv) -> float (nValid - i) * pv)
+
+            // running min from the right, capped at 1.0
+            let adjDesc =
+                raw
+                |> Array.rev
+                |> Array.scan (fun runningMin r -> min runningMin r) 1.0
+                |> Array.tail
+                |> Array.rev
+                |> Array.map (min 1.0)
+
+            // write adjusted values back into their original indices
+            Array.zip (sorted |> Array.map fst) adjDesc
+            |> Array.iter (fun (i, adj) -> result.[i] <- adj)
+
+            result
+
+
     /// Benjamini-Hochberg Correction (BH)
     /// 'projection' should return a tuple of any identifier and the pValues as float, when applied to 'rawP'
     /// This function applies the Benjamini-Hochberg multiple testing correcture and returns all False Discovery Rates to which the given p-values are still 

tests/FSharp.Stats.Tests/FSharp.Stats.Tests.fsproj

@@ -62,6 +62,9 @@
     <EmbeddedResource Include="data/SparseMatrixFormat1NoInfo.txt" CopyToOutputDirectory="PreserveNewest" />
     <EmbeddedResource Include="data/SparseMatrixFormat1WithInfo.txt" CopyToOutputDirectory="PreserveNewest" />
     <EmbeddedResource Include="data/TableFormat.txt" CopyToOutputDirectory="PreserveNewest" />
+    <EmbeddedResource Include="data/fwer_hochberg_results.csv" CopyToOutputDirectory="PreserveNewest" />
+    <EmbeddedResource Include="data/fwer_holm_results.csv" CopyToOutputDirectory="PreserveNewest" />
+    <EmbeddedResource Include="data/holmHochberg_Input_nan.csv" CopyToOutputDirectory="PreserveNewest" />
   </ItemGroup>
   <ItemGroup />
   <ItemGroup>

tests/FSharp.Stats.Tests/TestExtensions.fs

@@ -66,7 +66,6 @@
             | [|a;b|] -> a, float b
             | _ -> failwith "invalid csv format"
          )
-
 
     let comparisonMetricsEqualRounded (digits : int) (actual: ComparisonMetrics) (expected: ComparisonMetrics) message =
         let actual = 

tests/FSharp.Stats.Tests/Testing.fs

@@ -403,6 +403,88 @@ let pearsonTests =
     ]
 
 
+[<Tests>]
+let holmTests =
+
+
+    let largeSetnan        = 
+        Frame.ReadCsv(location = @"data/holmHochberg_Input_nan.csv",hasHeaders = true,separators = ",").GetColumn<float>("pValues")
+        |> Series.valuesAll
+        |> Array.ofSeq
+        |> Array.map (fun x -> if x.IsSome then x.Value else nan ) 
+ 
+    let largeSet        = 
+        largeSetnan |> Array.filter (fun x -> not (nan.Equals x))
+
+    let largeSet_Expectednan        = 
+        Frame.ReadCsv(location = @"data/fwer_holm_results.csv",hasHeaders = true,separators = ",").GetColumn<float>("pValues")
+        |> Series.valuesAll
+        |> Array.ofSeq
+        |> Array.map (fun x -> if x.IsSome then x.Value else nan ) 
+
+    let largeSet_Expected        = 
+        largeSet_Expectednan 
+        |> Array.filter (fun x -> not (nan.Equals x))
+
+    testList "Testing.MultipleTesting.holmFWER" [
+        
+        testCase "testHolmLarge" (fun () -> 
+            Expect.sequenceEqual 
+                (largeSet |> MultipleTesting.holmFWER |> Seq.map (fun x -> Math.Round(x,9))) 
+                (largeSet_Expected |> Seq.map (fun x -> Math.Round(x,9)))
+                "adjusted pValues should be equal to the reference implementation."
+        )
+
+        testCase "testHolmLargeNaN" (fun () -> 
+            TestExtensions.sequenceEqualRoundedNaN 9
+                (largeSetnan |> MultipleTesting.holmFWER |> Seq.ofArray) 
+                (largeSet_Expectednan |> Seq.ofArray)
+                "adjusted pValues should be equal to the reference implementation."
+        )
+
+    ]
+
+[<Tests>]
+let hochbergTests =
+
+
+    let largeSetnan        = 
+        Frame.ReadCsv(location = @"data/holmHochberg_Input_nan.csv",hasHeaders = true,separators = ",").GetColumn<float>("pValues")
+        |> Series.valuesAll
+        |> Array.ofSeq
+        |> Array.map (fun x -> if x.IsSome then x.Value else nan ) 
+ 
+    let largeSet        = 
+        largeSetnan |> Array.filter (fun x -> not (nan.Equals x))
+
+    let largeSet_Expectednan        = 
+        Frame.ReadCsv(location = @"data/fwer_hochberg_results.csv",hasHeaders = true,separators = ",").GetColumn<float>("pValues")
+        |> Series.valuesAll
+        |> Array.ofSeq
+        |> Array.map (fun x -> if x.IsSome then x.Value else nan ) 
+
+    let largeSet_Expected        = 
+        largeSet_Expectednan 
+        |> Array.filter (fun x -> not (nan.Equals x))
+
+    testList "Testing.MultipleTesting.hochbergFWER" [
+        
+        testCase "testHochbergLarge" (fun () -> 
+            Expect.sequenceEqual 
+                (largeSet |> MultipleTesting.hochbergFWER |> Seq.map (fun x -> Math.Round(x,9))) 
+                (largeSet_Expected |> Seq.map (fun x -> Math.Round(x,9)))
+                "adjusted pValues should be equal to the reference implementation."
+        )
+
+        testCase "testHochbergLargeNaN" (fun () -> 
+            TestExtensions.sequenceEqualRoundedNaN 9
+                (largeSetnan |> MultipleTesting.hochbergFWER |> Seq.ofArray) 
+                (largeSet_Expectednan |> Seq.ofArray)
+                "adjusted pValues should be equal to the reference implementation."
+        )
+
+    ]
+
 [<Tests>]
 let benjaminiHochbergTests =