added unconstrained fitting of continuous models with interactions

Author: davidsebfischer

diffxpy/testing/tests.py

@@ -595,6 +595,7 @@ def wald(
                 np.arange(design_loc.shape[-1])[design_loc.design_info.slice(x)]
                 for x in factor_loc_totest
             ])
+        print(col_indices)
         assert len(col_indices) > 0, "Could not find any matching columns!"
         if coef_to_test is not None:
             if len(factor_loc_totest) > 1:
@@ -1989,7 +1990,7 @@ def continuous_1d(
     for x in spline_basis.columns:
         sample_description[x] = spline_basis[x].values
 
-    # Add spline basis to continuous covariate list
+    # Add spline basis to continuous covariates list
     as_numeric.extend(new_coefs)
 
     if test.lower() == 'wald':
@@ -1999,23 +2000,34 @@ def continuous_1d(
         # Adjust factors / coefficients to test:
         # Note that the continuous covariate does not necessarily have to be tested,
         # it could also be a condition effect or similar.
-        # TODO handle interactions
         if continuous in factor_loc_totest:
             # Create reduced set of factors to test which does not contain continuous:
-            factor_loc_totest_new = [x for x in factor_loc_totest if x != continuous]
+            factor_loc_totest_intermediate = [x for x in factor_loc_totest if x != continuous]
             # Add spline basis terms in instead of continuous term:
-            factor_loc_totest_new.extend(new_coefs)
+            factor_loc_totest_intermediate.extend(new_coefs)
         else:
-            factor_loc_totest_new = factor_loc_totest
+            factor_loc_totest_intermediate = factor_loc_totest
+        # Replace continuous factor in interaction terms with new spline factors.
+        factor_loc_totest_final = []
+        for i, x in enumerate(factor_loc_totest_intermediate):
+            if len(x.split(":")) > 1:
+                if np.any([x == continuous for x in x.split(":")]):
+                    interaction_partner = [y for y in x.split(":") if y != continuous][0]
+                    for y in new_coefs:
+                        factor_loc_totest_final.append(y+":"+interaction_partner)
+                else:
+                    factor_loc_totest_final.append(x)
+            else:
+                factor_loc_totest_final.append(x)
 
         logging.getLogger("diffxpy").debug("model formulas assembled in de.test.continuos():")
-        logging.getLogger("diffxpy").debug("factor_loc_totest_new: " + ",".join(factor_loc_totest_new))
+        logging.getLogger("diffxpy").debug("factor_loc_totest_final: " + ",".join(factor_loc_totest_final))
         logging.getLogger("diffxpy").debug("formula_loc_new: " + formula_loc_new)
         logging.getLogger("diffxpy").debug("formula_scale_new: " + formula_scale_new)
 
         de_test = wald(
             data=data,
-            factor_loc_totest=factor_loc_totest_new,
+            factor_loc_totest=factor_loc_totest_final,
             coef_to_test=None,
             formula_loc=formula_loc_new,
             formula_scale=formula_scale_new,
@@ -2034,6 +2046,7 @@ def continuous_1d(
             dtype=dtype,
             **kwargs
         )
+        print(de_test.model_estim.input_data.loc_names)
         de_test = DifferentialExpressionTestWaldCont(
             de_test=de_test,
             noise_model=noise_model,

diffxpy/unit_test/test_continuous_null.py

@@ -39,6 +39,31 @@ def _fit_continuous(
         )
         return test
 
+    def _fit_continuous_interaction(
+        self,
+        sim,
+        sample_description,
+        constrained,
+        test,
+        spline_basis
+    ):
+        test = de.test.continuous_1d(
+            data=sim.input_data,
+            sample_description=sample_description,
+            gene_names=["gene" + str(i) for i in range(sim.input_data.num_features)],
+            formula_loc="~ 1 + continuous + batch + continuous:batch",
+            formula_scale="~ 1",
+            factor_loc_totest=["continuous", "continuous:batch"],
+            continuous="continuous",
+            size_factors="size_factors",
+            df=3,
+            spline_basis=spline_basis,
+            test=test,
+            quick_scale=True,
+            noise_model=self.noise_model
+        )
+        return test
+
     def _test_basic(
             self,
             ngenes: int,
@@ -67,6 +92,34 @@ def _test_basic(
         )
         return det
 
+    def _test_interaction(
+            self,
+            ngenes: int,
+            test: str,
+            constrained: bool,
+            spline_basis: str
+    ):
+        n_timepoints = 5
+        sim = Simulator(num_observations=n_timepoints*200, num_features=ngenes)
+        sim.generate_sample_description(num_batches=0, num_conditions=0)
+        sim.generate_params()
+        sim.generate_data()
+
+        random_sample_description = pd.DataFrame({
+            "continuous": np.asarray(np.random.randint(0, n_timepoints, size=sim.nobs), dtype=float)
+        })
+        random_sample_description["batch"] = [str(np.random.randint(0, 3))
+                                              for x in random_sample_description["continuous"]]
+        random_sample_description["size_factors"] = np.random.uniform(0.9, 1.1, sim.nobs)  # TODO put into simulation.
+        det = self._fit_continuous_interaction(
+            sim=sim,
+            sample_description=random_sample_description,
+            test=test,
+            constrained=constrained,
+            spline_basis=spline_basis,
+        )
+        return det
+
     def _test_null_model(
             self,
             ngenes: int,
@@ -82,6 +135,21 @@ def _test_null_model(
         )
         return self._eval(det=det)
 
+    def _test_null_model_interaction(
+            self,
+            ngenes: int,
+            test: str,
+            constrained: bool,
+            spline_basis: str
+    ):
+        det = self._test_interaction(
+            ngenes=ngenes,
+            test=test,
+            constrained=constrained,
+            spline_basis=spline_basis
+        )
+        return self._eval(det=det)
+
     def _eval(self, det):
         pval_h0 = stats.kstest(det.pval, 'uniform').pvalue
         logging.getLogger("diffxpy").info(
@@ -145,6 +213,15 @@ def _test_null_model_all_splines(
         for x in ["bs", "cr", "cc"]:
             self._test_null_model(ngenes=ngenes, test=test, constrained=constrained, spline_basis=x)
 
+    def _test_null_model_all_splines_interaction(
+            self,
+            ngenes: int,
+            test: str,
+            constrained: bool
+    ):
+        for x in ["bs", "cr", "cc"]:
+            self._test_null_model_interaction(ngenes=ngenes, test=test, constrained=constrained, spline_basis=x)
+
 
 class TestContinuousNb(_TestContinuous, unittest.TestCase):
 
@@ -203,7 +280,8 @@ def test_null_distribution_wald_unconstrained(self):
 
         self.noise_model = "nb"
         np.random.seed(1)
-        self._test_null_model_all_splines(ngenes=100, test="wald", constrained=False)
+        #self._test_null_model_all_splines(ngenes=100, test="wald", constrained=False)
+        self._test_null_model_all_splines_interaction(ngenes=100, test="wald", constrained=False)
         return True
 
     def test_null_distribution_wald_constrained(self):
@@ -223,6 +301,7 @@ def test_null_distribution_wald_constrained(self):
         self.noise_model = "nb"
         np.random.seed(1)
         self._test_null_model_all_splines(ngenes=100, test="wald", constrained=True)
+        # Interaction not supported yet.
         return True
 
     def _test_null_distribution_lrt(self):