Merge pull request #110 from theislab/dev

Dev v0.6.10

Author: davidsebfischer

diffxpy/fit/fit.py

@@ -17,7 +17,7 @@
 
 
 def model(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         formula_loc: Union[None, str] = None,
         formula_scale: Union[None, str] = "~1",
         as_numeric: Union[List[str], Tuple[str], str] = (),
@@ -226,7 +226,7 @@ def model(
 
 
 def residuals(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         formula_loc: Union[None, str] = None,
         formula_scale: Union[None, str] = "~1",
         as_numeric: Union[List[str], Tuple[str], str] = (),
@@ -400,7 +400,7 @@ def residuals(
 
 
 def partition(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         parts: Union[str, np.ndarray, list],
         gene_names: Union[np.ndarray, list] = None,
         sample_description: pd.DataFrame = None,
@@ -454,7 +454,7 @@ class _Partition:
 
     def __init__(
             self,
-            data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+            data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
             parts: Union[str, np.ndarray, list],
             gene_names: Union[np.ndarray, list] = None,
             sample_description: pd.DataFrame = None,
@@ -481,7 +481,7 @@ def __init__(
             same order as in data or string-type column identifier of size-factor containing
             column in sample description.
         """
-        if isinstance(data, glm.typing.InputDataBaseTyping):
+        if isinstance(data, glm.typing.InputDataBase):
             self.x = data.x
         elif isinstance(data, anndata.AnnData) or isinstance(data, Raw):
             self.x = data.X

diffxpy/stats/stats.py

@@ -211,6 +211,7 @@ def wald_test(
         if theta_mle.shape[0] != theta0.shape[0]:
             raise ValueError('stats.wald_test(): theta_mle and theta0 have to contain the same number of entries')
 
+    theta_sd = np.nextafter(0, np.inf, out=theta_sd, where=theta_sd < np.nextafter(0, np.inf))
     wald_statistic = np.abs(np.divide(theta_mle - theta0, theta_sd))
     pvals = 2 * (1 - scipy.stats.norm(loc=0, scale=1).cdf(wald_statistic))  # two-tailed test
     return pvals
@@ -313,7 +314,10 @@ def two_coef_z_test(
     if theta_mle0.shape[0] != theta_sd0.shape[0]:
         raise ValueError('stats.two_coef_z_test(): theta_mle0 and theta_sd0 have to contain the same number of entries')
 
-    z_statistic = np.abs((theta_mle0 - theta_mle1) / np.sqrt(np.square(theta_sd0) + np.square(theta_sd1)))
+    divisor = np.square(theta_sd0) + np.square(theta_sd1)
+    divisor = np.nextafter(0, np.inf, out=divisor, where=divisor < np.nextafter(0, np.inf))
+    divisor = np.sqrt(divisor)
+    z_statistic = np.abs((theta_mle0 - theta_mle1)) / divisor
     pvals = 2 * (1 - scipy.stats.norm(loc=0, scale=1).cdf(z_statistic))  # two-tailed test
     return pvals
 

diffxpy/testing/det.py

@@ -1544,7 +1544,7 @@ def __init__(
         super().__init__()
         if isinstance(data, anndata.AnnData) or isinstance(data, anndata.Raw):
             data = data.X
-        elif isinstance(data, glm.typing.InputDataBaseTyping):
+        elif isinstance(data, glm.typing.InputDataBase):
             data = data.x
         self._x = data
         self.sample_description = sample_description
@@ -1669,7 +1669,7 @@ def __init__(
         super().__init__()
         if isinstance(data, anndata.AnnData) or isinstance(data, anndata.Raw):
             data = data.X
-        elif isinstance(data, glm.typing.InputDataBaseTyping):
+        elif isinstance(data, glm.typing.InputDataBase):
             data = data.x
         self._x = data
         self.sample_description = sample_description
@@ -2102,11 +2102,13 @@ def __init__(
         self.grouping = grouping
         self.groups = list(np.asarray(groups))
 
-        # values of parameter estimates: coefficients x genes array with one coefficient per group
-        self._theta_mle = model_estim.par_link_loc
-        # standard deviation of estimates: coefficients x genes array with one coefficient per group
-        # theta_sd = sqrt(diagonal(fisher_inv))
-        self._theta_sd = np.sqrt(np.diagonal(model_estim.fisher_inv, axis1=-2, axis2=-1)).T
+        # Values of parameter estimates: coefficients x genes array with one coefficient per group
+        self._theta_mle = model_estim.a_var
+        # Standard deviation of estimates: coefficients x genes array with one coefficient per group
+        # Need .copy() here as nextafter needs mutabls copy.
+        theta_sd = np.diagonal(model_estim.fisher_inv, axis1=-2, axis2=-1).T.copy()
+        theta_sd = np.nextafter(0, np.inf, out=theta_sd, where=theta_sd < np.nextafter(0, np.inf))
+        self._theta_sd = np.sqrt(theta_sd)
         self._logfc = None
 
         # Call tests in constructor.
@@ -2307,11 +2309,13 @@ def __init__(
         else:
             self.groups = groups.tolist()
 
-        # values of parameter estimates: coefficients x genes array with one coefficient per group
+        # Values of parameter estimates: coefficients x genes array with one coefficient per group
         self._theta_mle = model_estim.a_var
-        # standard deviation of estimates: coefficients x genes array with one coefficient per group
-        # theta_sd = sqrt(diagonal(fisher_inv))
-        self._theta_sd = np.sqrt(np.diagonal(model_estim.fisher_inv, axis1=-2, axis2=-1)).T
+        # Standard deviation of estimates: coefficients x genes array with one coefficient per group
+        # Need .copy() here as nextafter needs mutabls copy.
+        theta_sd = np.diagonal(model_estim.fisher_inv, axis1=-2, axis2=-1).T.copy()
+        theta_sd = np.nextafter(0, np.inf, out=theta_sd, where=theta_sd < np.nextafter(0, np.inf))
+        self._theta_sd = np.sqrt(theta_sd)
 
     def _correction(self, pvals, method="fdr_bh") -> np.ndarray:
         """
@@ -2349,7 +2353,6 @@ def _test(self, **kwargs):
 
     def _test_pairs(self, groups0, groups1):
         num_features = self.model_estim.x.shape[1]
-
         pvals = np.tile(np.NaN, [len(groups0), len(groups1), num_features])
 
         for i, g0 in enumerate(groups0):

diffxpy/testing/tests.py

@@ -39,7 +39,7 @@ def _fit(
         quick_scale: bool = None,
         close_session=True,
         dtype="float64"
-) -> glm.typing.InputDataBaseTyping:
+) -> glm.typing.InputDataBase:
     """
     :param noise_model: str, noise model to use in model-based unit_test. Possible options:
 
@@ -186,7 +186,7 @@ def _fit(
 
 
 def lrt(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         full_formula_loc: str,
         reduced_formula_loc: str,
         full_formula_scale: str = "~1",
@@ -370,7 +370,7 @@ def lrt(
 
 
 def wald(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         factor_loc_totest: Union[str, List[str]] = None,
         coef_to_test: Union[str, List[str]] = None,
         formula_loc: Union[None, str] = None,
@@ -547,7 +547,7 @@ def wald(
     if isinstance(as_numeric, str):
         as_numeric = [as_numeric]
 
-    # # Parse input data formats:
+    # Parse input data formats:
     gene_names = parse_gene_names(data, gene_names)
     if dmat_loc is None and dmat_scale is None:
         sample_description = parse_sample_description(data, sample_description)
@@ -644,7 +644,7 @@ def wald(
 
 
 def t_test(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         grouping,
         gene_names: Union[np.ndarray, list] = None,
         sample_description: pd.DataFrame = None,
@@ -686,7 +686,7 @@ def t_test(
 
 
 def rank_test(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         grouping: Union[str, np.ndarray, list],
         gene_names: Union[np.ndarray, list] = None,
         sample_description: pd.DataFrame = None,
@@ -728,7 +728,7 @@ def rank_test(
 
 
 def two_sample(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         grouping: Union[str, np.ndarray, list],
         as_numeric: Union[List[str], Tuple[str], str] = (),
         test: str = "t-test",
@@ -819,8 +819,8 @@ def two_sample(
     :param kwargs: [Debugging] Additional arguments will be passed to the _fit method.
     """
     if test in ['t-test', 'rank'] and noise_model is not None:
-        raise ValueError('base.two_sample(): Do not specify `noise_model` if using test t-test or rank_test: ' +
-                         'The t-test is based on a gaussian noise model and wilcoxon is model free.')
+        raise Warning('two_sample(): Do not specify `noise_model` if using test t-test or rank_test: ' +
+                      'The t-test is based on a gaussian noise model and the rank sum test is model free.')
 
     gene_names = parse_gene_names(data, gene_names)
     grouping = parse_grouping(data, sample_description, grouping)
@@ -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,
@@ -883,16 +887,14 @@ def two_sample(
             data=data,
             gene_names=gene_names,
             grouping=grouping,
-            is_sig_zerovar=is_sig_zerovar,
-            dtype=dtype
+            is_sig_zerovar=is_sig_zerovar
         )
     elif test.lower() == 'rank':
         de_test = rank_test(
             data=data,
             gene_names=gene_names,
             grouping=grouping,
-            is_sig_zerovar=is_sig_zerovar,
-            dtype=dtype
+            is_sig_zerovar=is_sig_zerovar
         )
     else:
         raise ValueError('two_sample(): Parameter `test="%s"` not recognized.' % test)
@@ -901,19 +903,19 @@ def two_sample(
 
 
 def pairwise(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         grouping: Union[str, np.ndarray, list],
         as_numeric: Union[List[str], Tuple[str], str] = (),
-        test: str = 'z-test',
-        lazy: bool = False,
+        test: str = "z-test",
+        lazy: bool = True,
         gene_names: Union[np.ndarray, list] = None,
         sample_description: pd.DataFrame = None,
-        noise_model: str = None,
+        noise_model: str = "nb",
         size_factors: np.ndarray = None,
         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,
@@ -1036,6 +1038,8 @@ def pairwise(
             design_scale=dmat,
             gene_names=gene_names,
             size_factors=size_factors,
+            init_a="closed_form",
+            init_b="closed_form",
             batch_size=batch_size,
             training_strategy=training_strategy,
             quick_scale=quick_scale,
@@ -1058,6 +1062,10 @@ def pairwise(
                 correction_type=pval_correction
             )
     else:
+        if isinstance(data, anndata.AnnData) or isinstance(data, anndata.Raw):
+            data = data.X
+        elif isinstance(data, glm.typing.InputDataBase):
+            data = data.x
         groups = np.unique(grouping)
         pvals = np.tile(np.NaN, [len(groups), len(groups), data.shape[1]])
         pvals[np.eye(pvals.shape[0]).astype(bool)] = 0
@@ -1073,16 +1081,19 @@ def pairwise(
             for j, g2 in enumerate(groups[(i + 1):]):
                 j = j + i + 1
 
-                sel = (grouping == g1) | (grouping == g2)
+                idx = np.where(np.logical_or(
+                    grouping == g1,
+                    grouping == g2
+                ))[0]
                 de_test_temp = two_sample(
-                    data=data[sel],
-                    grouping=grouping[sel],
+                    data=data[idx, :],
+                    grouping=grouping[idx],
                     as_numeric=as_numeric,
                     test=test,
                     gene_names=gene_names,
-                    sample_description=sample_description.iloc[sel],
+                    sample_description=sample_description.iloc[idx, :],
                     noise_model=noise_model,
-                    size_factors=size_factors[sel] if size_factors is not None else None,
+                    size_factors=size_factors[idx] if size_factors is not None else None,
                     batch_size=batch_size,
                     training_strategy=training_strategy,
                     quick_scale=quick_scale,
@@ -1093,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
@@ -1112,7 +1123,7 @@ def pairwise(
 
 
 def versus_rest(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         grouping: Union[str, np.ndarray, list],
         as_numeric: Union[List[str], Tuple[str], str] = (),
         test: str = 'wald',
@@ -1274,7 +1285,7 @@ def versus_rest(
 
 
 def partition(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         parts: Union[str, np.ndarray, list],
         gene_names: Union[np.ndarray, list] = None,
         sample_description: pd.DataFrame = None
@@ -1317,7 +1328,7 @@ class _Partition:
 
     def __init__(
             self,
-            data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+            data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
             parts: Union[str, np.ndarray, list],
             gene_names: Union[np.ndarray, list] = None,
             sample_description: pd.DataFrame = None
@@ -1332,7 +1343,7 @@ def __init__(
         :param gene_names: optional list/array of gene names which will be used if `data` does not implicitly store these
         :param sample_description: optional pandas.DataFrame containing sample annotations
         """
-        if isinstance(data, glm.typing.InputDataBaseTyping):
+        if isinstance(data, glm.typing.InputDataBase):
             self.x = data.x
         elif isinstance(data, anndata.AnnData) or isinstance(data, Raw):
             self.x = data.X

diffxpy/testing/utils.py

@@ -18,13 +18,13 @@
 
 
 def parse_gene_names(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         gene_names: Union[list, np.ndarray, None]
 ):
     if gene_names is None:
         if anndata is not None and (isinstance(data, anndata.AnnData) or isinstance(data, Raw)):
             gene_names = data.var_names
-        elif isinstance(data, glm.typing.InputDataBaseTyping):
+        elif isinstance(data, glm.typing.InputDataBase):
             gene_names = data.features
         else:
             raise ValueError("Missing gene names")
@@ -33,7 +33,7 @@ def parse_gene_names(
 
 
 def parse_sample_description(
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         sample_description: Union[pd.DataFrame, None]
 ) -> pd.DataFrame:
     """
@@ -57,7 +57,7 @@ def parse_sample_description(
         assert data.X.shape[0] == sample_description.shape[0], \
             "data matrix and sample description must contain same number of cells: %i, %i" % \
             (data.X.shape[0], sample_description.shape[0])
-    elif isinstance(data, glm.typing.InputDataBaseTyping):
+    elif isinstance(data, glm.typing.InputDataBase):
         assert data.x.shape[0] == sample_description.shape[0], \
             "data matrix and sample description must contain same number of cells: %i, %i" % \
             (data.x.shape[0], sample_description.shape[0])
@@ -70,7 +70,7 @@ def parse_sample_description(
 
 def parse_size_factors(
         size_factors: Union[np.ndarray, pd.core.series.Series, np.ndarray],
-        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBaseTyping],
+        data: Union[anndata.AnnData, Raw, np.ndarray, scipy.sparse.csr_matrix, glm.typing.InputDataBase],
         sample_description: pd.DataFrame
 ) -> Union[np.ndarray, None]:
     """