Regresión lineal

El término regresión fue introducido por Galton en su libro “Natural inheritance” (1889) refiriéndose a la “ley de la regresión universal”: La primera forma de regresión lineal documentada fue el método de los mínimos cuadrados que fue publicada por Legendre en 1805, Gauss publicó un trabajo en donde desarrollaba de manera más profunda el método de los mínimos cuadrados,1​ y en dónde se incluía una versión del teorema de Gauss-Márkov.

https://www.youtube.com/watch?v=AAkV09ALjK8 En esta liga explico el metodo de los minimos cuadrados

Los métodos de regresión estudian la construcción de modelos para explicar o representar la dependencia entre una variable respuesta o dependiente (Y ) y la(s) variable(s) explicativa(s) o dependiente(s), X.

https://www.youtube.com/watch?v=SgPvRLxpdJE En esta liga explico paso a paso como se calcula cada uno de los coeficientes de Regresion y los estadisticos de ajuste. Le ayudara a entender mejor una regresion.

El modelo estadístico en regresión lineal simple se puede escribir de la siguiente forma:

\(Y= B0 + B1*x + e\)

En esta forma la variable respuesta Y se expresa como una suma de β0+β1Xi y un error aleatorio ϵi el cual tiene distribución N(0,σ2). El modelo en esta forma se puede expresar como sigue.

Las librerias

library(DescTools)
## Warning: package 'DescTools' was built under R version 4.1.3
library(tidyr)
library(pastecs)
library(PerformanceAnalytics)
## Warning: package 'zoo' was built under R version 4.1.3
library(nortest)
library(normtest)
library(correlation)
## Warning: package 'correlation' was built under R version 4.1.3
library(boot)
library(corrplot)
library(qgraph)

Los datos

data ("airquality")        # conjunto de datos 
Datos<-airquality
Datos<-Datos %>% drop_na() # eliminar celdas con NA
attach(Datos)

Lo tradicional

n<-length(Temp)
plot(Temp,Ozone, col="deepskyblue1", pch=16
     )   # plot entre temperatura y ozono
mod <- lm(Ozone ~ Temp, data=Datos)    # el modelo de regresion lineal simple
abline(mod,col="red")

# Podemos obtener el valor predicho de Y para cada valor de X
fitted <- predict(lm(Ozone ~ Temp))


# Con un bucle for y la funcion lines() podemos representar los RESIDUALES, ERRORES
for(i in 1:n)
{
  lines( c(Temp[i],Temp[i]), c(Ozone[i], fitted[i]), col="red")
}

Coeficientes del modelo

La columna Estimate indica el valor estimado de los parámetros o coeficientes del modelo; en regresión lineal simple, (Intercept) indica el punto donde la línea recta corta el eje Y (la ordenada al origen) mientras que el segundo estimado, coeficiente, indica la pendiente de la línea recta.

La columna Std. Error. Es el error estandar del coeficiente de regresion.

La columna t value indica el valor del estadístico t asociado a cada coeficiente.

La columna Pr(>|t|) es el valor p correspondiente de la prueba de hipótesis; la hipótesis a validar es que el coeficiente es igual a cero, es decir, no tiene impacto en el fenómeno. SI el valor p es menor a 0.05 los coeficientes son diferentes de cero (estadísticamente significativos) para un nivel de confianza del 95%.

Precision del modelo

Residual standard error:. Es el error tipico del modelo de regresion. Las unidades son en las mimas unidades que la variable y. Si Y es en ln, este valor es tambien en Ln.

Multiple R-squared:. Conocido como coeficiente de determinación. Indica la proporcion de la variable Y explicada por la varable X. Este valor aumenta si el numero de variables independientes aumenta, por lo que no es adecuado para comparar el ajuste entre modelos de diferente tamaño. Si se calcula la raiz de este valor, se tendra el coeficiente de correlacion multiple.

Adjusted R-squared:. Es el estadistico mas adecuado para determinar el ajuste del modelo, sobre todo cuando los modelos son de difrente tamaño. Tiene la misma interpretacion del anterior.

El valor F-statistic y p-value (valor p) corresponden a la prueba F del modelo de regresión que verifica la existencia de una relación lineal entre las variables.

p-value: Detemina si existe significancia estadistica entre la(s) variable(s) independiente(s) y dependiente. Un valor menor a 0.05, indica una relacion estadisticamente significativa.

summary(mod)
## 
## Call:
## lm(formula = Ozone ~ Temp, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.922 -17.459  -0.874  10.444 118.078 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -147.6461    18.7553  -7.872 2.76e-12 ***
## Temp           2.4391     0.2393  10.192  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.92 on 109 degrees of freedom
## Multiple R-squared:  0.488,  Adjusted R-squared:  0.4833 
## F-statistic: 103.9 on 1 and 109 DF,  p-value: < 2.2e-16
-147.6461+(2.4391*100)
## [1] 96.2639

Supuestos básicos a tener en cuenta al construir un modelo predictivo

Supuesto 1: Los parámetros del modelo de regresión lineal deben ser de naturaleza numérica y lineal.

Si los parámetros no son numéricos, como categóricos, utilice la codificación one-hot (python) o la codificación ficticia (R) para convertirlos a numéricos. ¿cómo la verifico?

Supuesto 2: La media de los residuos es cero.

Compruebe la media de los residuos. Si es cero (o muy cercano), entonces esta suposición se cumple para ese modelo. ¿cómo la verifico?

plot(Temp,Ozone, col="deepskyblue1")    # la plot entre las variables

mod <- lm(Ozone ~ Temp, data=Datos)     # el modelo lineal simple
mean(mod$residuals)                     # media de los residuales
## [1] 9.692377e-16

Supuesto 3: Homocedasticidad de residuos o igual varianza:

Este supuesto significa que la varianza alrededor de la línea de regresión es la misma para todos los valores de la variable predictora (X). ¿cómo verifico la homocedasticidad? ¿cómo la verifico?

H0:los errores tienen varianza constante.

H1:los errores no tienen varianza constante.

plot(mod, which = 1)  # la plot para visualizar la varianza (errores)

library(lmtest)       # la libreria para la prueba estadistica
bptest(mod)           # la prueba de Breusch-Pagan para 
## 
##  studentized Breusch-Pagan test
## 
## data:  mod
## BP = 1.5088, df = 1, p-value = 0.2193
library(car)
plot(mod, which = 3)  # tambien aqui se observa la distrobucion de la varianza

Se observa que el valor-P es mayor que el nivel de significancia usual de 5%, por lo tanto, hay evidencias para decir que se cumple la homocedasticidad de los errores.

Supuesto 4: No hay autocorrelación de residuos

La autocorrelación es especialmente aplicable para datos de series de tiempo. Cuando los residuos están autocorrelacionados, significa que el valor actual depende de los valores anteriores (históricos) y que existe un patrón definido e inexplicable en la variable Y que aparece en las perturbaciones. ¿cómo la verifico?

H0:los errores son independientes.

H1:los errores no son independientes.

acf(mod$residuals)  # plot la autocorrelacion de los residuales del modelo

library(lmtest)     # la libreria
dwtest(mod)         # aplica la prueba Durbin-Watson 
## 
##  Durbin-Watson test
## 
## data:  mod
## DW = 1.8644, p-value = 0.2123
## alternative hypothesis: true autocorrelation is greater than 0
plot(residuals(mod), pch=19, col="deepskyblue1", type = "l")  # plot de los residuales de otra forma

Se observa que el valor-P es mayor que el nivel de significancia usual de 5%, por lo tanto, hay evidencias para decir que los errores son independientes.

Supuesto 5: Los residuos y las variables X no deben estar correlacionados

¿cómo la verifico?

mod <- lm(Ozone ~ Temp, data=Datos)   # el modelo de regresion
cor.test(Temp, mod$residuals)         # la correlacion de Pearson
## 
##  Pearson's product-moment correlation
## 
## data:  Temp and mod$residuals
## t = 2.1232e-15, df = 109, p-value = 1
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1863929  0.1863929
## sample estimates:
##          cor 
## 2.033609e-16
plot(Temp, mod$residuals)             # graficacmente como se ve

Supuesto 6: El número de observaciones debe ser mayor que el número de predictores o X variables

Suposición 7: La variabilidad en los predictores o valores de X es positiva.

Lo que esto infiere es que la variación en los predictores no debería ser la misma (o incluso casi la misma). ¿cómo la verifico?

mod <- lm(Ozone ~ Temp, data=Datos)   # el modelo de regresion
var(Temp)
## [1] 90.82031

Supuesto 8: No hay multicolinealidad (correlacion) perfecta entre los predictores

Lo que esto significa es que no debería haber una relación lineal perfecta entre los predictores o las variables explicativas. ¿cómo la verifico?

Utilice el factor de inflación de varianza (VIF=1/(1-R2)). VIF es una métrica calculada para cada variable X que entra en un modelo lineal. Si el VIF de una variable es alto, significa que la información en esa variable ya está explicada por otras variables X presentes en el modelo dado, lo que significa que esa variable es más redundante. según algunas referencias, si el VIF es demasiado grande (más de 5 o 10), consideramos que la multicolinealidad existe. Entonces, baje el VIF (menos de 2) mejor.

mod <- lm(Ozone ~ Temp, data=Datos)   # el modelo de regresion
library(car)        # la liberia para el vif()
#vif(mod)           # No aplica por es un modelo lineal simple una variable
library(corrplot)   # plot de correlacion
corrplot(cor(Datos[, -1]), method="number")

Supuesto 9: La normalidad de los residuos.

Los residuos deben estar normalmente distribuidos. Esto se puede verificar visualmente usando el gráfico qqnorm().

# par(mfrow=c(2,2))  # use la linea para 4 plots juntas
mod <- lm(Ozone ~ Temp, data=Datos)   # el modelo de regresion
plot(mod, which = 2)                  # elegir solo la plot 2 del mod

hist(mod$residuals)                   # histograma de los residuales

boxplot(mod$residuals)                # boxplot de los residuales

shapiro.test(mod$residuals)           # prueba de normalidad
## 
##  Shapiro-Wilk normality test
## 
## data:  mod$residuals
## W = 0.88983, p-value = 1.522e-07

Comprueba automáticamente algunos supuestos antes mencionados

  • Global Stat: Rechazo de la hipótesis nula significa que no se cumple el supuesto de linealidad entre las variables.

  • Skewness: Rechazo de la hipótesis nula significa que no se cumple el supuesto de normalidad debido a una asimetría significativa.

  • Kurtosis: Rechazo de la hipótesis nula significa que no se cumple el supuesto de normalidad debido a una curtosis significativa.

  • Link Function: Rechazo de la hipótesis nula significa que el modelo de regresión lineal no es el más adecuado para el fenómeno estudiado.

  • Heteroscedasticy: Rechazo de la hipótesis nula significa que no se cumple el supuesto de homocedasticidad.

library(gvlma)
par(mfrow=c(2,2))  # draw 4 plots in same window
mod <- lm(Ozone ~ Temp, data=Datos)   # el modelo de regresion
gvlma::gvlma(x = mod)
## 
## Call:
## lm(formula = Ozone ~ Temp, data = Datos)
## 
## Coefficients:
## (Intercept)         Temp  
##    -147.646        2.439  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma::gvlma(x = mod) 
## 
##                        Value   p-value                   Decision
## Global Stat        1.924e+02 0.000e+00 Assumptions NOT satisfied!
## Skewness           4.941e+01 2.076e-12 Assumptions NOT satisfied!
## Kurtosis           1.312e+02 0.000e+00 Assumptions NOT satisfied!
## Link Function      1.185e+01 5.764e-04 Assumptions NOT satisfied!
## Heteroscedasticity 5.076e-03 9.432e-01    Assumptions acceptable.

Conclusion

El modelo generado no estadisticamente valido por no cumplir o pasar todos los supuestos de un modelo de regresion.