La prueba de hipótesis y el valor p

Llegamos al número más famoso y peor entendido de la estadística: el valor p 🏆

Y te lo voy a enseñar de una forma que casi nadie usa: calculándolo a mano, barajando cartas. Cuando lo veas así, ya no se te va a olvidar nunca, y las fórmulas de después van a ser un atajo y no una caja negra.

La pregunta, en su forma más honesta

En el capítulo 9 encontramos esto:

• Horeca cierra el 63
• 21% de sus ventas y Minimarket el 56
• 93%

Una diferencia de 6,28 puntos.

La pregunta es la de siempre: ¿eso es real o es ruido?

Y la forma de contestarla es preciosa. Consiste en suponer que no hay ninguna diferencia y ver qué tan raro sería, en ese mundo, encontrar los 6,28 puntos que encontramos.

A esa suposición se le llama hipótesis nula. Es el abogado del diablo: parte de que no pasa nada, y le toca a tus datos convencerla de lo contrario 👩‍⚖️

import pandas as pd
import numpy as np
from scipy import stats

URL = 'https://missyera.com/static/datasets/ventas-miss-yera.csv'


def carga_limpia(url):
    """La misma del capítulo 2."""
    v = pd.read_csv(url).drop_duplicates()
    v['ciudad'] = (v['ciudad'].str.strip().str.lower()
                   .str.normalize('NFKD')
                   .str.encode('ascii', 'ignore').str.decode('utf-8'))
    v['monto'] = pd.to_numeric(v['monto'].str.replace(',', '.'))
    for col in ['fecha', 'fecha_ultima_compra']:
        f = pd.to_datetime(v[col], format='%Y-%m-%d', errors='coerce')
        falta = f.isna() & v[col].notna()
        f[falta] = pd.to_datetime(v.loc[falta, col], format='%d/%m/%Y',
                                  errors='coerce')
        v[col] = f
    return v


v = carga_limpia(URL)

horeca = v.loc[v['segmento'] == 'Horeca', 'compro'].values
minimarket = v.loc[v['segmento'] == 'Minimarket', 'compro'].values
observada = horeca.mean() - minimarket.mean()

print('Horeca:     %.4f (%d ventas)' % (horeca.mean(), len(horeca)))
print('Minimarket: %.4f (%d ventas)' % (minimarket.mean(), len(minimarket)))
print('diferencia: %.6f' % observada)

Horeca:     0.6321 (742 ventas)
Minimarket: 0.5693 (801 ventas)
diferencia: 0.062787

Barajando las etiquetas

Aquí viene lo bonito 🃏

Si la hipótesis nula fuera cierta, o sea si el segmento no influyera en nada, entonces las etiquetas "Horeca" y "Minimarket" estarían puestas al azar sobre esas 1.543 ventas.

Así que vamos a hacer justo eso: mezclar las etiquetas mil veces y ver qué diferencias salen en un mundo donde no pasa nada.

generador = np.random.default_rng(7)
todas = np.concatenate([horeca, minimarket])
n_horeca = len(horeca)

diferencias = []
for _ in range(10000):
    generador.shuffle(todas)
    diferencias.append(todas[:n_horeca].mean() - todas[n_horeca:].mean())

diferencias = np.array(diferencias)

print('diferencias al azar: media %+.6f | desviacion %.6f'
      % (diferencias.mean(), diferencias.std()))
print('la mayor que salio por azar: %.6f' % np.abs(diferencias).max())
print('la nuestra:                  %.6f' % abs(observada))

diferencias al azar: media +0.000433 | desviacion 0.025086
la mayor que salio por azar: 0.092982
la nuestra:                  0.062787

Fíjate en lo que acabamos de fabricar: el mundo donde no pasa nada. En ese mundo las diferencias entre los dos grupos se mueven alrededor de cero con una desviación de 0,025 🌍

Y ahora el valor p, que es literalmente contar:

p = np.mean(np.abs(diferencias) >= abs(observada))

print('de 10000 mundos sin efecto, cuantos dieron una diferencia asi de grande o mas:')
print('   %d' % (np.abs(diferencias) >= abs(observada)).sum())
print('valor p = %.4f' % p)

de 10000 mundos sin efecto, cuantos dieron una diferencia asi de grande o mas:
   155
valor p = 0.0155

Ese es el valor p. Y ya está. 🎉

Se lee así: si el segmento no influyera en nada, solo el 1,55% de las veces vería una diferencia tan grande como la que vi.

No hay más misterio. Un valor p es la proporción de mundos-sin-efecto que producen algo tan llamativo como lo tuyo.

Y ahora la fórmula, que es el atajo

Lo que acabamos de hacer se llama prueba de permutación y tarda unos segundos. La prueba t hace lo mismo con una fórmula, en una línea:

resultado = stats.ttest_ind(horeca, minimarket, equal_var=False)

print('permutacion (barajando): p = %.4f' % p)
print('prueba t (formula):      p = %.4f' % resultado.pvalue)

permutacion (barajando): p = 0.0155
prueba t (formula):      p = 0.0118

0,0155 y 0,0118. Prácticamente lo mismo 🤝

La fórmula supone que las diferencias al azar siguen una campana, y por eso va más rápido. La permutación no supone nada: se fabrica el mundo sin efecto de verdad. Cuando las dos coinciden, como aquí, es que la suposición de la fórmula era razonable.

Yo uso la fórmula por comodidad y la permutación cuando el resultado es importante o los datos son raros. Y las dos veces sé lo que estoy calculando, que es de lo que iba este capítulo 😌

El contraste: cuando de verdad no hay nada

Hagamos lo mismo con dos ciudades, que llevamos nueve capítulos viendo que no separan nada:

trujillo = v.loc[v['ciudad'] == 'trujillo', 'compro'].values
piura = v.loc[v['ciudad'] == 'piura', 'compro'].values

print('trujillo %.4f | piura %.4f | diferencia %.4f'
      % (trujillo.mean(), piura.mean(), trujillo.mean() - piura.mean()))
print('prueba t: p = %.4f'
      % stats.ttest_ind(trujillo, piura, equal_var=False).pvalue)

trujillo 0.5924 | piura 0.5692 | diferencia 0.0232
prueba t: p = 0.4634

p = 0,4634. O sea que casi la mitad de los mundos sin efecto producirían una diferencia así 🤷

Y ojo con cómo se dice esto, porque es donde se resbala todo el mundo. No se dice "queda demostrado que las ciudades cierran igual". Se dice "no tengo evidencia de que sean distintas".

La diferencia importa: puede que Trujillo cierre mejor de verdad y que 471 ventas no alcancen para verlo. Ausencia de evidencia no es evidencia de ausencia 🕵️

El 0,05, que no tiene nada de mágico

Se corta en 0,05 por costumbre. Lo propuso Ronald Fisher en los años veinte como una guía razonable y se quedó para siempre.

Lo que sí es concreto es qué te cuesta ese corte. Vamos a verlo con datos donde sabemos que no hay ningún efecto porque los fabrico yo:

falsos = 0
for _ in range(1000):
    a = generador.normal(size=50)
    b = generador.normal(size=50)
    if stats.ttest_ind(a, b).pvalue < 0.05:
        falsos += 1

print('comparaciones sin ningun efecto: 1000')
print('salieron significativas:         %d' % falsos)
print('o sea el %.1f%%' % (falsos / 10))

comparaciones sin ningun efecto: 1000
salieron significativas:         57
o sea el 5.7%

El 5,7%, o sea el 5% prometido 🎯

Eso es lo que significa alfa = 0,05: estás aceptando equivocarte una de cada veinte veces cuando no hay nada. No es un defecto, es el precio que eliges pagar.

Y de ahí sale la tabla de los dos errores, que conviene tener clara:

Bajar alfa a 0,01 reduce las falsas alarmas y aumenta las veces que te pierdes algo real. No hay forma de reducir los dos a la vez sin conseguir más datos. Eso es el capítulo 12 💪

Las cuatro cosas que un valor p NO dice

Esta lista vale el capítulo entero. Cada una la he oído en una reunión 🙃

• ❌ "p = 0,0155 quiere decir que hay 98,45% de probabilidad de que la diferencia sea real." No. El p se calcula suponiendo que no hay diferencia; no puede decirte la probabilidad de que la haya.
• ❌ "p = 0,0155 quiere decir que hay 1,55% de probabilidad de que sea casualidad." Tampoco, y es la más común. Es la probabilidad de ver estos datos si fuera casualidad, que es al revés. Es el mismo error de dar la vuelta a una condicional del capítulo 7.
• ❌ "p chiquito quiere decir efecto grande." No. El p mezcla tamaño del efecto y cantidad de datos. Con muchísimas filas, un efecto ridículo sale con p diminuto.
• ❌ "p = 0,06 quiere decir que no hay nada." No. 0,049 y 0,051 son prácticamente el mismo resultado, y el corte en 0,05 es una convención.

Si te llevas una sola frase del capítulo, que sea esta: significativo no quiere decir importante 📌

El error del capítulo

Quiero comparar la satisfacción entre dos segmentos:

sa = v.loc[v['segmento'] == 'Horeca', 'satisfaccion']
sb = v.loc[v['segmento'] == 'Bodega', 'satisfaccion']

print(stats.ttest_ind(sa, sb))

TtestResult(statistic=np.float64(nan), pvalue=np.float64(nan), df=np.float64(nan))

nan por todas partes, sin error y sin aviso 🫥

La causa son los 231 nulos de satisfaccion. Cualquier cuenta que toque un nan devuelve nan, y como el resultado parece un resultado, se puede copiar a una diapositiva tal cual.

Un p de nan no es "no significativo": es "no calculé nada".

print(stats.ttest_ind(sa, sb, nan_policy='omit'))

TtestResult(statistic=np.float64(0.9507468810912626), pvalue=np.float64(0.341903692047498), df=np.float64(1345.0))

Ahora sí: p = 0,3419, o sea que la satisfacción no distingue esos dos segmentos. Que es lo que ya sospechábamos desde el capítulo 3.

Y el error que sí te frena, por si te lo preguntas:

stats.ttest_ind(sa.dropna(), sb.dropna(), alternative='mayor')

ValueError: `alternative` must be 'less', 'greater', or 'two-sided'.

Ese está en inglés y hay que escribirlo greater. Lo pongo porque la comparación es justa: el parámetro mal escrito te para, y los datos mal preparados no 😤

Practica 💪

may = v.loc[v['segmento'] == 'Mayorista', 'compro'].values
bod = v.loc[v['segmento'] == 'Bodega', 'compro'].values

r = stats.ttest_ind(may, bod, equal_var=False)
print('mayorista %.4f | bodega %.4f | diferencia %.4f'
      % (may.mean(), bod.mean(), may.mean() - bod.mean()))
print('t = %.4f' % r.statistic)
print('p = %.3e' % r.pvalue)

mayorista 0.7240 | bodega 0.3748 | diferencia 0.3492
t = 14.2705
p = 2.726e-43

for tamano in [100, 200, 400, 742]:
    a = generador.choice(horeca, min(tamano, len(horeca)), replace=False)
    b = generador.choice(minimarket, min(tamano, len(minimarket)), replace=False)
    r = stats.ttest_ind(a, b, equal_var=False)
    print('con %3d de cada uno -> diferencia %+.4f | p = %.4f'
          % (tamano, a.mean() - b.mean(), r.pvalue))

con 100 de cada uno -> diferencia +0.0000 | p = 1.0000
con 200 de cada uno -> diferencia +0.1050 | p = 0.0339
con 400 de cada uno -> diferencia +0.0650 | p = 0.0610
con 742 de cada uno -> diferencia +0.0580 | p = 0.0225

for referencia in [800, 750, 803]:
    r = stats.ttest_1samp(v['monto'], referencia)
    print('contra %d soles -> t = %+7.4f | p = %.4f' % (referencia, r.statistic, r.pvalue))

contra 800 soles -> t = +0.2740 | p = 0.7841
contra 750 soles -> t = +3.9158 | p = 0.0001
contra 803 soles -> t = +0.0554 | p = 0.9558

ma = v.loc[v['segmento'] == 'Horeca', 'monto'].values
mb = v.loc[v['segmento'] == 'Minimarket', 'monto'].values
obs = ma.mean() - mb.mean()

juntos = np.concatenate([ma, mb])
difs = []
for _ in range(10000):
    generador.shuffle(juntos)
    difs.append(juntos[:len(ma)].mean() - juntos[len(ma):].mean())
difs = np.array(difs)

print('diferencia observada: %.2f soles' % obs)
print('permutacion: p = %.4f' % np.mean(np.abs(difs) >= abs(obs)))
print('prueba t:    p = %.3e' % stats.ttest_ind(ma, mb, equal_var=False).pvalue)

diferencia observada: 340.26 soles
permutacion: p = 0.0000
prueba t:    p = 2.535e-143

for i in range(5):
    basura = generador.normal(size=len(v))
    grupo_a = basura[v['compro'] == 1]
    grupo_b = basura[v['compro'] == 0]
    r = stats.ttest_ind(grupo_a, grupo_b)
    marca = '  <-- SIGNIFICATIVA' if r.pvalue < 0.05 else ''
    print('columna de basura %d -> p = %.4f%s' % (i + 1, r.pvalue, marca))

columna de basura 1 -> p = 0.0700
columna de basura 2 -> p = 0.5952
columna de basura 3 -> p = 0.9215
columna de basura 4 -> p = 0.1748
columna de basura 5 -> p = 0.1566

print('comparando montos de Horeca contra Minimarket')
print('  medias:   %.2f contra %.2f' % (ma.mean(), mb.mean()))
print('  medianas: %.2f contra %.2f' % (np.median(ma), np.median(mb)))
print()
print('prueba t (medias):        p = %.3e'
      % stats.ttest_ind(ma, mb, equal_var=False).pvalue)
print('Mann-Whitney (posicion):  p = %.3e'
      % stats.mannwhitneyu(ma, mb).pvalue)

comparando montos de Horeca contra Minimarket
  medias:   755.18 contra 414.92
  medianas: 743.00 contra 413.74

prueba t (medias):        p = 2.535e-143
Mann-Whitney (posicion):  p = 1.022e-142

Lo que te llevas

• 🃏 Un valor p es contar: de 10.000 mundos sin efecto, ¿cuántos producen algo tan grande como lo tuyo? Aquí 155, o sea p = 0,0155.
• ⚡ La fórmula es el atajo de eso. La permutación dio 0,0155 y la prueba t 0,0118.
• 🚫 Un p alto no demuestra que no haya diferencia. Ausencia de evidencia no es evidencia de ausencia.
• 🎯 Alfa = 0,05 significa equivocarte 1 de cada 20 veces cuando no hay nada. Con 1.000 comparaciones vacías salieron 57 significativas.
• ❌ El p no es la probabilidad de que sea casualidad, ni el tamaño del efecto, ni una medida de importancia.
• 🎰 Con muestras chicas el p es un sorteo: el mismo efecto dio 0,2000, 0,3626, 0,0392 y 0,1263 según qué 100 filas tocaran.
• 👻 Un p de nan no es "no significativo", es "no calculé nada". Y llega sin avisar en cuanto hay nulos.

Qué viene ahora

Ya sabes qué es un valor p. En el capítulo 11 vemos las cuatro pruebas concretas que cubren el 90% de los casos reales, y sobre todo cómo elegir cuál toca, que es donde se atasca todo el mundo 🧰