Rank Korrelation mit Gewichten für Frequenzen, in Python

Question

Meine Daten sind ein Satz von n beobachtet Paaren zusammen mit ihren Frequenzen, das heißt zu jedem Paar (x _i, y _i) es einige entspricht k _i, die Anzahl der Male (x i , y i ) beobachtet wurde. Im Idealfall würde Ich mag sowohl Kendall-Tau und Spearman-rho für den Satz aller Kopien dieser Paare berechnen, die von k 1 + k 2 + ... + k _n Paaren besteht. Das Problem ist, dass k 1 + k 2 + ... + k _n, die Gesamtzahl der Beobachtungen, ist sehr groß und eine solche Datenstruktur wird nicht in den Speicher passen.

Natürlich dachte ich über die Häufigkeit der Zuweisung i -ten Paars, k i/(k + k 2 + ... + k _n), als Gewicht und Computing Rangkorrelation für die gewichtete Menge —, aber ich konnte keine Werkzeuge dafür finden. In den gewichteten Varianten der Rangkorrelation, die ich getroffen habe (z. B. scipy.stats.weightedtau), repräsentieren die Gewichte die Wichtigkeit von Rängen und nicht von Paaren, was für meine Ursache nicht relevant ist. Pearsons r scheint genau die Gewichtungsoption zu haben, die ich brauche, aber es entspricht nicht meinem Zweck, da x und y nirgends linear verwandt sind. Ich habe mich gefragt, ob ich einen Begriff für eine verallgemeinerte Korrelation bei gewichteten Datenpunkten verpasse.

Die einzige Idee, die ich bisher bekommen habe ist k zu verringern, k , ..., k _n von einem gemeinsamen Faktor c, so dass gestaffelte Anzahl an Kopien von i -ten Paars ist [k i /c] (hier [.] die Rundungs ​​ Operator ist, wie wir ganzzahlige Anzahl von Kopien jedes Paares haben müssen). Durch die Wahl c so dass [k/C] + [k/C] + ... + [k n/c] Paare können in den Speicher passen, könnten wir dann berechnen die Korrelationskoeffizienten tau und rho für die resultierende Menge. Allerdings k _i und k _j kann durch viele Größenordnungen unterscheiden, so c kann ich einige k _{signifikant groß sein} und damit k _i/c Rundung kann Informationsverlust verursachen.

UPD: Man kann, wie unten zusammen mit p-Wert auf einem Datensatz mit spezifizierten Frequenzgewichten Spearman-rho berechnen:

def frequency_pearsonr(data, frequencies): 
    """ 
    Calculates Pearson's r between columns (variables), given the 
    frequencies of the rows (observations). 

    :param data: 2-D array with data 
    :param frequencies: 1-D array with frequencies 
    :return: 2-D array with pairwise correlations, 
     2-D array with pairwise p-values 
    """ 
    df = frequencies.sum() - 2 
    Sigma = np.cov(data.T, fweights=frequencies) 
    sigma_diag = Sigma.diagonal() 
    Sigma_diag_pairwise_products = np.multiply.outer(sigma_diag, sigma_diag) 
    # Calculate matrix with pairwise correlations. 
    R = Sigma/np.sqrt(Sigma_diag_pairwise_products) 
    # Calculate matrix with pairwise t-statistics. Main diagonal should 
    # get 1/0 = inf. 
    with np.errstate(divide='ignore'): 
     T = R/np.sqrt((1 - R * R)/df) 
    # Calculate matrix with pairwise p-values. 
    P = 2 * stats.t.sf(np.abs(T), df) 

    return R, P 


def frequency_rank(data, frequencies): 
    """ 
    Ranks 1-D data array, given the frequency of each value. Same 
    values get same "averaged" ranks. Array with ranks is shaped to 
    match the input data array. 

    :param data: 1-D array with data 
    :param frequencies: 1-D array with frequencies 
    :return: 1-D array with ranks 
    """ 
    s = 0 
    ranks = np.empty_like(data) 
    # Compute rank for each unique value. 
    for value in sorted(set(data)): 
     index_grid = np.ix_(data == value) 
     # Find total frequency of the value. 
     frequency = frequencies[index_grid].sum() 
     ranks[index_grid] = s + 0.5 * (frequency + 1) 
     s += frequency  

    return ranks 


def frequency_spearmanrho(data, frequencies): 
    """ 
    Calculates Spearman's rho between columns (variables), given the 
    frequencies of the rows (observations). 

    :param data: 2-D array with data 
    :param frequencies: 1-D array with frequencies 
    :return: 2-D array with pairwise correlations, 
     2-D array with pairwise p-values 
    """ 
    # Rank the columns. 
    ranks = np.empty_like(data) 
    for i, data_column in enumerate(data.T): 
     ranks[:, i] = frequency_rank(data_column, frequencies) 
    # Compute Pearson's r correlation and p-values on the ranks. 
    return frequency_pearsonr(ranks, frequencies) 


# Columns are variables and rows are observations, whose frequencies 
# are specified. 
data_col1 = np.array([1, 0, 1, 0, 1]) 
data_col2 = np.array([.67, .25, .75, .2, .6]) 
data_col3 = np.array([.1, .3, .8, .3, .2]) 
data = np.array([data_col1, data_col2, data_col3]).T 
frequencies = np.array([2, 4, 1, 3, 2]) 

# Same data, but with observations (rows) actually repeated instead of 
# their frequencies being specified. 
expanded_data_col1 = np.array([1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1]) 
expanded_data_col2 = np.array([.67, .67, .25, .25, .25, .25, .75, .2, .2, .2, .6, .6]) 
expanded_data_col3 = np.array([.1, .1, .3, .3, .3, .3, .8, .3, .3, .3, .2, .2]) 
expanded_data = np.array([expanded_data_col1, expanded_data_col2, expanded_data_col3]).T 

# Compute Spearman's rho for data in both formats, and compare. 
frequency_Rho, frequency_P = frequency_spearmanrho(data, frequencies) 
Rho, P = stats.spearmanr(expanded_data) 
print(frequency_Rho - Rho) 
print(frequency_P - P) 

das spezielle Beispiel oben zeigt, dass beide Methoden gleich Korrelationen und gleiche p-Werte erzeugen:

[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00] 
[ 1.11022302e-16 0.00000000e+00 -5.55111512e-17] 
[ 0.00000000e+00 -5.55111512e-17 0.00000000e+00]] 
[[ 0.00000000e+00 -1.35525272e-19 4.16333634e-17] 
[ -9.21571847e-19 0.00000000e+00 -5.55111512e-17] 
[ 4.16333634e-17 -5.55111512e-17 0.00000000e+00]] 

Answer 1

Der von Paul vorgeschlagene Ansatz zur Berechnung von Kendalls Tau funktioniert. Sie müssen die Indizes des sortierten Arrays nicht als Ränge zuweisen, Indizes von unsortierten funktionieren genauso gut (wie im Beispiel mit gewichtetem Tau gezeigt). Die Gewichte müssen auch nicht normalisiert werden.

Regular (ungewichtet) Kendall-Tau (auf dem "erweiterte" Datensatz):

stats.kendalltau([0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1], 
       [.25, .25, .25, .25, .2, .2, .2, .667, .667, .75, .6, .6]) 
KendalltauResult(correlation=0.7977240352174656, pvalue=0.0034446936330652677) 

gewichteten tau Kendall (auf dem Datensatz mit dem Auftreten zählt als Gewicht):

stats.weightedtau([1, 0, 1, 0, 1], 
        [.667, .25, .75, .2, .6], 
        rank=False, 
        weigher=lambda r: [2, 4, 1, 3, 2][r], 
        additive=False) 
WeightedTauResult(correlation=0.7977240352174656, pvalue=nan) 

nun die Der p-Wert wird aufgrund der Spezifität der gewichteten Implementierung nicht berechnet. Wir könnten den p-Wert mit dem ursprünglich angebotenen Trick annähern, die Vorkommen zu reduzieren, aber ich würde andere Ansätze sehr schätzen. Entscheidungen über das Verhalten des Algorithmus auf der verfügbaren Speichermenge zu treffen, sieht für mich wie ein Schmerz aus.