From 44b050054bcf13ceccfb90cd4b0425e6b95d2244 Mon Sep 17 00:00:00 2001 From: Benjamin Pedigo Date: Tue, 16 May 2023 18:48:17 -0400 Subject: [PATCH] Added density and group connection tests (#1032) * Created new branch with files from jw-feature branch * Fix seaborn syntax * Increment bugfix version * Cleared many mypy errors * Corrected some typing and improved test * More fixes and edits * Fixing docs issues * Exclude 3.6.1 I added to line 35. I assumed the asterisk on 3.3 meant all versions of 3.3 ex (3.3.1, 3.3.2) are all excluded from the matplotlib. Is this correct? * Update setup.cfg * types need to be AdjacencyMatrix * bump scipy, which means deprecate python 3.7 * maybe fix labels error * add beartype * clean up docs and add a note * clean up line lengths and indentation in docstring * black * fix line lengths * add intersphinx * add another intersphinx * fix doc indentation * more docstring work * clarify number of nodes * add math for null and alt * remove unused imports * typo * add some math and notes and ref * add ref * disclaimers * add some warnings * try to fix changes from auto toc * fix isort * fix typing issues * Update README.md Removed outdated Zenodo DOI from README.md * Deleted bad kwargs from tests * some typing/tests * fix to non-unity probability ratio * fix code format * fix tests * fix citation * fix more citation * remove file * fix docs * clean up tests * docs * fix up tutorials * fix notebook * fix int issue * run black * black with updated * fix notebook math * fix imports * clean up tutorial --------- Co-authored-by: Jeremy Welland Co-authored-by: Dax Pryce Co-authored-by: hugwuoke <85888975+hugwuoke@users.noreply.github.com> Co-authored-by: Kartikeya Tripathi <96724863+ktwillcode@users.noreply.github.com> --- README.md | 2 +- docs/conf.py | 1 + docs/reference/reference/inference.rst | 4 + docs/sphinx-ext/toctree_filter.py | 11 +- docs/tutorials/index.rst | 2 + docs/tutorials/inference/density_test.ipynb | 194 ++++++++ .../inference/group_connection_test.ipynb | 238 +++++++++ graspologic/inference/__init__.py | 9 +- graspologic/inference/binomial.py | 84 ++++ graspologic/inference/density_test.py | 111 +++++ .../inference/group_connection_test.py | 471 ++++++++++++++++++ graspologic/inference/utils.py | 63 +++ graspologic/nominate/spectralVN.py | 6 +- mypy.ini | 3 + runtime.txt | 1 + setup.cfg | 3 +- tests/test_er_and_group_connection_tests.py | 85 ++++ 17 files changed, 1278 insertions(+), 10 deletions(-) create mode 100644 docs/tutorials/inference/density_test.ipynb create mode 100644 docs/tutorials/inference/group_connection_test.ipynb create mode 100644 graspologic/inference/binomial.py create mode 100644 graspologic/inference/density_test.py create mode 100644 graspologic/inference/group_connection_test.py create mode 100644 graspologic/inference/utils.py create mode 100644 runtime.txt create mode 100644 tests/test_er_and_group_connection_tests.py diff --git a/README.md b/README.md index 627ebe4b3..3e7c22a1f 100644 --- a/README.md +++ b/README.md @@ -7,7 +7,7 @@ [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) ## `graspologic` is a package for graph statistical algorithms. - + - [Overview](#overview) - [Documentation](#documentation) - [System Requirements](#system-requirements) diff --git a/docs/conf.py b/docs/conf.py index e4499c01c..8abfd50db 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -86,6 +86,7 @@ "scipy": ("https://docs.scipy.org/doc/scipy", None), "seaborn": ("https://seaborn.pydata.org", None), "sklearn": ("https://scikit-learn.org/dev", None), + "statsmodels": ("https://www.statsmodels.org/stable", None), } intersphinx_disabled_reftypes = [] diff --git a/docs/reference/reference/inference.rst b/docs/reference/reference/inference.rst index f73d41763..3383e70fd 100644 --- a/docs/reference/reference/inference.rst +++ b/docs/reference/reference/inference.rst @@ -6,6 +6,10 @@ Inference Two-graph hypothesis testing ---------------------------- +.. autofunction:: density_test + +.. autofunction:: group_connection_test + .. autofunction:: latent_position_test .. autofunction:: latent_distribution_test diff --git a/docs/sphinx-ext/toctree_filter.py b/docs/sphinx-ext/toctree_filter.py index e916fe4f3..34d6b5037 100644 --- a/docs/sphinx-ext/toctree_filter.py +++ b/docs/sphinx-ext/toctree_filter.py @@ -1,13 +1,15 @@ # Copied and modified from https://stackoverflow.com/questions/15001888/conditional-toctree-in-sphinx import re + from sphinx.directives.other import TocTree def setup(app): - app.add_config_value('toc_filter_exclude', [], 'html') - app.add_directive('toctree-filt', TocTreeFilt) - return {'version': '1.0.0'} + app.add_config_value("toc_filter_exclude", [], "html") + app.add_directive("toctree-filt", TocTreeFilt) + return {"version": "1.0.0"} + class TocTreeFilt(TocTree): """ @@ -21,7 +23,8 @@ class TocTreeFilt(TocTree): form `:secret:ultra-api` or `:draft:new-features` will be excuded from the final table of contents. Entries without a prefix are always included. """ - hasPat = re.compile('\s*(.*)$') + + hasPat = re.compile("\s*(.*)$") # Remove any entries in the content that we dont want and strip # out any filter prefixes that we want but obviously don't want the diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index dc97436e8..b74dcd7d5 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -80,6 +80,8 @@ are tutorials for robust statistical hypothesis testing on multiple graphs. :maxdepth: 1 :titlesonly: + inference/density_test + inference/group_connection_test inference/latent_position_test inference/latent_distribution_test diff --git a/docs/tutorials/inference/density_test.ipynb b/docs/tutorials/inference/density_test.ipynb new file mode 100644 index 000000000..7e99a2a98 --- /dev/null +++ b/docs/tutorials/inference/density_test.ipynb @@ -0,0 +1,194 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "9b546925", + "metadata": {}, + "source": [ + "# Testing Symmetry of Two Networks with the Density Test\n", + "The \"inference\" module of graspologic contains functions that enable quantitative comparison of two networks to assess whether they are statistically similar. This \"similarity\" can be assessed in a few different ways, depending on the details of the networks to be compared and the preferences of the user. \n", + "\n", + "The simplest test that can be performed is the density test, which is based upon the Erdos-Renyi model. Under this model, it is assumed that the probability of an edge between any two nodes of the network is some constant, p. To compare two networks, then, the question is whether the edge probability for the first network is different from the edge probability for the second network. This test can be performed easily with the inference module, and the procedure is described in greater detail below. " + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "75c9ea1b", + "metadata": {}, + "source": [ + "## The Erdos-Renyi (ER) model\n", + "The [**Erdos-Renyi (ER) model**\n", + "](https://en.wikipedia.org/wiki/Erd%C5%91s%E2%80%93R%C3%A9nyi_model)\n", + "is one of the simplest network models. This model treats\n", + "the probability of each potential edge in the network occuring to be the same. In\n", + "other words, all edges between any two nodes are equally likely.\n", + "\n", + "```{admonition} Math\n", + "Let $n$ be the number of nodes. We say that for all $(i, j), i \\neq j$, with $i$ and\n", + "$j$ both running\n", + "from $1 ... n$, the probability of the edge $(i, j)$ occuring is:\n", + "\n", + "$$ P[A_{ij} = 1] = p_{ij} = p $$\n", + "\n", + "Where $p$ is the the global connection probability.\n", + "\n", + "Each element of the adjacency matrix $A$ is then sampled independently according to a\n", + "[Bernoulli distribution](https://en.wikipedia.org/wiki/Bernoulli_distribution):\n", + "\n", + "$$ A_{ij} \\sim Bernoulli(p) $$\n", + "\n", + "For a network modeled as described above, we say it is distributed\n", + "\n", + "$$ A \\sim ER(n, p) $$\n", + "\n", + "```\n", + "\n", + "Thus, for this model, the only parameter of interest is the global connection\n", + "probability, $p$. This is sometimes also referred to as the **network density**." + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "27daff33", + "metadata": {}, + "source": [ + "## Testing under the ER model\n", + "In order to compare two networks $A^{(L)}$ and $A^{(R)}$ under this model, we\n", + "simply need to compute these network densities ($p^{(L)}$ and $p^{(R)}$), and then\n", + "run a statistical test to see if these densities are significantly different.\n", + "\n", + "```{admonition} Math\n", + "Under this\n", + "model, the total number of edges $m$ comes from a $Binomial(n(n-1), p)$ distribution,\n", + "where $n$ is the number of nodes. This is because the number of edges is the sum of\n", + "independent Bernoulli trials with the same probability. If $m^{(L)}$ is the number of\n", + "edges on the left\n", + "hemisphere, and $m^{(R)}$ is the number of edges on the right, then we have:\n", + "\n", + "$$m^{(L)} \\sim Binomial(n^{(L)}(n^{(L)} - 1), p^{(L)})$$\n", + "\n", + "and independently,\n", + "\n", + "$$m^{(R)} \\sim Binomial(n^{(R)}(n^{(R)} - 1), p^{(R)})$$\n", + "\n", + "To compare the two networks, we are just interested in a comparison of $p^{(L)}$ vs.\n", + "$p^{(R)}$. Formally, we are testing:\n", + "\n", + "$$H_0: p^{(L)} = p^{(R)}, \\quad H_a: p^{(L)} \\neq p^{(R)}$$\n", + "\n", + "Fortunately, the problem of testing for equal proportions is well studied.\n", + "Using graspologic.inference, we can conduct this comparison using either\n", + "Fisher's exact test or the chi-squared test by using method=\"fisher\" or \n", + "method = \"chi2\", respectively. In this example, we use Fisher's exact test.\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e9a74b29", + "metadata": { + "execution": { + "iopub.execute_input": "2022-04-19T19:39:17.453921Z", + "iopub.status.busy": "2022-04-19T19:39:17.453693Z", + "iopub.status.idle": "2022-04-19T19:39:24.698064Z", + "shell.execute_reply": "2022-04-19T19:39:24.697216Z" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from graspologic.inference.density_test import density_test\n", + "from graspologic.simulations import er_np\n", + "from graspologic.plot import heatmap\n", + "\n", + "np.random.seed(8888)\n", + "\n", + "%matplotlib inline" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "6895aedd", + "metadata": {}, + "source": [ + "# Performing the Density Test \n", + "\n", + "To illustrate the density test, we will first randomly generate two networks of known density to compare using the test. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a682bc4a", + "metadata": {}, + "outputs": [], + "source": [ + "A1 = er_np(500, 0.6)\n", + "A2 = er_np(400, 0.8)\n", + "heatmap(A1, title='Adjacency Matrix for Network 1')\n", + "heatmap(A2, title='Adjacency Matrix for Network 2')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "18611379", + "metadata": {}, + "source": [ + "Visibly, these networks have very different densities. We can statistically confirm this difference by conducting a density test." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44d3204c", + "metadata": {}, + "outputs": [], + "source": [ + "stat, pvalue, er_misc = density_test(A1,A2)\n", + "print(pvalue)\n" + ] + } + ], + "metadata": { + "jupytext": { + "cell_metadata_filter": "-all", + "main_language": "python", + "notebook_metadata_filter": "-all" + }, + "kernelspec": { + "display_name": "Python 3.9.5 ('venv')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.5" + }, + "vscode": { + "interpreter": { + "hash": "7b0fa2133086a4b245bc2fc6826174141053e0208e756a5fae09980b942619c9" + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorials/inference/group_connection_test.ipynb b/docs/tutorials/inference/group_connection_test.ipynb new file mode 100644 index 000000000..7abd05c0b --- /dev/null +++ b/docs/tutorials/inference/group_connection_test.ipynb @@ -0,0 +1,238 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "567eb99a", + "metadata": {}, + "source": [ + "# Group connection test\n", + "In this tutorial, we demonstrate the use of the ``graspologic.inference`` module to compare subgraph densities of two networks, both of which contain nodes belonging to some known set of families or groups. The number and identity of the families or groups must be the same in the two networks. The ``group_connection_test`` function can then be used to determine whether there are any statistical differences in the group-to-group connection densities of the two networks.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "78dfca69", + "metadata": { + "execution": { + "iopub.execute_input": "2022-04-13T20:42:53.628597Z", + "iopub.status.busy": "2022-04-13T20:42:53.627851Z", + "iopub.status.idle": "2022-04-13T20:42:59.829258Z", + "shell.execute_reply": "2022-04-13T20:42:59.829874Z" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from graspologic.inference.group_connection_test import group_connection_test\n", + "from graspologic.simulations import sbm\n", + "from graspologic.plot import heatmap\n", + "\n", + "np.random.seed(8888)\n", + "\n", + "%matplotlib inline" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "75aa28e9", + "metadata": {}, + "source": [ + "## The stochastic block model (SBM)\n", + "A [**stochastic block model (SBM)**\n", + "](https://en.wikipedia.org/wiki/Stochastic_block_model)\n", + "is a popular statistical model of networks. Put simply, this model treats the\n", + "probability of an edge occuring between node $i$ and node $j$ as purely a function of\n", + "the *communities* or *groups* that node $i$ and $j$ belong to. Therefore, this model\n", + "is parameterized by:\n", + "\n", + " 1. An assignment of each node in the network to a group. Note that this assignment\n", + " can be considered to be deterministic or random, depending on the specific\n", + " framing of the model one wants to use.\n", + " 2. A set of group-to-group connection probabilities\n", + "\n", + "Let $n$ be the number of nodes, and $K$ be the number of groups in an SBM. For a\n", + "network $A$ sampled from an SBM:\n", + "\n", + "$$ A \\sim SBM(B, \\tau)$$\n", + "\n", + "We say that for all $(i,j), i \\neq j$, with $i$ and $j$ both running\n", + "from $1 ... n$ the probability of edge $(i,j)$ occuring is:\n", + "\n", + "$$ P[A_{ij} = 1] = P_{ij} = B_{\\tau_i, \\tau_j} $$\n", + "\n", + "where $B \\in [0,1]^{K \\times K}$ is a matrix of group-to-group connection\n", + "probabilities and $\\tau \\in \\{1...K\\}^n$ is a vector of node-to-group assignments.\n", + "Note that here we are assuming $\\tau$ is a fixed vector of assignments, though other\n", + "formuations of the SBM allow these assignments to themselves come from a categorical\n", + "distribution.\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "c94ac2f1", + "metadata": {}, + "source": [ + "## Testing under the SBM model\n", + "Assuming this model, there are a few ways that one could test for differences between\n", + "two networks. In this example, we are interested in comparing the group-to-group\n", + "connection probability matrices, $B$, for the two graphs.\n", + "\n", + "We are interested in testing:\n", + "\n", + "$$ H_0: B^{(L)} = B^{(R)}, \\quad H_A: B^{(L)} \\neq B^{(R)} $$\n", + "\n", + "Rather than having to compare one proportion as in the density test, this test \n", + "requires comparing all $K^2$ probabilities between the SBM models for both\n", + "networks.\n", + "\n", + "\n", + "The hypothesis test above can be decomposed into $K^2$ indpendent hypotheses. \n", + "$B^{(L)}$ and $B^{(R)}$ are both $K \\times K$ matrices, where each element $b_{kl}$ \n", + "represents the probability of a connection from a node in group $k$ to one in group $l$. We\n", + "also know that group $k$ for the left network corresponds with group $k$ for the\n", + "right. In other words, the *groups* are matched. Thus, we are interested in testing,\n", + "for $k, l$ both running from $1...K$:\n", + "\n", + "$$ H_0: B_{kl}^{(L)} = B_{kl}^{(R)}, \\quad H_A: B_{kl}^{(L)} \\neq B_{kl}^{(R)}$$\n", + "\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "7dfedeec", + "metadata": {}, + "source": [ + "To perform the comparison, the user provides the adjacency matrices and label vectors for both networks. In this example, we will generate random networks with known properties to demonstrate the use of the ``group_connection_test`` function.\n", + "\n", + "First, we generate and plot an adjacency matrix representing a network with two groups and a specified array of group-to-group connection probabilities." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ecfa8db", + "metadata": { + "execution": { + "iopub.execute_input": "2022-04-13T20:42:59.842240Z", + "iopub.status.busy": "2022-04-13T20:42:59.841624Z", + "iopub.status.idle": "2022-04-13T20:43:04.453604Z", + "shell.execute_reply": "2022-04-13T20:43:04.454013Z" + }, + "lines_to_next_cell": 2, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "P = np.array([[0.8, 0.6],\n", + " [0.6, 0.8]])\n", + "csize = [50] * 2\n", + "A1, labels1 = sbm(csize, P, return_labels = True)\n", + "heatmap(A1, title='2-block SBM adjacency matrix')\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "84334ee7", + "metadata": {}, + "source": [ + "Next, we generate a second adjacency matrix for a second network. We will give this network a different number of nodes and a different connection probability matrix." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9220a838", + "metadata": {}, + "outputs": [], + "source": [ + "P = np.array([[0.865, 0.66],\n", + " [0.66, 0.865]])\n", + "csize = [60] * 2\n", + "A2, labels2 = sbm(csize, P, return_labels = True)\n", + "heatmap(A2, title='2-block SBM adjacency matrix')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "81e09cff", + "metadata": {}, + "source": [ + "Now, we can run ``group_connection_test`` to assess whether there is a statistical difference between these two networks. Of course, we expect that there will be as, by design, their group-to-group connection densities are different." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dddfad17", + "metadata": {}, + "outputs": [], + "source": [ + "stat, pvalue, misc = group_connection_test(A1, A2, labels1, labels2)\n", + "print(pvalue)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "935ebedb", + "metadata": {}, + "source": [ + "This extremely low p-value suggests that we should reject the null hypothesis and conclude that the two networks are statistically different under the stochastic block model. The individual p-values which compare each group-to-group connection are also extremely low:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9169bd1", + "metadata": {}, + "outputs": [], + "source": [ + "print(misc[\"corrected_pvalues\"])" + ] + } + ], + "metadata": { + "jupytext": { + "cell_metadata_filter": "-all", + "main_language": "python", + "notebook_metadata_filter": "-all" + }, + "kernelspec": { + "display_name": "Python 3.9.5 ('venv')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.5" + }, + "vscode": { + "interpreter": { + "hash": "7b0fa2133086a4b245bc2fc6826174141053e0208e756a5fae09980b942619c9" + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/graspologic/inference/__init__.py b/graspologic/inference/__init__.py index 9431c56aa..6980f542e 100644 --- a/graspologic/inference/__init__.py +++ b/graspologic/inference/__init__.py @@ -1,7 +1,14 @@ # Copyright (c) Microsoft Corporation and contributors. # Licensed under the MIT License. +from .density_test import density_test +from .group_connection_test import group_connection_test from .latent_distribution_test import latent_distribution_test from .latent_position_test import latent_position_test -__all__ = ["latent_position_test", "latent_distribution_test"] +__all__ = [ + "density_test", + "group_connection_test", + "latent_position_test", + "latent_distribution_test", +] diff --git a/graspologic/inference/binomial.py b/graspologic/inference/binomial.py new file mode 100644 index 000000000..5e4f19d3b --- /dev/null +++ b/graspologic/inference/binomial.py @@ -0,0 +1,84 @@ +from collections import namedtuple +from typing import Literal + +import numpy as np +from scipy.stats import chi2_contingency, fisher_exact +from statsmodels.stats.proportion import test_proportions_2indep + +BinomialResult = namedtuple("BinomialResult", ["stat", "pvalue"]) +BinomialTestMethod = Literal["score", "fisher", "chi2"] + + +def binom_2samp( + x1: int, + n1: int, + x2: int, + n2: int, + null_ratio: float = 1.0, + method: BinomialTestMethod = "score", +) -> BinomialResult: + """ + This function computes the likelihood that two binomial samples are drown from + identical underlying distributions. Null hypothesis is that the success probability + for each sample is identical (i.e. p1 = p2), and this function returns the + probability that the null hypothesis is accurate, under a variety of potential + statistical tests (default is score test). + + Parameters + ---------- + x1 : int + Success count for group 1 + n1 : int + The number of possible successes for group 1 + x2 : int + Success count for group 2 + n2 : int + The number of possible successes for group 2 + null_ratio : float, optional + Optional parameter for testing whether p1 is a fixed ratio larger or smaller + than p2, i.e. p1 = cp2, where c is the null_ratio. Default is 1.0. This + parameter can only be !=1 if the chosen statistical test is the score test. + method : str, optional + Defines the statistical test to be run in order to reject or fail to reject the + null hypothesis. By default, this is the score test (i.e. "score"). + + Returns + ------- + BinomialResult: namedtuple + This namedtuple contains the following data: + stat: float + Test statistic for the requested test. + pvalue: float + The p-value for the requested test. + + References + ------ + [1] Alan Agresti. Categorical data analysis. John Wiley & Sons, 3 edition, 2013. + + """ + if x1 == 0 or x2 == 0: + # logging.warn("One or more counts were 0, not running test and returning nan") + return BinomialResult(np.nan, np.nan) + if null_ratio != 1 and method != "score": + raise ValueError("Non-unity null odds only works with ``method=='score'``") + + cont_table = np.array([[x1, n1 - x1], [x2, n2 - x2]]) + if method == "fisher" and null_ratio == 1.0: + stat, pvalue = fisher_exact(cont_table, alternative="two-sided") + elif method == "chi2": + stat, pvalue, _, _ = chi2_contingency(cont_table) + elif method == "score": + stat, pvalue = test_proportions_2indep( + x1, + n1, + x2, + n2, + method="score", + compare="ratio", + alternative="two-sided", + value=null_ratio, + ) + else: + raise ValueError() + + return BinomialResult(stat, pvalue) diff --git a/graspologic/inference/density_test.py b/graspologic/inference/density_test.py new file mode 100644 index 000000000..1633df067 --- /dev/null +++ b/graspologic/inference/density_test.py @@ -0,0 +1,111 @@ +from collections import namedtuple + +import numpy as np +from beartype import beartype + +from ..types import AdjacencyMatrix +from .binomial import BinomialTestMethod +from .group_connection_test import group_connection_test + +DensityTestResult = namedtuple("DensityTestResult", ["stat", "pvalue", "misc"]) + + +@beartype +def density_test( + A1: AdjacencyMatrix, A2: AdjacencyMatrix, method: BinomialTestMethod = "fisher" +) -> DensityTestResult: + r""" + Compares two networks by testing whether the global connection probabilities + (densites) for the two networks are equal under an Erdos-Renyi model assumption. + + Parameters + ---------- + A1: np.array, int + The adjacency matrix for network 1. Will be treated as a binary network, + regardless of whether it was weighted. + A2: np.array, int + Adjacency matrix for network 2. Will be treated as a binary network, + regardless of whether it was weighted. + method: string, optional, default="fisher" + Specifies the statistical test to be used. The default option is "fisher", + which uses Fisher's exact test, but the user may also enter "chi2" to use a + chi-squared test. Fisher's exact test is more appropriate when the expected + number of edges are small. + + Returns + ------- + DensityTestResult: namedtuple + This named tuple returns the following data: + + stat: float + The statistic for the test specified by ``method``. + pvalue: float + The p-value for the test specified by ``method``. + misc: dict + Dictionary containing a number of computed statistics for the network + comparison performed: + + "probability1", float + The probability of an edge (density) in network 1 (:math:`p_1`). + "probability2", float + The probability of an edge (density) in network 2 (:math:`p_2`). + "observed1", pd.DataFrame + The total number of edge connections for network 1. + "observed2", pd.DataFrame + The total number of edge connections for network 2. + "possible1", pd.DataFrame + The total number of possible edges for network 1. + "possible2", pd.DataFrame + The total number of possible edges for network 1. + + + Notes + ----- + Under the Erdos-Renyi model, edges are generated independently with probability + :math:`p`. :math:`p` is also known as the network density. This function tests + whether the probability of an edge in network 1 (:math:`p_1`) is significantly + different from that in network 2 (:math:`p_2`), by assuming that both networks came + from an Erdos-Renyi model. In other words, the null hypothesis is + + .. math:: H_0: p_1 = p_2 + + And the alternative hypothesis is + + .. math:: H_A: p_1 \neq p_2 + + This test makes several assumptions about the data and test (which could easily be + loosened in future versions): + + We assume that the networks are directed. If the networks are undirected + (and the adjacency matrices are thus symmetric), then edges would be counted + twice, which would lead to an incorrect calculation of the edge probability. + We believe passing in the upper or lower triangle of the adjacency matrix + would solve this, but this has not been tested. + + We assume that the networks are loopless, that is we do not consider the + probability of an edge existing between a node and itself. This can be + weakened and made an option in future versions. + + We only implement the alternative hypothesis of "not equals" (two-sided); + future versions could implement the one-sided alternative hypotheses. + + References + ---------- + .. [1] Pedigo, B.D., Powell, M., Bridgeford, E.W., Winding, M., Priebe, C.E., + Vogelstein, J.T. "Generative network modeling reveals quantitative + definitions of bilateral symmetry exhibited by a whole insect brain + connectome," eLife (2023): e83739. + """ + stat, pvalue, misc = group_connection_test( + A1, + A2, + labels1=np.ones(A1.shape[0]), + labels2=np.ones(A2.shape[0]), + method=method, + ) + misc["probability1"] = misc["probabilities1"] + del misc["probabilities1"] + misc["probability2"] = misc["probabilities2"] + del misc["probabilities2"] + + return DensityTestResult(stat, pvalue, misc) diff --git a/graspologic/inference/group_connection_test.py b/graspologic/inference/group_connection_test.py new file mode 100644 index 000000000..e00b8cde8 --- /dev/null +++ b/graspologic/inference/group_connection_test.py @@ -0,0 +1,471 @@ +import warnings +from collections import namedtuple +from typing import Union + +import numpy as np +import pandas as pd +from beartype import beartype +from scipy.stats import combine_pvalues +from statsmodels.stats.multitest import multipletests + +from ..types import AdjacencyMatrix, List +from ..utils import import_graph, is_loopless, is_symmetric, is_unweighted, remove_loops +from .binomial import BinomialTestMethod, binom_2samp +from .utils import compute_density + +Labels = Union[np.ndarray, List] + +SBMResult = namedtuple( + "SBMResult", ["probabilities", "observed", "possible", "group_counts"] +) + +GroupTestResult = namedtuple("GroupTestResult", ["stat", "pvalue", "misc"]) + + +def fit_sbm(A: AdjacencyMatrix, labels: Labels, loops: bool = False) -> SBMResult: + """ + Fits a stochastic block model to data for a given network with known group + identities. Required inputs are the adjacency matrix for the + network and the group label for each node of the network. The number of labels must + equal the number of nodes of the network. + + For each possible group-to-group connection, e.g. group 1-to-group 1, + group 1-to-group 2, etc., this function computes the total number + of possible edges between the groups, the actual number of edges connecting the + groups, and the estimated probability of an edge + connecting each pair of groups (i.e. B_hat). The function also calculates and + returns the total number of nodes corresponding to each + provided label. + + Parameters + ---------- + A: np.array, int shape(n1,n1) + The adjacency matrix for the network at issue. Entries are either 1 + (edge present) or 0 (edge absent). This is a square matrix with + side length equal to the number of nodes in the network. + + labels: array-like, int shape(n1,1) + The group labels corresponding to each node in the network. This is a + one-dimensional array with a number of entries equal to the + number of nodes in the network. + + loops: boolean + This parameter instructs the function to either include or exclude self-loops + (i.e. connections between a node and itself) when + fitting the SBM model. This parameter is optional; default is false, meaning + self-loops will be excluded. + + Returns + ------- + SBMResult: namedtuple + This function returns a namedtuple with four key/value pairs: + ("probabilities",B_hat): + B_hat: array-like, float shape((number of unique labels)^2,1) + This variable stores the computed edge probabilities for all + possible group-to-group connections, computed as the ratio + of number of actual edges to number of possible edges. + ("observed",n_observed): + n_observed: dataframe + This variable stores the number of observed edges for each + group-to-group connection. Data is indexed as the number + of edges between each source group and each target group. + ("possibe",n_possible): + n_possible: dataframe + This variable stores the total number of possible edges for each + group-to-group connection. Indexing is identical to + the above. Network density for each group-to-group connection can + easily be determined as n_observed/n_possible. + ("group_counts",counts_labels): + counts_labels: pd.series + This variable stores the number of nodes belonging to each group + label. + + """ + + if not loops: + A = remove_loops(A) + + n = A.shape[0] + + node_to_comm_map = dict(zip(np.arange(n), labels)) + + # map edges to their incident communities + source_inds, target_inds = np.nonzero(A) + comm_mapper = np.vectorize(node_to_comm_map.get) + source_comm = comm_mapper(source_inds) + target_comm = comm_mapper(target_inds) + + # get the total number of possible edges for each community -> community cell + unique_labels, counts_labels = np.unique(labels, return_counts=True) + K = len(unique_labels) + + n_observed = ( + pd.crosstab( + source_comm, + target_comm, + dropna=False, + rownames=["source"], + colnames=["target"], + ) + .reindex(index=unique_labels, columns=unique_labels) + .fillna(0.0) + ) + + num_possible = np.outer(counts_labels, counts_labels) + + if not loops: + # then there would have been n fewer possible edges + num_possible[np.arange(K), np.arange(K)] = ( + num_possible[np.arange(K), np.arange(K)] - counts_labels + ) + + n_possible = pd.DataFrame( + data=num_possible, index=unique_labels, columns=unique_labels + ) + n_possible.index.name = "source" + n_possible.columns.name = "target" + + B_hat = np.divide(n_observed, n_possible) + + counts_labels = pd.Series(index=unique_labels, data=counts_labels) + + return SBMResult(B_hat, n_observed, n_possible, counts_labels) + + +def _make_adjacency_dataframe(data: AdjacencyMatrix, index: Labels) -> pd.DataFrame: + """ + Helper function to convert data with a given index into a dataframe data structure. + """ + df = pd.DataFrame(data=data, index=index.copy(), columns=index.copy()) + df.index.name = "source" + df.columns.name = "target" + return df + + +@beartype +def group_connection_test( + A1: AdjacencyMatrix, + A2: AdjacencyMatrix, + labels1: Labels, + labels2: Labels, + density_adjustment: Union[bool, float] = False, + method: BinomialTestMethod = "score", + combine_method: str = "tippett", + correct_method: str = "bonferroni", + alpha: float = 0.05, +) -> GroupTestResult: + r""" + Compares two networks by testing whether edge probabilities between groups are + significantly different for the two networks under a stochastic block model + assumption. + + This function requires the group labels in both networks to be known and to have the + same categories; although the exact number of nodes belonging to each group does not + need to be identical. Note that using group labels inferred from the data may + yield an invalid test. + + This function also permits the user to test whether one network's group connection + probabilities are a constant multiple of the other's (see ``density_adjustment`` + parameter). + + Parameters + ---------- + A1: np.array, shape(n1,n1) + The adjacency matrix for network 1. Will be treated as a binary network, + regardless of whether it was weighted. + A2 np.array, shape(n2,n2) + The adjacency matrix for network 2. Will be treated as a binary network, + regardless of whether it was weighted. + labels1: array-like, shape (n1,) + The group labels for each node in network 1. + labels2: array-like, shape (n2,) + The group labels for each node in network 2. + density_adjustment: boolean, optional + Whether to perform a density adjustment procedure. If ``True``, will test the + null hypothesis that the group-to-group connection probabilities of one network + are a constant multiple of those of the other network. Otherwise, no density + adjustment will be performed. + method: str, optional + Specifies the statistical test to be performed to compare each of the + group-to-group connection probabilities. By default, this performs + the score test (essentially equivalent to chi-squared test when + ``density_adjustment=False``), but the user may also enter "chi2" to perform the + chi-squared test, or "fisher" for Fisher's exact test. + combine_method: str, optional + Specifies the method for combining p-values (see Notes and [1]_ for more + details). Default is "tippett" for Tippett's method (recommended), but the user + can also enter any other method supported by + :func:`scipy.stats.combine_pvalues`. + correct_method: str, optional + Specifies the method for correcting for multiple comparisons. Default value is + "holm" to use the Holm-Bonferroni correction method, but + many others are possible (see :func:`statsmodels.stats.multitest.multipletests` + for more details and options). + alpha: float, optional + The significance threshold. By default, this is the conventional value of + 0.05 but any value on the interval :math:`[0,1]` can be entered. This only + affects the results in ``misc['rejections']``. + + Returns + ------- + GroupTestResult: namedtuple + A tuple containing the following data: + + stat: float + The statistic computed by the method chosen for combining + p-values (see ``combine_method``). + pvalue: float + The p-value for the overall network-to-network comparison using under a + stochastic block model assumption. Note that this is the p-value for the + comparison of the entire group-to-group connection matrices + (i.e., :math:`B_1` and :math:`B_2`). + misc: dict + A dictionary containing a number of statistics relating to the individual + group-to-group connection comparisons. + + "uncorrected_pvalues", pd.DataFrame + The p-values for each group-to-group connection comparison, before + correction for multiple comparisons. + "stats", pd.DataFrame + The test statistics for each of the group-to-group comparisons, + depending on ``method``. + "probabilities1", pd.DataFrame + This contains the B_hat values computed in fit_sbm above for + network 1, i.e. the hypothesized group connection density for + each group-to-group connection for network 1. + "probabilities2", pd.DataFrame + Same as above, but for network 2. + "observed1", pd.DataFrame + The total number of observed group-to-group edge connections for + network 1. + "observed2", pd.DataFrame + Same as above, but for network 2. + "possible1", pd.DataFrame + The total number of possible edges for each group-to-group pair in + network 1. + "possible2", pd.DataFrame + Same as above, but for network 2. + "group_counts1", pd.Series + Contains total number of nodes corresponding to each group label for + network 1. + "group_counts2", pd.Series + Same as above, for network 2 + "null_ratio", float + If the "density adjustment" parameter is set to "true", this + variable contains the null hypothesis for the quotient of + odds ratios for the group-to-group connection densities for the two + networks. In other words, it contains the hypothesized + factor by which network 1 is "more dense" or "less dense" than + network 2. If "density adjustment" is set to "false", this + simply returns a value of 1.0. + "n_tests", int + This variable contains the number of group-to-group comparisons + performed by the function. + "rejections", pd.DataFrame + Contains a square matrix of boolean variables. The side length of + the matrix is equal to the number of distinct group + labels. An entry in the matrix is "true" if the null hypothesis, + i.e. that the group-to-group connection density + corresponding to the row and column of the matrix is equal for both + networks (with or without a density adjustment factor), + is rejected. In simpler terms, an entry is only "true" if the + group-to-group density is statistically different between + the two networks for the connection from the group corresponding to + the row of the matrix to the group corresponding to the + column of the matrix. + "corrected_pvalues", pd.DataFrame + Contains the p-values for the group-to-group connection densities + after correction using the chosen correction_method. + + Notes + ----- + Under a stochastic block model assumption, the probability of observing an edge from + any node in group :math:`i` to any node in group :math:`j` is given by + :math:`B_{ij}`, where :math:`B` is a :math:`K \times K` matrix of connection + probabilities if there are :math:`K` groups. This test assumes that both networks + came from a stochastic block model with the same number of groups, and a fixed + assignment of nodes to groups. The null hypothesis is that the group-to-group + connection probabilities are the same + + .. math:: H_0: B_1 = B_2 + + The alternative hypothesis is that they are not the same + + .. math:: H_A: B_1 \neq B_2 + + Note that this alternative includes the case where even just one of these + group-to-group connection probabilities are different between the two networks. The + test is conducted by first comparing each group-to-group connection via its own + test, i.e., + + .. math:: H_0: {B_{1}}_{ij} = {B_{2}}_{ij} + + .. math:: H_A: {B_{1}}_{ij} \neq {B_{2}}_{ij} + + The p-values for each of these individual comparisons are stored in + ``misc['uncorrected_pvalues']``, and after multiple comparisons correction, in + ``misc['corrected_pvalues']``. The test statistic and p-value returned by this + test are for the overall comparison of the entire group-to-group connection + matrices. These are computed by appropriately combining the p-values for each of + the individual comparisons. For more details, see [1]_. + + When ``density_adjustment`` is set to ``True``, the null hypothesis is adjusted to + account for the fact that the group-to-group connection densities may be different + only up to a multiplicative factor which sets the densities of the two networks + the same in expectation. In other words, the null and alternative hypotheses are + adjusted to be + + .. math:: H_0: B_1 = c B_2 + + .. math:: H_A: B_1 \neq c B_2 + + where :math:`c` is a constant which sets the densities of the two networks the same. + + Note that in cases where one of the networks has no edges in a particular + group-to-group connection, it is nonsensical to run a statistical test for that + particular connection. In these cases, the p-values for that individual comparison + are set to ``np.nan``, and that test is not included in the overall test statistic + or multiple comparison correction. + + This test makes several assumptions about the data and test (which could easily be + loosened in future versions): + + We assume that the networks are directed. If the networks are undirected + (and the adjacency matrices are thus symmetric), then edges would be counted + twice, which would lead to an incorrect calculation of the edge probability. + We believe passing in the upper or lower triangle of the adjacency matrix + would solve this, but this has not been tested. + + We assume that the networks are loopless, that is we do not consider the + probability of an edge existing between a node and itself. This can be + weakened and made an option in future versions. + + We only implement the alternative hypothesis of "not equals" (two-sided); + future versions could implement the one-sided alternative hypotheses. + + References + ---------- + .. [1] Pedigo, B.D., Powell, M., Bridgeford, E.W., Winding, M., Priebe, C.E., + Vogelstein, J.T. "Generative network modeling reveals quantitative + definitions of bilateral symmetry exhibited by a whole insect brain + connectome," eLife (2023): e83739. + """ + + A1 = import_graph(A1) + A2 = import_graph(A2) + + if is_symmetric(A1) or is_symmetric(A2): + msg = ( + "This test assumes that the networks are directed, " + "but one or both adjacency matrices are symmetric." + ) + warnings.warn(msg) + if (not is_unweighted(A1)) or (not is_unweighted(A2)): + msg = ( + "This test assumes that the networks are unweighted, " + "but one or both adjacency matrices are weighted." + "Test will be run on the binarized version of these adjacency matrices." + ) + warnings.warn(msg) + if (not is_loopless(A1)) or (not is_loopless(A2)): + msg = ( + "This test assumes that the networks are loopless, " + "but one or both adjacency matrices have self-loops." + "Test will be run on the loopless version of these adjacency matrices." + ) + warnings.warn(msg) + + B1, n_observed1, n_possible1, group_counts1 = fit_sbm(A1, labels1) + B2, n_observed2, n_possible2, group_counts2 = fit_sbm(A2, labels2) + if not n_observed1.index.equals(n_observed2.index): + raise ValueError() + elif not n_observed1.columns.equals(n_observed2.columns): + raise ValueError() + elif not n_possible1.index.equals(n_possible2.index): + raise ValueError() + elif not n_observed1.columns.equals(n_observed2.columns): + raise ValueError() + + index = n_observed1.index.copy() + + if n_observed1.shape[0] != n_observed2.shape[0]: + raise ValueError() + + K = n_observed1.shape[0] + + uncorrected_pvalues_temp = np.empty( + (K, K), dtype=float + ) # had to make a new variable to keep mypy happy + uncorrected_pvalues = _make_adjacency_dataframe(uncorrected_pvalues_temp, index) + + stats_temp = np.empty((K, K), dtype=float) + stats = _make_adjacency_dataframe(stats_temp, index) + + if density_adjustment != False: # cause could be float + if density_adjustment == True: + density1 = compute_density(A1) + density2 = compute_density(A2) + adjustment_factor = density1 / density2 + else: + adjustment_factor = density_adjustment + else: + adjustment_factor = 1.0 + + for i in index: + for j in index: + curr_stat, curr_pvalue = binom_2samp( + n_observed1.loc[i, j], + n_possible1.loc[i, j], + n_observed2.loc[i, j], + n_possible2.loc[i, j], + method=method, + null_ratio=adjustment_factor, + ) + uncorrected_pvalues.loc[i, j] = curr_pvalue + stats.loc[i, j] = curr_stat + + misc = {} + misc["uncorrected_pvalues"] = uncorrected_pvalues + misc["stats"] = stats + misc["probabilities1"] = B1 + misc["probabilities2"] = B2 + misc["observed1"] = n_observed1 + misc["observed2"] = n_observed2 + misc["possible1"] = n_possible1 + misc["possible2"] = n_possible2 + misc["group_counts1"] = group_counts1 + misc["group_counts2"] = group_counts2 + misc["null_ratio"] = adjustment_factor + + run_pvalues = uncorrected_pvalues.values + indices = np.nonzero(~np.isnan(run_pvalues)) + run_pvalues = run_pvalues[indices] + n_tests = len(run_pvalues) + misc["n_tests"] = n_tests + + # correct for multiple comparisons + rejections_flat, corrected_pvalues_flat, _, _ = multipletests( + run_pvalues, + alpha, + method=correct_method, + is_sorted=False, + returnsorted=False, + ) + rejections = np.full((K, K), False, dtype=bool) + rejections[indices] = rejections_flat + rejections = _make_adjacency_dataframe(rejections, index) + misc["rejections"] = rejections + + corrected_pvalues = np.full((K, K), np.nan, dtype=float) + corrected_pvalues[indices] = corrected_pvalues_flat + corrected_pvalues = _make_adjacency_dataframe(corrected_pvalues, index) + misc["corrected_pvalues"] = corrected_pvalues + + # combine p-values (on the UNcorrected p-values) + if run_pvalues.min() == 0.0: + # TODO consider raising a new warning here + stat = np.inf + pvalue = 0.0 + else: + stat, pvalue = combine_pvalues(run_pvalues, method=combine_method) + return GroupTestResult(stat, pvalue, misc) diff --git a/graspologic/inference/utils.py b/graspologic/inference/utils.py new file mode 100644 index 000000000..6f4befbcd --- /dev/null +++ b/graspologic/inference/utils.py @@ -0,0 +1,63 @@ +import numpy as np + +from ..types import AdjacencyMatrix + + +def compute_density(adjacency: AdjacencyMatrix, loops: bool = False) -> float: + """ + For a given graph, this function computes the graph density, defined as the actual number of edges divided by the total possible number + of edges in the graph. + + Parameters + ---------- + adjacency: int, array shape (n_nodes,n_nodes) + The adjancy matrix for the graph. Edges are denoted by 1s while non-edges are denoted by 0s. + + loops: boolean + Optional variable to select whether to include self-loops (i.e. connections between a node and itself). Default is "false," meaning + such connections are ignored. + + Returns + ------- + n_edges/n_possible: float + The computed density, calculated as the total number of edges divided by the total number of possible edges. + + """ + n_edges = np.count_nonzero(adjacency) + n_nodes = adjacency.shape[0] + n_possible = n_nodes**2 + if not loops: + n_possible -= n_nodes + return n_edges / n_possible + + +def compute_density_adjustment( + adjacency1: AdjacencyMatrix, adjacency2: AdjacencyMatrix +) -> float: + """ + Computes the density adjustment to be used when testing the hypothesis that the density of one network is equal to a fixed parameter + times the density of a second network. This function first calls the compute_density function above to compute the densities of both + networks, then computes an odds ratio by calculating the odds of an edge in each network and taking the ratio of the results. + + Parameters + ---------- + adjacency1: int, array of size (n_nodes1,n_nodes1) + Adjacency matrix for the first graph. 1s represent edges while 0s represent the absence of an edge. The array is a square of side length + n_nodes1, where this corresponds to the number of nodes in graph 1. + + adjacency2: int, array of size (n_nodes2,n_nodes2) + Same as above, but for the second graph. + + Returns + --------- + odds_ratio: float + Computed as the ratio of the odds of an edge in graph 1 to the odds of an edge in graph 2. + + """ + density1 = compute_density(adjacency1) + density2 = compute_density(adjacency2) + # return density1 / density2 + odds1 = density1 / (1 - density1) + odds2 = density2 / (1 - density2) + odds_ratio = odds1 / odds2 + return odds_ratio diff --git a/graspologic/nominate/spectralVN.py b/graspologic/nominate/spectralVN.py index 239fa7886..16578f59e 100644 --- a/graspologic/nominate/spectralVN.py +++ b/graspologic/nominate/spectralVN.py @@ -111,7 +111,7 @@ def _check_x(self, X: np.ndarray) -> None: def _check_y(self, y: np.ndarray) -> None: # check y - if not np.issubdtype(y.dtype, np.integer): + if not np.issubdtype(y.dtype, int): raise TypeError("y must have dtype int") elif np.ndim(y) > 2 or (y.ndim == 2 and y.shape[1] > 1): raise IndexError("y must have shape (n) or (n, 1).") @@ -223,12 +223,12 @@ def predict(self, y: Union[List, np.ndarray]) -> Tuple[np.ndarray, np.ndarray]: y_checked: np.ndarray = check_array(y, ensure_2d=False) self._check_y(y_checked) y_checked = y_checked.reshape(-1) - y_vec = self.embedding_[y_checked.astype(np.int)] # type: ignore + y_vec = self.embedding_[y_checked.astype(int)] # type: ignore if not hasattr(self, "nearest_neighbors_"): raise ValueError("Fit must be called before predict.") distance_matrix, nomination_list = self.nearest_neighbors_.kneighbors(y_vec) # transpose for consistency with literature - return nomination_list.T.astype(np.int_), distance_matrix.T + return nomination_list.T.astype(int), distance_matrix.T def fit_predict( self, diff --git a/mypy.ini b/mypy.ini index 00412aa2a..09fa47789 100644 --- a/mypy.ini +++ b/mypy.ini @@ -56,5 +56,8 @@ ignore_missing_imports = True [mypy-sklearn.*] ignore_missing_imports = True +[mypy-statsmodels.*] +ignore_missing_imports = True + [mypy-umap] ignore_missing_imports = True diff --git a/runtime.txt b/runtime.txt new file mode 100644 index 000000000..cc1923a40 --- /dev/null +++ b/runtime.txt @@ -0,0 +1 @@ +3.8 diff --git a/setup.cfg b/setup.cfg index 3160a6621..0a21b45cd 100644 --- a/setup.cfg +++ b/setup.cfg @@ -38,7 +38,8 @@ install_requires = POT>=0.7.0 seaborn>= 0.11.0 scikit-learn>=0.22.0 - scipy>=1.4.0 + scipy>=1.9.0 + statsmodels>=0.13.2 typing-extensions>=4.4.0 umap-learn>=0.4.6 diff --git a/tests/test_er_and_group_connection_tests.py b/tests/test_er_and_group_connection_tests.py new file mode 100644 index 000000000..635bc4fbd --- /dev/null +++ b/tests/test_er_and_group_connection_tests.py @@ -0,0 +1,85 @@ +import unittest + +import numpy as np +from scipy.sparse import csr_array + +from graspologic.inference import density_test, group_connection_test +from graspologic.simulations import er_np, sbm + + +class TestGroupConnection(unittest.TestCase): + def test_gctest_works(self): + np.random.seed(8888) + B1 = np.array([[0.8, 0.6], [0.6, 0.8]]) + B2 = 0.8 * B1 + A1, labels1 = sbm([50, 50], B1, return_labels=True) + A2, labels2 = sbm([60, 60], B2, return_labels=True) + stat, pvalue, misc = group_connection_test( + A1, A2, labels1, labels2, density_adjustment=True + ) + self.assertTrue(pvalue > 0.05) + + def test_all_kwargs(self): + B1 = np.array([[0.4, 0.6], [0.6, 0.8]]) + B2 = np.array([[0.9, 0.4], [0.2, 0.865]]) + A1, labels1 = sbm([60, 60], B1, return_labels=True, directed=True) + A2, labels2 = sbm([50, 50], B2, return_labels=True, directed=True) + stat, pvalue, misc = group_connection_test( + A1, + A2, + labels1, + labels2, + combine_method="tippett", + method="score", + correct_method="Bonferroni", + density_adjustment=True, + ) + self.assertTrue(pvalue < 0.05) + self.assertTrue(misc["uncorrected_pvalues"].size == 4) + self.assertTrue(misc["probabilities1"].size == 4) + self.assertTrue(misc["probabilities2"].size == 4) + self.assertTrue(np.sum(misc["observed1"].to_numpy()) == np.count_nonzero(A1)) + self.assertTrue(np.sum(misc["observed2"].to_numpy()) == np.count_nonzero(A2)) + self.assertTrue(misc["null_ratio"] != 1.0) + self.assertTrue(misc["n_tests"] == 4) + self.assertTrue(misc["rejections"].to_numpy().size == 4) + self.assertTrue(misc["corrected_pvalues"].size == 4) + + def test_sparse(self): + B1 = np.array([[0.8, 0.6], [0.6, 0.8]]) + B2 = np.array([[0.87, 0.66], [0.66, 0.87]]) + A1, labels1 = sbm([50, 50], B1, return_labels=True) + A2, labels2 = sbm([60, 60], B2, return_labels=True) + sA1 = csr_array(A1) + sA2 = csr_array(A2) + + stat, pvalue, misc = group_connection_test(sA1, sA2, labels1, labels2) + self.assertTrue(pvalue <= 0.05) + + +class TestER(unittest.TestCase): + def test_er(self): + np.random.seed(234) + A1 = er_np(500, 0.6) + A2 = er_np(400, 0.8) + stat, pvalue, er_misc = density_test(A1, A2) + self.assertTrue(pvalue <= 0.05) + A3 = er_np(500, 0.8) + A4 = er_np(400, 0.8) + stat, pvalue, er_misc = density_test(A3, A4) + self.assertTrue(pvalue > 0.05) + + def test_all(self): + np.random.seed(234) + A1 = er_np(500, 0.6) + A2 = er_np(400, 0.8) + stat, pvalue, er_misc = density_test(A1, A2, method="chi2") + self.assertTrue(pvalue <= 0.05) + self.assertTrue(er_misc["probability1"].to_numpy() < 1.0) + self.assertTrue(er_misc["probability2"].to_numpy() < 1.0) + self.assertTrue(er_misc["observed1"].to_numpy() == np.count_nonzero(A1)) + self.assertTrue(er_misc["observed2"].to_numpy() == np.count_nonzero(A2)) + + +if __name__ == "__main__": + unittest.main()