Skip to content

Commit

Permalink
For poisson.
Browse files Browse the repository at this point in the history
  • Loading branch information
ilan-gold committed Aug 8, 2022
1 parent 002aa28 commit f3c69a6
Show file tree
Hide file tree
Showing 9 changed files with 220 additions and 76 deletions.
12 changes: 8 additions & 4 deletions diffxpy/testing/det.py
Original file line number Diff line number Diff line change
Expand Up @@ -468,17 +468,17 @@ class DifferentialExpressionTestLRT(_DifferentialExpressionTestSingle):

sample_description: pd.DataFrame
full_design_loc_info: patsy.design_info
full_estim: glm.train.numpy.glm_base.Estimator
full_estim: glm.train.base.BaseEstimatorGlm
reduced_design_loc_info: patsy.design_info
reduced_estim: glm.train.numpy.glm_base.Estimator
reduced_estim: glm.train.base.BaseEstimatorGlm

def __init__(
self,
sample_description: pd.DataFrame,
full_design_loc_info: patsy.design_info,
full_estim: glm.train.numpy.glm_base.Estimator,
full_estim: glm.train.base.BaseEstimatorGlm,
reduced_design_loc_info: patsy.design_info,
reduced_estim: glm.train.numpy.glm_base.Estimator
reduced_estim: glm.train.base.BaseEstimatorGlm
):
super().__init__()
self.sample_description = sample_description
Expand Down Expand Up @@ -1287,6 +1287,10 @@ def _assemble_gene_fits(
loc=loc,
scale=scale
)
elif self.noise_model == "poisson":
yhat = np.random.poisson(
lam=loc
)
else:
raise ValueError("noise model %s not yet supported for plot_gene_fits" % self.noise_model)

Expand Down
3 changes: 3 additions & 0 deletions diffxpy/testing/tests.py
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,9 @@ def _fit(
elif noise_model == "norm" or noise_model == "normal":
from batchglm.train.numpy.glm_norm import Estimator
from batchglm.models.glm_norm import Model
elif noise_model == "poisson":
from batchglm.train.numpy.glm_poisson import Estimator
from batchglm.models.glm_poisson import Model
else:
raise ValueError('noise_model="%s" not recognized.' % noise_model)
# Set default chunk size:
Expand Down
2 changes: 1 addition & 1 deletion diffxpy/unit_test/test_acc_glm_all_numpy_temp.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ def _test_full_model(self, noise_model):
_ = temp.summary()

def test(self):
for noise_model in ['norm', 'poisson', 'nb']:
for noise_model in ['poisson', 'norm', 'nb']:
self._test_full_model(noise_model)


Expand Down
21 changes: 13 additions & 8 deletions diffxpy/unit_test/test_backends.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import scipy.stats as stats
from batchglm.models.glm_nb import Model as NBModel
from batchglm.models.glm_norm import Model as NormModel
from batchglm.models.glm_poisson import Model as PoissonModel

import diffxpy.api as de

Expand Down Expand Up @@ -34,6 +35,9 @@ def _test_null_distribution_wald(
elif noise_model == "norm":
model = NormModel()
rand_fn_scale = lambda shape: np.random.uniform(1, 2, shape)
elif noise_model == "poisson":
model = PoissonModel()
rand_fn_scale = lambda shape: np.random.uniform(1, 2, shape)


model.generate_artificial_data(
Expand Down Expand Up @@ -69,13 +73,13 @@ def _test_null_distribution_wald(
return True


class TestSingleNullBackendsNb(_TestSingleNullBackends, unittest.TestCase):
class TestSingleNullBackends(_TestSingleNullBackends, unittest.TestCase):
"""
Negative binomial noise model unit tests that test whether a test generates uniformly
distributed p-values if data are sampled from the null model.
"""

def test_null_distribution_wald_nb_numpy(
def test_null_distribution_wald_numpy(
self,
n_cells: int = 2000,
n_genes: int = 200
Expand All @@ -91,12 +95,13 @@ def test_null_distribution_wald_nb_numpy(
logging.getLogger("diffxpy").setLevel(logging.WARNING)

np.random.seed(1)
_ = self._test_null_distribution_wald(
n_cells=n_cells,
n_genes=n_genes,
noise_model="nb",
backend="numpy"
)
for noise_model in ['poisson', 'nb', 'norm']:
_ = self._test_null_distribution_wald(
n_cells=n_cells,
n_genes=n_genes,
noise_model="nb",
backend="numpy"
)


if __name__ == '__main__':
Expand Down
25 changes: 25 additions & 0 deletions diffxpy/unit_test/test_continuous_de.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@

from batchglm.models.glm_nb import Model as NBModel
from batchglm.models.glm_norm import Model as NormModel
from batchglm.models.glm_poisson import Model as PoissonModel

class _TestContinuousDe:
noise_model: str
Expand All @@ -24,6 +25,10 @@ def _test_wald_de(
model = NormModel()
rand_fn_loc = lambda shape: np.random.uniform(500, 1000, shape)
rand_fn_scale = lambda shape: np.random.uniform(1, 2, shape)
elif self.noise_model == "poisson":
model = PoissonModel()
rand_fn_loc = lambda shape: np.random.uniform(2, 10, shape)
rand_fn_scale = lambda shape: np.random.uniform(1, 2, shape) # not used but called for uniformity
else:
raise ValueError("noise model %s not recognized" % self.noise_model)

Expand Down Expand Up @@ -142,6 +147,26 @@ def test_wald_de_norm(self):
self._test_wald_de_all_splines(ngenes=100, constrained=True)
return True

class TestContinuousDePoisson(_TestContinuousDe, unittest.TestCase):
"""
Normal noise model unit tests that tests false positive and false negative rates.
"""

def test_wald_de_poisson(self):
"""
:return:
"""
logging.getLogger("tensorflow").setLevel(logging.ERROR)
logging.getLogger("batchglm").setLevel(logging.WARNING)
logging.getLogger("diffxpy").setLevel(logging.WARNING)

self.noise_model = "poisson"
np.random.seed(1)
self._test_wald_de_all_splines(ngenes=100, constrained=False)
self._test_wald_de_all_splines(ngenes=100, constrained=True)
return True


if __name__ == '__main__':
unittest.main()
77 changes: 77 additions & 0 deletions diffxpy/unit_test/test_fit.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
import diffxpy.api as de
from batchglm.models.glm_nb import Model as NBModel
from batchglm.models.glm_norm import Model as NormModel
from batchglm.models.glm_poisson import Model as PoissonModel

class _TestFit:

Expand All @@ -31,6 +32,9 @@ def _test_model_fit(
elif noise_model == "norm":
model = NormModel()
rand_fn_scale = lambda shape: np.random.uniform(1, 2, shape)
elif noise_model == "poisson":
model = PoissonModel()
rand_fn_scale = lambda shape: np.random.uniform(1, 2, shape) # since it is called later but not used
else:
raise ValueError("noise model %s not recognized" % noise_model)

Expand Down Expand Up @@ -126,6 +130,8 @@ def _test_residuals_fit(
model = NBModel()
elif noise_model == "norm":
model = NormModel()
elif noise_model == "poisson":
model = PoissonModel()
else:
raise ValueError("noise model %s not recognized" % noise_model)

Expand Down Expand Up @@ -294,6 +300,77 @@ def test_residuals_fit(
noise_model="norm"
)

class TestFitNorm(_TestFit, unittest.TestCase):
"""
Normal noise model unit tests that tests whether model fit relay works.
"""

def test_model_fit(
self,
n_cells: int = 2000,
n_genes: int = 2
):
"""
Test if model fit for "norm" noise model works.
: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)

np.random.seed(1)
return self._test_model_fit(
n_cells=n_cells,
n_genes=n_genes,
noise_model="poisson"
)

def test_model_fit_partition(
self,
n_cells: int = 2000,
n_genes: int = 2
):
"""
Test if partitioned model fit for "norm" noise model works.
: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)

np.random.seed(1)
return self._test_model_fit_partition(
n_cells=n_cells,
n_genes=n_genes,
noise_model="poisson"
)

def test_residuals_fit(
self,
n_cells: int = 2000,
n_genes: int = 2
):
"""
Test if residual fit for "norm" noise model works.
: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)

np.random.seed(1)
return self._test_residuals_fit(
n_cells=n_cells,
n_genes=n_genes,
noise_model="poisson"
)


if __name__ == '__main__':
unittest.main()
48 changes: 46 additions & 2 deletions diffxpy/unit_test/test_pairwise_null.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import scipy.stats as stats
from batchglm.models.glm_nb import Model as NBModel
from batchglm.models.glm_norm import Model as NormModel
from batchglm.models.glm_poisson import Model as PoissonModel

import diffxpy.api as de

Expand All @@ -29,6 +30,10 @@ def _prepate_data(
rand_fn_loc = lambda shape: np.random.uniform(500, 1000, shape)
rand_fn_scale = lambda shape: np.random.uniform(1, 2, shape)
model = NormModel()
elif self.noise_model == "poisson":
rand_fn_loc = lambda shape: np.random.uniform(2, 10, shape)
rand_fn_scale = lambda shape: np.random.uniform(1, 2, shape) # not used b/c only loc param
model = PoissonModel()
else:
raise ValueError("noise model %s not recognized" % self.noise_model)

Expand Down Expand Up @@ -109,7 +114,7 @@ def test_null_distribution_ttest(self):
logging.getLogger("diffxpy").setLevel(logging.WARNING)

np.random.seed(1)
self.noise_model = None
self.noise_model = "norm"
self._test_null_distribution_basic(test="t-test", lazy=False)

def test_null_distribution_rank(self):
Expand All @@ -118,9 +123,48 @@ def test_null_distribution_rank(self):
logging.getLogger("diffxpy").setLevel(logging.WARNING)

np.random.seed(1)
self.noise_model = None
self.noise_model = "norm"
self._test_null_distribution_basic(test="rank", lazy=False)

class TestPairwiseNullPoisson(unittest.TestCase, _TestPairwiseNull):

def test_null_distribution_ztest(self):
logging.getLogger("tensorflow").setLevel(logging.ERROR)
logging.getLogger("batchglm").setLevel(logging.WARNING)
logging.getLogger("diffxpy").setLevel(logging.WARNING)

np.random.seed(1)
self.noise_model = "poisson"
self._test_null_distribution_basic(test="z-test", lazy=False, quick_scale=False)

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)

np.random.seed(1)
self.noise_model = "poisson"
self._test_null_distribution_basic(test="z-test", lazy=True, quick_scale=False)

def test_null_distribution_wald(self):
logging.getLogger("tensorflow").setLevel(logging.ERROR)
logging.getLogger("batchglm").setLevel(logging.WARNING)
logging.getLogger("diffxpy").setLevel(logging.WARNING)

np.random.seed(1)
self.noise_model = "poisson"
self._test_null_distribution_basic(test="wald", lazy=False, quick_scale=False)

def test_null_distribution_lrt(self):
logging.getLogger("tensorflow").setLevel(logging.ERROR)
logging.getLogger("batchglm").setLevel(logging.WARNING)
logging.getLogger("diffxpy").setLevel(logging.WARNING)

np.random.seed(1)
self.noise_model = "poisson"
self._test_null_distribution_basic(test="lrt", lazy=False, quick_scale=False)



class TestPairwiseNullNb(unittest.TestCase, _TestPairwiseNull):

Expand Down
Loading

0 comments on commit f3c69a6

Please sign in to comment.