fixed pairwise unit tests for null model

Author: davidsebfischer

diffxpy/testing/tests.py

@@ -848,6 +848,8 @@ def two_sample(
             sample_description=sample_description,
             noise_model=noise_model,
             size_factors=size_factors,
+            init_a="closed_form",
+            init_b="closed_form",
             batch_size=batch_size,
             training_strategy=training_strategy,
             quick_scale=quick_scale,
@@ -872,6 +874,8 @@ def two_sample(
             sample_description=sample_description,
             noise_model=noise_model,
             size_factors=size_factors,
+            init_a="closed_form",
+            init_b="closed_form",
             batch_size=batch_size,
             training_strategy=training_strategy,
             quick_scale=quick_scale,
@@ -911,7 +915,7 @@ def pairwise(
         batch_size: int = None,
         training_strategy: Union[str, List[Dict[str, object]], Callable] = "AUTO",
         is_sig_zerovar: bool = True,
-        quick_scale: bool = None,
+        quick_scale: bool = False,
         dtype="float64",
         pval_correction: str = "global",
         keep_full_test_objs: bool = False,
@@ -1082,7 +1086,7 @@ def pairwise(
                     grouping == g2
                 ))[0]
                 de_test_temp = two_sample(
-                    data=data[idx],
+                    data=data[idx, :],
                     grouping=grouping[idx],
                     as_numeric=as_numeric,
                     test=test,
@@ -1100,7 +1104,7 @@ def pairwise(
                 pvals[i, j] = de_test_temp.pval
                 pvals[j, i] = pvals[i, j]
                 logfc[i, j] = de_test_temp.log_fold_change()
-                logfc[j, i] = - logfc[i, j]
+                logfc[j, i] = -logfc[i, j]
                 if keep_full_test_objs:
                     tests[i, j] = de_test_temp
                     tests[j, i] = de_test_temp

diffxpy/unit_test/test_pairwise.py

@@ -4,265 +4,149 @@
 import pandas as pd
 import scipy.stats as stats
 
-from batchglm.api.models.glm_nb import Simulator
 import diffxpy.api as de
 
 
-class TestPairwiseNull(unittest.TestCase):
+class _TestPairwiseNull:
 
-    def test_null_distribution_ztest(self, n_cells: int = 2000, n_genes: int = 100, n_groups=2):
-        """
-        Test if de.wald() generates a uniform p-value distribution
-        if it is given data simulated based on the null model. Returns the p-value
-        of the two-side Kolmgorov-Smirnov test for equality of the observed 
-        p-value distriubution and a uniform distribution.
+    noise_model: str
 
-        :param n_cells: Number of cells to simulate (number of observations per test).
-        :param n_genes: Number of genes to simulate (number of tests).
-        """
-        logging.getLogger("tensorflow").setLevel(logging.ERROR)
-        logging.getLogger("batchglm").setLevel(logging.WARNING)
-        logging.getLogger("diffxpy").setLevel(logging.WARNING)
+    def _prepate_data(
+            self,
+            n_cells: int,
+            n_genes: int,
+            n_groups: int
+    ):
+        if self.noise_model == "nb":
+            from batchglm.api.models.glm_nb import Simulator
+            rand_fn_loc = lambda shape: np.random.uniform(0.1, 1, shape)
+            rand_fn_scale = lambda shape: np.random.uniform(0.5, 1, shape)
+        elif self.noise_model == "norm" or self.noise_model is None:
+            from batchglm.api.models.glm_norm import Simulator
+            rand_fn_loc = lambda shape: np.random.uniform(500, 1000, shape)
+            rand_fn_scale = lambda shape: np.random.uniform(1, 2, shape)
+        else:
+            raise ValueError("noise model %s not recognized" % self.noise_model)
 
         sim = Simulator(num_observations=n_cells, num_features=n_genes)
         sim.generate_sample_description(num_batches=0, num_conditions=0)
-        sim.generate()
+        sim.generate_params(
+            rand_fn_loc=rand_fn_loc,
+            rand_fn_scale=rand_fn_scale
+        )
+        sim.generate_data()
 
         random_sample_description = pd.DataFrame({
             "condition": np.random.randint(n_groups, size=sim.nobs)
         })
-
-        test = de.test.pairwise(
-            data=sim.x,
-            grouping="condition",
-            test="z-test",
-            noise_model="nb",
-            sample_description=random_sample_description,
-            dtype="float64"
-        )
-        summary = test.summary()
-
-        # Compare p-value distribution under null model against uniform distribution.
-        pval_h0 = stats.kstest(test.pval[~np.eye(test.pval.shape[0]).astype(bool)].flatten(), 'uniform').pvalue
-
-        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
-        assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)
-
-        return True
-
-    def test_null_distribution_z_lazy(self, n_cells: int = 2000, n_genes: int = 100):
+        return sim, random_sample_description
+
+    def _test_null_distribution_basic(
+            self,
+            test: str,
+            lazy: bool,
+            quick_scale: bool = False,
+            n_cells: int = 3000,
+            n_genes: int = 200,
+            n_groups: int = 3
+    ):
         """
-        Test if de.pairwise() generates a uniform p-value distribution for lazy z-tests
+        Test if de.wald() generates a uniform p-value distribution
         if it is given data simulated based on the null model. Returns the p-value
-        of the two-side Kolmgorov-Smirnov test for equality of the observed
+        of the two-side Kolmgorov-Smirnov test for equality of the observed 
         p-value distriubution and a uniform distribution.
 
         :param n_cells: Number of cells to simulate (number of observations per test).
         :param n_genes: Number of genes to simulate (number of tests).
         """
-        logging.getLogger("tensorflow").setLevel(logging.ERROR)
-        logging.getLogger("batchglm").setLevel(logging.WARNING)
-        logging.getLogger("diffxpy").setLevel(logging.WARNING)
-
-        sim = Simulator(num_observations=n_cells, num_features=n_genes)
-        sim.generate_sample_description(num_batches=0, num_conditions=0)
-        sim.generate()
-
-        random_sample_description = pd.DataFrame({
-            "condition": np.random.randint(4, size=sim.nobs)
-        })
-
+        sim, sample_description = self._prepate_data(
+            n_cells=n_cells,
+            n_genes=n_genes,
+            n_groups=n_groups
+        )
         test = de.test.pairwise(
-            data=sim.x,
+            data=sim.input_data,
+            sample_description=sample_description,
             grouping="condition",
-            test='z-test',
-            lazy=True,
-            noise_model="nb",
-            pval_correction="global",
-            quick_scale=True,
-            sample_description=random_sample_description,
-            dtype="float64"
+            test=test,
+            lazy=lazy,
+            quick_scale=quick_scale,
+            noise_model=self.noise_model
         )
+        _ = test.summary()
 
         # Compare p-value distribution under null model against uniform distribution.
-        pvals = test.pval_pairs(groups0=0, groups1=1)
-        pval_h0 = stats.kstest(pvals.flatten(), 'uniform').pvalue
+        if lazy:
+            pval_h0 = stats.kstest(test.pval_pairs(groups0=0, groups1=1).flatten(), 'uniform').pvalue
+        else:
+            pval_h0 = stats.kstest(test.pval[0, 1, :].flatten(), 'uniform').pvalue
 
         logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
         assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)
-
         return True
 
-    def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 100, n_groups=2):
-        """
-        Test if de.wald() generates a uniform p-value distribution
-        if it is given data simulated based on the null model. Returns the p-value
-        of the two-side Kolmgorov-Smirnov test for equality of the observed
-        p-value distriubution and a uniform distribution.
 
-        :param n_cells: Number of cells to simulate (number of observations per test).
-        :param n_genes: Number of genes to simulate (number of tests).
-        """
+class TestPairwiseNullStandard(unittest.TestCase, _TestPairwiseNull):
+
+    def test_null_distribution_ttest(self):
         logging.getLogger("tensorflow").setLevel(logging.ERROR)
         logging.getLogger("batchglm").setLevel(logging.WARNING)
         logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
-        sim = Simulator(num_observations=n_cells, num_features=n_genes)
-        sim.generate_sample_description(num_batches=0, num_conditions=0)
-        sim.generate()
-
-        random_sample_description = pd.DataFrame({
-            "condition": np.random.randint(n_groups, size=sim.nobs)
-        })
-
-        test = de.test.pairwise(
-            data=sim.x,
-            grouping="condition",
-            test="lrt",
-            noise_model="nb",
-            sample_description=random_sample_description,
-            dtype="float64"
-        )
-
-        # Compare p-value distribution under null model against uniform distribution.
-        pval_h0 = stats.kstest(test.pval[~np.eye(test.pval.shape[0]).astype(bool)].flatten(), 'uniform').pvalue
+        np.random.seed(1)
+        self.noise_model = None
+        self._test_null_distribution_basic(test="t-test", lazy=False)
 
-        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
-        assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)
-
-        return True
-
-    def test_null_distribution_ttest(self, n_cells: int = 2000, n_genes: int = 10000, n_groups=2):
-        """
-        Test if de.wald() generates a uniform p-value distribution
-        if it is given data simulated based on the null model. Returns the p-value
-        of the two-side Kolmgorov-Smirnov test for equality of the observed
-        p-value distriubution and a uniform distribution.
-
-        :param n_cells: Number of cells to simulate (number of observations per test).
-        :param n_genes: Number of genes to simulate (number of tests).
-        """
+    def test_null_distribution_rank(self):
         logging.getLogger("tensorflow").setLevel(logging.ERROR)
         logging.getLogger("batchglm").setLevel(logging.WARNING)
         logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
-        sim = Simulator(num_observations=n_cells, num_features=n_genes)
-        sim.generate_sample_description(num_batches=0, num_conditions=0)
-        sim.generate()
-
-        random_sample_description = pd.DataFrame({
-            "condition": np.random.randint(n_groups, size=sim.nobs)
-        })
-
-        test = de.test.pairwise(
-            data=sim.x,
-            grouping="condition",
-            test="t-test",
-            sample_description=random_sample_description,
-        )
-        summary = test.summary()
-
-        # Compare p-value distribution under null model against uniform distribution.
-        pval_h0 = stats.kstest(test.pval[~np.eye(test.pval.shape[0]).astype(bool)].flatten(), 'uniform').pvalue
-
-        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
-        assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)
+        np.random.seed(1)
+        self.noise_model = None
+        self._test_null_distribution_basic(test="rank", lazy=False)
 
-        return True
 
-    def test_null_distribution_wilcoxon(self, n_cells: int = 2000, n_genes: int = 10000, n_groups=2):
-        """
-        Test if de.wald() generates a uniform p-value distribution
-        if it is given data simulated based on the null model. Returns the p-value
-        of the two-side Kolmgorov-Smirnov test for equality of the observed
-        p-value distriubution and a uniform distribution.
+class TestPairwiseNullNb(unittest.TestCase, _TestPairwiseNull):
 
-        :param n_cells: Number of cells to simulate (number of observations per test).
-        :param n_genes: Number of genes to simulate (number of tests).
-        """
+    def test_null_distribution_ztest(self):
         logging.getLogger("tensorflow").setLevel(logging.ERROR)
         logging.getLogger("batchglm").setLevel(logging.WARNING)
         logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
-        sim = Simulator(num_observations=n_cells, num_features=n_genes)
-        sim.generate_sample_description(num_batches=0, num_conditions=0)
-        sim.generate()
-
-        random_sample_description = pd.DataFrame({
-            "condition": np.random.randint(n_groups, size=sim.nobs)
-        })
-
-        test = de.test.pairwise(
-            data=sim.x,
-            grouping="condition",
-            test="wilcoxon",
-            sample_description=random_sample_description,
-        )
-        summary = test.summary()
-
-        # Compare p-value distribution under null model against uniform distribution.
-        pval_h0 = stats.kstest(test.pval[~np.eye(test.pval.shape[0]).astype(bool)].flatten(), 'uniform').pvalue
+        np.random.seed(1)
+        self.noise_model = "nb"
+        self._test_null_distribution_basic(test="z-test", lazy=False, quick_scale=False)
+        self._test_null_distribution_basic(test="z-test", lazy=False, quick_scale=True)
 
-        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
-        assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)
-
-        return True
-
-
-class TestPairwiseDE(unittest.TestCase):
-
-    def test_ztest_de(self, n_cells: int = 2000, n_genes: int = 500):
-        """
-        Test if de.lrt() generates a uniform p-value distribution
-        if it is given data simulated based on the null model. Returns the p-value
-        of the two-side Kolmgorov-Smirnov test for equality of the observed
-        p-value distriubution and a uniform distribution.
-
-        :param n_cells: Number of cells to simulate (number of observations per test).
-        :param n_genes: Number of genes to simulate (number of tests).
-        """
+    def test_null_distribution_ztest_lazy(self):
         logging.getLogger("tensorflow").setLevel(logging.ERROR)
         logging.getLogger("batchglm").setLevel(logging.WARNING)
         logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
-        num_non_de = n_genes // 2
-        sim = Simulator(num_observations=n_cells, num_features=n_genes)
-        sim.generate_sample_description(num_batches=0, num_conditions=2)
-        # simulate: coefficients ~ log(N(1, 0.5)).
-        # re-sample if N(1, 0.5) <= 0
-        sim.generate_params(rand_fn=lambda shape: 1 + stats.truncnorm.rvs(-1 / 0.5, np.infty, scale=0.5, size=shape))
-        sim.params["a"][1, :num_non_de] = 0
-        sim.params["b"][1, :num_non_de] = 0
-        sim.params["isDE"] = ("features",), np.arange(n_genes) >= num_non_de
-        sim.generate_data()
+        np.random.seed(1)
+        self.noise_model = "nb"
+        self._test_null_distribution_basic(test="z-test", lazy=True, quick_scale=False)
+        self._test_null_distribution_basic(test="z-test", lazy=True, quick_scale=True)
 
-        sample_description = sim.sample_description
+    def test_null_distribution_wald(self):
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
-        test = de.test.pairwise(
-            data=sim.x,
-            grouping="condition",
-            test="z-test",
-            noise_model="nb",
-            sample_description=sample_description,
-        )
-        summary = test.summary()
+        np.random.seed(1)
+        self.noise_model = "nb"
+        self._test_null_distribution_basic(test="wald", lazy=False, quick_scale=False)
+        self._test_null_distribution_basic(test="wald", lazy=False, quick_scale=True)
 
-        frac_nonde_sig = np.mean(
-            np.sum(test.qval[~np.eye(test.pval.shape[0]).astype(bool), :num_non_de] < 0.05) /
-            (2 * num_non_de)
-        )
-        frac_de_sig = np.mean(
-            np.sum(test.qval[~np.eye(test.pval.shape[0]).astype(bool), num_non_de:] < 0.05) /
-            (2 * (n_genes - num_non_de))
-        )
-        logging.getLogger("diffxpy").info('fraction of non-DE genes with q-value < 0.05: %.1f%%' %
-                                          str(np.round(100. * frac_nonde_sig, 3)))
-        logging.getLogger("diffxpy").info('fraction of DE genes with q-value < 0.05: %.1f%%' %
-                                          str(np.round(100. * frac_de_sig, 3)))
+    def test_null_distribution_lrt(self):
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
-        assert frac_de_sig > 0.5, "too many DE"
-        assert frac_nonde_sig < 0.5, "too many non-DE"
-        return True
+        np.random.seed(1)
+        self.noise_model = "nb"
+        self._test_null_distribution_basic(test="lrt", lazy=False, quick_scale=False)
 
 
 if __name__ == '__main__':