Saltar al contenido

Coeficientes de correlación y valores p para todos los pares de filas de una matriz

Siéntete en la libertad de divulgar nuestra página y códigos con otro, necesitamos de tu ayuda para aumentar esta comunidad.

Solución:

Hoy me he encontrado con el mismo problema.

Después de media hora de buscar en Google, no puedo encontrar ningún código en la biblioteca numpy / scipy que pueda ayudarme a hacer esto.

Así que escribí mi propia versión de 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

La primera versión usa el resultado de np.corrcoef y luego calcula el valor p basado en los valores superiores del triángulo de la matriz corrcoef.

La segunda versión de bucle simplemente iterando sobre filas, hace pearsonr manualmente.

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)

La prueba pasó, son los mismos.

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()

El rendimiento de mi Macbook frente a la matriz 100x2500

corrcoef toma 0.06608104705810547 segundos bucles = 10

corrcoef_loop tarda 7.585600137710571 segundos bucles = 10

La forma más conveniente de hacerlo podría ser el método de construcción. .corr en pandas, para obtener 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

Para obtener valores de p usando la prueba 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

También puedes usar el scipy.stats.pearsonr que mencionaste en OP:

En [95]: # devuelve una lista de tuplas de herramientas de iteración de importación (r, p, index1, index2)
[ss.pearsonr(m[:,i],metro[:,j]) + (i, j) para i, j en itertools.product (range (n), range (n))]Fuera[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

Sort of hackish and possibly inefficient, but I think this could be what you're looking for:

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])) # valores p imprimir dist.squareform (dist.pdist (datos, lambda x, y: ss.pearsonr (x, y)[1]))

La pdist de Scipy es una función muy útil, que está destinada principalmente a encontrar distancias por pares entre observaciones en un espacio n-dimensional.

Pero permite "métricas de distancia" definidas por el usuario que se pueden llamar, que pueden explotarse para llevar a cabo cualquier tipo de operación por pares. El resultado se devuelve en una forma de matriz de distancia condensada, que se puede cambiar fácilmente a la forma de matriz cuadrada utilizando la función 'forma cuadrada' de Scipy.

Si para ti ha resultado de utilidad este post, sería de mucha ayuda si lo compartes con otros seniors de esta forma contrubuyes a difundir este contenido.

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada.