11. Inferencia estadística con tests no-paramétricos¶

11.1. Introducción¶

En esta lección revisaremos metodologías y tests estadísticos de la categoría no parámetrica, es decir que no presuponen una distribución en particular para las observaciones

Esto es contrario a las técnicas parámetricas que vimos en lecciones anteriores. Por ejemplo, un t-test parte del supuesto de que nuestros datos se distribuyen normal. Con este supuesto somos capaces de obtener una distribución para el estadístico de prueba ante la hipótesis nula. Esto a su vez nos permite calcular p-values y finalmente concluir si la hipótesis nula puede ser rechazada.

Advertencia

Si los supuestos que hicimos sobre la población son incorrectos entonces las conclusiones del método parámetrico pierden validez.

Sin embargo, en la práctica, los tests parámetricos funcionan bien ante “ligeras” desviaciones de los supuestos

Como recomendación general considere:

Comprobar los supuestos antes de hacer el análisis con métodos parámetricos
Si los supuestos no se cumplen una opción es usar inferencia no paramétrica

A continuación veremos tests estadísticos que nos permiten comprobar si nuestros datos siguen una distribución en particular. Luego veremos las alternativas no-parámetricas para los tests paramétricos más comunes. Finalmente veremos como usar la técnica de bootstrap resampling para estimar intervalos de confianza.

11.2. Evaluar la normalidad de una distribución¶

La normalidad es un supuesto en muchos procedimientos parámetricos, por ejemplo t-test, ANOVA y regresion lineal.

Tengamos en cuenta las características principales de la distribución normal:

Es unimodal con media igual a su moda
Es simétrica en torno a la media
Está concentrada en torno a la media (~68% a \(\pm \sigma\) y ~95% a \(\pm 2\sigma\))

Si tenemos suficientes muestras podriamos observar el histograma y comprobar visualmente si esto se cumple. A modo de ejemplo observe las siguientes distribuciones y discuta:

¿Cuáles distribuciones son no-normales?
¿Cuál o cuáles características de la distribución normal no se están cumpliendo en cada caso?

np.random.seed(1234) # Por reproducibilidad
data1 = ss.norm(loc=4, scale=2).rvs(1000) # Normal
data2 = ss.gamma(a=2, scale=1).rvs(1000) # Gamma
data3 = ss.uniform(loc=-2, scale=10).rvs(1000) #Uniforme

p = []
for i, data in enumerate([data1, data2, data3]):
    bins, edges = np.histogram(data, bins=20, range=(-4, 10), density=True)
    p.append(hv.Histogram((edges, bins), kdims='Datos', vdims='Densidad'))
hv.Layout(p).cols(3).opts(hv.opts.Histogram(width=250, height=250))

La visualización mediante histograma es útil como un paso exploratorio. Sin embargo podemos comprobar normalidad de forma más rigurosa utilizando un test de hipótesis

Se puede formular un test sencillo basándonos en que la distribución normal debiese tener una simetría igual a cero. La simentría se puede medir con el tercer momento estadístico: skewness

Con scipy podemos calcular la simetría usando la función scipy.stats.skew. Además scipy.stats.skewtest implemente un test en base a skewness donde la hipótesis nula es que la distribución corresponde a una normal

Por ejemplo para las tres distribuciones anteriores tenemos que:

for i, data in enumerate([data1, data2, data3]):
    display(ss.skewtest(data))

SkewtestResult(statistic=-0.3084039968473185, pvalue=0.757774941053838)

SkewtestResult(statistic=13.939252143051059, pvalue=3.6575254745841706e-44)

SkewtestResult(statistic=0.010891784869360446, pvalue=0.9913097848349443)

Para la muestra data2 rechazamos la hipótesis de que la simetría corresponda al de una distribución normal (p-value menor a 0.05)

Para las muestras data1 y data3 no hay suficiente evidencia para rechazar la hipótesis nula

Podemos complementar lo anterior con un test sobre la kurtosis, el cuarto momento de la distribución. En la definición por defecto de la kurtosis (Fisher) su valor debiese ser cero si la distribución subyacente es normal

Con scipy podemos calcular la kurtosis con la función scipy.stats.kurtosis. Además está implementado un test para la hipótesis de que la kurtosis corresponde al de una distribución normal: scipy.stats.kurtosistest

Para las tres muestras anteriores tenemos

for i, data in enumerate([data1, data2, data3]):
    display(ss.kurtosistest(data))

KurtosistestResult(statistic=0.2922804424326845, pvalue=0.7700722109486459)

KurtosistestResult(statistic=8.282243229464225, pvalue=1.208763870924613e-16)

KurtosistestResult(statistic=-25.310331594927383, pvalue=2.4583647703898043e-141)

Para la muestra data2 rechazamos las hipótesis de que la kurtosis y simetría correspondan al de una distribución normal

Para la muestra data3 rechazamos la hipótesis de que la kurtosis corresponde al de una distribución normal

Para la muestra data1 no hay suficiente evidencia para rechazar las hipótesis

Estos tests son simples pero particulares para la distribución normal. A continuación veremos formas más generales para comparar muestras con distribuciones teóricas.

11.3. Comparando distribuciones empíricas y teóricas¶

A continuación revisaremos algunos métodos para evaluar si la distribución empírica de nuestros datos sigue alguna distribución teórica en particular.

En los ejemplos usaremos una distribución teórica normal sin embargo puede usarse cualquier otra distribución bajo ciertas condiciones que se detallan más adelante.

11.3.1. Test de Kolmogorov-Smirnov (KS)¶

El test de KS es un test no-paramétrico que examina la distancia de las funciones de distribución acumulada (CDF) entre dos muestras o entre una muestra y una distribución teórica. Nos concentraremos en el caso de una muestra, es decir cuando queremos verificar si nuestros datos siguen una distribución teórica en particular.

El estadístico de test para KS es

\[ D_n = \sup_x |F_n(x) - F(x)| \]

donde \(F(x) = \int_{-\infty}^x f(x) dx = P(X<x)\) es la CDF teórica y

\[ F_n(x) = \frac{1}{N} \sum_{i=1}^N \mathbb{1}[X_i<x] \]

es la distribución acumulada empírica (ECDF).

Según la definición anterior \(D_n\) es el mayor de los errores absolutos entre la CDF teórica y empírica

Las hipótesis de la prueba KS de una muestra son:

\(\mathcal{H}_0\): Los datos siguen la distribución teórica especificada
\(\mathcal{H}_A\): Los datos no siguen la distribución teórica especificada

Si \(D_n\) es menor que el valor crítico no podemos rechazar \(H_0\)

Python

La librería scipy tiene implementado el test de KS. Sus argumentos son

scipy.stats.kstest(rvs, # Arreglo de muestras o función rvs de scipy.stats
                   cdf, # String o función cdf de scipy.stats para test de una muestra
                   args=(), # Argumentos que se pasan a rvs y cdf si son funciones
                   N=20, # Número de muestras generadas si rvs es una función
                   alternative='two-sided', # String que define la hipótesis alternativa
                   ...
                  )

La función retorna una estructura con el valor del \(D_n\) y el p-value

R

El paquete stats tiene implementado el test de KS. Sus argumentos son

ks.test(x, # Vector con datos numéricos
        y, # String o función cdf para test de una muestra
        alternative = "two.sided", # String que define la hipótesis alternativa
        exact = NULL # Indica si el p-value se calcula de forma exacta
       )

La función returna una objeto htest con los componentes statistic para \(D_n\) y p.value para el p-value

A continuación se muestra como utilizar scipy.stats.kstest para comparar data1 con una distribución normal estándar:

alpha = 0.05
theorical_cdf = scipy.stats.norm(loc=0, scale=1).cdf
statistic, pvalue = scipy.stats.kstest(data1, theorical_cdf)
pvalue < alpha

True

Rechazamos la hipótesis nula de que data1 tiene distribución \(\mathcal{N}(\mu=0, \sigma=1)\) con confianza de 95%

En cambio si hubieramos utilizado una distribución normal con los parámetros que generaron la muestra tendríamos:

theorical_cdf = scipy.stats.norm(loc=4, scale=2).cdf
statistic, pvalue = scipy.stats.kstest(data1, theorical_cdf)
pvalue < alpha

False

No hay suficiente evidencia para rechazar la hipótesis de que data1 tiene distribución \(\mathcal{N}(\mu=4,\sigma=2)\)

El test de Kolmogorov-Smirnov es ampliamente usado pero tiene algunas limitaciones

Advertencia

El test de KS sólo sirve para comparar distribuciones continuas y univariadas

En el caso de querer comparar distribuciones discretas la alternativa es el test no paramétrico de \(\chi^2\)

Advertencia

El test de KS es más sensible en el centro de la distribución que en las colas de la distribución

Advertencia

El test de KS requiere que la distribución teórica esté completamente especificada. Estimar los parámetros a partir de los datos (por ejemplo utilizando MLE) altera los valores críticos del estadístico y por ende no se puede confiar en los p-values obtenidos.

Algunas alternativas:

Si se requiere probar normalidad irrespecto, de cual sea la media y la varianza, se puede utilizar el test de Shapiro-Wilks o el test de Lilliefors
Si se requiere usar una CDF teórica con parámetros estimados a partir de los datos se podrían descartar los p-values del test y obtener p-values (e intervalo de confianza) empíricos para el estadístico de KS por medio de la técnica de bootstrap

En la próxima lección se verá el procedimiento general de bootstrap. A continuación: Veamos como utilizar los test de Shapiro-Wilks y Lilliefors:

for data in [data1, data2, data3]:
    test_statistic, pvalue = ss.shapiro(data)
    print(pvalue < alpha)

False
True
True

for data in [data1, data2, data3]:
    test_statistic, pvalue = sm.stats.lilliefors(data, dist='norm')
    print(pvalue < alpha)

False
True
True

Ambos test concuerdan en que las muestras data2 y data3 no son normales. En cambio para data1 no hay suficiente evidencia para rechazar la hipótesis de normalidad

11.3.2. Probability plots (PP)¶

También se pueden comparar distribuciones de forma gráfica utilizando probability plots

Se calculan los cuantiles (qq) o percentiles (pp) de dos muestras o de una muestra y una distribución teórica
Se grafica uno en función del otro
Mientras más se parezca el resultado a una linea recta más similares son las distribuciones

Podemos usar scipy.stats.probplot para comparar una muestra con una distribución teórica. A continuación se muestra la comparación entre los probability plot para las tres muestras anteriores versus una distribución teórica normal.

Nota

Dado que la prueba es basada en rangos no necesitamos entregar parámetros de localización (media) o escala (varianza). Cualquier parámetro de otro tipo (e.g. shape) se debe especificar.

p = []
for i, data in enumerate([data1, data2, data3]):
    (osm, osr), (w, b, r2) = ss.probplot(data, dist="norm", fit=True)
    x = np.arange(-3, 4)
    y = x*w + b
    p.append(hv.Curve((osm, osr), 'Theorical', 'Empirical', label='Data') * hv.Curve((x, y), label='Ideal'))
hv.Layout(p).cols(3).opts(hv.opts.Curve(width=250, height=250, alpha=0.75), 
                          hv.opts.Overlay(legend_position='top_left'))

11.4. Elementos de inferencia no paramétrica¶

Inferencia se refiere a

estimar parámetros o probar hipótesis sobre una población en base a muestras.

En inferencia no parámetrica usamos estadísticos cuya distribución no depende de supuestos sobre la distribución de la población.

Importante

No paramétrico \(\equiv\) No asumimos una distribución específica para la población

Sin embargo:

Advertencia

No paramétrico \(\neq\) Libre de supuestos o Libre de (hiper)parámetros

Por ejemplo un supuesto que muchos tests no parámetricos tienen en común es que las observaciones son independientes. Otro supuesto o requisito común es que se cumpla un cierto número mínimo de observaciones.

¿Por que no usar siempre tests no paramétricos?

En general los tests no parámetricos tienen un menor poder estadístico (sensibilidad) que sus contrapartes paramétricas (siempre y cuando sus supuestos se cumplan).

A continuación veremos algunos conceptos fundamentales de inferencia no parámetrica para luego revisar algunos tests no paramétricos que son alternativas a los tests vistos en lecciones anteriores.

11.4.1. Función cuantil y estadístico de orden¶

La CDF se define como

\[ P(X<x) = F(x) = p, \]

donde \(p\in [0, 1]\)

En algunos casos nos interesa saber que valor \(x\) corresponde a un valor dado \(p\). Esto se hace con la función cuantil o CDF inversa

\[ Q(p) = \inf_{x\in \mathbb{R}} p \leq F(x) \]

En scipy podemos calcular \(x\) a partir de \(p\) usando el atributo ppf (percent-point function) de una distribución de scipy.stats, por ejemplo

# Que valor de x corresponde a un 50% de la distribución normal estándar?
ss.norm.ppf(0.5)

0.0

# y a un 99%?
ss.norm.ppf(0.98)

2.0537489106318225

Sea \({X_1, X_2, X_3, \ldots, X_N}\) un conjunto de VAs i.i.d. y \({x_1, x_2, x_3, \ldots, x_N}\) una muestra aleatoria.

Si asumimos que la población tiene una CDF continua (no hay dos VA con el mismo valor) entonces podemos hacer un ordenamiento único:

\[ x_{(1)} < x_{(2)} < x_{(3)} < \ldots < x_{(N)} \]

que se llaman colectivamente el estadístico de orden de la muestra aleatoria. El estadístico específico \(x_{(r)}\) se llama el estadístico de orden r

11.4.2. Distribución acumulada empírica (ECDF)¶

Con el estadístico de orden se construye la ECDF:

\[\begin{split} F_N(x) = \begin{cases} 0 & x < x_{(1)} \\ \frac{i}{N} & x_{(i-1)} \leq x < x_{(i)}, i=2,3,\ldots,N \\ 1 & x_{(N)} < x\end{cases} \end{split}\]

que podemos obtener en Python utilizando statsmodels.distributions.empirical_distribution.ECDF. Esta función retorna un objeto que podemos evaluar para visualizar la ECDF en el rango que deseemos.

from statsmodels.distributions.empirical_distribution import ECDF

p = []
for i, data in enumerate([data1, data2, data3]):
    x_t = np.linspace(-3, 10, num=100)
    y_t = ECDF(data)(x_t)
    p.append(hv.Curve((x_t, y_t), 'x', 'Fn(x)', label=f'data{i}'))
hv.Layout(p).cols(3).opts(hv.opts.Curve(width=250, height=250, alpha=0.75))

Notar que el ECDF es un estimador insesgado

\[ \mathbb{E}[F_N(x)] = F(x) \]

y su varianza tiende a cero con \(N\)

\[ \text{Var}[F_N(x)] = \frac{1}{N} F(x) ( 1- F(x)) \]

Adicionalmente la ECDF:

es un estimador consistente de \(F(x)\) (converge en probabilidad)
(estandarizada) es un asintoticamente normal estándar

11.4.3. Función cuantil empírica¶

También podemos usar el estadístico de orden para construir la función cuantil empírica

\[\begin{split} Q_N(u) = \begin{cases} x_{(1)} & 0 < u \leq \frac{1}{N} \\ x_{(2)} & \frac{1}{N} < u \leq\frac{2}{N} \\ \vdots & \vdots \\ x_{(N)} & \frac{(N-1)}{N} < u \leq 1 \\ \end{cases} \end{split}\]

En Python podemos usar np.quantile y np.percentile para calcular cuantiles y percentiles de una muestra, respectivamente

data = np.random.randn(1000)
np.quantile(data, 0.5) # Valores entre 0 y 1 

0.008435585886555123

np.percentile(data, 50) # Valores entre 0 y 100

0.008435585886555123

11.4.4. Aplicaciones de los estadísticos de orden¶

Dos estadísticos muy usuales que vienen de los estadísticos de orden son

El rango: \(x_{(N)} - x_{(1)}\)
La mediana: \(x_{(n/2)}\) (n par)

Los estadísticos de orden suelen ser más robustos que sus contrapartes en la presencia de outliers

En Python podemos calcular la mediana con np.median

x = np.random.randn(10) # media cero y desviación estándar uno
x[9] = 100 # Inserto un outlier

# Le media es:
np.mean(x)

10.173076601431145

# En cambio la mediana:
np.median(x)

0.16098636624373827

es claramente menos afectada por el outlier.

Otro estadístico muy usado es el rango intercuartil, que corresponde a la diferencia entre cuantil 0.75 y 0.25

q75, q25 = np.quantile(x, [0.75, 0.25])
q75 - q25

0.9432496155282291

este es una alternativa robusta para la desviación estándar

np.std(x)

29.951992146865294

que claramente ha sido afectada por el outlier insertamos.

11.5. Tests no paramétricos¶

Las propiedades de la ECDF nos permiten construir distribuciones nulas sin necesidad de asumir una distribución para la población

Existen muchos tests no-paramétricos. Anteriormente revisamos el test de KS, Wilcoxon y Lilliefors. A continuación revisaremos dos tests no paramétricos que son alternativas al t-test.

11.5.1. Test de Mann-Whitney U (MWU)¶

Es una prueba no-paramétrica para comparar dos muestras independientes.

Ejemplo de muestras independientes: notas de niños y niñas en una prueba
Ejemplo de muestras dependientes: notas de un mismo curso en pruebas consecutivas

El objetivo del test es

probar si las muestras independientes provienen de una misma población en función de su tendencia central.

Esto es similar a lo que realiza el test parámetrico t-test en donde se comparan dos medias \(\mu_0\) y \(\mu_1\), asumiendo que la distribución subyacente es normal.

El test de Mann-Whitney U no supone distribución para la población y algunos lo interpretan como una comparación entre medianas.

La hipótesis nula del test de MWU es que no hay diferencia entre las distribuciones muestrales

El algoritmo del test es:

Ordenar las observaciones de ambas muestras juntas
Sumar los rangos de la primera muestra \(R_1\). \(R_2\) queda definida ya que \(R_1 + R_2 = \frac{N(N+1)}{2}\)
Se calcula el estadístico \(U_1 = R_1 - \frac{N_1(N_1+1)}{2}\), donde \(N_1\) es el número de observaciones de la primera muestra. Por simetría \(U_2\) queda inmediatamente especificado
Se usa \(U = \min(U_1, U_2)\) para calcular el p-value

Intuición

La prueba mezcla y ordena las observaciones de ambas muestras. Intutivamente, si las muestras son similares entonces sus observaciones deberían mezclarse en el ordenamiento. Si en cambio se clusterizan entonces las distribuciones son distintas

Este test está implementado en:

scipy.stats.mannwhitneyu(x, # Muestra 1
                         y, # Muestra 2
                         alternative='two-sided' # Las otras opciones son 'less' y 'greater'
                         method='auto'
                        )

La función retorna el estadístico \(U\) y el p-value.

Nota

Si se usa method='auto', el comportamiento del test difiere en base al número de muestras. Si el número de muestras es menor a 8 se retorna un p-value exacto, de lo contrario se retorna un p-value asintótico. Esta versión del test está disponible en scipy>1.7.0

11.5.2. Test de los rangos con signo de Wilcoxon¶

Si las muestras que se quieren comparar son dependientes no se puede usar Mann Whitney U. En su lugar se puede usar el test de rango con signo de Wilcoxon.

Este test funciona como alternativa no-parámetrica al t-test pareado en caso de que los datos violen el supuesto de normalidad.

Un test pareado se usa para probar si hubo diferencia significativa entre el “antes” y el “después” de aplicar cierto tratamiento o intervención.

Los supuestos de este test son:

Los datos son continuos
Los datos son pareados y vienen de la misma población
Los pares son independientes y se escogen aleatoriamente

Sea \(x_{11}, x_{12}, \ldots, x_{1N}\) y sus pares \(x_{21}, x_{22}, \ldots, x_{2N}\).

La hipótesis nula del test es que la mediana de las diferencias entre pares de observaciones es nula.

El algoritmo del test de Wilcoxon es:

Se calculan las diferencias \(z_i = x_{2i} - x_{1i}\)
Se ordenan los valores absolutos de las diferencias \(|z|_{(1)}, |z|_{(2)}, \ldots, |z|_{(N)}\) y se reserva el signo de cada diferencia
Se suman los rangos de las diferencias con signo positivo: \(W^{+}\)
Se suman los rangos de las diferencias con signo negativo: \(W^{-}\)
El estadístico de prueba es \(\min (W^{+}, W^{-})\) sobre el dataset reducido (se eliminan los \(z_{(i)}=0\))

Intuición

Si la hipótesis nula es cierta, esperariamos que \(W^{-}\) y \(W^{+}\) sean similares. Si esto no se cumple se esperaría que \(W^{+}\) fuera mayor que \(W^{-}\)

El test está implementado en:

scipy.stats.wilcoxon(x, # El primer set de muestras. Si y=None entonces x debe ser igual a la diferencia z
                     y=None, # El segundo set de muestras
                     alternative='two-sided', # Para un test de dos colas
                     mode='auto', 
                     ...
                    )

donde el argumento:

mode=approx indica si se ocupa una aproximación normal para calcular el p-value
mode=exact indica si se obtiene el valor de una tabla. Sólo valido para menos de 25 muestras
mode=auto Usa approx para más de 25 muestras y exact para menos de 25

11.5.3. Test de Kruskall-Wallis: Alternativa a one-way ANOVA¶

Su objetivo es probar que existen diferencias significativas en una variable dependiente continua en función de una variable independiente categórica (2 o más grupos)

Está implementado en scipy.stats.kruskal

11.5.4. Kendal \(\tau\): Alternativa al coeficiente de correlación de Pearson¶

Su objetivo es probar si dos variables son estadísticamente dependientes

Está implementando en scipy.stats.kendalltau

Herramientas estadísticas para la investigación

Inferencia estadística con tests no-paramétricos

Contenido

11. Inferencia estadística con tests no-paramétricos¶

11.1. Introducción¶

11.2. Evaluar la normalidad de una distribución¶

11.3. Comparando distribuciones empíricas y teóricas¶

11.3.1. Test de Kolmogorov-Smirnov (KS)¶

11.3.2. Probability plots (PP)¶

11.4. Elementos de inferencia no paramétrica¶

11.4.1. Función cuantil y estadístico de orden¶

11.4.2. Distribución acumulada empírica (ECDF)¶

11.4.3. Función cuantil empírica¶

11.4.4. Aplicaciones de los estadísticos de orden¶

11.5. Tests no paramétricos¶

11.5.1. Test de Mann-Whitney U (MWU)¶

11.5.2. Test de los rangos con signo de Wilcoxon¶

11.5.3. Test de Kruskall-Wallis: Alternativa a one-way ANOVA¶

11.5.4. Kendal \(\tau\): Alternativa al coeficiente de correlación de Pearson¶