Correlación Policórica
Es bastante común encontrar escalas ordinales en encuestas, tipo [“Bueno” – “Regular” – “Malo”]. Este tipo de variables son difíciles de modelar matemáticamente. Algunos convierten estas variables en un número, por ejemplo, en la escala mencionada arriba, se podría mapear esa lista a esta lista de números: [1, 0, -1]. Esto implica asumir que entre los puntos de la escala se tiene el mismo espacio.
El constructo que está midiendo la pregunta categórica ordinal puede ser abstracto o arbitrario. La correlación policórica ha sido usada principalmente en psicología, sociología o econometría, para constructos abstractos, como la preferencia o la satisfacción, la afinidad política, etc. Sin embargo, no necesariamente nos limitamos a estos casos, cualquier pareja de variables ordinales podría ser analizada de esta manera.
Si vamos a tratar estas variables como si fueran continuas, no podemos hacerlo descuidadamente. Pearson ya había pensado en esto hace más de un siglo y desarrollo la correlación tetracórica y luego se desarrolló la correlación policórica para tablas de contingencia con más de 2×2 celdas. El nombre viene de una serie infinita que utilizó Pearson para resolver este problema matemático. Más detalles en los artículos originales [3,4,5]. Si se asume que existe la misma distancia entre cada nivel de la variable ordinal, se corre el riesgo de obtener una distribución asimétrica o con mucha kurtosis, lo cual rompe nuestro requisito de tener una distribución normal al calcular una correlación.
Para desarrollar la correlación policórica de dos variables ordinales, se asume que cada variable tiene una variable continua latente con distribución normal y que las dos variables tienen una distribución normal bivariada con una correlación rho. Una estimación inicial de los límites entre cada categoría pueden hacerse utilizando su distribución empírica acumulada, y convirtiendolos a valores z (desviaciones standard normales) usando la función inversa de la función de distribución acumulada normal.
De acuerdo a Uebersax, el test de independencia X^2 es equivalente a probar que la correlación policórica es igual a 0. Así que no es necesario meternos a sacar la CP si sólo se necesita demostrar independencia. La correlación policórica, según veo, se usa para hacer PCA y análisis factorial de la manera correcta. También se me ocurre que la CP, al igual que los valores z asociados a cada nivel pueden servir para meter variables ordinales apropiadamente en redes neuronales.
Si no ha quedado claro qué significa que haya una distribución normal bivariada latente a dos variables ordinales, en la siguiente sección pongo una gráfica que visualiza esto elegantemente para una tabla de 2×2 celdas.
En general, en [1] se demuestra que la correlación de Pearson común aplicada en variables categóricas ordinales convertidas en números espaciados uniformemente tiene un sesgo hacia el cero, es decir que subestima la correlación, bajo la suposición que es correcto asumir que se tiene una dist. normal bivariada latente.
Estimación por ML de la correlación policórica entre dos variables categóricas ordinales¶
import numpy as np
%matplotlib inline
import matplotlib as mlp
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.base.model import GenericLikelihoodModel
El paquete de software estadístico para python statsmodels incluye una clase genérica para construir modelos de estimación por verosimilitud máxima (MLE). Basado en el trabajo de Olsson [2] he hecho una clase en python que calcula la correlación latente entre dos variables categóricas ordinales. En el artículo de Olsson se describe cuál es la Log-likelihood de los datos bajo el modelo descrito.
class PolychoricCorrMLE(GenericLikelihoodModel):
"""
Based on the paper Maximum Likelihood Estimation of the Polychoric Correlation Coefficient by Olsson U.
Tested it with nm and powell algorithms, obtaining good results with both.
"""
def __init__(self, endog, ** kwds):
"""
endog will be the contingency table.
exog not used
"""
if len(endog.shape) != 2 or endog.shape[0] < 2 or endog.shape[1] < 2:
raise "Matrix expected with more than 2 columns and more than 2 rows"
# Initialize thresholds based on inverse normal CDF on marginal probabilities
self.endog = endog
marginalColsF = np.sum(endog, axis= 0)
marginalRowsF = np.sum(endog, axis= 1)
marginalColsP = marginalColsF/np.sum(marginalColsF)
marginalRowsP = marginalRowsF/np.sum(marginalRowsF)
self.init_zCols = stats.norm.ppf(np.cumsum(np.hstack([[0.0],marginalColsP])) )
self.init_zRows = stats.norm.ppf(np.cumsum(np.hstack([[0.0],marginalRowsP])) )
super(PolychoricCorrMLE, self).__init__(endog, None, **kwds)
def nloglikeobs(self, params):
# Pay attention to these 3 lines, it shows how "params" array is organized:
rho = params[-1]
zCols =np.hstack([-np.inf, params[0:self.endog.shape[1]-1], np.inf])
zRows =np.hstack([-np.inf, params[self.endog.shape[0]-1:self.endog.shape[0]+self.endog.shape[1]-2], np.inf])
ll = []
pdebug = []
for i in range(self.endog.shape[1]):
for j in range(self.endog.shape[0]):
error, p_ij, inform = stats.mvn.mvndst([zCols[i], zRows[j]], [zCols[i+1],zRows[j+1]],
[0 if i == 0 else 1 if (i+1) == self.endog.shape[0] else 2,
0 if j == 0 else 1 if (j+1) == self.endog.shape[1] else 2], rho)
ll.append(-self.endog[j,i] * np.log(p_ij))
pdebug.append((p_ij, self.endog[j,i]))
return np.array(ll)
def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):
if start_params == None:
# Reasonable starting values
start_params = np.hstack([self.init_zCols[1:-1], self.init_zRows[1:-1], 0.0])
return super(PolychoricCorrMLE, self).fit(start_params=start_params,
maxiter=maxiter, maxfun=maxfun,
**kwds)
Probando el código¶
Este ejemplo fue tomado de [6] (http://john-uebersax.com/stat/tetra.htm#ex1).
Se tienen dos clasificadores que determinan la presencia o ausencia de esquizofrenia. Se hace la tabla de contingencia para estos dos clasificadores:
--------------------------
Rater 2
---------
Rater 1 - + Total
---------------------------
- 40 10 50
+ 20 30 50
--------------------------
Total 60 40 100
--------------------------
Uebersax obtuvo para este ejemplo los siguientes resultados, que indican que los dos clasificadores están fuertemente correlacionados:
rho 0.6071 (0.1152)
Rater 1 0.0000 (0.1253)
Rater 2 0.2533 (0.1268)
Para comprender el concepto de que esta tabla de contingencia puede ser reproducida bajo un modelo de una distribución normal bivariada, he construido la siguiente gráfica, que se basa en estos resultados, los límites latentes del Rater 1 y del Rater 2 y la correlación latente entre estas dos variables (rho).
mlp.rcParams['figure.figsize'] = 5, 5
x, y = np.mgrid[-5:5:.01, -5:5:.01]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x; pos[:, :, 1] = y
rv = stats.multivariate_normal([0.0, 0.0], [[1.0, 0.6071], [0.6071, 1.0]])
fig, ax = plt.subplots()
ax.contourf(x, y, rv.pdf(pos), zorder=0, cmap=plt.cm.coolwarm)
ax.xaxis.grid(True, zorder=10, color='#bebeff', linestyle='--', linewidth=1)
ax.yaxis.grid(True, zorder=10, color='#bebeff', linestyle='--', linewidth=1)
ax.axvline(0.2533, color= "white")
ax.axhline(0.0, color= "white")
ax.text(-1,-1,"40", color="black", fontsize=20)
ax.text(1,-1,"10", color="black", fontsize=20)
ax.text(-1,1,"20", color="black", fontsize=20)
ax.text(1,1,"30", color="black", fontsize=20)
Probamos reproducir los resultados anteriores con nuestro propio código:
t = PolychoricCorrMLE(np.array([[40,10], [20,30]]))
r = t.fit()
print("\n".join([ "%s\t%f\t(%f)"%(a,b,c) for a,b,c in zip(["Rater 2", "Rater 1", "rho"], r.params, r.bse)]))
Ha funcionado bastante bien 🙂
Las estimaciones y su error standard son bastante cercanos a los resultados obtenidos por Uebersax.
Ahora vamos a poner a prueba nuestro código con los datos del segundo ejemplo:
t = PolychoricCorrMLE(np.array([[58,52,1], [26,58,3],[8,12,9]]))
r = t.fit()
print("rho : ", r.params[-1], "(",r.bse[-1],")")
Conclusión¶
La correlación policórica o latente es la manera correcta de medir la relación entre dos variables ordinales. Espero que la imagen deje bien claro el significado del modelo normal bivariado.
También ha sido una muy buena práctica para aprender el uso de la clase GenericLikelihoodModel de statsmodels, aunque aún me falta aprender sobre los métodos de optimización y cuándo es necesario definir el gradiente y la matriz Hessian. Afortunadamente esto nos permite probar modelos sin hacer el arduo trabajo matemático que hizo, por ejemplo, Olsson en su artículo para desarrollar la estimación por MLE. En este caso ha bastado con comprender superficialmente el método y saber implementarlo en python.
Referencias¶
- Angeles G. Kolenikov S. (2014) The Use of Discrete Data in PCA: Theory, Simulations, and Applications to Socioeconomic Indices
- Olsson U. (1979) Maximum Likelihood Estimation of the Polychoric Correlation Coefficient.
- Pearson K. (1900) Mathematical Contributions to the Theory of Evolution. VII. On the Correlation of Characters not Quantitatively Measurable
- Pearson K. (1922) On Polychoric Coefficients of Correlation.
- Richie-Scott A. (1918) The Correlation Coefficient of a Polychoric Table.
- Uebersax J. (2015) Introduction to the Tetrachoric and Polychoric Correlation Coefficients