Cuánto sube por cada uno

Última herramienta del libro, y la que más voy a usar yo en un trabajo de verdad 📏

El capítulo 13 terminó diciendo que dos cosas se mueven juntas. Está bien y se queda corto, porque en una reunión nadie pregunta si algo se relaciona: preguntan cuánto.

Y una pregunta para ti antes de empezar: ¿cuántas veces has dicho "sí influye" sin poder decir cuánto ni con qué margen? Este capítulo es para no tener que volver a hacerlo 🎯

La regresión contesta eso con un número que se lee en soles. Vamos a escribirla nosotras, que son dos líneas.

La recta que menos se equivoca

La idea es la que ya te imaginas: pasar una recta por la nube de puntos. Lo que hay que decidir es cuál recta, porque caben infinitas.

El criterio se llama mínimos cuadrados y es este: de todas las rectas posibles, quédate con la que hace más chica la suma de los errores al cuadrado. Al cuadrado por dos razones, y las dos importan: así los errores por arriba y por abajo no se cancelan, y así un error grande pesa mucho más que dos chicos.

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

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

def carga_limpia(url):
    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(',', '.'))
    return v

v = carga_limpia(URL)
v = v[v['monto'] > 0]
x = v['unidades'].to_numpy(float)
y = v['monto'].to_numpy(float)
print('filas:', len(v))
print('correlacion entre unidades y monto:', round(float(np.corrcoef(x, y)[0, 1]), 4))

filas: 2979
correlacion entre unidades y monto: 0.0125

Empezamos por la pregunta más obvia de todas: ¿cuánto sube la venta por cada unidad de más que se lleva el cliente?

Y ya vemos algo: la correlación es 0,0125, o sea prácticamente cero. Guárdalo, que en un momento nos va a servir 🔖

Escrita a mano, en dos líneas

pendiente = ((x - x.mean()) * (y - y.mean())).sum() / ((x - x.mean()) ** 2).sum()
intercepto = y.mean() - pendiente * x.mean()
print(f'pendiente  : {pendiente:.4f}')
print(f'intercepto : {intercepto:.4f}')

pendiente  : 1.6175
intercepto : 793.6360

Eso es una regresión lineal entera. En serio, no hay más 😌

La primera línea se lee bien si la miras por partes. Arriba, cuánto se mueven x e y juntas respecto de sus promedios. Abajo, cuánto se mueve x sola. La división es "de todo lo que se mueve x, qué parte arrastra a y".

Y el intercepto sale de una condición bonita: la recta siempre pasa por el punto de los dos promedios. Sabiendo la pendiente y ese punto, la recta queda fija.

Comprobemos contra la librería, que es como se comprueba todo en este libro:

r = stats.linregress(x, y)
print(f'scipy dice : pendiente {r.slope:.4f}   intercepto {r.intercept:.4f}')
print(f'diferencia : {abs(r.slope - pendiente):.2e}')

scipy dice : pendiente 1.6175   intercepto 793.6360
diferencia : 2.44e-15

Idénticas hasta la quinceava cifra 🎉

Cómo se lee una pendiente

Aquí está lo que hace que esto valga la pena, y es una sola frase que conviene aprenderse:

Nuestra pendiente es 1,6175, y las unidades son soles por unidad. O sea: por cada unidad más que se lleva el cliente, la venta sube 1,62 soles.

Y el intercepto, 793,64, es lo que predice la recta para un cliente que se lleva cero unidades. Que no significa nada aquí, porque nadie compra cero, y es justo el aviso: el intercepto solo se interpreta si el cero está dentro del rango de tus datos.

Ahora bien. 1,62 soles por unidad suena razonable, se puede poner en una diapositiva y alguien lo va a repetir en una reunión 😬

El número que decide si eso significa algo

La pendiente sale siempre. Le das cualquier par de columnas y te devuelve un número, sin quejarse. Lo que hay que preguntarle es con cuánta precisión, y para eso está su intervalo, igual que en el capítulo 9.

ic = 1.96 * r.stderr
print(f'pendiente        : {r.slope:.4f} soles por unidad')
print(f'error estandar   : {r.stderr:.4f}')
print(f'intervalo del 95%: [{r.slope - ic:.4f}, {r.slope + ic:.4f}]')
print(f'valor p          : {r.pvalue:.4f}')

pendiente        : 1.6175 soles por unidad
error estandar   : 2.3648
intervalo del 95%: [-3.0176, 6.2525]
valor p          : 0.4941

Mira el intervalo: de -3,02 a 6,25 😳

Cruza el cero de lado a lado. O sea que con estos datos no podemos descartar que la pendiente de verdad sea negativa, ni que sea cero, ni que sea 6. El valor p de 0,4941 dice lo mismo con otro formato.

Traducido al idioma de la reunión: llevarse más unidades no está relacionado con gastar más en esta distribuidora. Y tiene sentido si lo piensas, porque una cosa es cuántas cajas te llevas y otra es de qué producto.

Quédate con esto, que es la lección del capítulo entero: el número siempre sale, y decidir si significa algo es el trabajo 🧐

Entonces, ¿dónde está la venta?

La regresión sirve para variables de categoría, no solo para números, y ahí es donde estos datos sí tienen algo que contar. Metemos el segmento:

import statsmodels.formula.api as smf

modelo = smf.ols('monto ~ C(segmento)', data=v).fit()
print(f'R2 = {modelo.rsquared:.4f}')
for nombre, coef in modelo.params.items():
    bajo, alto = modelo.conf_int().loc[nombre]
    print(f'{nombre:26} {coef:9.2f}   IC 95% [{bajo:8.2f}, {alto:8.2f}]')

R2 = 0.7530
Intercept                     181.65   IC 95% [  154.19,   209.11]
C(segmento)[T.Horeca]         579.99   IC 95% [  541.66,   618.33]
C(segmento)[T.Mayorista]     1693.75   IC 95% [ 1655.50,  1731.99]
C(segmento)[T.Minimarket]     237.71   IC 95% [  200.06,   275.36]

Ahora sí 💰

Se lee así, y es más fácil de lo que parece. El Intercept de 181,65 es la categoría que falta, que es Bodega: una bodega compra 181,65 soles de promedio. Todo lo demás se mide contra ella.

O sea que un mayorista compra 1.693,75 soles más que una bodega, y su intervalo va de 1.655,50 a 1.731,99, que no se acerca al cero ni de lejos. Eso sí se puede llevar a una reunión.

Y el R cuadrado de 0,7530 dice qué parte de la variación del monto queda explicada: el 75,3%. Compáralo con lo de antes, donde con unidades salía 0,0002, y ahí tienes de un vistazo cuál de las dos preguntas era la buena 📊

Por qué existe la regresión, de verdad

Todo lo que hicimos hasta aquí se podía hacer comparando grupos. La razón de ser de la regresión es otra: mirar varias variables a la vez y decir cuánto aporta cada una con las demás fijas.

Añadamos las unidades al modelo del segmento:

completo = smf.ols('monto ~ C(segmento) + unidades', data=v).fit()
print(f'R2 solo con segmento     : {modelo.rsquared:.4f}')
print(f'R2 anadiendo unidades    : {completo.rsquared:.4f}')
print(f'coeficiente de unidades  : {completo.params["unidades"]:.4f}')
print(f'su valor p               : {completo.pvalues["unidades"]:.4f}')

R2 solo con segmento     : 0.7530
R2 anadiendo unidades    : 0.7530
coeficiente de unidades  : 0.0365
su valor p               : 0.9752

El R cuadrado no se movió y el coeficiente de unidades se desplomó de 1,6175 a 0,0365, con un p de 0,9752 😅

Esa es la lectura correcta: sabiendo el segmento, las unidades no aportan absolutamente nada. Y es lo que ninguna comparación de grupos te puede decir, porque para saberlo hay que mirar las dos a la vez.

Este es el mecanismo que hay detrás de la frase "controlando por". Cuando alguien dice "controlando por el tamaño de la empresa", está diciendo exactamente esto: metí esa variable en el modelo y miro qué queda del resto.

Y también es el aviso de siempre, el del capítulo 15: un coeficiente no es una causa. Que ser mayorista pese 1.693,75 no significa que convertir una bodega en mayorista le vaya a subir la venta en 1.693,75. Significa que los mayoristas compran más, que ya lo sabíamos, ahora con número e intervalo.

Los residuos, o dónde se ve si la recta sirve

Un residuo es lo que le sobra a la predicción: lo real menos lo predicho. Y son la parte del diagnóstico que casi nadie mira 🔍

pred = modelo.fittedvalues
res = modelo.resid
print(f'media de los residuos: {res.mean():.6f}')
tercio = pd.qcut(pred, 3, labels=['bajo', 'medio', 'alto'])
print('desviacion del residuo por tercio de prediccion:')
print(res.groupby(tercio, observed=True).std().round(1).to_string())

media de los residuos: -0.000000
desviacion del residuo por tercio de prediccion:
bajo     109.0
medio    262.2
alto     676.0

La media es cero, y eso no es una buena noticia: siempre da cero, por construcción de los mínimos cuadrados. No comprueba nada.

Lo que sí dice algo es la tabla de abajo, y dice bastante: el error típico pasa de 109 soles donde la predicción es baja a 676 donde es alta. Seis veces más 😳

Eso se llama heterocedasticidad, que es una palabra horrible para una idea simple: el modelo se equivoca mucho más en unos sitios que en otros. Aquí es evidente y tiene explicación de negocio: una bodega compra entre 100 y 300 soles y un mayorista entre 500 y 4.000, así que hay muchísimo más espacio para fallar arriba.

La consecuencia práctica es concreta: los intervalos de arriba están mal calculados. Suponen que el error es parejo y no lo es. Se arreglan pidiendo errores robustos, que es una palabra de más en la misma línea:

# los mismos coeficientes, pero con intervalos que aguantan
robusto = smf.ols('monto ~ C(segmento)', data=v).fit(cov_type='HC3')

Los coeficientes no cambian, solo los intervalos. Y aquí los coeficientes están tan lejos del cero que la conclusión aguanta igual, pero eso hay que comprobarlo y no suponerlo 🧯

Los dos errores que sí dan error

El primero no da ninguno, que es peor:

print('con la columna descuento, que trae nulos:')
print(stats.linregress(v['descuento'], v['monto']).slope)

con la columna descuento, que trae nulos:
nan

nan. Ni un aviso. Un solo nulo en la columna y la pendiente entera se vuelve nada, y si eso pasa dentro de un informe automático el número sale en blanco y nadie sabe por qué 🫠

El segundo sí revienta, y menos mal:

iguales = [5.0] * 20
stats.linregress(iguales, y[:20])

ValueError: Cannot calculate a linear regression if all x values are identical

if all x values are identical. Si x no varía, la división de la primera línea que escribimos a mano tiene un cero abajo.

Y esa es la condición de fondo de toda la regresión, dicha por la librería: sin variación en x no hay nada que medir. Si todos tus clientes compran lo mismo, ninguna herramienta te va a decir qué pasa cuando compran distinto.

Cuándo NO usarla

Y una más, que no cabe en la tabla y es la que más veo: meter en el modelo una variable que se calcula con la respuesta. Eso es la trampa de este capítulo 👇

v['precio_unitario'] = v['monto'] / v['unidades']
trampa = smf.ols('monto ~ precio_unitario', data=v).fit()
print(f'R2 = {trampa.rsquared:.4f}')
print(f'valor p = {trampa.pvalues["precio_unitario"]:.3g}')

# R2 = 0.1944   valor p = 5.54e-142

Ejercicios

Seis, y el 5 es el que más enseña. Intenta antes de abrir 💛

al_derecho = stats.linregress(x, y).slope
al_reves = stats.linregress(y, x).slope
print(f'monto sobre unidades : {al_derecho:.4f}')
print(f'unidades sobre monto : {al_reves:.6f}')
print(f'una es la inversa de la otra?: {1 / al_derecho:.6f}')

r = stats.linregress(x, y)
print(f'correlacion al cuadrado : {r.rvalue ** 2:.8f}')
print(f'R cuadrado              : {r.rvalue ** 2:.8f}')
print(f'correlacion             : {r.rvalue:.4f}')

x2 = np.append(x, 500.0)
y2 = np.append(y, 50000.0)
print(f'pendiente con 2.979 filas : {stats.linregress(x, y).slope:.4f}')
print(f'pendiente con una mas     : {stats.linregress(x2, y2).slope:.4f}')

medias = v.groupby('segmento')['monto'].mean()
print(medias.round(2).to_string())
print()
print(f'mayorista menos bodega: {medias["Mayorista"] - medias["Bodega"]:.2f}')

v2 = v.copy()
v2['log_monto'] = np.log(v2['monto'])
m = smf.ols('log_monto ~ C(segmento)', data=v2).fit()
tercio = pd.qcut(m.fittedvalues, 3, labels=['bajo', 'medio', 'alto'])
print(m.resid.groupby(tercio, observed=True).std().round(3).to_string())

print(f'prediccion para 10.000 unidades: {r.intercept + r.slope * 10000:.2f} soles')
print(f'unidades que hay en los datos  : de {x.min():.0f} a {x.max():.0f}')

Lo que te llevas

• 📏 La pendiente es cuánto cambia y cuando x sube en uno, con unidades puestas.
• 🧐 Sale siempre. Lo que decide si significa algo es su intervalo, y el nuestro iba de -3,02 a 6,25.
• 💰 Con variables de categoría, el intercepto es la que falta y el resto se mide contra ella.
• 🔗 Meter más variables sirve para saber qué aporta cada una con las demás fijas, que es lo que quiere decir "controlando por".
• 🔍 Los residuos por tramo son el diagnóstico. Los nuestros pasan de 109 a 676 soles.
• 🚫 Y ninguna regresión demuestra una causa. Ninguna.

Si quieres la versión que predice en vez de la que explica, con validación y métricas, eso es el libro de machine learning desde cero 🤖

Y si el pandas de este capítulo te costó más que la estadística, eso se arregla en el libro de Python desde cero 🐍

Que tengas lindo día! 🌸