Calculating Spearman Correlation from scratch

Question

Over the week, I was assigned to create a function in a Jupyter notebook that satisfies the Spearman Correlation formula without the use of any packages e.g. Numpy, Pandas. Below are the rules of the game, followed by solution.

I am seeking any input on how this could be done more efficiently.

The Spearman rank correlation is used to evaluate if the relationship between two variables, X and Y is monotonic. The rank correlation measures how closely related the ordering of one variable to the other variable, with no regard to the actual values of the variables. To calculate the Spearman rank correlation for \$X\$ and \$Y\$, \$ρ_S(X, Y)\$:

• the values of \$X\$ and \$Y\$ are ranked in an ascending order respectively. The smallest value will be given a rank of 1, the next smaller value will be given a rank of 2, etc. If there are multiple entries with the same value, then each entry is ranked equal to the average position.
• the rank difference for the i-th pair of \$x_i\$ and \$y_i\$ is then calculated and, $$ \rho_S(X, Y) = 1 - \frac{6\sum_{i=1}^{N} d_i^2}{N(N^2-1)} $$ where \$N\$ is the total pairs of data points and \$d_i\$ is the rank difference for the i-th pair of data.

Create a function calculate_spearman_rank_correlation(X, Y) that returns you the value of the rank correlation, given two sets of data X and Y.

My solution:

def calculate_spearman_rank_correlation(X,Y):
    '''
    #The Spearman rank correlation is used to evaluate if the relationship between two variables, 
    X and Y is monotonic. The rank correlation measures how closely related the ordering of one
    variable to the other variable, with no regard to the actual values of the variables. 
    '''
    if len(X) != len(Y): #Making sure length is the same
        return("Please check your dataset again for input errors.") 
    else:
        # for X
        xranks = [sorted(X).index(x) + 1 for x in X]
        dupesx = [x for n, x in enumerate(xranks) if x in xranks[:n]]
        if len(dupesx) == 0: 
            print(xranks)
        else:
            
            dupenumberx = dupesx[0]
            count_ax = xranks.count(dupesx[0])
        
            sumx = 0
            tx = 0
            while tx < count_ax:
                sumx = sumx + dupenumberx + 1*tx
                tx=tx+1
                sumx2 = sumx/count_ax
                for i in range(len(xranks)):
                    if xranks[i] == dupenumberx:
                        xranks[i] = sumx2
                        
        # for Y
        yranks = [sorted(Y).index(y) + 1 for y in Y]
        dupesy = [y for n, y in enumerate(yranks) if y in yranks[:n]]
        if len(dupesy) == 0: 
            print(yranks)
        else:
            
            dupenumbery = dupesy[0]
            count_ay = yranks.count(dupesy[0])
        
            sumy = 0
            ty = 0
            while ty < count_ay:
                sumy = sumy + dupenumbery + 1*ty
                ty=ty+1
                sumy2 = sumy/count_ay
                for i in range(len(yranks)):
                    if yranks[i] == dupenumbery:
                        yranks[i] = sumy2
        

        # calculate rank difference for the i-th pair of data
        sub_list = []
        zip_object = zip(xranks, yranks)
        for xranks_i, yranks_i in zip_object:
            sub_list.append(xranks_i - yranks_i)
            
        # inplementing the formula for SRC
        sub_list2 = [i **2 for i in sub_list]
        di = sum(sub_list2)
        N = len(X)
        SRC = (1 - (6 * (di))/(N*((N**2)-1)))
        SRC_1dp = round(SRC, 3)
    return "The Spearman rank correlation is", SRC_1dp

The test sets:

# test case a
Xa = [10, 20, 30, 40, 50]
Ya = [2, 4, 6, 8, 10]
calculate_spearman_rank_correlation(Xa,Ya)

Output:

[1, 2, 3, 4, 5]
[1, 2, 3, 4, 5]
('The Spearman rank correlation is', 1.0)

---

test:

# test case b
Xb= [10, 20, 20, 50, 40]
Yb= [2, 10, 6, 4, 8]
calculate_spearman_rank_correlation(Xb,Yb)

output:

[1, 5, 3, 2, 4]
('The Spearman rank correlation is', -0.45)

---

test:

# test case d
Xd= [50, 50, 50, 30, 30, 30]
Yd= [2, 10, 8, 6, 10, 10]
calculate_spearman_rank_correlation(Xd,Yd)

output:

 ('The Spearman rank correlation is', 0.883)

Answer 1

@RichardNeumann correctly points out that your implementation does not adhere to the specification you were given. You've written test cases - good! But the missing step is that you haven't validated your results against a known implementation of the algorithm, such as that at Social Science Statistics.

Do not return error message strings; raise an exception.

You should factor out a rank calculation function rather than copying and pasting for X and Y.

Aside from the rank calculation, you don't have to do any of the statistics yourself. You say you don't want to use any "packages": fine; the built-in module statistics is not a package.

Write unit tests via assert.

Suggested

from collections import Counter
from math import isclose
from statistics import stdev, covariance
from typing import Sequence, Iterator


def rank_series(values: Sequence[float]) -> Iterator[float]:
    rank_by_value = {}
    index = 0
    for x, n in sorted(Counter(values).items()):
        rank_by_value[x] = index + (1 + n)/2
        index += n
    for x in values:
        yield rank_by_value[x]


def calculate_spearman(x: Sequence[float], y: Sequence[float]) -> float:
    """
    The Spearman rank correlation is used to evaluate if the relationship between two variables,
    X and Y is monotonic. The rank correlation measures how closely related the ordering of one
    variable to the other variable, with no regard to the actual values of the variables.
    """
    if len(x) != len(y):
        raise ValueError(f'X length {len(x)} does not match Y length {len(y)}')

    x_ranks = tuple(rank_series(x))
    y_ranks = tuple(rank_series(y))
    covxy = covariance(x_ranks, y_ranks)
    stdx = stdev(x_ranks)
    stdy = stdev(y_ranks)
    r = covxy / stdx / stdy
    return r


def test(x: Sequence[float], y: Sequence[float], expected: float) -> None:
    r = calculate_spearman(x, y)
    assert isclose(r, expected, abs_tol=0, rel_tol=1e-6)


def test_all() -> None:
    test(
        (10, 20, 20, 50, 40),
        (2, 10, 6, 4, 8),
        0.20519567041703082,
    )
    test(
        (10, 20, 30, 40, 50),
        (2, 4, 6, 8, 10),
        1,
    )
    test(
        (50, 50, 50, 30, 30, 30),
        (2, 10, 8, 6, 10, 10),
        -0.3110855084191276,
    )


if __name__ == '__main__':
    test_all()