<jupyter_start><jupyter_text>Introduction<jupyter_code>from IPython.display import YouTubeVideo

YouTubeVideo(id="P-0CJpO3spg", width="100%")<jupyter_output><empty_output><jupyter_text>In this chapter, we are going to take a look at how to perform statistical inference on graphs.  Statistics refresherBefore we can proceed with statistical inference on graphs,we must first refresh ourselves with some ideas from the world of statistics.Otherwise, the methods that we will end up usingmay seem a tad _weird_, and hence difficult to follow along.To review statistical ideas,let's set up a few statements and explore what they mean.  We are concerned with models of randomnessAs with all things statistics, we are concerned with models of randomness.Here, probability distributions give us a way to think about random eventsand how to assign credibility points to them.  In an abstract fashion...The supremely abstract way of thinking about a probability distributionis that it is the space of all possibilities of "stuff"with different credibility points _distributed_ amongst each possible "thing".  More concretely: the coin flipA more concrete example is to consider the coin flip.Here, the space of all possibilities of "stuff" is the set of "heads" and "tails".If we have a fair coin, then we have 0.5 credibility points _distributed_to each of "heads" and "tails".  Another example: dice rollsAnother concrete example is to consider the six-sided dice.Here, the space of all possibilities of "stuff" is the set of numbers in the range $[1, 6]$.If we have a fair dice, then we have 1/6 credibility points assignedto each of the numbers.(Unfair dice will have an unequal _distribution_ of credibility points across each face.)  A graph-based example: social networksIf we receive an undirected social network graph with 5 nodes and 6 edges,we have to keep in mind that this graph with 6 edgeswas merely one of $15 \choose 6$ ways to construct 5 node, 6 edge graphs.(15 comes up because there are 15 edges that can be constructed in a 5-node undirected graph.)  Hypothesis TestingA commonplace task in statistical inferencesis calculating the probability of observing a value or something more extremeunder an assumed "null" model of reality.This is what we commonly call "hypothesis testing",and where the oft-misunderstood term "p-value" shows up.  Hypothesis testing in coin flips, by simulationAs an example, hypothesis testing in coin flips follows this logic:- I observe that 8 out of 10 coin tosses give me heads, giving me a probability of heads $p=0.8$ (a summary statistic).- Under a "null distribution" of a fair coin, I simulate the distribution of probability of heads (the summary statistic) that I would get from 10 coin tosses.- Finally, I use that distribution to calculate the probability of observing $p=0.8$ or more extreme.  Hypothesis testing in graphsThe same protocol applies when we perform hypothesis testing on graphs.Firstly, we calculate a _summary statistic_ that describes our graph.Secondly, we propose a _null graph model_, and calculate our summary statistic under simulated versions of that null graph model.Thirdly, we look at the probability of observing the summary statistic value that we calculated in step 1 or more extreme, under the assumed graph null model distribution.  Stochastic graph creation modelsSince we are going to be dealing with models of randomness in graphs,let's take a look at some examples.  Erdos-Renyi (a.k.a. "binomial") graphOn easy one to study is the Erdos-Renyi graph, also known as the "binomial" graph.The data generation story here is that we instantiate an undirected graph with $n$ nodes,giving $\frac{n^2 - n}{2}$ possible edges.Each edge has a probability $p$ of being created.<jupyter_code>import networkx as nx


G_er = nx.erdos_renyi_graph(n=30, p=0.2)
nx.draw(G_er)<jupyter_output><empty_output><jupyter_text>You can verify that there's approximately 20% of $\frac{30^2 - 30}{2} = 435$.<jupyter_code>len(G_er.edges())
len(G_er.edges()) / 435<jupyter_output><empty_output><jupyter_text>We can also look at the degree distribution:<jupyter_code>import pandas as pd
from nams.functions import ecdf
import matplotlib.pyplot as plt

x, y = ecdf(pd.Series(dict(nx.degree(G_er))))
plt.scatter(x, y)<jupyter_output><empty_output><jupyter_text>Barabasi-Albert GraphThe data generating story of this graph generator is essentially that nodes that have lots of edges preferentially get new edges attached onto them. This is what we call a "preferential attachment" process.<jupyter_code>G_ba = nx.barabasi_albert_graph(n=30, m=3)
nx.draw(G_ba)
len(G_ba.edges())<jupyter_output><empty_output><jupyter_text>And the degree distribution:<jupyter_code>x, y = ecdf(pd.Series(dict(nx.degree(G_ba))))
plt.scatter(x, y)<jupyter_output><empty_output><jupyter_text>You can see that even though the number of edges between the two graphs are similar,their degree distribution is wildly different.  Load DataFor this notebook, we are going to look at a protein-protein interaction network,and test the hypothesis that this network was _not_ generated by the data generating processdescribed by an Erdos-Renyi graph.Let's load a [protein-protein interaction network dataset](http://konect.uni-koblenz.de/networks/moreno_propro).> This undirected network contains protein interactions contained in yeast.> Research showed that proteins with a high degree> were more important for the surivial of the yeast than others.> A node represents a protein and an edge represents a metabolic interaction between two proteins. > The network contains loops.<jupyter_code>from nams import load_data as cf

G = cf.load_propro_network()
for n, d in G.nodes(data=True):
    G.nodes[n]["degree"] = G.degree(n)<jupyter_output><empty_output><jupyter_text>As is always the case, let's make sure we know some basic stats of the graph.<jupyter_code>len(G.nodes())
len(G.edges())<jupyter_output><empty_output><jupyter_text>Let's also examine the degree distribution of the graph.<jupyter_code>x, y = ecdf(pd.Series(dict(nx.degree(G))))
plt.scatter(x, y)<jupyter_output><empty_output><jupyter_text>Finally, we should visualize the graph to get a feel for it.<jupyter_code>import nxviz as nv
from nxviz import annotate

nv.circos(G, sort_by="degree", node_color_by="degree", node_aes_kwargs={"size_scale": 10})
annotate.node_colormapping(G, color_by="degree")<jupyter_output><empty_output><jupyter_text>One thing we might infer from this visualizationis that the vast majority of nodes have a very small degree,while a very small number of nodes have a high degree.That would prompt us to think:what process could be responsible for generating this graph?  Inferring Graph Generating ModelGiven a graph dataset, how do we identify which data generating model provides the best fit?One way to do this is to compare characteristics of a graph generating model against the characteristics of the graph.The logic here is that if we have a good graph generating model for the data,we should, in theory, observe the observed graph's characteristicsin the graphs generated by the graph generating model.  Comparison of degree distributionLet's compare the degree distribution between the data, a few Erdos-Renyi graphs, and a few Barabasi-Albert graphs.  Comparison with Barabasi-Albert graphs<jupyter_code>from ipywidgets import interact, IntSlider

m = IntSlider(value=2, min=1, max=10)

@interact(m=m)
def compare_barabasi_albert_graph(m):
    fig, ax = plt.subplots()
    G_ba = nx.barabasi_albert_graph(n=len(G.nodes()), m=m)
    x, y = ecdf(pd.Series(dict(nx.degree(G_ba))))
    ax.scatter(x, y, label="Barabasi-Albert Graph")
    
    x, y = ecdf(pd.Series(dict(nx.degree(G))))
    ax.scatter(x, y, label="Protein Interaction Network")
    ax.legend()<jupyter_output><empty_output><jupyter_text>Comparison with Erdos-Renyi graphs<jupyter_code>from ipywidgets import FloatSlider
p = FloatSlider(value=0.6, min=0, max=0.1, step=0.001)

@interact(p=p)
def compare_erdos_renyi_graph(p):
    fig, ax = plt.subplots()
    G_er = nx.erdos_renyi_graph(n=len(G.nodes()), p=p)
    x, y = ecdf(pd.Series(dict(nx.degree(G_er))))
    ax.scatter(x, y, label="Erdos-Renyi Graph")
    
    x, y = ecdf(pd.Series(dict(nx.degree(G))))
    ax.scatter(x, y, label="Protein Interaction Network")
    ax.legend()
    ax.set_title(f"p={p}")<jupyter_output><empty_output><jupyter_text>Given the degree distribution only, which model do you think better describes the generation of a protein-protein interaction network?  Quantitative Model ComparisonEach time we plug in a value of $m$ for the Barabasi-Albert graph model, we are using one of many possible Barabasi-Albert graph models, each with a different $m$.Similarly, each time we choose a different $p$ for the Erdos-Renyi model, we are using one of many possible Erdos-Renyi graph models, each with a different $p$.To quantitatively compare degree distributions, we can use the [Wasserstein distance][wasd] between the data.Let's see how to implement this.[wasd]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wasserstein_distance.html<jupyter_code>from scipy.stats import wasserstein_distance

def erdos_renyi_degdist(n, p):
    """Return a Pandas series of degree distribution of an Erdos-Renyi graph."""
    G = nx.erdos_renyi_graph(n=n, p=p)
    return pd.Series(dict(nx.degree(G)))

def barabasi_albert_degdist(n, m):
    """Return a Pandas series of degree distribution of an Barabasi-Albert graph."""
    G = nx.barabasi_albert_graph(n=n, m=m)
    return pd.Series(dict(nx.degree(G)))
deg = pd.Series(dict(nx.degree(G)))

er_deg = erdos_renyi_degdist(n=len(G.nodes()), p=0.001)
ba_deg = barabasi_albert_degdist(n=len(G.nodes()), m=1)
wasserstein_distance(deg, er_deg), wasserstein_distance(deg, ba_deg)<jupyter_output><empty_output><jupyter_text>Notice that because the graphs are instantiated in a non-deterministic fashion, re-running the cell above will give you different values for each new graph generated.Let's now plot the wasserstein distance to our graph data for the two particular Erdos-Renyi and Barabasi-Albert graph models shown above.<jupyter_code>import matplotlib.pyplot as plt
from tqdm.autonotebook import tqdm

er_dist = []
ba_dist = []
for _ in tqdm(range(100)):
    er_deg = erdos_renyi_degdist(n=len(G.nodes()), p=0.001)
    er_dist.append(wasserstein_distance(deg, er_deg))
    

    ba_deg = barabasi_albert_degdist(n=len(G.nodes()), m=1)
    ba_dist.append(wasserstein_distance(deg, ba_deg))

# er_degs = [erdos_renyi_degdist(n=len(G.nodes()), p=0.001) for _ in range(100)]
import seaborn as sns
import janitor


data = (
    pd.DataFrame(
        {
            "Erdos-Renyi": er_dist,
            "Barabasi-Albert": ba_dist,
        }
    )
    .melt(value_vars=["Erdos-Renyi", "Barabasi-Albert"])
    .rename_columns({"variable": "Graph Model", "value": "Wasserstein Distance"})
)
sns.swarmplot(data=data, x="Graph Model", y="Wasserstein Distance")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08472727272727272, 'fraction_numerical': 0.008363636363636363, 'mean_word_length': 5.407105416423995, 'pattern_counts': {'":': 5, '<': 57, '<?xml version=': 0, '>': 62, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '31112', 'license': 'MIT', 'n_tokens_mistral': 3471, 'n_tokens_neox': 3280, 'n_words': 1264, 'path': 'notebooks/04-advanced/03-stats.ipynb', 'repo_name': 'khanin-th/Network-Analysis-Made-Simple', 'size': '11000'}