Para llevar acabo un análisis vialble y confiable con un modelo de regresión lineal no depende de solo ajustar la recta. Los errores que presenta el modelo hay que análisarlos, viendo si cumplen con estos supuestos: De normalidad, homocedasticidad y Independecia. Usando estos métodos de supuestos podremos ver si los modelos de regresión simple y múltiple cumplen con los supuestos para ver si matematicamente estos modelos son válidos.
Estaré descargando datos de Gapminder para poder realizar los análisis primero del modelo de regresión lineal simple y luego regresión lineal multiple.
library(readr)
supuestos <- read_csv("supuestos.csv")
library (tidyverse)
library (ggplot2)
library(broom)
library (lmtest)
Para análisar los supuestos en el modelo de regresión lineal simple estarems tomando en cuenta estas variables:
Variable dependiente: expectativa de vida femenina (lifExpFem)
Variable independiente: médicos por 1000 habitantes (doctor)
Ahora se hace un data frame para poder ver la relación en los datos explicandos la expectativa de vida en relación a la cantidad de médicos por cada 1000 personas. Y depués veré un modelo ajustado para ver la relación entre ambas variables.
data1 <- data.frame(
y1 = supuestos$lifExpFem,
x1 = supuestos$doctor
)
Y depués veré un modelo ajustado para ver la relación entre ambas variables:
modelo1 <- lm(y1 ~ x1, data = data1)
summary (modelo1)
##
## Call:
## lm(formula = y1 ~ x1, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.053 -5.513 1.614 6.222 14.067
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 61.6726 0.8568 71.98 <2e-16 ***
## x1 5.3042 0.4314 12.29 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.771 on 170 degrees of freedom
## (22 observations deleted due to missingness)
## Multiple R-squared: 0.4706, Adjusted R-squared: 0.4675
## F-statistic: 151.1 on 1 and 170 DF, p-value: < 2.2e-16
df1 <- data.frame(
yhat = fitted.values(modelo1),
res = rstandard(modelo1))
El supuesto de normalidad lo que nos ayuda a estudiar si los errores del modelo poseen una distribución normal. Esto garantiza la validez exacta de pruebas e intervalos en muestras pequeñas. Se va a la gráfica de dispersión para ver si hay algo diferente que se pueda identificar.
ggplot(df1, aes(sample = res)) +
stat_qq(color = "lightskyblue", alpha = 0.7) +
stat_qq_line(linewidth = 1, color = "black", linetype = "solid") +
labs(x = "Cuantiles teóricos", y = "Cuantiles muestrales") +
theme_minimal(base_size = 14)
Viendo esta gráfica es evidente decir que hay patrones raros en ella, vemos una desviación en ambas colas de la linea, y el centro de la gráfica se ve un movimiento raro como un poco de desviación de la linea.
Visiaulmente podemos asumir que no cumple con el supuesto, pero necesitamos tomar otras medidas más certeras para confimar si la conlusión es correcta. Voy a usar el Shapiro Test para estudiar el p-Value ver si rechazamos o aceptamos la hipotesis núla.
H0: Datos son normales
Ha: Los datos no son normales
Shapiro Test
shapiro.test(df1$res)
##
## Shapiro-Wilk normality test
##
## data: df1$res
## W = 0.9619, p-value = 0.000122
Al ver que el p-Value es aproximadamente 0 rechazamos la hipóstesis nula, dejandonos saber que los datos no son normales así que este modelo no es válido.
El supuesto de Homoscendasticidad nos ayuda a ver si la varianza se ha mantenido constante a lo largo de todos los niveles de los valores ajustados. Viendo las gráficas hay que evaluar como se distribuye la variabilidad de los datos, velando que la nube de puntos debe permanecer constante alrededor de 0 (sin ningún patrón).
ggplot(df1, aes(x = yhat, y = res)) +
geom_point(alpha = 0.7, color = "blue") +
geom_hline(yintercept = 0, linetype = "dashed", color = "orange4") +
labs(x = "Valores ajustados", y = "Residuales estandarizados") +
theme_minimal(base_size = 14)
Mirando la gráfica se puede ver un patrón en la nube de puntos, viendo al principio de la gráficas estan ascendiendo y rapidamente la dirección de los puntos cae. Así que por ese patrón se puede decir que hay una homoscendasticidad pobre o débil.
Como mencioné anteriormente hay que corroborar que nuestra observación es correcta hay que hacer una prueba de Hipotesis pero esta vez con el Breusch-Pagan Test.
La prueba de esta hipotesis es:
H0: La varianza es constante
Ha: La varianza no es constante
Breusch-Pagan Test:
bptest(modelo1)
##
## studentized Breusch-Pagan test
##
## data: modelo1
## BP = 6.7478, df = 1, p-value = 0.009386
Viendo el resultado del p-Value de este modelo podemos ver que que es mejor muy cercano a 0 así que rechazamos la hipótesis núla y confirmando nuestra conclución que que la varianza de este modelo no es constane. Así que al iguall que el supuesto anterior, este supuesto no se cumple.
El supuesto de Independencia muestra si los errores no están correlacionados entre sí. Es para ver que los puntos no vayan a depender de uno y del otro. Para poder hacer esta gráfica hay que hacer una columna para enumerar los residuales estandarizados para después ponerlos en la gráfica.
df1.5 <- data.frame(
res = rstandard(modelo1)
) %>%
mutate(orden = 1:length(res))
ggplot(data = df1.5, aes(x = orden, y = res)) +
geom_point(alpha = 0.7, color = "aquamarine4") +
geom_hline(yintercept = 0, linetype = "dashed", color = "orange4") +
labs(x = "Orden/tiempo", y = "Residuales estandarizados") +
theme_minimal(base_size = 14)
Viendo la gráfica los puntos se ven bastantes independientes están expandidos. Esto muestra que los errores son independientes entre sí. Pero para corroborar se estará haciendo una prueba de Hipótesis. Usando el Durbin-Watson Test
H0: No hay correlación entre los residuales.
Ha: Hay correlación entre los reciduales.
Durbin-Watson Test
dwtest (modelo1)
##
## Durbin-Watson test
##
## data: modelo1
## DW = 1.9794, p-value = 0.4393
## alternative hypothesis: true autocorrelation is greater than 0
Viendo el resultado del p-Value podemos ver que es mayor que el alpha standard de 0.05. Así que aceptamos la hipotesis núla y podemos decir que no hay correlación entre los residuales. Son independientes entre sí. Así que este modelo cumple con este supuesto
Para la regresión multiple el proceso es bastante similar. En este caso estaré analizando:
Variable dependiente: tasa de fecundidad (tfr)
Variables independientes: uso de anticonceptivos (contracep) y promedio de años de educación femenina (yearSchF)
data2 <- data.frame(
y2 = supuestos$tfr,
x1 = supuestos$contracep,
x2 = supuestos$yearSchF
)
Ahora se estará formando el modelo ajustado en el que se incluyan las variables mecionados.
modelo2 <- lm(y2 ~ x1 + x2, data = data2)
summary(modelo2)
##
## Call:
## lm(formula = y2 ~ x1 + x2, data = data2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.88868 -0.52076 0.06251 0.50355 2.22631
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.950649 0.169315 41.051 < 2e-16 ***
## x1 -0.042085 0.004648 -9.054 3.25e-15 ***
## x2 -0.194993 0.030539 -6.385 3.44e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7979 on 119 degrees of freedom
## (72 observations deleted due to missingness)
## Multiple R-squared: 0.8018, Adjusted R-squared: 0.7985
## F-statistic: 240.7 on 2 and 119 DF, p-value: < 2.2e-16
df2 <- data.frame(
yhat = fitted.values(modelo2),
res = rstandard(modelo2))
Este modelo ajustado no se puede considerar válido hasta que se verifiqué que cumple con los tres tipos de supuestos.
ggplot(df2, aes(sample = res)) +
stat_qq(color = "blue1") +
stat_qq_line(linewidth = 1, color = "black", linetype = "solid") +
labs(x = "Cuantiles teóricos", y = "Cuantiles muestrales") +
theme_minimal(base_size = 14)
Cuando vemos esta gráfica lineal se puede observar que los residuales de
la gráfica estan bastante cerca de la linea. Al estos residuales np
desviarse tanto de la linea podemos concluir que los errores de este
modelo son normales. Para corroborar esta conclusión debemos hacer la
prueba de hipótesis.
H0: Datos son normales
Ha: Los datos no son normales
Shapiro Test:
shapiro.test(df2$res)
##
## Shapiro-Wilk normality test
##
## data: df2$res
## W = 0.99231, p-value = 0.7397
Viendo el resultado del p-Value (0.73) es mayor que alpha (0.05), podemos concluir que la distribución de los datos son normales. Así que este modelos cumple con el supuesto.
ggplot(df2, aes(x = yhat, y = res)) +
geom_point(alpha = 0.7, color = "red") +
geom_hline(yintercept = 0, linetype = "dashed", color = "orange4") +
labs(x = "Valores ajustados", y = "Residuales estandarizados") +
theme_minimal(base_size = 14)
Viendo esta gráfica es difícil identificar algún patrón en la gráfica
auqnue se ve que al inicio de la gráfica los puntos estan más
concentrados en la linea y poco a poco se van dispersando. No puedo
tomar una conclusión solo viendo la gráfica así que haré una prueba de
hipótesis para aclarar la duda.
H0: La varianza es constante
Ha: La varianza no es constante
Breusch-Pagan Test:
bptest(modelo2)
##
## studentized Breusch-Pagan test
##
## data: modelo2
## BP = 4.5645, df = 2, p-value = 0.1021
Viendo el valor del p-Value (0.10) podemos aceptar la hipótesis nula y decir que el modelo tiene una varianza constante. Así que este modelo cumple con el supuesto.
df2.5 <- data.frame(
res = rstandard(modelo2)) %>%
mutate(orden = 1:length(res))
ggplot(df2.5, aes(x = orden, y = res)) +
geom_point(alpha = 0.7, color = "lightskyblue") +
geom_hline(yintercept = 0, linetype = "dashed", color = "orange4") +
labs(x = "Orden/tiempo", y = "Residuales estandarizados") +
theme_minimal(base_size = 14)
No se puede logra captar ningún patrón que muestre una dependencia entre
los puntos, todos estan dispersos alrededor del 0. Para poder podrar esa
conclusión estaré haciend una prueba de hipótesis nuevamente.
H0: No hay correlación entre los residuales.
Ha: Hay correlación entre los reciduales.
Durbin-Watson Test
dwtest(modelo2)
##
## Durbin-Watson test
##
## data: modelo2
## DW = 1.9599, p-value = 0.4015
## alternative hypothesis: true autocorrelation is greater than 0
El p-Value nos ayuda a aceptar la hipótesis nula, que como pudimos ver en la gráfica podemos decir que no hay correlación entre los residuales, haciendo que este modelo cumpla con el supuesto de independencia.
Se puede concluir que el modelo de regresión simple lineal no cumple con los supuesto de normalidad o homocedasticidad. Así que podemos decir que los resultados del modelo de regresión simple no son válidos porque no cumple con los tres supuesto. Por el otro lado, el modelo de Regresión múltiple lineal si cumple con cada uno de los supestos mostrando que sus resultados de los errores son válidos y se puede continuar usando el modelo para un análisis.