sábado, 18 de agosto de 2012

Regresión: selección de variables. Stepwise, Forward, Backward

¿Cómo realizar una selección de variables predictoras para un modelo de regresión lineal múltiple?


A menudo nos encontramos estudiando un gran número de variables explicativas (regresoras, predictoras, causales o independientes) y nos preguntamos cuál de ellas incluir en un modelo de regresión múltiple tanto por problemas de multicolinealidad como por su significación en la regresión.


Estrategias de regresión

Los métodos más comunes de selección de variables son: el de pasos sucesivos (stepwise), el de introducción progresiva (forward) y el de eliminación progresiva (backward).

  • Eliminación hacia atrás (Backward Stepwise Regression). Se introducen todas las variables en la ecuación y después se van excluyendo una tras otra. En cada etapa se elimina la variable menos influyente según el contraste individual (de la t o de la F).
  • 􏰐  Selección hacia adelante (Fordward Stepwise Regression). Las variables se introducen secuencialmente en el modelo. La primera variable que se introduce es la de mayor correlación (+ o -) con la variable dependiente. Dicha variable se introducirá en la ecuación solo si cumple el criterio de entrada. A continuación se
    considerará la variable independiente cuya correlación parcial sea la mayor y que no esté en la ecuación. El procedimiento termina cuando ya no quedan variables que cumplan el criterio de entrada.

  • 􏰐  Pasos sucesivos (Stepwise Regression). Este método es una combinación de los procedimientos anteriores. En cada paso se introduce la variable independiente que no se encuentre ya en la ecuación y que tenga la probabilidad para F más pequeña. Las variables ya introducidas en la ecuación de regresión pueden ser eliminadas del modelo. El método termina cuando ya no hay más variables candidatas a ser incluidas o eliminadas. 
    Fuente: http://eio.usc.es/eipc1/BASE/BASEMASTER/FORMULARIOS-PHP-DPTO/MATERIALES/Mat_50140129_RegresionMultiple.pdf


Criterios de elección de variables

Cuando se dispone de muchas variables explicativas potenciales, las estrategias de regresión anteriores definen normalmente un subconjunto posible de modelos y el problema es seleccionar entre ellos. Para ellos se puede utilizar distintos estadísticos, como el criterio de información de Akaike AIC, el citerior de información Bayesiana BIC, el coeficiente de determinación R2, el coeficiente de determinación corregido R2aj, la varianza residual o el estadístico Cp de Mallows.


Implementación en R-project

El programa R-project, al igual que otros programas (SPSS, SAS, StatSoft Statistica, etc.), cuenta con librerías específicas para la selección de variables explicativas. Entre ellas encontramos la función stepAIC del paquete MASS, que permite seleccionar el método mediante el argumento direction = c("both", "backward", "forward").

También está la función leaps del paquete leaps, que permite realizar todas las regresiones con subconjuntos de variables (package).

La función eleaps del paquete subselect permite realizar la selección de variables mediante ocho criterios diferentes (vignettes).

Sin embargo, algunos investigadores plantean que la selección del modelo mediante un criterio automático no es recomendable (e.g. Peña 1998). La práctica, dicen, muestra que los buenos modelos suelen coincidir con cualquier criterio razonable y la selección final del modelo es mejor hacerla en base a su adecuación lógica a la realidad que describe. Si esto no es suficiente, recomiendan utilizar criterios de validez externa.
Por ello, para explicar el funcionamiento de los métodos de selección de variables predictoras, en esta entrada no utilizaré paquetes sino los pasos específicos para los tres métodos que he comentado.


Funciones de correlación parcial

Primero, es necesario especificar las funciones de correlación parcial tanto de primer orden como de segundo y tercer orden. (NOTA: podemos encontrar librerías específicas para estas funciones en R-project).

#correlación parcial de primer orden
pcor = function(v1, v2, v3)
{
c12 = cor(v1, v2)
c23 = cor(v2, v3)
c13 = cor(v1, v3)
 
partial = (c12-(c13*c23))/(sqrt(1-(c13^2)) * sqrt(1-(c23^2)))
 
return(partial)
}
 
#correlación parcial de segundo orden
pcor2 = function(v1, v2, v3, v4)
{
c123 = pcor(v1, v2, v3)
c143 = pcor(v1, v4, v3)
c243 = pcor(v2, v4, v3)
 
partial = (c123-(c143*c243))/(sqrt(1-(c143^2)) * sqrt(1-(c243^2)))
 
return(partial)
}
 
#correlación parcial de tercer orden
pcor3 = function(v1, v2, v3, v4, v5) #r12|345=(r12|34-r15|34*r25|34)/sqrt((1-r15|34^2)*(1-r25|34^2))
{
c1234 = pcor2(v1, v2, v3, v4)
c1534 = pcor2(v1, v5, v3, v4)
c2534 = pcor2(v2, v5, v3, v4)
 
partial = (c1234-(c1534*c2534))/(sqrt(1-(c1534^2)) * sqrt(1-(c2534^2)))
 
return(partial)
}
 


Algoritmos de selección de variables en modelos de regresión.

(próximamente)


Ejemplo en R.

Segundo, utilizaremos un conjunto de datos bien conocido para estos procedimientos. Se trata de un estudio que analiza la relación entre la composición de un cemento tipo Portland y el calor desprendido durante la fase de fraguado. La muestra está formada por 13 cementos. Y es la cantidad de calor desprendido (cals/gr), mientras que las variables X1 , X2 , X3 y X4 representan el contenido ( %) de cuatro ingredientes activos.

library(MASS)
data(cement)
datos=apply(cement,2,as.numeric)
str(datos)
attach(datos)
 
 
Utilizaremos los tres métodos más comunes de selección de variables predictivas para determinar qué variables serán incluidas en el modelo final.


#### FORWARD
N=13
p=1 #el parámetro del intercepto
 
#mayor correlación: x4
cor(cement)
summary.aov(lm(y~x4)) #Fexp=22.8
Fexp=(cor(x4, y)^2*(N-p-1))/(1-cor(x4, y)^2)
Fteo=qf(.95,1,(N-p-1))
Fexp #[1] 24.87111
Fteo #[1] 4.747225
#como Fexp>Fteo, entra x4
p=p+1
 
#mayor correlación parcial de primer orden: x1
pcor(x1, y, x4)
pcor(x2, y, x4)
pcor(x3, y, x4)
Fexp=(pcor(x1, y, x4)^2*(N-p-1))/(1-pcor(x1, y, x4)^2)
Fteo=qf(.95,1,(N-p-1))
Fexp #[1] 108.2239
Fteo #[1] 4.964603
#como Fexp>Fteo, entra x1
p=p+1
 
#mayor correlación parcial de segundo orden: x2
pcor2(x2,y,x4,x1) # 0.5986053
pcor2(x3,y,x4,x1) #-0.5657105
Fexp=(pcor2(x2,y,x4,x1)^2*(N-p-1))/(1-pcor2(x2,y,x4,x1)^2)
Fteo=qf(.95,1,(N-p-1))
Fexp #[1] 5.025865
Fteo #[1] 5.117355
#Como Fexp<Fteo, no entra x2
 
#FIN: el modelo final incluye a x1 y x4


El modelo final incluirá las variables X1 y X4.



 
#### BACKWARD
N=13
p=4 #todos los parámetros
 
#menor correlación parcial de tercer orden: x3
pcor3(x1,y,x2,x3,x4) #0.5929326
pcor3(x2,y,x1,x3,x4) #0.2418094
pcor3(x3,y,x1,x2,x4) #0.04768649
pcor3(x4,y,x1,x2,x3) #-0.07164829
 
Fexp=(pcor3(x3,y,x1,x2,x4)^2*(N-p-1))/(1-pcor3(x3,y,x1,x2,x4)^2)
Fteo=qf(.95,1,(N-p-1))
Fexp #[1] 0.01823347
Fteo #[1] 5.317655
#como Fexp<Fteo, sale x3
p=p-1
 
#menor correlación parcial de segundo orden: x4 (menor r)
pcor2(x1,y,x2,x4) #0.972002
pcor2(x2,y,x1,x4) #0.5986053
pcor2(x4,y,x1,x2) #-0.4141492
 
Fexp=(pcor2(x4,y,x1,x2)^2*(N-p-1))/(1-pcor2(x4,y,x1,x2)^2)
Fteo=qf(.95,1,(N-p-1))
Fexp #[1] 1.863262
Fteo #[1] 5.117355
#como Fexp<Fteo, sale x4
p=p-1
 
#menor correlación parcial de primer orden: x1 (menor r)
pcor(x1,y,x2) #0.9675285
pcor(x2,y,x1) #0.9768575
 
Fexp=(pcor(x1,y,x2)^2*(N-p-1))/(1-pcor(x1,y,x2)^2)
Fteo=qf(.95,1,(N-p-1))
Fexp #[1] 146.5227
Fteo #[1] 4.964603
#como Fexp>Fteo,no sale x1
 
#FIN: el modelo final incluye a x1 y x2
 
 
El modelo final incluirá las variables X1 y X2.


#### STEPWISE
N=13
p=0 #sin el parámetro del intercepto
 
#mayor correlación: x4
cor(cement)
summary.aov(lm(y~x4)) #Fexp=22.8
Fexp=(cor(x4, y)^2*(N-p-1))/(1-cor(x4, y)^2) #22.79852
Fteo=qf(.95,1,(N-p-1)) #4.844336
Fexp #[1] 24.87111
Fteo #[1] 4.747225
#como Fexp>Fteo, entra x4
p=p+1
 
#mayor correlación parcial de primer orden: x1
pcor(x1, y, x4)
pcor(x2, y, x4)
pcor(x3, y, x4)
Fexp=(pcor(x1, y, x4)^2*(N-p-1))/(1-pcor(x1, y, x4)^2) #108.2239
Fteo=qf(.95,1,(N-p-1)) # 4.964603
Fexp #[1] 119.0463
Fteo #[1] 4.844336
#como Fexp>Fteo, entra x1
p=p+1
 
#pruebo la significación de las variables introducidas en pasos anteriores: x4
Fexp=(pcor(x4, y, x1)^2*(N-p-1))/(1-pcor(x4, y, x1)^2) #159.2952
Fteo=qf(.95,1,(N-p-1)) # 4.964603
Fexp #[1] 159.2952
Fteo #[1] 4.964603
#como Fexp>Fteo, sigo sin eliminar x4
 
#mayor correlación parcial de segundo orden: x2
pcor2(x2,y,x4,x1) # 0.5986053
pcor2(x3,y,x4,x1) #-0.5657105
Fexp=(pcor2(x2,y,x4,x1)^2*(N-p-1))/(1-pcor2(x2,y,x4,x1)^2) #5.584294
Fteo=qf(.95,1,(N-p-1)) # 4.964603
Fexp #[1] 5.584294
Fteo #[1] 4.964603
# como Fexp>Fteo entra x2
p=p+1
 
#pruebo la significación de las variables introducidas en pasos anteriores: x4 y x1
Fexp=(pcor2(x4, y, x1,x2)^2*(N-p-1))/(1-pcor2(x4, y, x1,x2)^2) #2.070292
Fteo=qf(.95,1,(N-p-1)) # 4.964603
Fexp #[1] 1.863262
Fteo #[1] 5.117355
#como Fexp<Fteo, debo eliminar x4
 
Fexp=(pcor2(x1, y, x2,x4)^2*(N-p-1))/(1-pcor2(x1, y,x2,x4)^2) #171.1196
Fteo=qf(.95,1,(N-p-1)) # 4.964603
Fexp #[1] 154.0076
Fteo #[1] 5.117355
#como Fexp>Fteo, se queda x1
 
#solo queda x3, correlación parcial de segundo orden: x3
pcor2(x3,y,x2,x1) # 0.4112643
Fexp=(pcor2(x3,y,x2,x1)^2*(N-p-1))/(1-pcor2(x3,y,x2,x1)^2)
Fteo=qf(.95,1,(N-p-1))
Fexp #[1] 1.832128
Fteo #[1] 5.117355
#como Fexp<Fteo, no entra x3.
 
#FIN: modelo final con x1 y x2



El modelo final incluirá las variables X1 y X2.

(NOTA: También hay investigadores que fijan de antemano el criterio de inclusión o exclusión de antemano, el Fteórico).

Espero que este material les sea de ayuda. Saludos.