C09.4 Regresión lineal (Guía para hacer una buena regresión lineal).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"logo.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pasos para una regresión lineal.\n",
    "\n",
    "## Planteamiento del modelo\n",
    "\n",
    "***\n",
    "**Paso 1. Establecer el modelo**\n",
    "\n",
    "Función|Descripción\n",
    ":--:|:--\n",
    "``resultado_regresion <- lm(nombre_variable_predictora~nombre_variable_1+nombre_variable_2+...+nombre_variable_k,data=tabla)``| Construye el modelo\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Paso 2. Establecer $R^2$.**\n",
    "\n",
    "El valor de $R^2*100$ representa el porcentaje de la variabilidad total que es explicada por el modelo lineal.\n",
    "\n",
    "Función|Descripción|Decisión\n",
    ":--:|:--|:--\n",
    "``summary(resultado_regresion)``|Muestra un resumen de la regresión|\n",
    "``summary(resultado_regresion)$r.squared``|Muestra el $R^2$| Si es $>0.7$ se considera que vale la pena continuar\n",
    "``summary(resultado_regresion)$adj.r.squared``|Muestra el $R^2$ ajustado| A mayor valor, mejor es el modelo\n",
    "\n",
    "***\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Paso 3. Elegir el modelo mas simple (solo en el caso de regresión múltiple).**\n",
    "\n",
    "Basado en el Coeficiente de información de Akaike, se elige el modelo con el menor AIC\n",
    "\n",
    "Función|Descripción|Decisión\n",
    ":--:|:--|:--\n",
    "``step(resultado_regresion)``|Indica el mejor modelo|El nuevo modelo es modelo_simplificado que es el *call* del resultado del ``step``\n",
    "\n",
    "Hacemos ``summary(modelo_simplificado)$r.squared`` para ver la nueva $R^2$. \n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Diagnósticos de regresión\n",
    "***\n",
    "\n",
    "**Paso 4. Verificamos la normalidad de los residuos**\n",
    "\n",
    "Se aplica la prueba de gaussianidad a los residuos. Es decir, a ``modelo_simplificado$residuals``\n",
    "\n",
    "Paqutería|Prueba|Abreviatura|Tamaño de<br> muestra|Desventaja|Estadístico|**R**\n",
    ":--:|:--|:--:|:--:|:--|:--:|:--\n",
    "``fBasics``|Kolmogorov-Smirnov-Lilliefors|KSL|$5\\le n$|No detecta bien las colas|D|``lillieTest()``\n",
    "``fBasics``|Anderson-Darling|AD|$7\\le n\\le1000$|No trabaja bien en muestras grandes|A|``adTest()``\n",
    "``fBasics``|Shapiro-Wilks|SW|En **R** admite $\\le5000$|Sensible a repeticiones|W|``shapiroTest()``\n",
    "``fBasics``|D'Agostino|D|$20\\le n$||Omnibus|``dagoTest()``\n",
    "\n",
    "Regla de decisión:\n",
    "\n",
    "$p$-valor|Decisión|Significado\n",
    ":--|:--|:--\n",
    "Pequeño|Rechazar $H_0$|Hay buena probabilidad de que **NO es gaussiana**\n",
    "Grande|Rechazar $H_1$|Hay buena probabilidad de que **SÍ es gaussiana**\n",
    "\n",
    "***\n",
    "\n",
    "**Paso 5. Homocedasticidad de los residuos: Test de Breuch-Pagan**\n",
    "\n",
    "Paquetería|Función\n",
    ":--:|:--|\n",
    "``lmtest``|``bptest(modelo_simplificado)``\n",
    "\n",
    "Regla de decisión \n",
    "\n",
    "$p$-valor|Decisión|Significado\n",
    ":--|:--|:--\n",
    "Pequeño|Rechazar $H_0$|Hay buena probabilidad de que los residuos **NO sean homocedásticos**\n",
    "Grande|Rechazar $H_1$|Hay buena probabilidad de que los residuos **SÍ sean homocedásticos**\n",
    "\n",
    "**Método gráfico**\n",
    "\n",
    "``ggplot(data = datos, aes(modelo_simplificado$fitted.values, modelo_simplificado$residuals)) +\n",
    "    geom_point() +\n",
    "    geom_hline(yintercept = 0) +\n",
    "    theme_bw()``\n",
    "\n",
    "donde ``datos`` es tu tabla original. Deberemos observar un **cielo estrellado** centrado en la recta $y=0$.\n",
    "\n",
    "***\n",
    "\n",
    "**Paso 6. Residuos no correlacionados (independientes): Test de Durbin-Watson**\n",
    "\n",
    "Paqutería|Función\n",
    ":--:|:--|\n",
    "``car``|``dwt(modelo_simplificado,alternative=\"two.sided\")``\n",
    "\n",
    "Regla de decisión\n",
    "\n",
    "$p$-valor|Decisión|Significado\n",
    ":--|:--|:--\n",
    "Pequeño|Rechazar $H_0$|Hay buena probabilidad de que los residuos **SÍ tengan autocorrelación**\n",
    "Grande|Rechazar $H_1$|Hay buena probabilidad de que los residuos **NO tengan autocorrelación**\n",
    "\n",
    "***\n",
    "\n",
    "**Paso 7. Aditividad del modelo: Test de Tukey**\n",
    "\n",
    "Paquetería|Función|Descripción\n",
    ":--:|:--:|:--\n",
    "``car``|``residualPlots(modelo_simplificado,plot=FALSE)``|Da el valor del estadístico de Tukey y su $p$-valor\n",
    "``car``|``residualPlots(modelo_simplificado,plot=TRUE)``|Devuelve los gráficos de los residuos frente a las variables regresoras y frente a los valores estimados junto con una curva de color azul indicando su tendencia\n",
    "\n",
    "Regla de decisión\n",
    "\n",
    "$p$-valor|Decisión|Significado\n",
    ":--|:--|:--\n",
    "Pequeño|Rechazar $H_0$|Hay buena probabilidad de que el modelo **NO es aditivo**\n",
    "Grande|Rechazar $H_1$|Hay buena probabilidad de que el modelo **SÍ es aditivo**\n",
    "\n",
    "Si optamos por ver los gráficos, entonces para concluir satisfactoriamente todos ellos deben parecer cielos estrellados.\n",
    "\n",
    "***\n",
    "\n",
    "**Paso 8. Linealidad del modelo**\n",
    "\n",
    "Paquetería|Función\n",
    ":--:|:--:\n",
    "``car``|``crPlots(modelo_simplificado)``\n",
    "\n",
    "Los residuos parciales se dibujan contra los valores de  $x_j$  y se hace su recta de regresión.\n",
    "\n",
    "Si esta no se ajusta a la curva dada por una regresión no paramétrica suave (las variables independientes no están predeterminadas y se construyen con los datos), el modelo no es lineal.\n",
    "\n",
    "***\n",
    "\n",
    "**Paso 9. Observaciones anómalas**\n",
    "\n",
    "Sea $n$ el número de renglones de tu tabla y $k$ el número de variables predictoras.\n",
    "\n",
    "* **Leverages.** son observaciones con un valor anómalo de las variables de control. No tienen por qué afectar los coeficientes de la regresión.\n",
    "\n",
    "* **Outliers de regresión** son observaciones que tienen un valor anómalo de la variable $Y$, condicionado a los valores de sus variables independientes $X_i$. Tendrán un residuo muy alto pero no pueden afectar demasiado a los coeficientes de la regresión.\n",
    "\n",
    "* **Observaciones influyentes** son aquellas que tienen un leverage alto; son outliers de regresión y afectan fuertemente a la regresión.\n",
    "\n",
    "Anomalía|Paquetería|Función|Descripción|Uso\n",
    ":--:|:--:|:--:|:--|:--:\n",
    "Leverage (apalancamiento)||``hatvalues(modelo_simplificado)``|Calcula los valores hat de todas las observaciones|``which(hatvalues(modelo_simplificado)>2(k+1)/n)`` \n",
    "Outlier|``car``|``outlierTest(modelo_simplificado)``||Nos da los outliers\n",
    "Influeyente|``car``|``cooks.distance(modelo_simplificado)``|Calcula las distancias de Cook de todas las observaciones|``which(cooks.distance(modelo_simplificado)>4/(n-k-2))``\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Intervalos de confianza y significancia\n",
    "\n",
    "Siempre es buena idea incluir los intervalos de confianza para los coeficientes $\\beta_i$ y su significancia.\n",
    "***\n",
    "***\n",
    "\n",
    "**Paso 10. Intervalos de confianza de los coeficientes $\\beta_i$**\n",
    "\n",
    "Función|Descripción\n",
    ":--|:--:\n",
    "``confint(modelo_simplificado)``|Muestra los intervalos de confianza de los $\\beta$'s\n",
    "***\n",
    "\n",
    "**Paso 11. Significación de los coeficientes $\\beta_i$**\n",
    "\n",
    "Función|Descripción|$p$-valor|Significado\n",
    ":--:|:--|:--:|:--\n",
    "``summary(modelo_simplificado)``|¿Algún coeficiente es diferente de cero?|$<0.05$|Vale la pena hacer la regresión\n",
    "``summary(modelo_simplificado)``|¿Qué coeficientes no son cero?|$<0.05$|El coeficiente respectivo tiene efectos significativos\n",
    "\n",
    "Una vez que tienes los efectos significativos, los enuncias de acuerdo a la siguiente tabla:\n",
    "\n",
    "Símbolo|Significado\n",
    ":--:|:--\n",
    "``*``|la variable $X_i$ tiene efectos significativos\n",
    "``**``|la variable $X_i$ tiene efectos muy significativos\n",
    "``***``|la variable $X_i$ tiene efectos fuertemente significativos\n",
    "``.``|la variable $X_i$ probablemente tiene efectos significativos\n",
    "`` ``|la variable $X_i$ no tiene efectos significativos\n",
    "\n",
    "***\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Predicciones\n",
    "***\n",
    "***\n",
    "\n",
    "**Paso 12.**\n",
    "\n",
    "Supongamos que llegan los nuevos $j$ datos que no fueron parte de las observaciones de nuestra tabla original:\n",
    "\n",
    "Individuo|$X_1$|$X_2$|...|$X_k$\n",
    ":--:|:--:|:--:|:--:|:--:\n",
    "1|$x_1^{(1)}$|$x_1^{(2)}$|...|$x_1^{(k)}$\n",
    "2|$x_2^{(1)}$|$x_2^{(2)}$|...|$x_2^{(k)}$\n",
    "...|...|...|...|...\n",
    "$j$|$x_j^{(1)}$|$x_j^{(2)}$|...|$x_j^{(k)}$\n",
    "\n",
    "Aquí, $x_{r}^{(s)}$ representa el valor de la característica $X_s$ del nuevo individuo $r$. \n",
    "\n",
    "Llama ``newdata`` a la tabla anterior sin la columna individuo. Es decir, la columna 1 son los valores de la característica $X_1$ en los nuevos individuos; la columna 2 son los valores de la característica $X_2$ en los nuevos individuos; etcétera. **Es importante que los nombres de las columnas coincidan con los nombres de las características de la tabla original con la cual hiciste el modelo**.\n",
    "\n",
    "Función|Descripción\n",
    ":--:|:--\n",
    "``predict.lm(modelo_simplificado,newdata,interval = \"confidence\",level = nivel.confianza)``|Te indica el valor promedio y el intervalo de confianza de todos los individuos $r$ cuyos características son $(x_r^{(1)},x_r^{(2)},...,x_r^{(k)})$\n",
    "``predict.lm(modelo_simplificado,newdata,interval = \"prediction\",level = nivel.confianza)``|Te indica el valor calculado por el modelo y el intervalo de confianza de un individuo $r$ cuyos características son $(x_r^{(1)},x_r^{(2)},...,x_r^{(k)})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "***\n",
    "\n",
    "## Ejemplo\n",
    "\n",
    "Se postula que la altura actual de un bebé tiene una relación lineal con su edad actual (en días), su altura al nacer (en cm), su peso al nacer (en Kg)  y el aumento, en tanto por ciento, de su peso actual respecto de su peso al nacer. \n",
    "\n",
    "En una muestra de 9 niños, se obtuvieron los siguientes resultados.\n",
    "\n",
    "|Peso_actual|Edad_actual|Altura_nacer|Peso_nacer|Porcentaje_cambio_peso|\n",
    "|:--:|:--:|:--:|:--:|:--:|\n",
    "|57.5|78|48.2|2.75|29.5|\n",
    "|52.8|69|45.5|2.15|26.3|\n",
    "|61.3|77|46.3|4.41|32.2|\n",
    "|67|88|49|5.52|36.5|\n",
    "|53.5|67|43|3.21|27.2|\n",
    "|62.7|80|48|4.32|27.7|\n",
    "|56.2|74|48|2.31|28.3|\n",
    "|68.5|94|53|4.3|30.3|\n",
    "|69.2|102|58|3.71|28.7|\n",
    "\n",
    "La tabla la encuentras en nuestro repositorio con el nombre de [bebes.csv](https://github.com/scidatmath2020/Inferencia-Estad-stica-2022/blob/main/data/bebes.csv)\n",
    "\n",
    "**¿Vale la pena estudiar un modelo lineal para calcular el peso actual? Si la respuesta es sí, estudiar completamente el modelo lineal**\n",
    "\n",
    "***\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problema final\n",
    "\n",
    "Se probó un modelo de simulación para el flujo máximo de agua de las cuencas hidrográficas comparando el flujo máximo medido de 10 tormentas con predicciones del flujo máximo obtenido del modelo de simulación. $Q_o$ y $Q_p$ son los flujos máximos observados y pronosticados, respectivamente. Se registraron cuatro variables independientes:\n",
    "\n",
    "* **area.** Área de la cuenca (m^2),\n",
    "\n",
    "* **pendiente.** Pendiente promedio de la cuenca (en porcentaje),\n",
    "\n",
    "* **absorbencia.** Índice de absorbencia superficial (0 = absorbencia completa, 100 = sin absorbencia), y\n",
    "\n",
    "* **intensidad.** Intensidad de pico de lluvia calculada en intervalos de media hora.\n",
    "\n",
    "|$Q_o$ | $Q_p$| area| pendiente|absorbencia |intensidad|\n",
    "|:--:|:--:|:--:|:--:|:--:|:--:|\n",
    "|28|              32|              .03|             3.0|             70|              .6|\n",
    "|112  |           142     |        .03  |           3.0    |         80 |             1.8|\n",
    "|398  |           502      |       .13  |           6.5    |         65 |             2.0|\n",
    "|772  |           790   |          1.00   |         15.0  |          60   |           .4|\n",
    "|2294  |          3075|            1.00  |          15.0|            65|              2.3|\n",
    "|2484   |         3230  |          3.00    |        7.0     |        67 |             1.0|\n",
    "|2586    |        3535   |         5.00 |           6.0 |            62|              .9|\n",
    "|3024   |         4265 |           7.00    |        6.5     |        56      |        1.1|\n",
    "|4179   |         6529   |         7.00    |        6.5     |        56|              1.4|\n",
    "|710    |         935  |           7.00   |         6.5   |          56   |           .7|\n",
    "\n",
    "\n",
    "La tabla la encuentras en [flujos_pronosticados.csv](https://github.com/scidatmath2020/Inferencia-Estad-stica-2022/blob/main/data/flujos_pronosticados.csv)\n",
    "\n",
    "\n",
    "Tomemos $Y=\\ln\\left(\\frac{Q_o}{Q_p}\\right)$ como variable dependiente y consideramos la regresión de Y como función de area, pendiente, absorbencia e intensidad. \n",
    "\n",
    "**Realiza un análisis completo de regresión lineal para explicar $Y$ a través de las otras variables**"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}