Repository: networkx/networkx
Pull request: Fix numeric and degree assortativity coefficient calculation
It fixes #4776 and wrong assortativity coefficient calculation in some cases that follow PR #4851.

A cause of the issue is in `numeric_ac` in [correlation.py](https://github.com/networkx/networkx/blob/main/networkx/algorithms/assortativity/correlation.py). Before my fixes in `numeric_ac`
```python
x = np.arange(nx)
y = np.arange(ny)
```
we create supports for random variables x and y. We assume that node values are integer and coincide with mixing matrix indices — 0 is in the first row/column, 1 is in the second row/column and so on. But after PR #4851, node values can also be negative or float. In such cases we need to know a mixing matrix indices for each unique pair of values. After my fixes, supports of x and y contain all node values, not only integers. The same for degree assortativity coefficient calculation.

Note that the PR #4851 version also correctly calculates assortativity coefficient in the case of *equal intervals* between support values due to normalization (in fact, this is a correlation between x and y, not covariance). For example, it correctly calculates a degree assortativity coefficient for float node degrees `[0, 0.5, 1, 1.5]` or an assortativity coefficient for mixed node values `[-1, 0, 1, 2]`.  But in #4776 we have a graph with node values `[-1, 1, 2]` and then the old coefficient is incorrect. I added tests that check cases with non-equal intervals between support values — one for degree and one for node values.

Problem statement:
numeric assortativity coefficient only works for nonnegative numeric attributes
<!-- If you have a general question about NetworkX, please use the discussions tab to create a new discussion -->

<!--- Provide a general summary of the issue in the Title above -->


### Current Behavior
<!--- Tell us what happens instead of the expected behavior -->

If the range of numeric attributes ranges from negative to positive (and not all negative), `nx.numeric_mixing_matrix` effectively disregards it since it is generated using `mapping = {x: x for x in range(m + 1)}`.

### Expected Behavior
<!--- Tell us what should happen -->
The negative integers should be included in generating the matrix. Alternatively:
1. just use the Pearson correlation coefficient as given in Peel et al, eq. S4-S5 (Leto Peel, Jean-Charles Delvenne, Renaud Lambiotte, Multiscale mixing patterns in networks, PNAS Apr 2018, 115 (16) 4057-4062; DOI: 10.1073/pnas.1713019115). This has the added benefit that the scalar quantities need not be integers.
2. require that the integers be nonnegative
3. increase the size of the matrix so it includes the negative integers as well.



### Steps to Reproduce
<!--- Provide a minimal example that reproduces the bug -->

```
G = nx.Graph()
G.add_nodes_from([1,2], a=-1)
G.add_nodes_from([3], a=1)
G.add_nodes_from([4], a=2)
G.add_edges_from([(3,4),(1,2),(1,3)])
# numeric_assortativity_coefficient gives -1.0
print(nx.numeric_assortativity_coefficient(G, 'a'))
```

### Environment
<!--- Please provide details about your local environment -->
Python version: 3.8
NetworkX version: 2.5


### Additional context
<!--- Add any other context about the problem here, screenshots, etc. -->

Hints from issue discussion:
I don't think `mapping = {x: x for x in range(m + 1)}` affects or is affected by the attributes of nodes.
This constructs a mapping of nodes to integer indexes in the matrix. In this case it is mapping from the indexes created earlier to the indexes used here -- so it is the identity mapping.

The numeric attributes can range with both positive and negative values.  Also, your example gives coefficient -1, which seems to indicate that it is not ignoring the negative values of attribute `a`. What value would you expect?

Can you clarify your description? 
I actually discovered this error using a larger dataset, comparing the Pearson correlation coefficient (directly from scipy) with the one obtained using `numeric_assortativity_coefficient`. I tried taking the absolute value of the attribute and the discrepancy disappeared.

I think the mapping is the culprit. `numeric_mixing_matrix` takes the dictionary `attribute_mixing_dict`which has _attribute values_ (not node labels) as the keys and the values being another dictionary, with the keys being the neighbor attribute values and the dictionary values being the counts of how many times these attribute values occur in an edge.

Therefore, in `numeric_mixing_matrix`, the value of `m` is the maximum of the attribute values (not the node labels).

I've annotated the snippet below to point out what I think is wrong with the code. I hope that helps.

```
def numeric_mixing_matrix(G, attribute, nodes=None, normalized=True):
    """Returns numeric mixing matrix for attribute.

    The attribute must be an integer.

    Parameters
    ----------
    G : graph
       NetworkX graph object.

    attribute : string
       Node attribute key.  The corresponding attribute must be an integer.

    nodes: list or iterable (optional)
        Build the matrix only with nodes in container. The default is all nodes.

    normalized : bool (default=True)
       Return counts if False or probabilities if True.

    Returns
    -------
    m: numpy array
       Counts, or joint, probability of occurrence of node attribute pairs.
    """
    d = attribute_mixing_dict(G, attribute, nodes) # {attr: {attr: count}}
    s = set(d.keys())  # effectively a set of attribute _values_
    for k, v in d.items():
        s.update(v.keys())
    m = max(s) # the maximum of attribute _values_
    mapping = {x: x for x in range(m + 1)} # begins from zero to m, so if there are negative values, they won't be included
    a = dict_to_numpy_array(d, mapping=mapping)
    if normalized:
        a = a / a.sum()
    return a
```
Sorry -- my mistake -- I didn't understand the `mixing _dict` data structure.

So, would a fix be to replace `range(m + 1)` with `range(min(s), max(s)+1)`?

Also, what value would you expect for your example instead of -1?
Peel et al 2018 published a paper in PNAS about assortativity, and they did include the formula for scalar attributes in the Supplementary Information. Based on this, the expected answer is `r=4/9`.

The advantage of Peel's 2018 formula is that it doesn't require integer values for the attribute values.

Leto Peel, Jean-Charles Delvenne, Renaud Lambiotte, Multiscale mixing patterns in networks, PNAS Apr 2018, 115 (16) 4057-4062; DOI: 10.1073/pnas.1713019115

Commit message: Reading node degrees once

Patch to review:
diff --git a/networkx/algorithms/assortativity/correlation.py b/networkx/algorithms/assortativity/correlation.py
index b885a22ee1b..7fe1b2324a8 100644
--- a/networkx/algorithms/assortativity/correlation.py
+++ b/networkx/algorithms/assortativity/correlation.py
@@ -73,8 +73,12 @@ def degree_assortativity_coefficient(G, x="out", y="in", weight=None, nodes=None
     .. [2] Foster, J.G., Foster, D.V., Grassberger, P. & Paczuski, M.
        Edge direction and the structure of networks, PNAS 107, 10815-20 (2010).
     """
-    M = degree_mixing_matrix(G, x=x, y=y, nodes=nodes, weight=weight)
-    return numeric_ac(M)
+    if nodes is None:
+        nodes = G.nodes
+    degrees = set([d for n, d in G.degree(nodes, weight=weight)])
+    mapping = {d: i for i, d, in enumerate(degrees)}
+    M = degree_mixing_matrix(G, x=x, y=y, nodes=nodes, weight=weight, mapping=mapping)
+    return numeric_ac(M, mapping=mapping)
 
 
 def degree_pearson_correlation_coefficient(G, x="out", y="in", weight=None, nodes=None):
@@ -223,8 +227,12 @@ def numeric_assortativity_coefficient(G, attribute, nodes=None):
     .. [1] M. E. J. Newman, Mixing patterns in networks
            Physical Review E, 67 026126, 2003
     """
-    a = numeric_mixing_matrix(G, attribute, nodes)
-    return numeric_ac(a)
+    if nodes is None:
+        nodes = G.nodes
+    vals = set(G.nodes[n][attribute] for n in nodes)
+    mapping = {d: i for i, d, in enumerate(vals)}
+    M = attribute_mixing_matrix(G, attribute, nodes, mapping)
+    return numeric_ac(M, mapping)
 
 
 def attribute_ac(M):
@@ -254,7 +262,7 @@ def attribute_ac(M):
     return r
 
 
-def numeric_ac(M):
+def numeric_ac(M, mapping):
     # M is a numpy matrix or array
     # numeric assortativity coefficient, pearsonr
     import numpy as np
@@ -262,12 +270,13 @@ def numeric_ac(M):
     if M.sum() != 1.0:
         M = M / float(M.sum())
     nx, ny = M.shape  # nx=ny
-    x = np.arange(nx)
-    y = np.arange(ny)
+    x = np.array(list(mapping.keys()))
+    y = x  # x and y have the same support
+    idx = list(mapping.values())
     a = M.sum(axis=0)
     b = M.sum(axis=1)
-    vara = (a * x ** 2).sum() - ((a * x).sum()) ** 2
-    varb = (b * x ** 2).sum() - ((b * x).sum()) ** 2
+    vara = (a[idx] * x ** 2).sum() - ((a[idx] * x).sum()) ** 2
+    varb = (b[idx] * x ** 2).sum() - ((b[idx] * x).sum()) ** 2
     xy = np.outer(x, y)
-    ab = np.outer(a, b)
+    ab = np.outer(a[idx], b[idx])
     return (xy * (M - ab)).sum() / np.sqrt(vara * varb)
diff --git a/networkx/algorithms/assortativity/tests/base_test.py b/networkx/algorithms/assortativity/tests/base_test.py
index 5b3712310ee..ee5b126c071 100644
--- a/networkx/algorithms/assortativity/tests/base_test.py
+++ b/networkx/algorithms/assortativity/tests/base_test.py
@@ -52,6 +52,10 @@ def setup_class(cls):
         cls.W = nx.Graph()
         cls.W.add_edges_from([(0, 3), (1, 3), (2, 3)], weight=0.5)
         cls.W.add_edge(0, 2, weight=1)
+        S1 = nx.star_graph(4)
+        S2 = nx.star_graph(4)
+        cls.DS = nx.disjoint_union(S1, S2)
+        cls.DS.add_edge(4, 5)
 
 
 class BaseTestNumericMixing:
@@ -70,3 +74,10 @@ def setup_class(cls):
         F.add_edge(0, 2, weight=1)
         nx.set_node_attributes(F, dict(F.degree(weight="weight")), "margin")
         cls.F = F
+
+        M = nx.Graph()
+        M.add_nodes_from([1, 2], margin=-1)
+        M.add_nodes_from([3], margin=1)
+        M.add_nodes_from([4], margin=2)
+        M.add_edges_from([(3, 4), (1, 2), (1, 3)])
+        cls.M = M
diff --git a/networkx/algorithms/assortativity/tests/test_correlation.py b/networkx/algorithms/assortativity/tests/test_correlation.py
index 763f0a41571..8da72652f17 100644
--- a/networkx/algorithms/assortativity/tests/test_correlation.py
+++ b/networkx/algorithms/assortativity/tests/test_correlation.py
@@ -42,6 +42,10 @@ def test_degree_assortativity_weighted(self):
         r = nx.degree_assortativity_coefficient(self.W, weight="weight")
         np.testing.assert_almost_equal(r, -0.1429, decimal=4)
 
+    def test_degree_assortativity_double_star(self):
+        r = nx.degree_assortativity_coefficient(self.DS)
+        np.testing.assert_almost_equal(r, -0.9339, decimal=4)
+
 
 class TestAttributeMixingCorrelation(BaseTestAttributeMixing):
     def test_attribute_assortativity_undirected(self):
@@ -91,3 +95,7 @@ def test_numeric_assortativity_negative(self):
     def test_numeric_assortativity_float(self):
         r = nx.numeric_assortativity_coefficient(self.F, "margin")
         np.testing.assert_almost_equal(r, -0.1429, decimal=4)
+
+    def test_numeric_assortativity_mixed(self):
+        r = nx.numeric_assortativity_coefficient(self.M, "margin")
+        np.testing.assert_almost_equal(r, 0.4340, decimal=4)


Code review comments:
diff_hunk: @@ -254,20 +262,21 @@ def attribute_ac(M):
     return r
 
 
-def numeric_ac(M):
+def numeric_ac(M, mapping):
     # M is a numpy matrix or array
     # numeric assortativity coefficient, pearsonr
     import numpy as np
 
     if M.sum() != 1.0:
         M = M / float(M.sum())
     nx, ny = M.shape  # nx=ny
-    x = np.arange(nx)
-    y = np.arange(ny)
+    x = np.array(list(mapping.keys()))
+    y = x  # x and y have the same support
+    idx = list(mapping.values())
     a = M.sum(axis=0)
     b = M.sum(axis=1)
-    vara = (a * x ** 2).sum() - ((a * x).sum()) ** 2
-    varb = (b * x ** 2).sum() - ((b * x).sum()) ** 2
+    vara = (a[idx] * x ** 2).sum() - ((a[idx] * x).sum()) ** 2
+    varb = (b[idx] * x ** 2).sum() - ((b[idx] * x).sum()) ** 2
path: networkx/algorithms/assortativity/correlation.py
text: @user1:
I know the values are the same, but could you change “x” to “y” in “varb”?
This also makes me wonder if the mixing matrix for a directed graph ought to have in_degree for rows and out_degree for columns.  But I think that is beyond this issue.  It might be a topic of a whole paper on mixing over graphs (if it hasn’t been written already).

@author:
Ok, I will change it.
In the current implementation, a mixing matrix is squared and supports (sets of possible values) for in_degree and out_degree are the same. But there is another way to represent a mixing matrix in a non-squared view. That is as you said, in_degree is for rows and out_degree is for columns (or vice versa). Both are correct, and do not affect assortativity coefficient. The disadvantage of a non-squared view is that it requires separate mappings for rows and columns, say two variables `out_mapping` and `in_mapping` instead of a single `mapping`. The advantage of a non-squared view is that a matrix is much smaller in some limit cases, say in a graph with many trees. For example, a graph

<img src='https://user-images.githubusercontent.com/12996098/123511014-80841180-d687-11eb-92d7-4b8de61bad3f.png' width=400>

has a 3x3 degree mixing matrix 
```
mapping={1: 0, 3: 1, 4: 2}
```
```
[[0.   0.   0.  ]
 [0.43 0.   0.  ]
 [0.57 0.   0.  ]]
```
or a 2x1 degree mixing matrix
```
out_mapping = {3: 0, 4: 1}
in_mapping = {1: 0}
```
```
[[0.43]
 [0.57]]
```
As we see, if we increase the number of such trees with different number of leafs, the non-squared mixing matrix will still be just a one-column matrix.

@user1:
Nice description -- Thanks! 

I think this is ready to merge.  Anything else you know of? 

@author:
Yes, it is ready to merge