Correlacion lineal con python

Correlación lineal con Python

Joaquín Amat Rodrigo
Marzo, 2020

Más sobre ciencia de datos: cienciadedatos.net

Introducción


La correlación lineal es un método estadístico que permite cuantificar la relación lineal existente entre dos variables. Existen varios estadísticos, llamados coeficientes de correlación lineal, desarrollados con el objetivo de medir este tipo de asociación, algunos de los más empleados son Pearson, Spearman y Kendall.

Con frecuencia, los estudios de correlación lineal preceden a análisis más complejos, como la creación de modelos de regresión. Primero, se analiza si las variables están correlacionadas y, en caso de estarlo, se procede a generar modelos.

Es importante destacar que, la existencia de correlación entre dos variables, no implica necesariamente causalidad. La asociación observada puede deberse a un tercer factor (confounder).

Correlación lineal con Python

Existen múltiples implementaciones que permiten calcular correlaciones lineales en Python, tres de las más utilizadas están disponibles en las librerías: SciPy, Pandas y Pingouin. A lo largo de este documento se muestran ejemplos de cómo utilizarlas.

Coeficientes de correlación lineal


Covarianza


Para estudiar la relación lineal existente entre dos variables continuas es necesario disponer de parámetros que permitan cuantificar dicha relación. Uno de estos parámetros es la covarianza, que mide el grado de variación conjunta de dos variables aleatorias.

$$\text{Covarianza muestral} = Cov(X,Y)=\dfrac {\sum _{i=1}^{n}\left( x_{i}-\overline {x}\right) \left( y_{i}-\overline {y}\right) }{N-1}$$

donde $\overline{x}$ e $\overline{y}$ son la media de cada variable, y $x_i$ e $y_i$ son el valor de las variables para la observación $i$.

Valores positivos indican que las dos variables cambian en la misma dirección y, valores negativos, que lo hacen en direcciones opuestas.

La principal limitación de la covarianza es que, su magnitud, depende de las escalas en que se miden las variables estudiadas. Esto implica que no puede utilizarse para comparar el grado de asociación entre pares de variables medidas en distintas escalas. Una forma de evitar esta limitación y poder hacer comparaciones consiste en estandarizar la covarianza, generando lo que se conoce como coeficientes de correlación.

Coeficientes de correlación lineal


Los coeficientes de correlación lineal son estadísticos que cuantifican la asociación lineal entre dos variables numéricas. Existen diferentes tipos, de entre los que destacan el Pearson, Rho de Spearman y Tau de Kendall. Todos ellos comparten que:

  • Su valor está comprendido en el rango [+1 , -1]. Siendo +1 una correlación positiva perfecta y -1 una correlación negativa perfecta.

  • Se emplean como medida de la fuerza de asociación entre dos variables (tamaño del efecto):

    • 0: asociación nula.

    • 0.1: asociación pequeña.

    • 0.3: asociación mediana.

    • 0.5: asociación moderada.

    • 0.7: asociación alta.

    • 0.9: asociación muy alta.

Desde el punto de vista práctico, las principales diferencias entre estos tres coeficientes son:

  • La correlación de Pearson funciona bien con variables cuantitativas que tienen una distribución normal o próxima a la normal. Es más sensible a los valores extremos que las otras dos alternativas.

  • La correlación de Spearman se emplea con variables cuantitativas (continuas o discretas). En lugar de utilizar directamente el valor de cada variable, los datos son ordenados y reemplazados por su respectivo orden ranking. Es un método no paramétrico muy utilizado cuando no se satisface la condición de normalidad necesaria para aplicar la correlación de Pearson.

  • La correlación de Kendall es otra alternativa no paramétrica que, al igual que la correlación de Spearman, utiliza la ordenación de las observaciones ranking. Es recomendable cuando se dispone de pocos datos y muchos de ellos ocupan la misma posición en el rango, es decir, cuando hay muchas ligaduras.

Significancia estadística


Además del valor obtenido para el coeficiente de correlación, es necesario calcular su significancia estadística. Por muy cercano que sea el valor del coeficiente de correlación a $+1$ o $-1$, si no es significativo, no se dispone de evidencias suficiente para afirmar que existe una correlación real, ya que el valor observado podría deberse a simple aleatoriedad.

El test paramétrico de significancia estadística empleado para el coeficiente de correlación es el t-test. Donde el estadístico t se obtiene acorde a la ecuación:

$$t=\dfrac {r\sqrt{N-2}}{\sqrt{1-r^{2}}}$$

$r$ es el calor del coeficiente de correlación y N es el número de observaciones disponibles. Los grados de libertad se calculan como $df= N-2$.

En este test, se considera como hipótesis nula ($H_0$) que las variables son independientes (coeficiente de correlación poblacional = 0), y como hipótesis alternativa ($H_a$), que sí existe relación (coeficiente de correlación poblacional $\neq$ 0)

También se puede calcular la significancia de un coeficiente de correlación recurriendo a métodos no paramétricos como el bootstrapping.

Tamaño de efecto


La correlación lineal, además del valor del coeficiente de correlación y de sus significancia, también tiene un tamaño de efecto asociado conocido como coeficiente de determinación $R^2$.

$R^2$ se interpreta como la cantidad de varianza de $Y$ explicada por $X$. En el caso del coeficiente de Pearson y el de Spearman, $R^2$ se obtiene elevando al cuadrado el coeficiente de correlación. En el caso de Kendall no se puede calcular de este modo.

Coeficiente de Pearson


El coeficiente de correlación de Pearson es la covarianza estandarizada.

$$\rho=\dfrac {Cov(X,Y)}{\sigma_{x}\sigma_{y}}$$

La anterior ecuación se corresponde con el coeficiente de Pearson poblacional ($\rho$). En la práctica, raramente se tiene acceso a toda la población, por lo que su valor se estima a partir de una muestra mediante el coeficiente de Pearson muestral ($r$).

$$r_{xy}=\dfrac {\sum _{i=1}^{n}\left( x_{i}-\overline {x}\right) \left( y_{i}-\overline {y}\right) } {\sqrt {\sum _{i=1}^{n}\left( x_{i}-\overline {x}\right) ^{2}\sum _{i=1}^{n}}\overline {\left( y_{i}-\overline {y}\right) ^{2}}}$$

Condiciones

Las condiciones que se deben de cumplir para que el coeficiente de correlación de Pearson sea válido son:

  • La relación que se quiere estudiar es de tipo lineal (de lo contrario, el coeficiente de Pearson no la puede detectar).

  • Las dos variables deben de ser numéricas.

  • Normalidad: ambas variables se tienen que distribuir de forma normal. En la práctica, se suele considerar válido aun cuando se alejan moderadamente de la normalidad.

  • Homocedasticidad: la varianza de $Y$ debe ser constante a lo largo de la variable $X$. Esto se puede contrastar si en un scatterplot los valores de $Y$ mantienen la misma dispersión en las distintas zonas de la variable $X$.

Características

  • Toma valores entre [-1, +1], siendo +1 una correlación lineal positiva perfecta y -1 una correlación lineal negativa perfecta.

  • Es independiente de las escalas en las que se midan las variables.

  • No varía si se aplican transformaciones a las variables.

  • No tiene en consideración que las variables sean dependientes o independientes.

  • El coeficiente de correlación de Pearson no equivale a la pendiente de la recta de regresión.

  • Es sensible a outliers, por lo que se recomienda, en caso de poder justificarlos, excluirlos antes de realizar el cálculo.

Coeficiente de Spearman (Spearman's rho)


El coeficiente de Spearman es el equivalente al coeficiente de Pearson pero con una previa transformación de los datos a ranking. Esto significa que, en lugar de utilizar directamente el valor de cada variable, los datos son ordenados y reemplazados por su respectivo orden (primero, segundo, tercero...). Se emplea como alternativa no paramétrica al coeficiente de Pearson cuando los valores son ordinales, o bien, cuando los valores son continuos pero no satisfacen la condición de normalidad.

$$r_{s}=1-\frac{6\sum d_{i}^{2}}{n (n^{2}-1)},$$

Siendo $d_{i}$ la distancia entre los rangos de cada observación ($x_{i} - y_{i}$) y n el número de observaciones.

El coeficiente de Spearman requiere que la relación entre las variables sea monótona, es decir, que cuando una variable crece la otra también lo hace o cuando una crece la otra decrece (que la tendencia sea constante).

Al trabajar con el orden de las observaciones en lugar de su valor real, tiene la característica de ser menos sensible que Pearson a valores extremos.

Coeficiente Tau de Kendall


La correlación de Kendall es un método no paramétrico que, al igual que la correlación de Spearman, utiliza la ordenación de las observaciones ranking.

Es otra alternativa al coeficiente de correlación de Pearson cuando no se cumple la condición de normalidad. Suele utilizarse en lugar del coeficiente de Spearman cuando el número de observaciones es pequeño o los valores se acumulan en una región por lo que el número de ligaduras al generar el ranking es alto.

$$\tau= \frac{C-D}{\frac{1}{2}n(n-1)}$$

siendo $C$ el número de pares concordantes, aquellos en los que el rango de la segunda variable es mayor que el rango de la primera variable. $D$ el número de pares discordantes, cuando el rango de la segunda es igual o menor que el rango de la primera variable.

Jackknife correlation


El coeficiente de correlación de Pearson resulta efectivo en ámbitos muy diversos. Sin embargo, tiene la desventaja de no ser robusto frente a outliers a pesar de que se cumpla la condición de normalidad. Si dos variables tienen un pico o un valle común en una única observación, por ejemplo por un error de lectura, la correlación va a estar dominada por este registro a pesar de que entre las dos variables no haya correlación real alguna. Lo mismo puede ocurrir en la dirección opuesta. Si dos variables están altamente correlacionadas excepto para una observación en la que los valores son muy dispares, entonces la correlación existente quedará enmascarada. Una forma de evitarlo es recurrir a la Jackknife correlation, que consiste en calcular todos los posibles coeficientes de correlación entre dos variables si se excluye cada vez una de las observaciones. El promedio de todas las Jackknife correlations calculadas atenúa en cierta medida el efecto del outlier.

$$\bar{\theta}_{(A,B)} = \text{Promedio Jackknife correlation (A,B)} = \frac{1}{n} \sum_{i=1}^n \hat r_i$$

donde $n$ es el número de observaciones y $\hat r_i$ es el coeficiente de correlación de Pearson estimado entre las variables $A$ y $B$, habiendo excluido la observación $i$.

Además del promedio, se puede estimar su error estándar ($SE$) y así obtener intervalos de confianza para la Jackknife correlation y su correspondiente p-value.

$$SE = \sqrt{\frac{n-1}{n} \sum_{i=1}^n (\hat r_i - \bar{\theta})^2}$$

Intervalo de confianza del 95% $(Z=1.96)$

$$\text{Promedio Jackknife correlation (A,B)} \pm \ 1.96 *SE$$$$\bar{\theta} - 1.96 \sqrt{\frac{n-1}{n} \sum_{i=1}^n (\hat r_i - \bar{\theta})^2} \ , \ \ \bar{\theta}+ 1.96 \sqrt{\frac{n-1}{n} \sum_{i=1}^n (\hat r_i - \bar{\theta})^2}$$

P-value para la hipótesis nula de que $\bar{\theta} = 0$:

$$Z_{calculada} = \frac{\bar{\theta} - H_0}{SE}= \frac{\bar{\theta} - 0}{\sqrt{\frac{n-1}{n} \sum_{i=1}^n (\hat r_i - \bar{\theta})^2}}$$

$$p_{value} = P(Z > Z_{calculada})$$

Cuando se emplea este método es conveniente calcular la diferencia entre el valor de correlación obtenido por Jackknife correlation $(\bar{\theta})$ y el que se obtiene si se emplean todas las observaciones $(\bar{r})$. A esta diferencia se le conoce como Bias. Su magnitud es un indicativo de cuanto está influenciada la estimación de la correlación entre dos variables debido a un valor atípico u outlier.

$$Bias = (n-1)*(\bar{\theta} - \hat{r})$$

Si se calcula la diferencia entre cada correlación $(\hat r_i)$ estimada en el proceso de Jackknife y el valor de correlación $(\hat r)$ obtenido si se emplean todas las observaciones, se puede identificar que observaciones son más influyentes.

Cuando el estudio requiere minimizar al máximo la presencia de falsos positivos, a pesar de que se incremente la de falsos negativos, se puede seleccionar como valor de correlación el menor de entre todos los calculados en el proceso de Jackknife.

$$\text{Correlacion}= min \{ \hat r_1, \hat r_2,..., \hat r_n \}$$

A pesar de que el método de Jackknife permite aumentar la robustez de la correlación de Pearson, si los outliers son muy extremos su influencia seguirá siendo notable. Siempre es conveniente una representación gráfica de los datos para poder identificar si hay valores atípicos y eliminarlos. Otras alternativas robustas son la correlación de Spearman o el método de Bootstrapping.

Correlación parcial


Como se ha explicado, la correlación estudia la relación (lineal o monotónica) existente entre dos variables. Puede ocurrir que, la relación que muestran dos variables, se deba a una tercera variable que influye sobre las otras dos. A este fenómeno se le conoce como confounding. Por ejemplo, si se correlaciona el tamaño del pie de una persona con su inteligencia, se encuentra una correlación positiva alta. Sin embargo, dicha relación se debe a una tercera variable que está relacionada con las otras dos, la edad. La correlación parcial permite estudiar la relación lineal entre dos variables bloqueando el efecto de una tercera (o más) variables. Si el valor de correlación de dos variables es distinto al valor de correlación parcial de esas mismas dos variables, cuando se controla una tercera, significa que la tercera variable influye en las otras dos. La función en partial_corr() del paquete pingouin permite estudiar correlaciones parciales.

Ejemplo correlación lineal


Un estudio pretende analizar si existe una correlación lineal positiva entre la altura y el peso de las personas. Los datos utilizados en este ejemplo se han obtenido del libro Statistical Rethinking by Richard McElreath. El set de datos contiene información recogida por Nancy Howell a finales de la década de 1960 sobre el pueblo !Kung San, que viven en el desierto de Kalahari entre Botsuana, Namibia y Angola.

Librerías


Las librerías utilizadas en este documento son:

In [30]:
# Tratamiento de datos
# ==============================================================================
import pandas as pd
import numpy as np
from sklearn.datasets import load_diabetes

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from matplotlib import style
import seaborn as sns

# Preprocesado y análisis
# ==============================================================================
import statsmodels.api as sm
import pingouin as pg
from scipy import stats
from scipy.stats import pearsonr

# Configuración matplotlib
# ==============================================================================
plt.style.use('ggplot')

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')

Datos

In [31]:
url = ('https://raw.githubusercontent.com/JoaquinAmatRodrigo/' +
       'Estadistica-machine-learning-python/master/data/Howell1.csv')
datos = pd.read_csv(url)

# Se utilizan únicamente información de individuos mayores de 18 años.
datos = datos[datos.age > 18]

datos.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 346 entries, 0 to 543
Data columns (total 4 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   height  346 non-null    float64
 1   weight  346 non-null    float64
 2   age     346 non-null    float64
 3   male    346 non-null    int64  
dtypes: float64(3), int64(1)
memory usage: 13.5 KB

Análisis gráfico


En primer lugar se representan las dos variables mediante un diagrama de dispersión (scatterplot) para intuir si existe relación lineal o monotónica. Si no la hay, no tiene sentido calcular este tipo de correlaciones.

In [32]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(1, 1, figsize=(6,4))
ax.scatter(x=datos.height, y=datos.weight, alpha= 0.8)
ax.set_xlabel('Altura')
ax.set_ylabel('Peso');

El diagrama de dispersión parece indicar una relación lineal positiva entre ambas variables.

Para poder elegir el coeficiente de correlación adecuado, se tiene que analizar el tipo de variables y la distribución que presentan. En este caso, ambas variables son cuantitativas continuas y pueden ordenarse para convertirlas en un ranking, por lo que, a priori, los tres coeficientes podrían aplicarse. La elección se hará en función de la distribución que presenten las observaciones: normalidad, homocedasticidad y presencia de outliers.

Normalidad

In [33]:
# Gráfico distribución variables
# ==============================================================================
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))

axs[0].hist(x=datos.height, bins=20, color="#3182bd", alpha=0.5)
axs[0].plot(datos.height, np.full_like(datos.height, -0.01), '|k', markeredgewidth=1)
axs[0].set_title('Distribución altura (height)')
axs[0].set_xlabel('height')
axs[0].set_ylabel('counts')

axs[1].hist(x=datos.weight, bins=20, color="#3182bd", alpha=0.5)
axs[1].plot(datos.weight, np.full_like(datos.weight, -0.01), '|k', markeredgewidth=1)
axs[1].set_title('Distribución peso (weight)')
axs[1].set_xlabel('weight')
axs[1].set_ylabel('counts')


plt.tight_layout();
In [34]:
# Gráfico Q-Q
# ==============================================================================
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))

sm.qqplot(
    datos.height,
    fit   = True,
    line  = 'q',
    alpha = 0.4,
    lw    = 2,
    ax    = axs[0]
)
axs[0].set_title('Gráfico Q-Q height', fontsize = 10, fontweight = "bold")
axs[0].tick_params(labelsize = 7)

sm.qqplot(
    datos.height,
    fit   = True,
    line  = 'q',
    alpha = 0.4,
    lw    = 2,
    ax    = axs[1]
)
axs[1].set_title('Gráfico Q-Q height', fontsize = 10, fontweight = "bold")
axs[1].tick_params(labelsize = 7)

Además del estudio gráfico, se recurre a dos test estadísticos que contrasten la normalidad de los datos: Shapiro-Wilk test y D'Agostino's K-squared test. Este último es el que incluye el summary de statsmodels bajo el nombre de Omnibus.

En ambos test, la hipótesis nula considera que los datos siguen una distribución normal, por lo tanto, si el p-value no es inferior al nivel de referencia alpha seleccionado, no hay evidencias para descartar que los datos se distribuyen de forma normal.

In [35]:
# Normalidad de los residuos Shapiro-Wilk test
# ==============================================================================
shapiro_test = stats.shapiro(datos.height)
print(f"Variable height: {shapiro_test}")
shapiro_test = stats.shapiro(datos.weight)
print(f"Variable weight: {shapiro_test}")
Variable height: ShapiroResult(statistic=0.9910696148872375, pvalue=0.03439393267035484)
Variable weight: ShapiroResult(statistic=0.9911819696426392, pvalue=0.03673496097326279)
In [36]:
# Normalidad de los residuos D'Agostino's K-squared test
# ==============================================================================
k2, p_value = stats.normaltest(datos.height)
print(f"Variable height: Estadítico = {k2}, p-value = {p_value}")
k2, p_value = stats.normaltest(datos.weight)
print(f"Variable weight: Estadítico = {k2}, p-value = {p_value}")
Variable height: Estadítico = 7.210790495766356, p-value = 0.02717670115638557
Variable weight: Estadítico = 8.402628478646022, p-value = 0.014975881988445145

El análisis gráfico y los test estadísticos muestran evidencias de que no se puede asumir normalidad en ninguna de las dos variables. Siendo estrictos, este hecho excluye la posibilidad de utilizar el coeficiente de Pearson, dejando como alternativas el de Spearman o Kendall. Sin embargo, dado que la distribución no se aleja mucho de la normalidad y de que el coeficiente de Pearson tiene cierta robustez, a fines prácticos sí que se podría utilizar siempre y cuando se tenga en cuenta este hecho y se comunique en los resultados. Otra posibilidad es tratar de transformar las variables para mejorar su distribución, por ejemplo, aplicando el logaritmo.

In [37]:
# Transformación logarítmica de los datos
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 4))

sm.qqplot(
    np.log(datos.height),
    fit   = True,
    line  = 'q',
    alpha = 0.4,
    lw    = 2,
    ax    = ax
)
ax.set_title('Gráfico Q-Q log(height)', fontsize = 13)
ax.tick_params(labelsize = 7)


shapiro_test = stats.shapiro(np.log(datos.height))
print(f"Variable height: {shapiro_test}")
Variable height: ShapiroResult(statistic=0.9922674894332886, pvalue=0.06951289623975754)

La trasformación logarítmica de la variable altura (height) consigue una distribución más próxima a la normal.

Homocedasticidad


La homocedasticidad implica que la varianza se mantiene constante. Puede analizarse de forma gráfica representando las observaciones en un diagrama de dispersión y viendo si mantiene una homogeneidad en su dispersión a lo largo del eje X. Una forma cónica es un claro indicativo de falta de homocedasticidad. Dos test estadísticos utilizados para contrastar la homocedasticidad son: test de Goldfeld-Quandt y el de Breusch-Pagan.

Tal como muestra el diagrama de dispersión generado al inicio del ejercicio, no se aprecia ningún patrón cónico y la dispersión es constante.

Coeficientes correlación


Debido a la falta de normalidad, los resultados generados por Pearson no son del todo precisos. Sin embargo, dado que la desviación de la normalidad es leve y no se aprecian outliers, con fines ilustrativos, se procede a calcular los tres tipos de coeficientes.

De nuevo recordar que, cuando alguna de las condiciones asumidas por un modelo o test estadístico no se cumplen, no significa que obligatoriamente se tenga que descartar, pero hay que ser consciente de las implicaciones que tiene y reportarlo siempre en los resultados.

Pandas

Pandas permite calcular la correlación de dos Series (columnas de un DataFrame). El cálculo se hace por pares, eliminando automáticamente aquellos con valores NA/null. Una limitación de Pandas es que no calcula la significancia estadística.

In [38]:
# Cálculo de correlación con Pandas
# ==============================================================================
print('Correlación Pearson: ', datos['weight'].corr(datos['height'], method='pearson'))
print('Correlación spearman: ', datos['weight'].corr(datos['height'], method='spearman'))
print('Correlación kendall: ', datos['weight'].corr(datos['height'], method='kendall'))
Correlación Pearson:  0.7528177220327662
Correlación spearman:  0.7510966609219974
Correlación kendall:  0.5639709660523899

Scypy.stats

La implementación de Scypy.stats sí permite calcular la significancia estadística además del coeficiente de correlación. La función stats.pearsonr(), devuelve un error si alguna de las observaciones contienen valores NA/null. Las funciones stats.spearmanr() y stats.kendalltau() sí permiten excluirlos de forma automática si se indica nan_policy='omit'.

In [39]:
# Cálculo de correlación y significancia con Scipy
# ==============================================================================
r, p = stats.pearsonr(datos['weight'], datos['height'])
print(f"Correlación Pearson: r={r}, p-value={p}")

r, p = stats.spearmanr(datos['weight'], datos['height'])
print(f"Correlación Spearman: r={r}, p-value={p}")

r, p = stats.kendalltau(datos['weight'], datos['height'])
print(f"Correlación Pearson: r={r}, p-value={p}")
Correlación Pearson: r=0.752817722032767, p-value=1.8941037794176005e-64
Correlación Spearman: r=0.7510966609219974, p-value=5.2882247217804375e-64
Correlación Pearson: r=0.5639709660523899, p-value=3.162649137764771e-54

Pingouin

La librería Pingouin tiene una de las implementaciones más completas. Con la función corr() se obtiene, además del coeficiente de correlación, su significancia, intervalo de confianza y poder estadístico entre otros.

In [40]:
# Cálculo de correlación, significancia e intervalos con pingouin
# ==============================================================================
display(pg.corr(datos['weight'], datos['height'], method='pearson'))
display(pg.corr(datos['weight'], datos['height'], method='spearman'))
display(pg.corr(datos['weight'], datos['height'], method='kendall'))
n r CI95% r2 adj_r2 p-val BF10 power
pearson 346 0.752818 [0.7, 0.8] 0.566735 0.564208 1.894104e-64 8.84e+60 1.0
n r CI95% r2 adj_r2 p-val power
spearman 346 0.751097 [0.7, 0.79] 0.564146 0.561605 5.288225e-64 1.0
n r CI95% r2 adj_r2 p-val power
kendall 346 0.563971 [0.49, 0.63] 0.318063 0.314087 3.162649e-54 1.0

Conclusión


Los test estadísticos muestran una correlación lineal entre moderada y alta, con claras evidencias estadísticas de que la relación observada no se debe al azar ($pvalue \approx 0$).

Ejemplo Jackknife correlation


Un equipo de investigadores quiere estudiar si existe correlación en la presencia de dos sustancias (A y B) en el agua de los ríos. Para ello han realizado una serie de mediciones en las que se cuantifica la concentración de las dos sustancias en 10 muestras independientes de agua. Se sospecha que el instrumento de lectura sufre alguna avería que provoca que algunas lecturas se disparen, por esta razón se quiere emplear un método de correlación robusto. El objetivo de este ejemplo es ilustrar el método de Jackknife, por lo que se asume que se cumplen las condiciones para la correlación de Pearson.

Datos

In [41]:
# Datos simulados de dos variables A y B
a = np.array([12,9,6,7,2,5,4,0,1,8])
b = np.array([3,5,1,9,5,3,7,2,10,5])

# Se introduce un outlier
a[5] = 20
b[5] = 16
In [42]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(1, 1, figsize=(6,4))
ax.plot(a, label='A')
ax.plot(b, label='B')
ax.set_xlabel('ID muestra')
ax.set_ylabel('Concentración')
ax.set_title('Concentración sustancias A y B en las muestras')
ax.legend();

Correlación de Pearson


Se procede a calcular la correlación de Pearson con y sin el outlier.

In [43]:
# Correlación con outlier
r, p = stats.pearsonr(a, b)
print(f"Correlación Pearson con outlier: r={r}, p-value={p}")
Correlación Pearson con outlier: r=0.517373115168915, p-value=0.12563522982639538
In [44]:
# Correlación sin outlier
r, p = stats.pearsonr(np.delete(a, 5), np.delete(b, 5))
print(f"Correlación Pearson sin outlier: r={r}, p-value={p}")
Correlación Pearson sin outlier: r=-0.18420184544326057, p-value=0.6351961086690547

Se confirma que, La observación número 5, tiene una gran influencia en el resultado de la correlación, siendo de 0.52 si está presente y de -0.18 si se excluye.

Jackknife Pearson correlation

In [45]:
# Función Jackknife correlation
# ==============================================================================

def correlacion_jackknife(x, y):
    '''
    Esta función aplica el método de Jackknife para el cálculo del coeficiente
    de correlación de Pearson.
    
    
    Parameters
    ----------
    x : 1D np.ndarray, pd.Series 
        Variable X.
        
    y : 1D np.ndarray, pd.Series
        Variable y.     

    Returns 
    -------
    correlaciones: 1D np.ndarray
        Valor de correlación para cada iteración de Jackknife
    '''
    
    n = len(x)
    valores_jackknife = np.full(shape=n, fill_value=np.nan, dtype=float)
    
    for i in range(n):
        # Loop para excluir cada observación y calcular la correlación
        r = stats.pearsonr(np.delete(x, i), np.delete(y, i))[0]
        valores_jackknife[i] = r

    promedio_jackknife = np.nanmean(valores_jackknife)
    standar_error = np.sqrt(((n - 1) / n) * \
                    np.nansum((valores_jackknife - promedio_jackknife) ** 2))
    bias = (n - 1) * (promedio_jackknife - stats.pearsonr(x, y)[0])
    
    resultados = {
        'valores_jackknife' : valores_jackknife,
        'promedio'          : promedio_jackknife,
        'se'                : standar_error,
        'bias'              : bias
    }
    
    return resultados
In [46]:
correlacion = correlacion_jackknife(x=a, y=b)
print(f"Correlación jackknife: {correlacion['promedio']}")
print(f"Error estándar: {correlacion['se']}")
print(f"Error bias: {correlacion['bias']}")
print(f"Valores_jackknife: {correlacion['valores_jackknife']}")
Correlación jackknife: 0.4781856690553491
Error estándar: 0.6915097651492329
Error bias: -0.35268701502209354
Valores_jackknife: [ 0.64719522  0.53705998  0.54597653  0.52827609  0.51215595 -0.18420185
  0.53554935  0.44125573  0.69065085  0.52793883]

Para identificar si existe algún valor muy influyente, se puede graficar el cambio que produce en el coeficiente de regresión la exclusión de cada observación.

In [47]:
variacion_corr = correlacion['valores_jackknife'] - stats.pearsonr(a, b)[0]
fig, ax = plt.subplots(1, 1, figsize=(6,4))
ax.plot(variacion_corr)
ax.set_xlabel('ID muestra')
ax.set_ylabel('Variación correlación');

El método Jackknife correlation solo ha sido capaz de amortiguar una pequeña parte de la influencia del outlier, sin embargo, sí ha permitido identificar qué observación está afectando en mayor medida.

Ejemplo matriz de correlaciones


Cuando se dispone de múltiples variables numéricas, por ejemplo en problemas de modelado estadístico y machine learning, es conveniente estudiar el grado de correlación entre las variables disponibles.

Una forma de hacerlo es mediante matrices de correlación, en las que se muestra el coeficiente de correlación para cada par de variables.

In [48]:
# Datos
# ==============================================================================
url = ('https://raw.githubusercontent.com/JoaquinAmatRodrigo/' +
       'Estadistica-machine-learning-python/master/data/SaratogaHouses.csv')
datos = pd.read_csv(url, sep=",")

# Se renombran las columnas para que sean más descriptivas
datos.columns = ["precio", "metros_totales", "antiguedad", "precio_terreno",
                 "metros_habitables", "universitarios", "dormitorios", 
                 "chimenea", "banyos", "habitaciones", "calefaccion",
                 "consumo_calefacion", "desague", "vistas_lago",
                 "nueva_construccion", "aire_acondicionado"]
       
# Variables numéricas
datos = datos.select_dtypes(include=['float64', 'int'])
In [49]:
# Matriz de correlación
# ==============================================================================
corr_matrix = datos.corr(method='pearson')
corr_matrix
Out[49]:
precio metros_totales antiguedad precio_terreno metros_habitables universitarios dormitorios chimenea banyos habitaciones
precio 1.000000 0.158333 -0.188793 0.581266 0.712390 0.200119 0.400349 0.376786 0.597250 0.531170
metros_totales 0.158333 1.000000 -0.016352 0.059222 0.163450 -0.033148 0.113982 0.085226 0.084823 0.137604
antiguedad -0.188793 -0.016352 1.000000 -0.021818 -0.174242 -0.037785 0.027125 -0.172022 -0.361897 -0.082264
precio_terreno 0.581266 0.059222 -0.021818 1.000000 0.423441 0.228427 0.202449 0.211727 0.297498 0.298865
metros_habitables 0.712390 0.163450 -0.174242 0.423441 1.000000 0.209981 0.656196 0.473788 0.718564 0.733666
universitarios 0.200119 -0.033148 -0.037785 0.228427 0.209981 1.000000 0.162919 0.246626 0.179541 0.157068
dormitorios 0.400349 0.113982 0.027125 0.202449 0.656196 0.162919 1.000000 0.284475 0.458033 0.671863
chimenea 0.376786 0.085226 -0.172022 0.211727 0.473788 0.246626 0.284475 1.000000 0.436234 0.319894
banyos 0.597250 0.084823 -0.361897 0.297498 0.718564 0.179541 0.458033 0.436234 1.000000 0.517585
habitaciones 0.531170 0.137604 -0.082264 0.298865 0.733666 0.157068 0.671863 0.319894 0.517585 1.000000

Las matrices de correlación tienen el inconveniente de tener un tamaño notable cuando se dispone de muchas variables. Para facilitar la identificación de pares de variables con correlaciones altas, es conveniente convertirlas en formato de tabla larga (tidy).

In [50]:
def tidy_corr_matrix(corr_mat):
    '''
    Función para convertir una matriz de correlación de pandas en formato tidy.
    '''
    corr_mat = corr_mat.stack().reset_index()
    corr_mat.columns = ['variable_1','variable_2','r']
    corr_mat = corr_mat.loc[corr_mat['variable_1'] != corr_mat['variable_2'], :]
    corr_mat['abs_r'] = np.abs(corr_mat['r'])
    corr_mat = corr_mat.sort_values('abs_r', ascending=False)
    
    return(corr_mat)

tidy_corr_matrix(corr_matrix).head(10)
Out[50]:
variable_1 variable_2 r abs_r
94 habitaciones metros_habitables 0.733666 0.733666
49 metros_habitables habitaciones 0.733666 0.733666
84 banyos metros_habitables 0.718564 0.718564
48 metros_habitables banyos 0.718564 0.718564
4 precio metros_habitables 0.712390 0.712390
40 metros_habitables precio 0.712390 0.712390
69 dormitorios habitaciones 0.671863 0.671863
96 habitaciones dormitorios 0.671863 0.671863
46 metros_habitables dormitorios 0.656196 0.656196
64 dormitorios metros_habitables 0.656196 0.656196
In [51]:
# Heatmap matriz de correlaciones
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))

sns.heatmap(
    corr_matrix,
    annot     = True,
    cbar      = False,
    annot_kws = {"size": 8},
    vmin      = -1,
    vmax      = 1,
    center    = 0,
    cmap      = sns.diverging_palette(20, 220, n=200),
    square    = True,
    ax        = ax
)

ax.set_xticklabels(
    ax.get_xticklabels(),
    rotation = 45,
    horizontalalignment = 'right',
)

ax.tick_params(labelsize = 10)

Ejemplo correlación parcial


Se quiere estudiar la relación entre las variables precio y peso de los automóviles. Se sospecha que esta relación podría estar influenciada por la variable potencia del motor, ya que a mayor peso del vehículo se requiere mayor potencia y, a su vez, motores más potentes son más caros.

Datos

Para este ejemplo se utiliza el set de datos Cars93 disponible en el paquete de R MASS.

In [52]:
url = ('https://raw.githubusercontent.com/JoaquinAmatRodrigo/' +
       'Estadistica-machine-learning-python/master/data/Cars93.csv')
datos = pd.read_csv(url)
datos['log_Price'] = np.log(datos.Price)
In [53]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(1, 1, figsize=(6,4))
ax.scatter(x=datos.Weight, y=datos.log_Price, alpha= 0.8)
ax.set_xlabel('Peso')
ax.set_ylabel('Log(Precio)');

El gráfico muestra una clara relación lineal entre el peso de un coche y el logaritmo de su precio.

Correlación

In [54]:
# Cálculo de correlación lineal
# ==============================================================================
pg.corr(x=datos['Weight'], y=datos['log_Price'], method='pearson')
Out[54]:
n r CI95% r2 adj_r2 p-val BF10 power
pearson 93 0.763544 [0.66, 0.84] 0.582999 0.573733 5.640674e-19 1.069e+16 1.0
In [55]:
# Cálculo de correlación lineal parcial
# ==============================================================================
pg.partial_corr(data=datos, x='Weight', y='log_Price', covar='Horsepower', method='pearson')
Out[55]:
n r CI95% r2 adj_r2 p-val BF10 power
pearson 93 0.404741 [0.22, 0.56] 0.163816 0.145234 0.000057 374.629 0.983505

Conclusión


La correlación entre el peso y el logaritmo del precio es alta (r=0.764) y significativa ($pvalue \approx 0$). Sin embargo, cuando se estudia su relación bloqueando la variable potencia de motor, a pesar de que la relación sigue siendo significativa, pasa a ser baja (r=0.4047).

Se puede afirmar que, la relación lineal existente entre el peso y el logaritmo del precio, está influenciada por el efecto de la variable potencia de motor.

Información de sesión

In [56]:
from sinfo import sinfo
sinfo()
-----
ipykernel   5.4.3
matplotlib  3.3.2
numpy       1.19.2
pandas      1.0.1
pingouin    0.3.10
scipy       1.6.0
seaborn     0.10.1
sinfo       0.3.1
sklearn     0.24.1
statsmodels 0.11.1
-----
IPython             7.20.0
jupyter_client      6.1.11
jupyter_core        4.7.1
notebook            6.2.0
-----
Python 3.7.9 (default, Aug 31 2020, 12:42:55) [GCC 7.3.0]
Linux-5.4.0-1038-aws-x86_64-with-debian-buster-sid
2 logical CPU cores, x86_64
-----
Session information updated at 2021-03-10 14:02

Bibliografía


Linear Models with R by Julian J.Faraway

OpenIntro Statistics: Fourth Edition by David Diez, Mine Çetinkaya-Rundel, Christopher Barr

Introduction to Machine Learning with Python: A Guide for Data Scientists

Points of Significance: Association, correlation and causation. Naomi Altman & Martin Krzywinski Nature Methods

¿Cómo citar este documento?

Correlacion lineal con Python por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/pystats05-correlacion-lineal-python.html DOI


¿Te ha gustado el artículo? Tu ayuda es importante

Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊


Creative Commons Licence
Este contenido, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.

Se permite:

  • Compartir: copiar y redistribuir el material en cualquier medio o formato.

  • Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

  • Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.

  • NoComercial: No puedes utilizar el material para fines comerciales.

  • CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.