{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cuando la pregunta no es predecir\n",
    "\n",
    "Coeficientes con su valor p y su intervalo, odds que se dicen en voz alta, y la columna que el boosting adora y la regresión no distingue del ruido.\n",
    "\n",
    "Cuaderno de soluciones del capítulo 22 de **Machine learning desde cero**, de Miss Yera.\n",
    "\n",
    "Corre de arriba abajo. Si lo abres en Google Colab no necesitas instalar nada.\n",
    "\n",
    "Capítulo completo: https://missyera.com/guias/machine-learning-desde-cero/regresion-con-statsmodels/\n",
    "\n",
    "Este es el cuaderno de **soluciones**. Trae el código de cada ejercicio, la\n",
    "explicación de la trampa y la respuesta del quiz. Si vienes del cuaderno de\n",
    "práctica sin haberlo intentado, vuelve 🙂"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Todo el libro hasta aquí ha ido de **predecir**: dame una fila y\n",
    "te digo la probabilidad. Este capítulo va de la otra pregunta, la que también\n",
    "llega y que se contesta con otras herramientas 🔬\n",
    "\n",
    "Y te la hago a ti primero: **¿cuántas veces has contestado \"sí, influye\" sin poder decir cuánto ni con qué margen?** Este capítulo es para no tener que volver a hacerlo 🔬\n",
    "\n",
    "\"¿Ser Mayorista influye de verdad en que cierre la venta, o es que los\n",
    "mayoristas casualmente compran por WhatsApp y lo que influye es el canal?\"\n",
    "\n",
    "Eso no es una predicción. Es una pregunta sobre **si el efecto\n",
    "existe**, y scikit-learn no está hecho para contestarla. Le pides los\n",
    "coeficientes y te da trece números pelados, sin errores estándar, sin\n",
    "intervalos y sin valores p.\n",
    "\n",
    "Para eso está statsmodels, que es la librería que usaría alguien de\n",
    "estadística y que en este mundo se nombra poco 💜"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Los mismos coeficientes, con lo que les falta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Empecemos por lo importante: **no es otro modelo**. Es la misma\n",
    "regresión logística del capítulo 13, ajustada igual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "URL = 'https://missyera.com/static/datasets/ventas-miss-yera.csv'\n",
    "\n",
    "def carga_limpia(url):\n",
    "    v = pd.read_csv(url).drop_duplicates()\n",
    "    v['ciudad'] = (v['ciudad'].str.strip().str.lower()\n",
    "                   .str.normalize('NFKD')\n",
    "                   .str.encode('ascii', 'ignore').str.decode('utf-8'))\n",
    "    v['monto'] = pd.to_numeric(v['monto'].str.replace(',', '.'))\n",
    "    for col in ['fecha', 'fecha_ultima_compra']:\n",
    "        f = pd.to_datetime(v[col], format='%Y-%m-%d', errors='coerce')\n",
    "        falta = f.isna() & v[col].notna()\n",
    "        f[falta] = pd.to_datetime(v.loc[falta, col], format='%d/%m/%Y', errors='coerce')\n",
    "        v[col] = f\n",
    "    return v\n",
    "\n",
    "ventas = carga_limpia(URL)\n",
    "ventas['sin_compra_previa'] = ventas['fecha_ultima_compra'].isna().astype(int)\n",
    "ventas['precio_unitario'] = ventas['monto'] / ventas['unidades']\n",
    "datos = ventas.dropna(subset=['satisfaccion']).copy()\n",
    "\n",
    "FORMULA = ('compro ~ satisfaccion + monto + sin_compra_previa '\n",
    "           '+ C(segmento) + C(canal)')\n",
    "modelo = smf.logit(FORMULA, data=datos).fit(disp=0)\n",
    "\n",
    "print('statsmodels', sm.__version__)\n",
    "print('filas usadas:', int(modelo.nobs))\n",
    "print(modelo.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esa tabla es lo que statsmodels añade, y son cuatro columnas 🔎\n",
    "\n",
    "- 📏 **`std err`**: cuánto se movería ese coeficiente\n",
    "si repitieras el estudio con otra muestra. Es la barra de error.\n",
    "\n",
    "- ➗ **`z`**: el coeficiente dividido entre esa barra,\n",
    "o sea cuántas veces cabe el ruido dentro del efecto.\n",
    "\n",
    "- 🎲 **`P>|z|`**: la probabilidad de ver algo así\n",
    "de grande si el efecto real fuera cero.\n",
    "\n",
    "- 📐 **`[0.025 0.975]`**: el intervalo donde\n",
    "razonablemente está el valor de verdad."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "z=β^ee(β^),IC95%=exp(β^±1,96ee)\n",
    "\n",
    "el coeficiente dividido entre su error estándar dice cuántas veces cabe el ruido dentro del efecto, y el intervalo se lee en odds al pasarlo por la exponencial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si esas palabras te suenan a chino, están todas explicadas de cero en el\n",
    "[libro de estadística](https://missyera.com/guias/estadistica-desde-cero/), que es el que\n",
    "va antes que este 📐"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lo mismo con sklearn, para que no haya duda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "X = pd.get_dummies(\n",
    "    datos[['satisfaccion', 'monto', 'sin_compra_previa', 'segmento', 'canal']],\n",
    "    drop_first=True).astype(float)\n",
    "sk = LogisticRegression(max_iter=5000, C=1e9).fit(X, datos['compro'])\n",
    "\n",
    "print(pd.Series(sk.coef_[0], index=X.columns).round(4).to_string())\n",
    "print()\n",
    "print('satisfaccion  statsmodels', round(modelo.params['satisfaccion'], 4),\n",
    "      '| sklearn', round(float(sk.coef_[0][list(X.columns).index('satisfaccion')]), 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Idénticos hasta la cuarta cifra 🤝\n",
    "\n",
    "Ese `C=1e9` es lo único que hay que saber para que coincidan:\n",
    "scikit-learn regulariza por defecto y statsmodels no, así que hay que apagar la\n",
    "regularización poniendo la `C` altísima.\n",
    "\n",
    "Y ojo, que ese detalle explica una cosa que confunde a mucha gente: si sacas\n",
    "coeficientes de scikit-learn sin tocar la `C`, están encogidos hacia\n",
    "cero a propósito. Sirven para predecir. Para contarlos en una reunión, no."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Los odds, que es como se dice en voz alta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un coeficiente de 1,2799 no se puede decir en una reunión. Pasado por la\n",
    "exponencial, sí."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "log(p1−p)=b+∑jwjxj\n",
    "\n",
    "la misma fórmula al revés: cada peso dice cuánto mueve el logaritmo de la ventaja, y por eso un peso de 0,7 se lee como multiplicar la ventaja por dos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tabla = pd.DataFrame({'odds': np.exp(modelo.params), 'p': modelo.pvalues})\n",
    "intervalo = np.exp(modelo.conf_int())\n",
    "tabla['ic_bajo'], tabla['ic_alto'] = intervalo[0], intervalo[1]\n",
    "print(tabla.round(4).to_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora sí se puede leer en voz alta 🗣️\n",
    "\n",
    "- 🏭 Ser **Mayorista multiplica por 3,60** las probabilidades a\n",
    "favor de que cierre, comparado con Bodega. Y el intervalo va de 2,33 a 5,55, así\n",
    "que aunque el número exacto no lo sepamos, que multiplica al menos por dos sí.\n",
    "\n",
    "- ⭐ Cada punto de **satisfacción multiplica por 1,33**.\n",
    "\n",
    "- 🆕 Ser **cliente nuevo multiplica por 0,41**, o sea que las\n",
    "reduce a menos de la mitad.\n",
    "\n",
    "- 📱 **WhatsApp multiplica por 2,34** contra Marketplace, que es\n",
    "la categoría de referencia.\n",
    "\n",
    "Ese \"comparado con Bodega\" y ese \"contra Marketplace\" no son un detalle: en\n",
    "una variable categórica **todos los coeficientes se leen contra la\n",
    "categoría que falta**, que es la primera por orden alfabético. Si te\n",
    "olvidas de decirlo, la frase queda mal 🏷️"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## La que no se distingue de nada"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hay una fila de esa tabla que no se parece a las demás."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cruza = tabla[(tabla['ic_bajo'] < 1) & (tabla['ic_alto'] > 1)]\n",
    "print(cruza.round(4).to_string() if len(cruza) else 'ninguna')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`monto`, con p de 0,1561 y un intervalo que va de 0,9999 a 1,0004,\n",
    "o sea que cruza el 1 😐\n",
    "\n",
    "Cruzar el 1 en odds es lo mismo que cruzar el cero en el coeficiente: con\n",
    "estos datos **no se puede distinguir de no tener efecto**. Y fíjate\n",
    "en cómo lo he dicho, porque la frase importa: no es \"el monto no influye\". Es\n",
    "\"con estas 2.769 filas no lo puedo separar del ruido\".\n",
    "\n",
    "Y ahora el detalle que a mí me parece lo mejor de todo el libro 🤯\n",
    "\n",
    "En el capítulo 21, `monto` era la **columna\n",
    "número uno** del LightGBM, con 0,2342 de SHAP. Aquí no llega a\n",
    "distinguirse de cero.\n",
    "\n",
    "No se contradicen. Dicen cosas distintas sobre la misma columna:\n",
    "\n",
    "- 📈 La regresión pregunta si el monto empuja **en línea recta**:\n",
    "más soles, más probabilidad. Y en línea recta no dice nada.\n",
    "\n",
    "- 🪜 El árbol puede trocear el monto donde quiera y quedarse con \"por debajo\n",
    "de 300 soles pasa una cosa y por encima de 1.500 otra\". Eso sí lo encuentra.\n",
    "\n",
    "La lección práctica: **un valor p alto no significa que la columna sea\n",
    "inútil para predecir**. Significa que su efecto no es una recta. Tirar\n",
    "columnas por su valor p es de los errores que más modelos empeora 🚮"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lo que de verdad se dice en la reunión"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los odds son correctos y siguen sin ser lo que la gente entiende. Nadie\n",
    "piensa en \"multiplicar las probabilidades a favor\". La gente piensa en puntos\n",
    "porcentuales."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "marginal = modelo.get_margeff(at='overall')\n",
    "print(marginal.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esto ya es castellano 🎯\n",
    "\n",
    "- ⭐ Cada punto de satisfacción sube la probabilidad de cierre\n",
    "**6,05 puntos**.\n",
    "\n",
    "- 🏭 Ser Mayorista en vez de Bodega la sube **27,00 puntos**.\n",
    "\n",
    "- 🆕 Ser cliente nuevo la baja **18,99 puntos**.\n",
    "\n",
    "El efecto marginal es la pendiente media: cuánto se mueve la probabilidad si\n",
    "esa variable sube una unidad y todo lo demás se queda quieto.\n",
    "\n",
    "Y digo \"media\" a propósito, porque en una logística la pendiente no es la\n",
    "misma en todos lados: donde la probabilidad ya es 0,95 casi no se puede subir\n",
    "más. El número de arriba es el promedio de todas las filas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Más datos afilan lo real y no arreglan lo que no está"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La demostración que enseña qué es de verdad un valor p 🔬"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in (500, 1000, 2000, len(datos)):\n",
    "    trozo = datos.sample(n=min(n, len(datos)), random_state=7)\n",
    "    m = smf.logit(FORMULA, data=trozo).fit(disp=0)\n",
    "    print(f'{len(trozo):5d} filas -> monto coef {m.params[\"monto\"]:+.6f} '\n",
    "          f'p {m.pvalues[\"monto\"]:.4f}   satisfaccion p {m.pvalues[\"satisfaccion\"]:.2e}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mira las dos columnas por separado, que cuentan historias opuestas.\n",
    "\n",
    "**La satisfacción** va de p 5,44e-06 a 4,12e-22. Cada vez más\n",
    "seguro. Eso es lo que hace un efecto que existe: con más datos la barra de error\n",
    "se encoge y el efecto se ve más claro.\n",
    "\n",
    "**El monto** pasea. Su coeficiente empieza en −0,000358, se va a\n",
    "−0,000017, sube a +0,000120 y acaba en +0,000159. **Cambia de\n",
    "signo** por el camino. Y su p va de 0,21 a 0,93 a 0,36 a 0,16, sin\n",
    "acercarse nunca a nada.\n",
    "\n",
    "Un efecto que existe se afila con datos. Uno que no existe se pasea, y ese\n",
    "paseo es exactamente lo que mide el valor p 🎲"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## El error: cuando el modelo se niega a ajustar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos a meterle la trampa del capítulo 12: una\n",
    "columna que *es* la respuesta."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Esto revienta a propósito.** Se ejecuta dentro de un `try` para que puedas seguir con \"ejecutar todo\" y aun así ver la queja."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    con_soplon = datos.copy()\n",
    "    con_soplon['soplon'] = con_soplon['compro']\n",
    "    smf.logit('compro ~ soplon + satisfaccion', data=con_soplon).fit(disp=0)\n",
    "except Exception as e:\n",
    "    print(f'{type(e).__name__}: {e}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Y la queja que tiene que salir es esta:\n",
    "\n",
    "```\n",
    "LinAlgError: Singular matrix\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esto es precioso y no lo esperaba la primera vez 🤩\n",
    "\n",
    "Se llama **separación perfecta**: existe una columna que separa\n",
    "las dos clases sin errores, así que el coeficiente que mejor ajusta es infinito.\n",
    "Y como no hay número que sirva, el álgebra se rompe y sale un\n",
    "`LinAlgError`.\n",
    "\n",
    "Compara eso con lo que pasó en el capítulo de fuga: scikit-learn cogió la\n",
    "columna soplona, entrenó tan contento y devolvió 0,9994 de AUC. Nadie avisó de\n",
    "nada.\n",
    "\n",
    "Statsmodels **se niega a ajustar**. Y aunque sea molesto, es la\n",
    "mejor detección de fuga que hay: si tu regresión no converge, mira si alguna\n",
    "columna sabe demasiado antes de tocar nada más 🚨"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Y el silencio que sí duele"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solo_satisfaccion = smf.logit('compro ~ satisfaccion', data=ventas).fit(disp=0)\n",
    "print('filas del archivo :', len(ventas))\n",
    "print('filas que usó     :', int(solo_satisfaccion.nobs))\n",
    "print('descartadas       :', len(ventas) - int(solo_satisfaccion.nobs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "231 filas fuera y ni un aviso 😨\n",
    "\n",
    "Statsmodels borra las filas con nulos y sigue. Le pasé el archivo entero y\n",
    "ajustó con 2.769, o sea que el 7,7% de las ventas no participó en las\n",
    "conclusiones y en ningún sitio del `summary()` lo dice a gritos.\n",
    "\n",
    "Por eso en este capítulo el `dropna` va escrito arriba, a mano.\n",
    "Que la decisión de qué filas entran la tomes tú y no la librería, y que quede\n",
    "por escrito cuántas eran 📋\n",
    "\n",
    "Y acuérdate del capítulo 9: esas filas no son un\n",
    "descuido, son la columna `sin_satisfaccion`. Si el hueco significa\n",
    "algo, descartar las filas te lo tira a la basura."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ejercicios"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. ¿Se pisan las columnas entre ellas?\n",
    "\n",
    "Si dos columnas dicen casi lo mismo, sus coeficientes se\n",
    "vuelven locos y sus intervalos se ensanchan. Eso se mide con el VIF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "\n",
    "M = datos[['satisfaccion', 'monto', 'unidades', 'precio_unitario',\n",
    "           'sin_compra_previa']].astype(float)\n",
    "M = sm.add_constant(M)\n",
    "vif = pd.Series([variance_inflation_factor(M.values, i) for i in range(M.shape[1])],\n",
    "                index=M.columns)\n",
    "print(vif.round(2).to_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "const                12.05\n",
    "satisfaccion          1.00\n",
    "monto                 1.33\n",
    "unidades              1.28\n",
    "precio_unitario       1.60\n",
    "sin_compra_previa     1.00\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La regla de andar por casa es que por encima de 5 hay que mirarlo y por\n",
    "encima de 10 hay problema. Aquí lo más alto es 1,60 🎉\n",
    "\n",
    "Y eso que `precio_unitario` se calcula dividiendo\n",
    "`monto` entre `unidades`, que suena a estar contando lo\n",
    "mismo tres veces. Pues no: una división no es una combinación lineal, y el VIF\n",
    "solo ve relaciones lineales.\n",
    "\n",
    "La constante sale en 12,05 y esa no cuenta. Su VIF alto solo dice que las\n",
    "columnas no están centradas en cero, que es lo normal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Meter las tres y ver si algo se mueve\n",
    "\n",
    "El VIF dice que no hay problema. Compruébalo de la forma\n",
    "que de verdad importa: añadiendo las columnas y mirando si los demás\n",
    "coeficientes se mueven."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ampliado = smf.logit(FORMULA + ' + unidades + precio_unitario',\n",
    "                     data=datos).fit(disp=0)\n",
    "comparar = pd.DataFrame({'sin': modelo.params, 'con': ampliado.params,\n",
    "                         'p_sin': modelo.pvalues, 'p_con': ampliado.pvalues})\n",
    "print(comparar.round(4).to_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "                              sin     con   p_sin   p_con\n",
    "C(canal)[T.Tienda]         0.6429  0.6414  0.0000  0.0000\n",
    "C(canal)[T.Web]            0.4781  0.4781  0.0000  0.0000\n",
    "C(canal)[T.WhatsApp]       0.8502  0.8470  0.0000  0.0000\n",
    "C(segmento)[T.Horeca]      1.0444  1.0438  0.0000  0.0000\n",
    "C(segmento)[T.Mayorista]   1.2799  1.2829  0.0000  0.0000\n",
    "C(segmento)[T.Minimarket]  0.7989  0.7931  0.0000  0.0000\n",
    "Intercept                 -1.7627 -1.8732  0.0000  0.0000\n",
    "monto                      0.0002  0.0001  0.1561  0.4575\n",
    "precio_unitario               NaN  0.0005     NaN  0.0618\n",
    "satisfaccion               0.2868  0.2881  0.0000  0.0000\n",
    "sin_compra_previa         -0.9003 -0.8961  0.0000  0.0000\n",
    "unidades                      NaN  0.0094     NaN  0.2492\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los coeficientes que ya estaban casi no se enteran: Tienda pasa de 0,6429 a\n",
    "0,6414 y la satisfacción de 0,2868 a 0,2881. Eso es lo que se espera cuando no\n",
    "hay colinealidad 👍\n",
    "\n",
    "Pero mira `monto`: su p pasa de 0,1561 a 0,4575 al entrar\n",
    "`precio_unitario`. Las dos comparten información y se reparten el\n",
    "poco efecto que había, así que las dos quedan más borrosas.\n",
    "\n",
    "Es la versión suave del problema: no rompe nada, pero si tu pregunta era\n",
    "justo sobre el monto, acabas de contestarla peor por añadir una columna 🤷‍♀️"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. La categoría de referencia, cambiada a propósito\n",
    "\n",
    "Todos los coeficientes se leen contra Bodega porque es la\n",
    "primera del abecedario. Ponla contra Mayorista y mira qué cambia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "otra = smf.logit(\n",
    "    'compro ~ satisfaccion + monto + sin_compra_previa '\n",
    "    \"+ C(segmento, Treatment(reference='Mayorista')) + C(canal)\",\n",
    "    data=datos).fit(disp=0)\n",
    "\n",
    "print('contra Bodega:')\n",
    "print(np.exp(modelo.params.filter(like='segmento')).round(4).to_string())\n",
    "print('\\ncontra Mayorista:')\n",
    "print(np.exp(otra.params.filter(like='segmento')).round(4).to_string())\n",
    "print('\\nsatisfaccion en los dos:', round(modelo.params['satisfaccion'], 4),\n",
    "      'y', round(otra.params['satisfaccion'], 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "contra Bodega:\n",
    "C(segmento)[T.Horeca]        2.8417\n",
    "C(segmento)[T.Mayorista]     3.5961\n",
    "C(segmento)[T.Minimarket]    2.2232\n",
    "\n",
    "contra Mayorista:\n",
    "C(segmento, Treatment(reference='Mayorista'))[T.Bodega]        0.2781\n",
    "C(segmento, Treatment(reference='Mayorista'))[T.Horeca]        0.7902\n",
    "C(segmento, Treatment(reference='Mayorista'))[T.Minimarket]    0.6182\n",
    "\n",
    "satisfaccion en los dos: 0.2868 y 0.2868\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los odds del segmento cambian todos y los de la satisfacción no se mueven ni\n",
    "una milésima. Es el mismo modelo contado desde otro sitio 🔄\n",
    "\n",
    "Elige la referencia que haga la frase útil. Si la reunión va de si vale la\n",
    "pena ir a por bodegas, la referencia es Bodega. Si va de si conviene dejar de\n",
    "atenderlas, ponla contra el segmento fuerte y que se vea la distancia."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Un modelo dentro de otro\n",
    "\n",
    "¿Aporta algo el canal, sabiendo ya el segmento? Se compara\n",
    "el modelo con y sin él, con una prueba de razón de verosimilitud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "sin_canal = smf.logit('compro ~ satisfaccion + monto + sin_compra_previa '\n",
    "                      '+ C(segmento)', data=datos).fit(disp=0)\n",
    "\n",
    "estadistico = 2 * (modelo.llf - sin_canal.llf)\n",
    "grados = int(modelo.df_model - sin_canal.df_model)\n",
    "p = stats.chi2.sf(estadistico, grados)\n",
    "\n",
    "print('log-verosimilitud sin canal:', round(sin_canal.llf, 2))\n",
    "print('log-verosimilitud con canal:', round(modelo.llf, 2))\n",
    "print(f'chi2 = {estadistico:.2f} con {grados} grados de libertad, p = {p:.3e}')\n",
    "print('pseudo R2:', round(sin_canal.prsquared, 4), '->', round(modelo.prsquared, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "log-verosimilitud sin canal: -1718.59\n",
    "log-verosimilitud con canal: -1689.7\n",
    "chi2 = 57.77 con 3 grados de libertad, p = 1.758e-12\n",
    "pseudo R2: 0.0901 -> 0.1054\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esta es la forma correcta de preguntar \"¿este bloque de columnas aporta?\",\n",
    "y sirve para cualquier grupo: las tres columnas del canal a la vez, no una por\n",
    "una 🧱\n",
    "\n",
    "Es lo mismo que hace la validación cruzada del capítulo\n",
    "16 pero contestando otra cosa: allí se pregunta si\n",
    "predice mejor, aquí si el ajuste mejora más de lo que mejoraría por azar al\n",
    "añadir parámetros."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Los residuos, para encontrar filas raras\n",
    "\n",
    "Un residuo grande es una fila que el modelo no explica.\n",
    "Míralas, que suelen contar algo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "residuos = modelo.resid_pearson\n",
    "extremos = datos.assign(residuo=residuos.values,\n",
    "                        ajustado=modelo.fittedvalues.values)\n",
    "extremos = extremos.reindex(extremos['residuo'].abs().sort_values(ascending=False).index)\n",
    "\n",
    "print(extremos[['segmento', 'canal', 'satisfaccion', 'monto', 'compro', 'residuo']]\n",
    "      .head(6).round(3).to_string(index=False))\n",
    "print()\n",
    "print('residuos por encima de |3|:', int((np.abs(residuos) > 3).sum()), 'de', len(residuos))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    " segmento    canal  satisfaccion   monto  compro  residuo\n",
    "Mayorista WhatsApp           5.0 3440.83       0   -3.233\n",
    "Mayorista   Tienda           5.0 1991.81       0   -2.599\n",
    "Mayorista   Tienda           5.0 1980.53       0   -2.596\n",
    "Mayorista   Tienda           5.0 1796.06       0   -2.559\n",
    "Mayorista WhatsApp           4.0 2135.48       0   -2.526\n",
    "Mayorista   Tienda           5.0 1588.35       0   -2.517\n",
    "\n",
    "residuos por encima de |3|: 1 de 2769\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Las seis de arriba son la misma película: **mayoristas con satisfacción\n",
    "5 y montos altos que no compraron**. Y solo una pasa de |3| en las 2.769,\n",
    "así que el modelo no está roto: son casos que el archivo no sabe explicar 🕵️‍♀️\n",
    "\n",
    "Ahora fíjate en la cuarta fila, la de 1.796,06 soles. Es **exactamente\n",
    "la misma** que salió en el capítulo 21 como la predicción\n",
    "más segura y equivocada del LightGBM.\n",
    "\n",
    "Dos modelos distintos, dos formas distintas de buscar y la misma venta\n",
    "arriba. Eso ya no es casualidad: es que a esa venta le falta una columna que\n",
    "nadie está midiendo, y ninguna herramienta la va a inventar 🔍"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Predecir con statsmodels, y por qué no lo harías\n",
    "\n",
    "Se puede, claro. Mide su AUC contra la logística de\n",
    "scikit-learn del resto del libro."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "entrena, examen = train_test_split(datos, test_size=0.25, random_state=42,\n",
    "                                   stratify=datos['compro'])\n",
    "m = smf.logit(FORMULA, data=entrena).fit(disp=0)\n",
    "p = m.predict(examen)\n",
    "\n",
    "print('AUC de statsmodels:', round(roc_auc_score(examen['compro'], p), 4))\n",
    "print('filas de examen   :', len(examen))\n",
    "print()\n",
    "print('lo que NO trae de fabrica:')\n",
    "for cosa in ('pipeline', 'validacion cruzada', 'imputacion',\n",
    "             'one-hot con handle_unknown', 'predict_proba'):\n",
    "    print('  -', cosa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "AUC de statsmodels: 0.7343\n",
    "filas de examen   : 693\n",
    "\n",
    "lo que NO trae de fabrica:\n",
    "  - pipeline\n",
    "  - validacion cruzada\n",
    "  - imputacion\n",
    "  - one-hot con handle_unknown\n",
    "  - predict_proba\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El AUC es del mismo orden que el del resto del libro, porque es el mismo\n",
    "modelo. Lo que no trae es toda la maquinaria de producción 🏭\n",
    "\n",
    "Y esa es la regla con la que yo trabajo: **statsmodels para entender y\n",
    "para el informe, scikit-learn para el pipeline que se despliega**. No\n",
    "compiten. Se usan en momentos distintos del mismo proyecto."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. La misma pregunta sin regresión, para comprobar\n",
    "\n",
    "Antes de creerte un coeficiente, mira el número crudo. Si\n",
    "no se parecen, algo pasa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cruda = datos.groupby('segmento')['compro'].agg(['mean', 'count'])\n",
    "cruda['odds'] = cruda['mean'] / (1 - cruda['mean'])\n",
    "cruda['odds_vs_bodega'] = cruda['odds'] / cruda.loc['Bodega', 'odds']\n",
    "\n",
    "print(cruda.round(4).to_string())\n",
    "print()\n",
    "print('lo que dice el modelo, descontando canal, satisfaccion y lo demas:')\n",
    "print(np.exp(modelo.params.filter(like='segmento')).round(4).to_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "              mean  count    odds  odds_vs_bodega\n",
    "segmento                                         \n",
    "Bodega      0.3740    655  0.5976          1.0000\n",
    "Horeca      0.6315    692  1.7137          2.8679\n",
    "Mayorista   0.7149    684  2.5077          4.1965\n",
    "Minimarket  0.5678    738  1.3135          2.1981\n",
    "\n",
    "lo que dice el modelo, descontando canal, satisfaccion y lo demas:\n",
    "C(segmento)[T.Horeca]        2.8417\n",
    "C(segmento)[T.Mayorista]     3.5961\n",
    "C(segmento)[T.Minimarket]    2.2232\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los dos números contestan preguntas distintas y por eso no coinciden 📊\n",
    "\n",
    "La tabla cruda dice cuánto cierran los mayoristas *tal como están*,\n",
    "con su canal y su satisfacción incluidos. El coeficiente dice cuánto cierran de\n",
    "más **a igualdad de todo lo demás**.\n",
    "\n",
    "Y sale más chico: en crudo el mayorista multiplica por 4,20 y descontando lo\n",
    "demás por 3,60. Una sexta parte de su ventaja **no era suya**, era\n",
    "del canal por el que compra y de la satisfacción que reporta.\n",
    "\n",
    "Eso cambia la decisión, y no poco: significa que empujar el canal puede\n",
    "servirle también a los demás segmentos, y que \"los mayoristas cierran más\" es\n",
    "verdad pero no por lo que parecía 💡\n",
    "\n",
    "Cuando el coeficiente y el número crudo se separan mucho más que esto, la\n",
    "señal es más seria: algo de lo que estás midiendo se explica por otra cosa.\n",
    "Ahí es donde empieza la conversación sobre causalidad, que no es este libro."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comprueba que lo tienes\n",
    "\n",
    "Un coeficiente sale con p de 0,16 y su intervalo de odds va de 0,9999 a 1,0004. ¿Qué escribes en el informe?\n",
    "\n",
    "a) Que con estos datos no se distingue de no tener efecto\n",
    "\n",
    "b) Que esa variable no influye en la compra\n",
    "\n",
    "c) Que influye poco, porque el odds es casi 1\n",
    "\n",
    "d) Que hace falta quitarla del modelo\n",
    "\n",
    "---\n",
    "\n",
    "**La correcta es la a.**\n",
    "\n",
    "*b)* No poder distinguirlo de cero no es lo mismo que demostrar que es cero.\n",
    "\n",
    "*c)* El tamaño que mides no significa nada si el intervalo cruza el 1.\n",
    "\n",
    "*d)* Quitarla es una decisión de modelado y esta es una conclusión sobre los datos. No son lo mismo.\n",
    "\n",
    "Un intervalo que cruza el 1 dice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Y el resumen que sale demasiado bonito"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### La trampa\n",
    "\n",
    "Depuración de variables de manual: se quita la menos significativa, se vuelve a ajustar, y así hasta que todas bajan de 0,05. Al final el resumen sale precioso.\n",
    "\n",
    "```\n",
    "while True:\n",
    "    m = sm.Logit(y, X).fit(disp=0)\n",
    "    peor = m.pvalues.idxmax()\n",
    "    if m.pvalues[peor] < 0.05:\n",
    "        break\n",
    "    X = X.drop(columns=[peor])\n",
    "```\n",
    "\n",
    "**Qué está mal**\n",
    "\n",
    "Los valores p del último modelo **ya no significan lo que dicen**. Un valor p contesta \"si esta columna no sirviera para nada, qué probabilidad habría de ver algo así de fuerte por azar\", y esa pregunta supone que miraste el modelo una vez. Aquí lo miraste veinte, quedándote cada vez con lo que mejor se veía 🎣\n",
    "\n",
    "Lo que sale al final es la combinación de columnas que mejor le va a estas filas concretas, con unos p calculados como si nada de eso hubiera pasado. Todos salen por debajo de 0,05 **por construcción**: es el criterio con el que paraste.\n",
    "\n",
    "Si el objetivo es decidir, las columnas se eligen **antes**, por criterio de negocio, y se reporta el modelo entero con los p que salgan, feos incluidos. Un coeficiente que no se distingue de cero también es un resultado, y de los útiles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lo que te llevas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 🧮 statsmodels no es otro modelo: es la misma regresión logística con lo\n",
    "que le falta a scikit-learn. Error estándar, z, valor p e intervalo.\n",
    "\n",
    "- 🔧 Para que los coeficientes coincidan hay que apagar la regularización de\n",
    "scikit-learn con `C=1e9`. Sin eso están encogidos a propósito y no\n",
    "se cuentan en una reunión.\n",
    "\n",
    "- 🗣️ Los odds se dicen en voz alta: Mayorista multiplica por 3,60 (entre 2,33\n",
    "y 5,55), cada punto de satisfacción por 1,33 y ser cliente nuevo por 0,41.\n",
    "\n",
    "- 🏷️ Todo coeficiente categórico se lee contra la categoría que falta, que es\n",
    "la primera del abecedario. Sin decirlo, la frase queda mal.\n",
    "\n",
    "- 🎯 El efecto marginal es lo que de verdad entiende la gente: cada punto de\n",
    "satisfacción sube 6,05 puntos de probabilidad y ser cliente nuevo baja 18,99.\n",
    "\n",
    "- 😐 `monto` tiene p de 0,1561 y su intervalo cruza el 1. Eso es\n",
    "\"con estos datos no lo distingo del ruido\", no es \"no influye\".\n",
    "\n",
    "- 🤯 Y esa misma columna era la número uno del boosting en el capítulo\n",
    "21. No se contradicen: la regresión pregunta si empuja en línea\n",
    "recta y el árbol puede trocearla. Tirar columnas por su valor p empeora\n",
    "modelos.\n",
    "\n",
    "- 🔬 Con más filas, la satisfacción pasa de p 5,44e-06 a 4,12e-22 y el monto\n",
    "se pasea hasta cambiar de signo. Un efecto real se afila con datos; uno que no\n",
    "está, no.\n",
    "\n",
    "- 🚨 Una columna que es la respuesta hace que la regresión ni ajuste:\n",
    "separación perfecta y `LinAlgError`. Es la mejor detección de fuga\n",
    "que hay, y scikit-learn en el mismo caso devolvía 0,9994 de AUC tan contento.\n",
    "\n",
    "- 🕳️ Y el silencio que sí duele: statsmodels borra las filas con nulos sin\n",
    "avisar. Le pasé 3.000 y ajustó con 2.769.\n",
    "\n",
    "Y si de todo el capítulo te llevas una sola frase, que sea esta:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un coeficiente sin intervalo no es un hallazgo. Es una cifra con suerte."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Intervalos, error estándar y qué contesta de verdad un valor p: todo eso viene del [libro de estadística desde cero](https://missyera.com/guias/estadistica-desde-cero/) 📐"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Qué viene ahora"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ya sabemos predecir, explicar el modelo entero, explicar una fila y decir si\n",
    "un efecto existe. Falta la pregunta más incómoda de todas 😬\n",
    "\n",
    "El capítulo 24 abre el modelo por grupos: ¿acierta\n",
    "igual con las bodegas que con los mayoristas? Y después, cuánta confianza merece\n",
    "cada predicción por separado.\n",
    "\n",
    "Lo que verás allí:\n",
    "\n",
    "- ⚖️ Las métricas del capítulo 14 calculadas dentro de cada\n",
    "segmento, que es donde salen las diferencias que el promedio esconde.\n",
    "\n",
    "- 🎲 Cuánto de esa diferencia es real y cuánto es la lotería de la partición,\n",
    "medido con diez sorteos.\n",
    "\n",
    "- 📏 Calibración: si el modelo dice 70%, ¿de cada diez cierran siete?\n",
    "\n",
    "- 🛡️ Y la predicción conforme, que en vez de un número da un conjunto y\n",
    "promete cubrir la respuesta el 90% de las veces."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Ese era el capítulo 22 de **Machine learning desde cero**. El texto completo, con las salidas de cada bloque, está en https://missyera.com/guias/machine-learning-desde-cero/regresion-con-statsmodels/\n",
    "\n",
    "Que tengas lindo día! 🌸"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}