La regresión lineal, busca encontrar una relación funcional entre 2 variables, una de las cuales es la Dependiente (Y) y una independiente (X). Para que la f(x) que se encuentre funcione siempre y con el menor error posible, es necesario que se cumplan ciertos requisitos llamados supuestos.

El modelo teórico tiene la siguiente forma

\[Y = \beta_0+\beta_1*x+e\] y el modelo que debemos encontrar es

\[\hat{Y} = \hat{\beta_0}+\hat{\beta_1}*x\]

Este procedimiento contempla múltiples análisis numéricos y gráficos que procederé a mostrar.

¿qué librerías se requieren?:

Librerias

library(tidyverse)
library(PerformanceAnalytics)
library(AER)
library(car)
library(lmtest)
library(MASS)
library(nortest)
library(moments)
library(ppcor)
library(ggfortify)
library(lawstat)

Cargando datos

Idealmente trabaja con el archivo de datos en la misma carpeta, donde grabas el script de R, abre el script desde esa misma carpeta para que se cargue R, asi no tendras que escribir PATH largos.

Supuestos (CALIDAD DEL MODELO)

Lo que aconsejo en primer lugar y para este caso de modelo es buscar si realmente hay una asociación lineal, para ello, junto con el análisis de robustez se puede analizar dicha relación.

El concepto de calidad del modelo se basa en que si hay algún supuesto que no se cumpla, uno puede estimar que el modelo no será de tan buena calidad, ya que algunos de los requisitos puede generar efecto falso positivo.

- Robustez de Datos

¿Qué tan confiables son los datos, que se utilizaran para crear el modelo?. Se sabe que la Presencia de Valores Atípicos, influye de varias formas a la relación funcional, tanto en el nivel de correlación, como en la dirección de la pendiente.

chart.Correlation(datos)

Como se observa en el gráfico, existen 2 valores en la esquina inferior derecha, lo que provoca un cambio en la pendiente.

Si consideramos dichos valores, una regresión lineal, no sería una opción a considerar para encontrar la relación funcional.

Debemos conocer el contexto de los datos, para poder generar preguntas sobre los datos. En este caso, la pregunta que podemos hacernos es ¿Los valores sobre 650 es posible que obtengan menos de 82?

Si la respuesta es positiva, entonces la relación no es lineal, mientras que si la respuesta en negativa, entonces lo podemos considerar valores atípicos.

Para asegurarnos la presencia de valores atípicos, lo podemos realizar mediante el análisis de residuales, para ello se debe crear el modelo, de la siguiente forma.

mod1<-lm(y~x,data=datos)

Luego se puede realizar con mayor confiabilidad, una gráfica que muestra los valores que tienen mayor influencia sobre el modelo,

plot(mod1,5)

se puede apreciar que el valor 1 esta en la región que indica que tiene mayor incidencia sobre el cambio de pendiente, mientras que el valor 11 esta cerca, al buscar los valore, estos coinciden con lo comentado anteriormente.

Por tanto, para estos datos, es conveniente eliminar los valores atípicos o corregirlos si es que tiene la posibilidad de hacer trazabilidad (buscar tratamiento de valores atípicos)

dado que se tiene claro que los datos atípicos están por sobre el valor 650, considerare dicho valor como cota.

datos1<-datos %>% 
  filter(x<650)

Lo que hicimos fue crear un nuevo dataframe, con 18 datos, Entonces haremos nuevamente el gráfico de correlaciones para ver que ocurrió.

chart.Correlation(datos1)

Como pueden ver, el valor de la correlación cambio significativamente, luego podemos ver si el modelo tiene algunos valores que pudiesen afectar negativamente al modelo ( rehacer el proceso )

mod2<-lm(y~x,datos1)

plot(mod2,5)

Se ve que ningún valor que en la zonas demarcadas por las lineas segmentadas, de hecho ni aparecen dichas zonas. Luego asumimos que no existen valores de alta influencia.

NOTA: para trabajar con Hipótesis, es necesario considerar que existen 2 tipos en estadística y que son complementarias \(H_0\) y \(H_1\) y la idea es aceptar o rechazar \(H_0\) y esto depende de lo que requiere el usuario. Para este caso necesitamos aprobar cada uno de los \(H_0\), pero solo se hace donde se apliquen TEST ESTADISTICOS

Rechazar \(H_0\) ssi el p-value es menor que un nivel de significancia \(\alpha\) que tiene un valor definido por acuerdo de 0,05 pudiendo ser cambiado, entre 0,001 a 0,1.

Este valor \(\alpha\) se entiende como el porcentaje de veces que puede equivocarse al tomar la decisión de rechazar.

- Media del Error o Residuales es 0 o E(e)=0

La Hipótesis es \[H_0: E(e)=0\]

mean(mod2$residuals)
## [1] 6.79072e-17

como se puede apreciar, el valor tiene mas de 15 ceros antes del 6, por tanto tiende a 0.

- Varianza Minima de los residuales

Hipótesis es \[ H_0: min(\sigma^2_{e})\]

La explicación es teórica. Tenemos un modelo \[ Y = \beta_0+\beta_1*x+e\] Como los residuales(e) son aleatorios, es imposible estimarlos por tanto, se espera que los parámetros \(\beta_i\) sean los mejores. ¿Cómo saberlo?, esto está asociado a la forma en que se obtienen los estimadores de los parámetros \(\hat{\beta_i}\) Existen varios métodos para obtener los estimadores, como por ejemplo Mínimos Cuadrados Ordinarios (MCO), Máximo Verosímil, de los Momentos, Bayesiano entre otros. El que se usa para este modelo es el de MCO, el cual, entre todos los métodos asegura varianza mínima, ya que lo que hace es obtener distancias entre \(Y - \hat{Y}\) en linea recta, haciendo que los errores tiendan a la menor variación.

- Normalidad de los errores

Hipótesis es \[ H_0: X\sim N(0,\sigma^2_{e})\] para lograr lo anterior se pueden usar varias funciones que se denominan \(Prueba\) \(de\) \(Normalidad\), pero se debe tener cuidado, dado que tienen condiciones de aplicación.

shapiro.test(mod2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod2$residuals
## W = 0.94019, p-value = 0.292

Los siguientes asociados a la librería NORTEST

ad.test(mod2$residuals)
## 
##  Anderson-Darling normality test
## 
## data:  mod2$residuals
## A = 0.4019, p-value = 0.3228
cvm.test(mod2$residuals)
## 
##  Cramer-von Mises normality test
## 
## data:  mod2$residuals
## W = 0.063548, p-value = 0.3222
pearson.test(mod2$residuals)
## 
##  Pearson chi-square normality test
## 
## data:  mod2$residuals
## P = 4.5556, p-value = 0.336
lillie.test(mod2$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  mod2$residuals
## D = 0.15911, p-value = 0.2657
sf.test(mod2$residuals)
## 
##  Shapiro-Francia normality test
## 
## data:  mod2$residuals
## W = 0.95676, p-value = 0.4604

La siguiente función esta asociada a la librería MOMENTS

jarque.test(mod2$residuals)
## 
##  Jarque-Bera Normality Test
## 
## data:  mod2$residuals
## JB = 1.1894, p-value = 0.5517
## alternative hypothesis: greater
agostino.test(mod2$residuals)
## 
##  D'Agostino skewness test
## 
## data:  mod2$residuals
## skew = -0.06629, z = -0.14344, p-value = 0.8859
## alternative hypothesis: data have a skewness

Como se puede apreciar en cada uno de los p-value, tiene un valor mayor a 0,05 (valor por defecto de \(\alpha\)), por tanto se No se puede rechazar \(H_0\), indicando que los residuales tendrían un comportamiento aproximadamente de una Distribución Normal o Gaussiano.

NOTA: solo escoja un procedimiento, en ocasiones hay ambigüedades entre los métodos.

- Independencia de los errores

Hipótesis es \(H_0: R(e_i,e_j)\) es 0.

esta supuesto, busca que los datos se hayan obtenido de forma aleatoria, tratando de interpretar lo mejor posible al contexto de donde se obtiene los datos. Para esto, se busca que No exista relación entre los errores.

Si bien, se puede aplicar un test, lo que es de mayor relevancia, es que el proceso de obtención de la muestra, el diseño experimental, hayan incorporado el proceso de aleatoriedad,

lawstat::runs.test(mod2$residuals)
## 
##  Runs Test - Two sided
## 
## data:  mod2$residuals
## Standardized Runs Statistic = 0, p-value = 1

NOTA: el análisis de los p-value se efectúa de la misma forma que lo explicado en la nota anterior

Como se puede apreciar el p-value, tiene un valor mayor a 0,05 (valor por defecto de \(\alpha\)), por tanto se No se puede rechazar \(H_0\), indicando que los residuales tendrían un comportamiento aleatorio indicando independencia entre los errores

- Homocedasticidad de los Errores

Hipótesis es \(H_0: V(e_i)\) es constante

que la varianza de los errores sea constante tiene que ver con que al estimar \(Y_i\) para cada uno de los valores de X, las estimaciones tengan el mismo error probable.

en la imagen, el gráfico de la izquierda presenta un ejemplo de Heterocedasticidad, el de la derecha, es Homocedástico. Para entenderlo mejor, fíjense en las lineas, mientras estas sean paralelas con cualquier pendiente, entonces será de varianza constante.

plot(mod2,3)

Por lo general, a muchos les cuesta detectarlo, por esto pueden aplicar el test de Breush-Pagan.

bptest(mod2)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod2
## BP = 0.71707, df = 1, p-value = 0.3971

Como se puede apreciar el p-value, tiene un valor mayor a 0,05 (valor por defecto de \(\alpha\)), por tanto se No se puede rechazar \(H_0\), indicando que los residuales serían Homocedásticos

- No Autocorrelación

Hipótesis es \(H_0: R(e_t,e_{t-1})=0\)

La autocorrelación implica que el estado actual de una observación podría depender de un estado anterior de la misma observación. Esto necesariamente implica de que las medidas son tomadas a través del tiempo. Por tanto, existiría otra variable a considerar en este modelo, dejando de ser Simple para pasar a Múltiple.

en este caso se aplica el test de Durbin Watson.

dwtest(mod2)
## 
##  Durbin-Watson test
## 
## data:  mod2
## DW = 2.491, p-value = 0.8545
## alternative hypothesis: true autocorrelation is greater than 0

Como se puede apreciar el p-value, tiene un valor mayor a 0,05 (valor por defecto de \(\alpha\)), por tanto se No se puede rechazar \(H_0\), indicando que los residuales tendrían un comportamiento que no dependería de estados temporales anteriores.

Evaluación de la Calidad del Modelo

La calidad de los datos, solo asegura que el modelo obtenido siendo de buen o mal ajuste, esta correctamente calculado. Por ende no puedo dudar de sus resultados.

Lo primero a evaluar es la existencia de los parámetros \(\beta_i\), si alguno no existe, tienen diferente impacto.

Nuestro modelo puede tener la forma \[\hat{Y} = \hat{\beta_0}+\hat{\beta_1}*x\]

Para tener claridad debemos analizar la siguiente tabla y reconocer los siguientes valores

- Parámetros del modelo

La hipótesis que se evalua es \(H_0: \beta_i=0\)

La no existencia de \(\beta_0\), significa que es un modelo que pasa por el origen, por tanto, se debe realizar un nuevo modelo, quedando de la siguiente forma

\[\hat{Y} = \hat{\beta_1}*x\]

Si no existe \(\beta_1\), entonces estaríamos en el caso de que la variable independiente No tiene influencia sobre la variable dependiente.

NOTA: el análisis del modelo de regresión múltiple es distinto, para llegar a la misma conclusión, deben asegurarse muchas cosas previas.

En este ejemplo al observar los valores de Pr(t), se ve que los valores son menores que \(\alpha\) (0,05 para este caso), por tanto se rechaza \(H_0\), esto implica que los \(\beta_i\) existen, por ende el potencial modelo es: \[\hat{Y}=60.0487 + 0.0526*X\]

- Correlación

hipótesis \(H_0: R(x,y)=0\), es decir, no existe relación entre las variables.

habiendo obtenido un modelo completo, es decir, con ambos parámetros es necesario corroborar la calidad de la estimación a través de \(R\) y \(R^2\), dichos parámetros al presentar un valor por 0,8 y 0,6 respectivamente, se considera que existe alguna relación CONFIABLE o ALTA o FUERTE.

algunos comentarios, cuando el valor de R, esta en el intervalo (-0,2 : 0,2), se considera que tiene una nula relación (aunque algunos más pragmáticos consideran el rango entre -0.3 a 0.3). Desde 0.3 a 0.6 varían entre débil a media relación.

cor.test(predict(mod2),datos1$y,data=datos1,method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  predict(mod2) and datos1$y
## t = 5.9522, df = 16, p-value = 2.028e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5928520 0.9346767
## sample estimates:
##       cor 
## 0.8299927

El Valor del p-value es menor que \(\alpha\) (0,05), por tanto se rechaza \(H_0\), implicando de que la Relación entre X e Y, existe. No obstante, no dice el grado de fuerza de la relación, para ello vea el valor ubicado al final de la tabla con el valor 0.89299927, este es R.

Por tanto, las variables X, Y presentan una correlación de aproximadamente el 89,3%, es decir, 9 de cada 10 datos, presentan la misma tendencia.

- Modelo

Adicionalmente, debe leerse la tabla anterior buscando el valor \(R^2\) (marcados en rojo).

de esos 2 valores, tenemos \(R^2=0.6889\) y \(R^2_{adj}=0.6694\) (ajustado). Al existir una diferencia se entiende que hay un efecto tamaño de la muestra, es decir, para el comportamiento de los datos de la muestra, la cantidad de unidades a muestrear debería ser mayor. Por lo anterior, se aconseja, cada vez que haya una diferencia entre \(R^2\) y \(R^2_{adj}\), lea el \(R^2_{adj}\) y siempre este debe ser menor a \(R^2\).

Luego,

hipótesis \(H_0: Modelo = 0\) es decir, No existe Este procedimiento reafirma lo encontrado anteriormente, pero adicionalmente entrega un valor que dice relación con el Supuesto de Varianza del error mínima.

acá debes fijarte en el valor encerrado en circulo y el valor encerrado en un cuadrado.

summary(aov(mod2))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## x            1  173.6   173.6   35.43 2.03e-05 ***
## Residuals   16   78.4     4.9                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El Valor del p-value es menor que \(\alpha\) (0,05), por tanto se rechaza \(H_0\), indicando la existencia de un modelo que tiene como \(V(e) = 4.9\).

Luego si encuentra otro modelo, idealmente usar el que tenga varianza minima, asociado al resto de valores analizados.

Estabilidad del Modelo

por cuanto debe seguir usando el mismo modelo. En general, debería estar probando si existe una modificación en el comportamiento. En este aspecto, estudie los intervalos de confianza-

confint(mod2)
##                   2.5 %      97.5 %
## (Intercept) 50.50429521 69.59319120
## x            0.03388029  0.07136356

Entonces de acuerdo a esto, deberiamos esperar que los datos futuros que definen el modelo, generen nuevos \(\beta_0\) y \(\beta_1\) con comportamiento aleatorio en el rango de los valores mostrados. Si percibe de algún cambio en la tendencia, desarrolle un nuevo modelo.

Nota de Cierre

No siempre encontrará un modelo que explique la relación entre los datos. ¿eso es malo?, si y solo si, hizo mal algún procedimiento.

el modelo puede no tener ajuste como si tenerlos, pero puede ser de mala calidad, si es que no se cumple alguno(s) de los supuestos.