El análisis de regresión es un enfoque estadístico ampliamente utilizado que busca identificar relaciones entre variables. La idea es agrupar datos relevantes para tomar mejores decisiones. La regresión paso a paso (stepwise) es la construcción iterativa paso a paso de un modelo de regresión que implica la selección automática de variables independientes. La disponibilidad de paquetes de software estadístico hace posible la regresión stepwise, incluso en modelos con cientos de variables.
El objetivo de la regresión stepwise es, a través de una serie de pruebas (pruebas F, pruebas t) encontrar un conjunto de variables independientes que influyan significativamente en la variable dependiente. Esto se hace con las computadoras a través de la iteración, que es el proceso de llegar a resultados o decisiones pasando por rondas repetidas o ciclos de análisis. Realizar pruebas automáticamente con la ayuda de paquetes de software estadístico tiene la ventaja de ahorrar tiempo para el individuo.
Stepwise regression se puede lograr probando una variable independiente a la vez e irincluyéndola en el modelo de regresión si es estadísticamente significativa o incluyendo todas las variables independientes potenciales en el modelo y eliminando aquellas que no son estadísticamente significativas. Algunos usan una combinación de ambos métodos y, por lo tanto, existen tres enfoques para la regresión por pasos:
Forward selection comienza sin variables en el modelo, prueba cada variable a medida que se agrega al modelo, luego mantiene las que se consideran estadísticamente más significativas, repitiendo el proceso hasta que los resultados sean óptimos.
Backward elimination comienza con un conjunto de variables independientes, eliminando una a la vez, luego probando para ver si la variable eliminada es estadísticamente significativa.
Bidirectional elimination es una combinación de los dos primeros métodos que prueba qué variables deben incluirse o excluirse.
rm(list=ls(all=TRUE)) # Borramos todos los objetos y bases de datos existentes
getwd() # Ver en que directorio de trabajo estamos
## [1] "G:/UAAAN/MATERIAS/2019/EPIDOMETRIA"
library(pacman)
p_load(corrplot, ecospat, ggplot2, dplyr, tidyr, sampler, reshape, tidyverse, knitr, psych)
Los datos que se utilizarán son de las librerías de r, estos son: data(“trees”). Estos datos vienen como: “Girth”, “Height”, y “Volume”, es decir; es la circunferencia de árbol a 1.3 m, la altura total del árbol, y el volumen. Son 31 observaciones y 3 variables de Black Cherry (Prunus serótina). Note que las unidades son en pulgadas, pies y pies cubicos, respectivamente. Para eso será necesario convertirlas al sistema métrico decimal (S.M.D.), para mejor entendimiento. Ademas re-nombraremos las variables.
data("trees") # Datos de la libreria
attach(trees) # Adjuntamos los datos como dataframe
str(trees) # Vemos la estructura de los datos
## 'data.frame': 31 obs. of 3 variables:
## $ Girth : num 8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
## $ Height: num 70 65 63 72 81 83 66 75 80 75 ...
## $ Volume: num 10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...
Como se indicó anteriormente, los datos vienen en pulgadas, pies y pies cúbicos. Para trabajar en nuestro sistema métrico, es necesario hacer algunas conversiones. Aquí hacemos todo esto. Además, para mayor facilidad, renombraremos las variables.
D<-Girth / (pi) # Convertir de "inches" a "cm" y llamar "Girth" como "D"
H<-(Height*0.3048) # Convertir de "ft" a "m" y llamar "Height" como "H"
V<-(Volume*0.0283168) # Convertir de "cubic ft" a "m3" y llamar "Volume" como "V"
Densidad<-((0.49+0.67)/2)*1000 # Promediar la densidad y transformarla de g/cm3 a kg/m3
B<-Densidad*V # Transformar el volumen (m3) a Biomasa (kg)
Datos<-data.frame(D, H, V, B)
attach(Datos) #La base de datos prelimina
## The following objects are masked _by_ .GlobalEnv:
##
## B, D, H, V
En este apartado, el investigador hace transformaciones o genera variables a partir de las existentes. estas transformaciones pueden ser en la variable dependiente, sobre las variables independientes, o sobre ambas. Lo más común es transformar las variables utilizando el ln, log, elevar al cuadrado, multiplicar, o dividir las variables.
Hagamos algunas trasformaciones para generar un modelo. Suponga que deseamos generar un modelo para estimar biomasa (B).
head(Datos) # Ver con que datos contamos y sus nombres
## D H V B
## 1 2.641972 21.3360 0.2916630 169.1646
## 2 2.737465 19.8120 0.2916630 169.1646
## 3 2.801127 19.2024 0.2888314 167.5222
## 4 3.342254 21.9456 0.4643955 269.3494
## 5 3.405916 24.6888 0.5323558 308.7664
## 6 3.437747 25.2984 0.5578410 323.5478
D2<-D*D # Generr una variable (diametro al cuadrado)
#D3<-D^3 # Esta no la voy a incluir en el ejemplo
H2<-H*H # Generr una variable (H al cuadrado)
#H3<-H^3 # Esta no la voy a incluir tampoco
DH<-D*H # .... ASI SUCESIVAMENTE
D2H<-D2*H
DH2<-D*H2
lnD<-log(D)
Una vez generadas las variables, hagamos un “daframe”, solo con las variables de interés. Esto es importante (solo las de interés) porque se van a incluir todas al mismo tiempo, las que existan en el “daframe”. Llamaremos al “daframe” Set1.
Set1<-data.frame(B, D, H, D2, H2, DH, D2H, lnD ); Set1
## B D H D2 H2 DH D2H lnD
## 1 169.1646 2.641972 21.3360 6.980016 455.2249 56.36912 148.9256 0.9715256
## 2 169.1646 2.737465 19.8120 7.493715 392.5153 54.23466 148.4655 1.0070323
## 3 167.5222 2.801127 19.2024 7.846312 368.7322 53.78836 150.6680 1.0300218
## 4 269.3494 3.342254 21.9456 11.170660 481.6094 73.34777 245.1468 1.2066454
## 5 308.7664 3.405916 24.6888 11.600262 609.5368 84.08797 286.3966 1.2255139
## 6 323.5478 3.437747 25.2984 11.818103 640.0090 86.96949 298.9791 1.2348162
## 7 256.2104 3.501409 20.1168 12.259863 404.6856 70.43714 246.6292 1.2531654
## 8 298.9121 3.501409 22.8600 12.259863 522.5796 80.04220 280.2605 1.2531654
## 9 371.1766 3.533240 24.3840 12.483783 594.5795 86.15452 304.4046 1.2622152
## 10 326.8325 3.565071 22.8600 12.709729 522.5796 81.49752 290.5444 1.2711839
## 11 397.4546 3.596902 24.0792 12.937702 579.8079 86.61052 311.5295 1.2800728
## 12 344.8986 3.628733 23.1648 13.167701 536.6080 84.05887 305.0272 1.2888835
## 13 351.4681 3.628733 23.1648 13.167701 536.6080 84.05887 305.0272 1.2888835
## 14 349.8257 3.724226 21.0312 13.869857 442.3114 78.32493 291.6997 1.3148590
## 15 313.6935 3.819719 22.8600 14.590250 522.5796 87.31877 333.5331 1.3401768
## 16 364.6071 4.106198 22.5552 16.860858 508.7370 92.61611 380.3000 1.4124974
## 17 555.1225 4.106198 25.9080 16.860858 671.2245 106.38337 436.8311 1.4124974
## 18 450.0106 4.233521 26.2128 17.922704 687.1109 110.97245 469.8043 1.4430341
## 19 422.0902 4.360845 21.6408 19.016973 468.3242 94.37218 411.5425 1.4726659
## 20 408.9512 4.392676 19.5072 19.295606 380.5309 85.68882 376.4032 1.4799387
## 21 566.6192 4.456338 23.7744 19.858952 565.2221 105.94677 472.1347 1.4943274
## 22 520.6327 4.520000 24.3840 20.430403 594.5795 110.21569 498.1750 1.5085121
## 23 596.1819 4.615493 22.5552 21.302779 508.7370 104.10338 480.4884 1.5294188
## 24 629.0294 5.092958 21.9456 25.938223 481.6094 111.76802 569.2299 1.6278588
## 25 699.6515 5.188451 23.4696 26.920025 550.8221 121.77087 631.8022 1.6464352
## 26 909.8754 5.506761 24.6888 30.324417 609.5368 135.95532 748.6735 1.7059766
## 27 914.8025 5.570423 24.9936 31.029612 624.6800 139.22492 775.5417 1.7174710
## 28 957.5043 5.697747 24.3840 32.464320 594.5795 138.93386 791.6100 1.7400708
## 29 845.8228 5.729578 24.3840 32.828064 594.5795 139.71003 800.4795 1.7456419
## 30 837.6109 5.729578 24.3840 32.828064 594.5795 139.71003 800.4795 1.7456419
## 31 1264.6283 6.557184 26.5176 42.996657 703.1831 173.88077 1140.1682 1.8805612
Antes de todo hagamos un análisis de correlación, para ver de forma general, si existe correlación o no entre las variables generadas (independientes) y la variable dependiente, es decir, de la que deseamos predecir o estimar (Biomasa, B).
Corr<-cor(Set1, method = "pearson") # Hacer la correlacion
round(Corr, digits = 3)
## B D H D2 H2 DH D2H lnD
## B 1.000 0.967 0.598 0.979 0.603 0.977 0.989 0.940
## D 0.967 1.000 0.519 0.993 0.520 0.972 0.980 0.992
## H 0.598 0.519 1.000 0.508 0.999 0.697 0.596 0.530
## D2 0.979 0.993 0.508 1.000 0.512 0.971 0.992 0.970
## H2 0.603 0.520 0.999 0.512 1.000 0.699 0.601 0.529
## DH 0.977 0.972 0.697 0.971 0.699 1.000 0.987 0.961
## D2H 0.989 0.980 0.596 0.992 0.601 0.987 1.000 0.953
## lnD 0.940 0.992 0.530 0.970 0.529 0.961 0.953 1.000
Basta con ver la segunda línea de los resultados.
Aquí se observa que la Biomasa (B), se correlaciona mas con el Diámetro normal (B), el cual llega a tener un valor de 0.967, mientras que la menor correlación de B, ocurre con la variable altura (H), cuyo valor es de solo 0.598. No obstante, esto es correlaciona múltiple, en https://rpubs.com/jorge_mendez/602205, puede ver cómo hacer la correlación parcial, y determinar la influencia de una tercera variable. .
Ahora veamos de forma gráfica como se ve esta correlación.
library(psych)
pairs.panels(x = Set1, ellipses = FALSE, lm = TRUE, method = "pearson", hist.col = "#CD1076") # Hacer la correlacion en figuras
En la diagonal de la Figura superior, puede ver el histograma de frecuencias de cada una de las variables. Debajo de la diagonal, la figura de dispersión de cada variable con las demás. Y por encima de la diagonal, se observan los valores del coeficiente de correlación.
Los resultados demuestran que la Biomasa (B), no parece tener una distribución normal, es posible que, al hacer la regresión, no se cumpla con el supuesto de normalidad. La variable H2, parece ser la más normalmente distribuida. De la misma forma, la correlación más fuerte de Biomasa (B) es con diámetro (D). Esto es solo otra forma de ver la correlaciona entre las variables.
library(qgraph)
## Registered S3 methods overwritten by 'huge':
## method from
## plot.roc pROC
## plot.sim BDgraph
## print.roc pROC
## print.sim BDgraph
Otra_Corr=cor(Set1)
qgraph(Otra_Corr, shape="circle", posCol="darkgreen", negCol="darkred", layout="spring", vsize=10)
Lo mismo puede ver en la figura anterior, mientras más cercano este los círculos (las variables) y entre más gruesa sea la línea que las une, la correlación es más alta y viceversa.
Ahora procederemos a generar el modelo. Lo primero que tenemos que hacer es estructurar el modelo y correrlo sin ninguna variable independiente. A este objeto le llamaremos ”regvacia”. Recuerde que el objeto ”Set1”, contiene solo las variables de interés. La variable Volumen (V) ya no fue incluida en este set. Recuerde también que este es un modelo de regresión lineal múltiple.
regvacia<-lm(formula = B~1, Set1); summary(regvacia) # B, es la variable dependiente (Biomasa)
##
## Call:
## lm(formula = B ~ 1, data = Set1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -328.00 -176.90 -98.07 117.09 769.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 495.52 48.49 10.22 2.75e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 270 on 30 degrees of freedom
Como puede observar, después de correr al modelo, en los resultados aparece solo con el coeficiente beta cero o intercepta, cuyo valor es de 495.52, pero no aparece el coeficiente ni ninguna variable independiente.
En este apartado, haremos un nuevo modelo pero ahora, se incluirán todas las variables independientes que existen en objeto ”Set1”. A este nuevo objeto el llamaremos, ”regcompleta”.
regcompleta<-lm(formula = B~.,Set1); # B, es la variable dependiente (Biomasa)
summary(regcompleta)
##
## Call:
## lm(formula = B ~ ., data = Set1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -80.693 -15.976 1.042 31.643 67.493
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1393.4269 3869.9804 -0.360 0.722
## D 1466.2000 3570.8365 0.411 0.685
## H 68.5103 173.3075 0.395 0.696
## D2 -119.1389 319.0984 -0.373 0.712
## H2 -0.7985 2.4081 -0.332 0.743
## DH -16.6705 71.4788 -0.233 0.818
## D2H 3.4673 8.7350 0.397 0.695
## lnD -2199.4890 4049.3601 -0.543 0.592
##
## Residual standard error: 45.12 on 23 degrees of freedom
## Multiple R-squared: 0.9786, Adjusted R-squared: 0.9721
## F-statistic: 150.1 on 7 and 23 DF, p-value: < 2.2e-16
De este procedimiento se observa que ya existe el modelo con todas (siete) las variables independientes. Cada una tiene su propio coeficiente. Puede ver también, que ninguno de estos coeficientes es estadísticamente significativo, puesto que el valor de Pr(>|t|), supera por mucho el valor de 0.05. Esto No es importante, aun no es el modelo definitivo.. Lo que si puede ver es que como se explicó en clase, entre mayor sea el número de variables independientes en el modelo, mayor será la R cuadrada, esta llega a ser de 0.9786, para este moldeo. Pero, el modelo no es estadísticamente bueno.
Finalmente, una vez tenido lo anterior, solo resta decirle al programa que utilice ambas regresiones, y que use el procedimiento Forward para generar el modelo. Es decir, el procedimiento, comienza sin variables en el modelo ”regvacia”, después toma una a una de las variables del objeto ”regcompleta”, y las agrega al modelo, luego mantiene las que se consideran estadísticamente más significativas, repitiendo el proceso hasta tener el modelo solo con las variables estadísticamente significativas.
regforw<-step(regvacia, scope = list(lower=regvacia, upper=regcompleta),direction = "forward"); summary(regforw)
## Start: AIC=348.08
## B ~ 1
##
## Df Sum of Sq RSS AIC
## + D2H 1 2137913 48617 232.09
## + D2 1 2097700 88830 250.78
## + DH 1 2085993 100537 254.61
## + D 1 2045105 141425 265.19
## + lnD 1 1932293 254237 283.37
## + H2 1 794559 1391971 336.08
## + H 1 782565 1403965 336.35
## <none> 2186530 348.08
##
## Step: AIC=232.09
## B ~ D2H
##
## Df Sum of Sq RSS AIC
## <none> 48617 232.09
## + H 1 253.30 48363 233.93
## + H2 1 231.98 48385 233.94
## + D2 1 208.02 48409 233.96
## + DH 1 123.14 48494 234.01
## + D 1 118.27 48498 234.01
## + lnD 1 102.79 48514 234.02
##
## Call:
## lm(formula = B ~ D2H, data = Set1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -75.87 -18.07 -2.72 28.66 68.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.88901 15.82519 -0.309 0.76
## D2H 1.12976 0.03164 35.711 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.94 on 29 degrees of freedom
## Multiple R-squared: 0.9778, Adjusted R-squared: 0.977
## F-statistic: 1275 on 1 and 29 DF, p-value: < 2.2e-16
De ese procedimiento observe lo siguiente: lm(formula = B ~ D2H, data = Set1) Esto indica que el modelo generado es: B= B0+B1D2H, es decir: Biomasa = -4.88901+ 1.12976* D2H. De las siete variables que se generaron, solo la variable D2H fue estadísticamente significativa. Ninguna otra, ni combinación de ellas pueden predecir la biomasa de esta especie.
Todos los estadísticos generados ya son correctos. El error de estimación es de 40.94 kg. La Adjusted R-squared es de 0.977; es decir esta variable (D2H) explica el 97.7 % de la biomasa aérea de esta especie. Es un buen modelo.
Observe detenidamente que el procedimiento se realizó para generar este modelo fue con forward, si desea utilizar backward o ambos both solo cambie la palabra en la línea correspondiente. El resultado será el mismo.
regstep<-step(regvacia, scope = list(lower=regvacia, upper=regcompleta),direction = "both"); summary(regstep)
## Start: AIC=348.08
## B ~ 1
##
## Df Sum of Sq RSS AIC
## + D2H 1 2137913 48617 232.09
## + D2 1 2097700 88830 250.78
## + DH 1 2085993 100537 254.61
## + D 1 2045105 141425 265.19
## + lnD 1 1932293 254237 283.37
## + H2 1 794559 1391971 336.08
## + H 1 782565 1403965 336.35
## <none> 2186530 348.08
##
## Step: AIC=232.09
## B ~ D2H
##
## Df Sum of Sq RSS AIC
## <none> 48617 232.09
## + H 1 253 48363 233.93
## + H2 1 232 48385 233.94
## + D2 1 208 48409 233.96
## + DH 1 123 48494 234.01
## + D 1 118 48498 234.01
## + lnD 1 103 48514 234.02
## - D2H 1 2137913 2186530 348.08
##
## Call:
## lm(formula = B ~ D2H, data = Set1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -75.87 -18.07 -2.72 28.66 68.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.88901 15.82519 -0.309 0.76
## D2H 1.12976 0.03164 35.711 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.94 on 29 degrees of freedom
## Multiple R-squared: 0.9778, Adjusted R-squared: 0.977
## F-statistic: 1275 on 1 and 29 DF, p-value: < 2.2e-16
El bootstrap es un mecanismo propio de la estadística y la econometría que se centra en el remuestreo de datos dentro de una muestra aleatoria o al azar. Su principal uso es hallar una aproximación a la distribución de la variable analizada.
Este proceso también es conocido en el argot estadístico como bootstrapping y es fruto dentro de los estudios en el campo del muestreo estadístico por el matemático Bradley Efron a finales de los años 70.
La principal utilidad del empleo del bootstrap es reducir el sesgo dentro de análisis o, en otras palabras, aproximar la varianza gracias a la realización de remuestreos aleatorios de la muestra inicial y no de la población. De este modo, se hace más sencillo la construcción de modelos estadísticos mediante la creación de intervalos de confianza y contrastes de hipótesis.
Aunque pueda parecer una práctica muy compleja a priori, el procedimiento en que se basa el bootstrapping es simplemente la creación de un gran número de muestras reposicionando los datos tomando como referencia una muestra poblacional inicial.
Para este ejemplo, utilizaremos la librería, “bootStepAIC”. El procedimiento es similar al anterior. Hay que indicar cuál es la variable dependiente (B) y cuál es el set de variables independientes que debe tomar para generar el modelo, en este caso, todas las variables contenidas en el objeto ”Set1”.
library(pacman) # La libreria
p_load(bootStepAIC,MASS ) # Cargasr las librerias
lmFit <-lm(formula = B~.,Set1) # El modelo, B = variable dependiente.
boot.stepAIC(regstep, Set1)
##
## Summary of Bootstrapping the 'stepAIC()' procedure for
##
## Call:
## lm(formula = B ~ D2H, data = Set1)
##
## Bootstrap samples: 100
## Direction: backward
## Penalty: 2 * df
##
## Covariates selected
## (%)
## D2H 100
##
## Coefficients Sign
## + (%) - (%)
## D2H 100 0
##
## Stat Significance
## (%)
## D2H 100
##
##
## The stepAIC() for the original data-set gave
##
## Call:
## lm(formula = B ~ D2H, data = Set1)
##
## Coefficients:
## (Intercept) D2H
## -4.889 1.130
##
##
## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## B ~ D2H
##
## Final Model:
## B ~ D2H
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 29 48616.73 232.0898
boot.stepAIC(regstep, Set1, B = 100, alpha = 0.05, direction = "forward", # Hacer 100 remuestreos
k = 2, verbose = F) # Usar el procedimiento "forward"
##
## Summary of Bootstrapping the 'stepAIC()' procedure for
##
## Call:
## lm(formula = B ~ D2H, data = Set1)
##
## Bootstrap samples: 100
## Direction: forward
## Penalty: 2 * df
##
## Covariates selected
## (%)
## D2H 100
##
## Coefficients Sign
## + (%) - (%)
## D2H 100 0
##
## Stat Significance
## (%)
## D2H 100
##
##
## The stepAIC() for the original data-set gave
##
## Call:
## lm(formula = B ~ D2H, data = Set1)
##
## Coefficients:
## (Intercept) D2H
## -4.889 1.130
##
##
## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## B ~ D2H
##
## Final Model:
## B ~ D2H
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 29 48616.73 232.0898
# Pruebde cambiar "T" por "F", vera todas las iteraciones que
# Hace este prcedimiento
Como puede ver, los resultados son igual (VER DONDE DICE Final Model:).Se genera el mismo modelo, pero este procedimiento es mejor que los anteriores ya que de la muestra, se hace “n” re muestreos, dando mayor validez al modelo generado. Este procedimiento se hizo para generae el modelo, pero “bootstrap” se puede utilizar para calcular los Intervalos de Confianza (I.C.) de los coeficientes de regresión, de la R2, del Error, del coeficiente de correlacion, etc. Es decir, de cualquier estadístico.
MUY IMPORTANTE UNA VEZ GENERADO EL MODELO, PRUEBE QUE CUMPLA CON TODOS LOS SUPUESTOS DE UN MODELO DE REGRESIÓN LINEAL MÚLTIPLE. https://rpubs.com/jorge_mendez/602205
De datos de las librerías de r, data(“trees”) realice los siguiente:
1.- Obtenga una muestra aleatoria de 25 datos. (https://rpubs.com/jorge_mendez/601482)
2.- Genere el número de variables que considere adecuadas (mínimo siete) y haga los siguientes cuatro sets, para generar un modelo de cada set para estimar Biomasa. Nota Transformar significa aplicar ln)
Set | Variable dependiete (y) | Variables independienets (x) | |
---|---|---|---|
Uno | Sin trasnformar | Sin trasnformar | |
Dos | Transformada con ln | Transformadas con ln | |
Tres | Sin trasformar | Transformadas con ln | |
Cuatro | Sin trasformar | Transformadas con ln y no transformadas |
3.- De los cuatro sets, seleccione el mejor modelo y pruebe todos los supuestos de un modelo de regresión multiple https://rpubs.com/jorge_mendez/602205