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

Tengo una matriz de data con m filas yn columnas. np.corrcoef para calcular los coeficientes de correlación entre todos los pares de filas usando np.corrcoef :

 import numpy as np data = np.array([[0, 1, -1], [0, -1, 1]]) np.corrcoef(data) 

Ahora también me gustaría echar un vistazo a los valores p de estos coeficientes. np.corrcoef no proporciona estos; scipy.stats.pearsonr hace. Sin embargo, scipy.stats.pearsonr no acepta una matriz en la entrada.

¿Hay una forma rápida de calcular tanto el coeficiente como el valor p para todos los pares de filas (que llegan, por ejemplo, a matrices de dos por metros , una con coeficientes de correlación y otra con los valores de p correspondientes) sin tener que pasar manualmente? todas las parejas?

Me he encontrado con el mismo problema hoy.

Después de media hora de búsqueda 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)] = pf 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 de triángulo superior de la matriz de corrcoef.

La segunda versión de bucle solo iterando sobre filas, hacer pearsonr manualmente.

 def test_corrcoef(): a = np.array([ [1, 2, 3, 4], [1, 3, 1, 4], [8, 3, 8, 5]]) r1, p1 = corrcoef(a) r2, p2 = corrcoef_loop(a) assert np.allclose(r1, r2) assert np.allclose(p1, p2) 

La prueba pasada, son las mismas.

 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 en mi Macbook contra la matriz 100×2500

corrcoef toma 0,06608104705810547 segundos bucles = 10

corrcoef_loop toma 7.585600137710571 segundos bucles = 10

La forma más conveniente de hacerlo podría ser el método de .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 los 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:

 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 

Tipo de piratería y posiblemente ineficiente, pero creo que esto podría ser lo que estás buscando:

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

El pdist de Scipy es una función muy útil, que está diseñada principalmente para encontrar distancias por pares entre las observaciones en el espacio n-dimensional.

Pero permite que el usuario pueda definir “métricas de distancia”, que pueden ser explotadas para realizar cualquier tipo de operación por pares. El resultado se devuelve en 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 no tiene que usar el coeficiente de correlación de Pearson , puede usar el coeficiente de correlación de spearman , ya que devuelve la matriz de correlación y los valores p (tenga en cuenta que el primero requiere que sus datos se distribuyan normalmente, mientras que la correlación de spearman no es -Medida paramétrica, por lo tanto no asumiendo la distribución normal de sus datos). Un código de ejemplo:

 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