Коэффициенты корреляции и значения p для всех пар строк матрицы

Сегодня я столкнулся с той же проблемой.

После получаса поиска в Google я не могу найти какой-либо код в библиотеке numpy / scipy, который может помочь мне в этом.

Итак, я написал свою собственную версию corrcoef

import numpy as np
from scipy.stats import pearsonr, betai

def corrcoef(matrix):
    r = np.corrcoef(matrix)
    rf = r[np.triu_indices(r.shape[0], 1)]
    df = matrix.shape[1] - 2
    ts = rf * rf * (df / (1 - rf * rf))
    pf = betai(0.5 * df, 0.5, df / (df + ts))
    p = np.zeros(shape=r.shape)
    p[np.triu_indices(p.shape[0], 1)] = pf
    p[np.tril_indices(p.shape[0], -1)] = p.T[np.tril_indices(p.shape[0], -1)]
    p[np.diag_indices(p.shape[0])] = np.ones(p.shape[0])
    return r, p

def corrcoef_loop(matrix):
    rows, cols = matrix.shape[0], matrix.shape[1]
    r = np.ones(shape=(rows, rows))
    p = np.ones(shape=(rows, rows))
    for i in range(rows):
        for j in range(i+1, rows):
            r_, p_ = pearsonr(matrix[i], matrix[j])
            r[i, j] = r[j, i] = r_
            p[i, j] = p[j, i] = p_
    return r, p

Первая версия использует результат np.corrcoef, а затем вычисляет p-значение на основе верхних треугольных значений матрицы corrcoef.

Вторая версия цикла просто перебирает строки и выполняет pearsonr вручную.

def test_corrcoef():
    a = np.array([
        [1, 2, 3, 4],
        [1, 3, 1, 4],
        [8, 3, 8, 5],
        [2, 3, 2, 1]])

    r1, p1 = corrcoef(a)
    r2, p2 = corrcoef_loop(a)

    assert np.allclose(r1, r2)
    assert np.allclose(p1, p2)

Тест пройден, они такие же.

def test_timing():
    import time
    a = np.random.randn(100, 2500)

    def timing(func, *args, **kwargs):
        t0 = time.time()
        loops = 10
        for _ in range(loops):
            func(*args, **kwargs)
        print('{} takes {} seconds loops={}'.format(
            func.__name__, time.time() - t0, loops))

    timing(corrcoef, a)
    timing(corrcoef_loop, a)


if __name__ == '__main__':
    test_corrcoef()
    test_timing()

Производительность на моем Macbook против матрицы 100x2500

corrcoef занимает 0,06608104705810547 секундных циклов = 10

corrcoef_loop занимает 7,585600137710571 секунд цикла = 10

Самым изощренным способом сделать это может быть метод сборки .corr в pandas, чтобы получить r:

In [79]:

import pandas as pd
m=np.random.random((6,6))
df=pd.DataFrame(m)
print df.corr()
          0         1         2         3         4         5
0  1.000000 -0.282780  0.455210 -0.377936 -0.850840  0.190545
1 -0.282780  1.000000 -0.747979 -0.461637  0.270770  0.008815
2  0.455210 -0.747979  1.000000 -0.137078 -0.683991  0.557390
3 -0.377936 -0.461637 -0.137078  1.000000  0.511070 -0.801614
4 -0.850840  0.270770 -0.683991  0.511070  1.000000 -0.499247
5  0.190545  0.008815  0.557390 -0.801614 -0.499247  1.000000

Чтобы получить значения p с помощью t-теста:

In [84]:

n=6
r=df.corr()
t=r*np.sqrt((n-2)/(1-r*r))

import scipy.stats as ss
ss.t.cdf(t, n-2)
Out[84]:
array([[ 1.        ,  0.2935682 ,  0.817826  ,  0.23004382,  0.01585695,
         0.64117917],
       [ 0.2935682 ,  1.        ,  0.04363408,  0.17836685,  0.69811422,
         0.50661121],
       [ 0.817826  ,  0.04363408,  1.        ,  0.39783538,  0.06700715,
         0.8747497 ],
       [ 0.23004382,  0.17836685,  0.39783538,  1.        ,  0.84993082,
         0.02756579],
       [ 0.01585695,  0.69811422,  0.06700715,  0.84993082,  1.        ,
         0.15667393],
       [ 0.64117917,  0.50661121,  0.8747497 ,  0.02756579,  0.15667393,
         1.        ]])
In [85]:

ss.pearsonr(m[:,0], m[:,1])
Out[85]:
(-0.28277983892175751, 0.58713640696703184)
In [86]:
#be careful about the difference of 1-tail test and 2-tail test:
0.58713640696703184/2
Out[86]:
0.2935682034835159 #the value in ss.t.cdf(t, n-2) [0,1] cell

Также вы можете просто использовать scipy.stats.pearsonr, который вы упомянули в OP:

In [95]:
#returns a list of tuples of (r, p, index1, index2)
import itertools
[ss.pearsonr(m[:,i],m[:,j])+(i, j) for i, j in itertools.product(range(n), range(n))]
Out[95]:
[(1.0, 0.0, 0, 0),
 (-0.28277983892175751, 0.58713640696703184, 0, 1),
 (0.45521036266021014, 0.36434799921123057, 0, 2),
 (-0.3779357902414715, 0.46008763115463419, 0, 3),
 (-0.85083961671703368, 0.031713908656676448, 0, 4),
 (0.19054495489542525, 0.71764166168348287, 0, 5),
 (-0.28277983892175751, 0.58713640696703184, 1, 0),
 (1.0, 0.0, 1, 1),
#etc, etc

Вроде хакерский и, возможно, неэффективный, но я думаю, что это может быть то, что вы ищете:

import scipy.spatial.distance as dist

import scipy.stats as ss

# Pearson's correlation coefficients
print dist.squareform(dist.pdist(data, lambda x, y: ss.pearsonr(x, y)[0]))    

# p-values
print dist.squareform(dist.pdist(data, lambda x, y: ss.pearsonr(x, y)[1]))

pdist Scipy - это очень полезная функция, которая в первую очередь предназначена для нахождения попарных расстояний между наблюдениями в n-мерном пространстве.

Но он позволяет определять вызываемые пользователем «метрики расстояния», которые можно использовать для выполнения любых парных операций. Результат возвращается в форме сжатой матрицы расстояний, которую можно легко преобразовать в форму квадратной матрицы, используя Функция Scipy 'squareform'.

Если вам не нужно использовать pearson, коэффициент корреляции, вы можете использовать коэффициент корреляции Спирмена, поскольку он возвращает как матрицу корреляции, так и p-значения (обратите внимание, что первое требует, чтобы ваши данные были нормально распределены, тогда как корреляция Спирмена является непараметрической мерой, поэтому не предполагает нормальное распределение ваших данных). Пример кода:

from scipy import stats
import numpy as np

data = np.array([[0, 1, -1], [0, -1, 1], [0, 1, -1]])
print 'np.corrcoef:', np.corrcoef(data)
cor, pval = stats.spearmanr(data.T)
print 'stats.spearmanr - cor:\n', cor
print 'stats.spearmanr - pval\n', pval

это точно такая же производительность, как и у corrcoef в MATLAB:

чтобы эта функция работала, вам нужно будет установить pandas, а также scipy.

# Compute correlation correfficients matrix and p-value matrix
# Similar function as corrcoef in MATLAB
# dframe: pandas dataframe
def corrcoef(dframe):

    fmatrix = dframe.values
    rows, cols = fmatrix.shape

    r = np.ones((cols, cols), dtype=float)
    p = np.ones((cols, cols), dtype=float)

    for i in range(cols):
        for j in range(cols):
            if i == j:
                r_, p_ = 1., 1.
            else:
                r_, p_ = pearsonr(fmatrix[:,i], fmatrix[:,j])

            r[j][i] = r_
            p[j][i] = p_

    return r, p

Вот минимальная версия ответа @CT Zhu. Нам не нужен pandas, так как корреляцию можно вычислить непосредственно из numpy, что должно быть быстрее, поскольку нам не нужен этап преобразования во фрейм данных.

import numpy as np
import scipy.stats as ss

def corr_significance_two_sided(cc, nData):
    # We will divide by 0 if correlation is exactly 1, but that is no problem
    # We would simply set the test statistic to be infinity if it evaluates to NAN
    with np.errstate(divide='ignore'):
        t = -np.abs(cc) * np.sqrt((nData - 2) / (1 - cc**2))
        t[t == np.nan] = np.inf
        return ss.t.cdf(t, nData - 2) * 2  # multiply by two to get two-sided p-value

x = np.random.uniform(0, 1, (8, 1000))
cc = np.corrcoef(x)
pVal = corr_significance_two_sided(cc, 1000)