Fix numerical issues inside p-value computation

indirect_gwas/src/utils.cpp

@@ -7,10 +7,10 @@
 
 #include <boost/math/distributions/students_t.hpp>
 
-float compute_log_p_value(float t_statistic, unsigned int degrees_of_freedom)
+double compute_log_p_value(double t_statistic, unsigned int degrees_of_freedom)
 {
     boost::math::students_t dist(degrees_of_freedom);
-    float p_value = boost::math::cdf(dist, -std::abs(t_statistic));
+    double p_value = boost::math::cdf(dist, -std::abs(t_statistic));
     return -1. * (std::log(p_value) + std::log(2)) / std::log(10);
 }
 

indirect_gwas/src/utils.hpp

@@ -39,6 +39,6 @@ struct ResultsChunk
     Eigen::VectorXf sample_size;
 };
 
-float compute_log_p_value(float t_statistic, unsigned int degrees_of_freedom);
+double compute_log_p_value(double t_statistic, unsigned int degrees_of_freedom);
 
 std::string get_formatted_time();

indirect_gwas/src/wrapper.cpp

@@ -84,6 +84,7 @@ PYBIND11_MODULE(_igwas, m)
         py::arg("n_covariates"),
         py::arg("chunksize"),
         py::arg("single_file_output"));
-    m.def("compute_pvalue", &compute_log_p_value, "Compute negative log10 p-value");
+    m.def("compute_pvalue", &compute_log_p_value, "Compute negative log10 p-value",
+          py::arg("t_statistic"), py::arg("degrees_of_freedom"));
     m.def("head", &head, "Print the first n lines of a file", py::arg("filename"), py::arg("n_lines"));
 }

indirect_gwas/test/test_pvalues.py

@@ -0,0 +1,31 @@
+import numpy as np
+import pytest
+import scipy.stats
+
+from indirect_gwas._igwas import compute_pvalue
+
+
+@pytest.mark.parametrize(
+    "t_stat, df",
+    [
+        (0.0, 1),
+        (1.0, 1),
+        (1.0, 2),
+        (1e-4, 100_000),
+        (1e-8, 100_000),
+        (10, 100_000),
+        (-10, 100_000),
+        (20, 100_000),
+        (50, 100_000),
+        (100, 100_000),
+        (500, 100_000),
+        (1000, 100_000),
+    ],
+)
+def test_compute_pvalue(t_stat, df):
+    """
+    Test that the p-value is computed correctly. All p-values are negative log10.
+    """
+    cpp_pvalue = compute_pvalue(t_stat, df)
+    py_pvalue = -(scipy.stats.t(df=df).logsf(abs(t_stat)) + np.log(2)) / np.log(10)
+    assert cpp_pvalue == pytest.approx(py_pvalue, rel=1e-6)