lunes, 1 de mayo de 2017

Aventura de Regresión Lineal y el método de eliminación hacia atrás paso a paso

En esta entrada, haremos el ejercicio de regresión lineal planteado en el semestre 2012-1, cuyo trabajo práctico es igual al de este semestre (2017-1). Se hace uso del R y mi librería estUNA para construir modelos de regresión lineal mediante el método de eliminación hacia atrás. En este proceso, enfatizamos la importancia de realizar un análisis de residuos, entre otras cosas para sugerirnos posibles transformaciones de las variables con las que podamos mejorar los modelos de regresión.

En el semestre 2012-1, los trabajos de estadística se basaron en una data correspondiente a 150 observaciones de algunas variables relacionadas al tema de la nutrición infantil. Específicamente, la muestra de 150 niños con edades comprendida entre 5 y 13 años recopilaba información sobre las siguientes variables:

  • X1 Índice de masa corporal o IMC, calculado según la siguiente fórmula: IMC=peso/estatura^2.
  • X2 Edad en años.
  • X3 Edad ósea en años.
  • X4 Sexo.
  • X5 Peso en kilogramos.
  • X6 Estatura en metros.
  • X7 Índice de actividad física.
  • X8 Calorías consumidas por día en kcal.
  • X9 Nivel de colesterol en la sangre medido en mg/dl.
  • X10 Nivel de glicemia en la sangre medido en mg/dl..

Para evaluar los objetivos referentes a la regresión lineal múltiple de las asignaturas 746 y 738/748, se planteaban dos modelos lineales con todas las variables que se debían depurar mediante el método de eliminación paso a paso hasta obtener modelos con todas sus variables significativas. El primer modelo consideraba la variable X1 como variable independiente, según la siguiente fórmula:

\[X_1=b_0+b_1\cdot X_2+ b_2\cdot X_4+b_3\cdot X_8+b_4\cdot X_9+b_5\cdot X_{10}+\varepsilon\]

El enunciado pedía aplicar el procedimiento de eliminación de variables paso a paso (eliminación hacia atrás) para obtener el "mejor" modelo en cada caso. Seguidamente, había que realizar un análisis de residuos sobre el modelo final obtenido y por supuesto, justificar cada uno de los pasos. A continuación les ofrezco una breve relatoria sobre los resultados de esta experiencia.

Todo el análisis fue realizado en R (ver la salida del script al final de esta entrada). En R (cómo en SPSS y otros programas de estadística) existen librerías que incluyen funciones para realizar la eliminación hacia atrás o la inclusión hacia adelante de forma automatizada, obteniendo el modelo final mediante una sola llamada a la función. Yo no utilicé tales funciones porque requieren la instalación de paquetes que no vienen incluidos por defecto en la instalación base de R. Además, el uso de estos procedimientos automatizados requiere de mucha discreción. De hecho, los procedimientos como el de la eliminación paso a paso, que se basan en ir eliminando la variable con el menor estadístico F en la tabla ANOVA o el menor estadístico T en las pruebas parciales de los coeficientes son bastante controversiales, por las siguientes razones:

  1. Estos procedimientos se justifican cuando hay una grán cantidad de variables predictoras y no existe un cuerpo de teoría que le permita al investigador tener una idea a priori sobre cuales predictores incluir en el modelo. Este no es el caso aquí: son pocas variables predictoras y si el estudiante hizo la tarea investigando sobre el tema del sobre peso y el índice de masa corporal, debería tener una vaga idea sobre cuales variables inciden sobre el IMC o los niveles de colesterol en la sangre (X9).
  2. Hay mucha controversia entre los especialistas sobre los criterios de inclusión/exclusión de variables en estos algoritmos. Por ejemplo, los valores de los estadísticos F parciales (o los estadísticos t de cada coeficiente) pueden depender del orden de inclusión de las variables en el modelo, con el agravante de que cada secuencia de pruebas es condicionada por las secuencias previas y la resultante exclusión/inclusión. Por último, esto tiende a crear modelos demasiado pequeños (muy pocas variables).

Por mi parte no me gustan los procedimientos automáticos de ajuste de modelos- hay demasiados factores a considerar. Ejemplo de ello es el ajuste que vamos a realizar de los modelos que piden en el enunciado. Para el modelo 1, que inicialmente es X1 ~ X2 + X4 + X9 + X8 + X10 + intercepto. Se dieron los siguientes pasos:

  1. El R2 ajustado inicialmente fue de 29,11%. Se eliminó la variable X10.
  2. Se eliminó la variable X8.
  3. Se eliminó la variable X9.
  4. Se eliminó la variable X4 (fíjese que el p-valor del estadístico t no es mucho mayor que 5%). De hecho, mi "intuición" me dice que esta variable no se debería eliminar. X4 es el sexo y si creo que es significativa sobre el IMC...
  5. Tenemos el modelo X1 ~ X2 + intercepto, con un R2 ajustado de apenas 27,47%.
  6. En este momento se hace el análisis de residuos del modelo final, observando que en la gráfica de residuos versus predictores, la gráfica de dispersión de X2 versus residuos (ver fig. 1) presenta una forma característica de U que indica que hay alguna relación no lineal entre esa variable y la variable X1. Por cierto, en mi librería estUNA, estas gráficas incluyen una línea de tendencia verde (estimador LOWESS de la tendencia en una gráfica de dispersión) para guiar al estudiante a apreciar este tipo de patrones. Esto no lo hace Excel.
  7. Guiado por otra "intuición" (¿vé porqué no me gustan los procedimientos automáticos?), me pregunto porqué diablos no hice un análisis de residuos con el modelo inicial (el que incluye todas las variables) desde el principio. Esto reveló que en efecto, sigue habiendo un indicativo de relación no lineal entre X2 y X1 que este modelo no capturó (con una gráfica de residuos versus predictores muy similar a la gráfica de la fig. 1). No observó otras relaciones no lineales con otras variables (las desviaciones de horizontalidad de la línea verde en las otras gráficas se debe a la poca presencia de valores de las variables en el eje X o a valores atípicos con mucha influencia sobre la regresión). Es entonces cuando decido incluir un término cuadrdático para X2 a mi modelo. Al incluir X2^2, mi R2 ajustado sube a 33,72% y la variable X2^2 es significativa. Vamos bien.
  8. De golpe, se me aviene a la mente otra "intuición"- ¿porqué no volver a incluir la variable X4 (sexo)? Lo hago así y sube el R2 ajustado a 36,11%. La variable X4 también resulta significativa en este modelo (p-valor igual a 1%). Observe que aún cuando X4 no era significativa según la eliminación paso a paso, al incluir una nueva variable - X2^2 - sí resultó significativa! ¿Se dá cuenta porqué hay controversia sobre estos métodos de ajuste de modelos?
  9. Mi interpretación de este modelo es la siguiente: El IMC depende de la edad y el sexo. Durante la adolescencia, los varones suelen "estirarse" en estatura haciendo que el IMC baje, pero con la edad ganan peso y el IMC sube nuevamente. Este fenómeno se diferencia según el sexo- las hembras se estiran a una edad distinta. El efecto de "estirarse" durante la adolescencia explica porqué la relación entre X2 y el IMC es la de una parábola convexa (U hacia arriba). Sin embargo, creo que este modelo sólo sería efectivo para este rango de edades. Para edades mucho mayores que los 13 años, el modelo predice un IMC muy por encima de los esperados.

Fig. 1 - Análisis de residuos. Observe el patrón
de U en el diagrama de dispersión.

Para el modelo 2 se siguieron pasos similares a lo arriba expuesto, obteniendo que el mejor modelo para X9 (colesterol) es aquél que incluye un polinomio de grado dos con la variable X8 (calorías consumidas). Aún así, el R2 de este modelo es demasiado bajo como para ser útil.

Por cierto, en el modelo 1 por ejemplo, los residuos parecen conformarse a los supuestos teóricos de los modelos de regresión lineal: los residuos no tienen desviaciones importantes de normalidad (Fig. 2) y parecen ser heteroscedásticos e independientes entre si y las variables del modelo (Fig. 3).



Fig. 2 - Gráficas de ajuste de normalidad de los residuos
y plot cuantíl-cuantíl para verificar la normalidad.


Fig. 3 - Análisis de residuos final para el modelo 1.
Verificación de homocedasticidad e independencia
de los residuos y las variables de predicción.


La moraleja- hay que hacer el análisis de residuos desde el principio para detectar no-linealidades entre las variables o cualquier problema con los modelos (revise mi artículo titulado "Análisis de Residuos para Dummies"). Otra moraleja- déjese guiar por su intuición: la estadística tiene más de arte que de procedimientos mecánicos. La estadística es una ciencia oculta.

R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R es un software libre y viene sin GARANTIA ALGUNA.
Usted puede redistribuirlo bajo ciertas circunstancias.
Escriba 'license()' o 'licence()' para detalles de distribucion.

R es un proyecto colaborativo con muchos contribuyentes.
Escriba 'contributors()' para obtener más información y
'citation()' para saber cómo citar R o paquetes de R en publicaciones.

Escriba 'demo()' para demostraciones, 'help()' para el sistema on-line de ayuda,
o 'help.start()' para abrir el sistema de ayuda HTML con su navegador.
Escriba 'q()' para salir de R.

> load("estUNA")
> #Utilizaremos el dataset del semestre 2012-1, que es el mismo
> #que el del semestre 2017-1
> attach(d20121)
> png("trabajo%02d.png")
> #Genera la variable X1
> X1 <- round(X5/X6^2,digits=2)
> #----------------------------------------------------
> #MODELO 1 - ELIMINACIÓN HACIA ATRAS PASO A PASO
> #----------------------------------------------------
> modelo1 <- regresion.lineal(X1~X2+X4+X9+X8+X10)
> resumen(modelo1)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo1 
  Marco de datos : variables globales 
  Formula        : X1 ~ X2 + X4 + X9 + X8 + X10 

Estimacion de los coeficientes poblacionales

                Estimacion   Error Est. Estadistico T    p-valor
[Intercepto]  9.9393823481 2.0222851574     4.9149262 2.3869e-06
X2            0.3609888788 0.0491676987     7.3419926 1.4257e-11
X41           0.4504463538 0.2315323235     1.9455009   0.053663
X9            0.0124392038 0.0080671047     1.5419663   0.125276
X8           -0.0003852448 0.0002983862    -1.2910948   0.198740
X10           0.0099836186 0.0135616497     0.7361655   0.462828

Prueba F global

  Valor F : 13.24102   gl. num: 5   gl. den : 144 
  p-valor : 1.2941e-10 

Coeficientes de determinacion

   R^2 : 0.3149548   R^2 ajustado : 0.2911685 

Residuos

  Minimo  : -3.207713 
  Mediana : 0.004358484 
  Maximo  : 3.229979 

  Desv. estandar residual:  1.367149 

--------------------------------------------------------------
> graficar(modelo1)
> modelo1 <- regresion.lineal(X1~X2+X4+X9+X8)
> resumen(modelo1)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo1 
  Marco de datos : variables globales 
  Formula        : X1 ~ X2 + X4 + X9 + X8 

Estimacion de los coeficientes poblacionales

               Estimacion   Error Est. Estadistico T    p-valor
[Intercepto] 10.817027359 1.6309212992      6.632464 6.0987e-10
X2            0.353884004 0.0481349612      7.351912 1.3192e-11
X41           0.446365447 0.2311000578      1.931481   0.055374
X9            0.012749762 0.0080433321      1.585134   0.115114
X8           -0.000399444 0.0002972914     -1.343611   0.181173

Prueba F global

  Valor F : 16.46782   gl. num: 4   gl. den : 145 
  p-valor : 3.8168e-11 

Coeficientes de determinacion

   R^2 : 0.3123767   R^2 ajustado : 0.2934078 

Residuos

  Minimo  : -3.313198 
  Mediana : -0.007141708 
  Maximo  : 3.304436 

  Desv. estandar residual:  1.364988 

--------------------------------------------------------------
> modelo1 <- regresion.lineal(X1~X2+X4+X9)
> resumen(modelo1)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo1 
  Marco de datos : variables globales 
  Formula        : X1 ~ X2 + X4 + X9 

Estimacion de los coeficientes poblacionales

               Estimacion Error Est. Estadistico T    p-valor
[Intercepto] 10.726946123 1.63403053      6.564716 8.5259e-10
X2            0.359219626 0.04810299      7.467719 6.8329e-12
X41           0.457919817 0.23157602      1.977406   0.049879
X9            0.009890414 0.00777806      1.271579   0.205545

Prueba F global

  Valor F : 21.23818   gl. num: 3   gl. den : 146 
  p-valor : 1.791e-11 

Coeficientes de determinacion

   R^2 : 0.3038156   R^2 ajustado : 0.2895104 

Residuos

  Minimo  : -3.4445 
  Mediana : 0.01699545 
  Maximo  : 3.436306 

  Desv. estandar residual:  1.368747 

--------------------------------------------------------------
> modelo1 <- regresion.lineal(X1~X2+X4)
> resumen(modelo1)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo1 
  Marco de datos : variables globales 
  Formula        : X1 ~ X2 + X4 

Estimacion de los coeficientes poblacionales

             Estimacion Error Est. Estadistico T    p-valor
[Intercepto] 12.7293771 0.43701846     29.127779 < 2.22e-16
X2            0.3592637 0.04820381      7.453016 7.2307e-12
X41           0.4277187 0.23083769      1.852898   0.065902

Prueba F global

  Valor F : 30.91906   gl. num: 2   gl. den : 147 
  p-valor : 6.1916e-12 

Coeficientes de determinacion

   R^2 : 0.2961055   R^2 ajustado : 0.2865287 

Residuos

  Minimo  : -3.231319 
  Mediana : 0.01775983 
  Maximo  : 3.431739 

  Desv. estandar residual:  1.371616 

--------------------------------------------------------------
> modelo1 <- regresion.lineal(X1~X2)
> resumen(modelo1)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo1 
  Marco de datos : variables globales 
  Formula        : X1 ~ X2 

Estimacion de los coeficientes poblacionales

             Estimacion Error Est. Estadistico T    p-valor
[Intercepto] 12.8272191 0.43736813     29.328198 < 2.22e-16
X2            0.3670029 0.04841566      7.580251 3.4923e-12

Prueba F global

  Valor F : 57.4602   gl. num: 1   gl. den : 148 
  p-valor : 3.4923e-12 

Coeficientes de determinacion

   R^2 : 0.2796658   R^2 ajustado : 0.2747987 

Residuos

  Minimo  : -2.966373 
  Mediana : -0.03913403 
  Maximo  : 3.668746 

  Desv. estandar residual:  1.382845 

--------------------------------------------------------------
> modelo1 <- regresion.lineal(X1~X2+I(X2^2))
> resumen(modelo1)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo1 
  Marco de datos : variables globales 
  Formula        : X1 ~ X2 + I(X2^2) 

Estimacion de los coeficientes poblacionales

             Estimacion Error Est. Estadistico T    p-valor
[Intercepto] 19.1813771 1.69562250     11.312292 < 2.22e-16
X2           -1.1762778 0.40178676     -2.927617 0.00395964
X2^2          0.0871831 0.02254666      3.866785 0.00016511

Prueba F global

  Valor F : 38.91451   gl. num: 2   gl. den : 147 
  p-valor : 2.7332e-14 

Coeficientes de determinacion

   R^2 : 0.3461698   R^2 ajustado : 0.3372742 

Residuos

  Minimo  : -2.509407 
  Mediana : -0.1023154 
  Maximo  : 3.559675 

  Desv. estandar residual:  1.321938 

--------------------------------------------------------------
> modelo1 <- regresion.lineal(X1~X2+I(X2^2)+X4)
> resumen(modelo1)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo1 
  Marco de datos : variables globales 
  Formula        : X1 ~ X2 + I(X2^2) + X4 

Estimacion de los coeficientes poblacionales

             Estimacion Error Est. Estadistico T    p-valor
[Intercepto] 19.6459116 1.67480152     11.730292 < 2.22e-16
X2           -1.3305022 0.39911327     -3.333646  0.0010864
X2^2          0.0953210 0.02236676      4.261726  3.621e-05
X41           0.5620888 0.22070200      2.546822  0.0119065

Prueba F global

  Valor F : 29.07335   gl. num: 3   gl. den : 146 
  p-valor : 8.4809e-15 

Coeficientes de determinacion

   R^2 : 0.3739818   R^2 ajustado : 0.3611184 

Residuos

  Minimo  : -2.814932 
  Mediana : -0.1555243 
  Maximo  : 3.214668 

  Desv. estandar residual:  1.29794 

--------------------------------------------------------------
> graficar(modelo1)
> #----------------------------------------------------
> #MODELO 2 - ELIMINACIÓN HACIA ATRAS PASO A PASO
> #----------------------------------------------------
> modelo2 <- regresion.lineal(X9~X1+X2+X4+X8+X10)
> resumen(modelo2)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo2 
  Marco de datos : variables globales 
  Formula        : X9 ~ X1 + X2 + X4 + X8 + X10 

Estimacion de los coeficientes poblacionales

               Estimacion   Error Est. Estadistico T    p-valor
[Intercepto] 161.97869371 17.864219175     9.0672137 8.1917e-16
X1             1.30581723  0.846851974     1.5419663 0.12527648
X2            -0.27761681  0.590117211    -0.4704435 0.63874995
X41           -3.06137905  2.389623457    -1.2811136 0.20221308
X8             0.01022012  0.002954539     3.4591235 0.00071311
X10            0.07344650  0.139076249     0.5281024 0.59824116

Prueba F global

  Valor F : 3.103189   gl. num: 5   gl. den : 144 
  p-valor : 0.010855 

Coeficientes de determinacion

   R^2 : 0.09726893   R^2 ajustado : 0.06592411 

Residuos

  Minimo  : -36.61831 
  Mediana : -0.4628569 
  Maximo  : 37.4289 

  Desv. estandar residual:  14.0075 

--------------------------------------------------------------
> modelo2 <- regresion.lineal(X9~X1+X4+X8+X10)
> resumen(modelo2)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo2 
  Marco de datos : variables globales 
  Formula        : X9 ~ X1 + X4 + X8 + X10 

Estimacion de los coeficientes poblacionales

               Estimacion  Error Est. Estadistico T    p-valor
[Intercepto] 161.67709284 17.80471072      9.080580 7.2407e-16
X1             1.09803148  0.72060593      1.523761 0.12974637
X41           -3.04447122  2.38292885     -1.277617 0.20342613
X8             0.01026969  0.00294472      3.487493 0.00064569
X10            0.08660736  0.13586725      0.637441 0.52484285

Prueba F global

  Valor F : 3.844302   gl. num: 4   gl. den : 145 
  p-valor : 0.0053307 

Coeficientes de determinacion

   R^2 : 0.09588151   R^2 ajustado : 0.07094031 

Residuos

  Minimo  : -35.96905 
  Mediana : -0.4268156 
  Maximo  : 37.82659 

  Desv. estandar residual:  13.96984 

--------------------------------------------------------------
> modelo2 <- regresion.lineal(X9~X1+X4+X8)
> resumen(modelo2)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo2 
  Marco de datos : variables globales 
  Formula        : X9 ~ X1 + X4 + X8 

Estimacion de los coeficientes poblacionales

               Estimacion   Error Est. Estadistico T    p-valor
[Intercepto] 169.56149016 12.780959448     13.266726 < 2.22e-16
X1             1.07744252  0.718416600      1.499746 0.13583948
X41           -3.09919146  2.376535563     -1.304080 0.19425870
X8             0.01019272  0.002936255      3.471333 0.00068139

Prueba F global

  Valor F : 5.010667   gl. num: 3   gl. den : 146 
  p-valor : 0.0024586 

Coeficientes de determinacion

   R^2 : 0.09334791   R^2 ajustado : 0.07471807 

Residuos

  Minimo  : -35.83157 
  Mediana : -0.1794051 
  Maximo  : 39.26263 

  Desv. estandar residual:  13.94141 

--------------------------------------------------------------
> modelo2 <- regresion.lineal(X9~X1+X8)
> resumen(modelo2)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo2 
  Marco de datos : variables globales 
  Formula        : X9 ~ X1 + X8 

Estimacion de los coeficientes poblacionales

               Estimacion   Error Est. Estadistico T    p-valor
[Intercepto] 170.56986487 12.787913387     13.338366 < 2.22e-16
X1             0.92115246  0.710035602      1.297333 0.19654904
X8             0.01038678  0.002939462      3.533565 0.00054835

Prueba F global

  Valor F : 6.63407   gl. num: 2   gl. den : 147 
  p-valor : 0.001744 

Coeficientes de determinacion

   R^2 : 0.08278713   R^2 ajustado : 0.07030804 

Residuos

  Minimo  : -37.58161 
  Mediana : 0.4793699 
  Maximo  : 40.21525 

  Desv. estandar residual:  13.97459 

--------------------------------------------------------------
> modelo2 <- regresion.lineal(X9~X8+I(X8^2))
> resumen(modelo2)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo2 
  Marco de datos : variables globales 
  Formula        : X9 ~ X8 + I(X8^2) 

Estimacion de los coeficientes poblacionales

                Estimacion   Error Est. Estadistico T    p-valor
[Intercepto]  1.522697e+02 1.647125e+01      9.244572 2.5171e-16
X8            5.414282e-02 2.091370e-02      2.588868   0.010596
X8^2         -1.359345e-05 6.368918e-06     -2.134342   0.034474

Prueba F global

  Valor F : 8.182145   gl. num: 2   gl. den : 147 
  p-valor : 0.00042734 

Coeficientes de determinacion

   R^2 : 0.1001705   R^2 ajustado : 0.08792796 

Residuos

  Minimo  : -39.82377 
  Mediana : 0.7860429 
  Maximo  : 42.5595 

  Desv. estandar residual:  13.84153 

--------------------------------------------------------------
> graphics.off()
> 
Como citar esta entrada

Romero, J. (Noviembre, 2012). Aventura de regresión lineal y el método de eliminación hacia atrás paso a paso. [Entrada de blog]. Recuperado desde https://unamatematicaseltigre.blogspot.com/2012/11/aventura-de-regresion-lineal-y-el.html.


Si te gustó o te pareció útil este contenido, compártelo en las redes sociales y dale tu voto positivo en el botón "me gusta" de G+, para que otros puedan encontrar el contenido también.
Abajo encontrarás un enlace para descargar una versión pdf imprimible de este artículo más el código fuente en R pero esencialmente es la misma información que estoy compartiendo contigo arriba. El enlace de llevará a un gateway de pago en Bitcoin que es muy seguro y anónimo (no recogeré tu información confidencial sensible). El monto es modesto - 0,0004 BTC o aproximadamente $0,56 al cambio de hoy. Considéralo como una donación y no una transacción comercial que me ayudará a seguir produciendo contenido y a sobrevivir en este hueco infernal en que se ha convertido Venezuela. !Muchas gracias de antemano!