Introducción

Un modelo de regresión no conlleva meramente ajustar una recta. Para asegurar que la interpretación sea válida, los errores deben cumplir ciertos supuestos (De normalidad, homocedasticidad e independencia a un tiempo u orden). El objetivo aquí es analizar los supuestos de dos modelos planteados, uno siendo de Regresión Lineal Simple y el otro Múltiple, más lograr identificar en base a los supuestos si cada modelo es útil y válido. Si estos modelos no cumplen con los supuestos, no podemos seguir adelante con el análisis, ya que matemáticamente el modelo no es válido.

Para llevar a cabo este ejercicio, trabajaremos con datos descargados de la base de datos Gapminder. Estos datos contienen informacion específica sobre cada país del mundo. En el primer ejercicio, utilizaremos una regresión lineal simple y para el segundo ejercicio, una lineal múltiple.

Como siempre, iniciaremos insertando nuestra base de datos para que el sistema los guarde. De paso, corremos todas las librerias de Rstudio necesarias para hacer nuestro ejercicio.

# Base de Datos
library(readr)
supuestos <- read_csv("supuestos.csv")

# Librerías
library(tidyverse)
library(ggplot2)
library(broom)
library(lmtest)

Una vez terminado esto, iniciamos el ejercicio.

Regresión Lineal Simple

Para este modelo tendremos como variable dependiente “expectativa de vida femenina” y como variable independiente “médicos por 1000 habitantes”.

Por lo tanto, creamos un data frame que contenga columnas donde se reflejen los datos sobre la expectativa de vida femenina y los médicos por 1000 habitantes de cada país.

data1 <- data.frame(
  y1 = supuestos$lifExpFem,
  x1 = supuestos$doctor
)

Luego, creamos un modelo ajustado que nos permita analizar la relacion 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))

Supuesto de Normalidad

El supuesto de normalidad lo que nos perimte analizar es si los errores poseen una distribución normal, es decir, que no hayan sido alterados y contengan sesgos que arrastren y dañen los datos. Para esto, utilizamos una gráfica de dispersión que nos ayude a visualizar algún patrón extraño.

## Gráfica
ggplot(df1, aes(sample = res)) +
  stat_qq(color = "red", 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)

De tan solo verla, nos damos cuenta que, incluso en el centro, los datos parecen no estar bien fijados. Además de que los datos se desvían considerablemente de la recta en las colas.

Al no siempre ser suficiente con mirar las gráficas, realizaremos una prueba Shapiro-Wilk que analiza si los residuales son normales o no. Esta es una prueba de hipótesis realizada de la siguente manera:

H0: Los datos son normales

Ha: Los datos no son normales

## Prueba Shapiro-Wilk
shapiro.test(df1$res)
## 
##  Shapiro-Wilk normality test
## 
## data:  df1$res
## W = 0.9619, p-value = 0.000122

Al tener un pvalue aproximadamente cero, esto nos indica que debemos rechazar la hipótesis nula, lo que nos significa que los datos efectivamente no son normales y que el modelo no es válido.

Supuesto de Homoscedasticidad

El supuesto de homoscedasticidad lo que nos dice es que tan homogénea es la variabilidad de los residuales del modelo. Es decir, que la varianza se haya mantenido constante a lo largo de todos los niveles de los valores ajustados. Una varianza heterogénea puede resultar en eroores al calcular pvalues e intervalos de confianza.

Para visualizar esto, utilizaremos una gráfica que nos permita ver como se distribuye la variabilidad de los datos.

## Gráfica
ggplot(df1, aes(x = yhat, y = res)) +
  geom_point(alpha = 0.7, color = "magenta3") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "orange4") +
  labs(x = "Valores ajustados", y = "Residuales estandarizados") +
  theme_minimal(base_size = 14)

En esta gráfica se observa como la variabilidad de los datos ascienden abruptamente al principio y luego decaen poco a poco. Esto muestra una muy baja homoscedasticidad.

Para estar seguros de esta afirmacion, realizamos la prueba Breusch-Pagan para calcular si la variabilidad de los datos realmente es homogenea o no. Esta prueba tiene las siguientes hipótesis:

H0: La varianza es constante

Ha: La varianza no es constante

#Prueba Breusch-Pagan
bptest(modelo1)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo1
## BP = 6.7478, df = 1, p-value = 0.009386

Al igual que con la varianza, el pvalue es cercano a cero por lo que rechazamos la hiótesis nula y confirmamos que la varianza no es constante y no se cumple el supuesto de homoscedasticidad.

Supuesto de Independencia

El supuesto de independencia nos enseña que tan independientes son los errores entre sí. Es decir, velar que un error no dependa del anterior ni del siguiente, siendo esto de máxima importancia en datos temporales. Si lo hacen, pueden volver las pruebas de hiótesis muy optimistas y malestima los errores estándar.

Al igual que en los casos anteriores, realizaremos una gráfica para analizar visualmente si los datos siguen un patrón determinado o si realmente son independientes entre sí. Antes de, debemos crear una columna en la que enumeremos los residuales estandarizados pra luego insertarlos en la gráfica.

## 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 = "hotpink2") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "orange4") +
  labs(x = "Orden/tiempo", y = "Residuales estandarizados") +
  theme_minimal(base_size = 14)

Al mirar la gráfica, no se logra visualizar ningún patrón en particular y todos los puntos parecen estar graficados aleatoriamente. Esto, de momento, nos indica que los errores son independientes entre sí.

Para confirmarlo, realizaremos la prueba Durbin-Watson que analiza la independencia de los residuales estandarizados. Esta prueba contiene una prueba de hipótesis de la siguiente manera:

H0: No hay autocorrelación entre los residuales.

Ha: Hay correlación entre los residuales.

#Prueba Durbin-Watson
dwtest(modelo1)
## 
##  Durbin-Watson test
## 
## data:  modelo1
## DW = 1.9794, p-value = 0.4393
## alternative hypothesis: true autocorrelation is greater than 0

Debido a que el pvalue es mayor que cualquier valor razonable de alpha, concluimos que no se puede rechazar la hipótesis nula y por tanto los residuales son independientes.

Regresión Lineal Múltiple

Para la regresión lineal múltiple haremos un ejercicio parecido. Analizaremos la relación entre la tasa de fecundidad(tfr) con el uso de anticonceptivos(contracep) y promedio de años de educación femenina(yearSchF).

Por tanto, creamos un data frame que incluya nuestras variables respuesta y explicativas.

data2 <- data.frame(
  y2 = supuestos$tfr,
  x1 = supuestos$contracep,
  x2 = supuestos$yearSchF
)

Una vez listo, creamos nuestro modelo ajustado en el que se incluyan las variables anteriormente mencionadas.

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))

Como hemos puntualizado anteriormente, este modelo ajustado no puede ser considerado válido sin antes estar seguros de que cumple con los tres supuestos.

Supuesto de Normalidad

##Grafica
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)

Comenzamos analizando la gráfica y observamos que los residuales están bastante alineados a la línea, desviándose poco en las esquinas. Lo que esto nos puede indicar es que los errores del modelo son normales.

## Prueba Shapiro-Wilk
shapiro.test(df2$res)
## 
##  Shapiro-Wilk normality test
## 
## data:  df2$res
## W = 0.99231, p-value = 0.7397

Una vez realizada la prueba Shapiro-Wilk, confirmamos que los errores si tienen una distribución normal. Esto lo sabemos debido a que el pvalue es mayor a cualquier valor de alpha razonable. Esto implica que, al menos por la normalidad de sus residuales, este modelo es válido.

Supuesto de Homoscedasticidad

## Gráfica
ggplot(df2, aes(x = yhat, y = res)) +
  geom_point(alpha = 0.7, color = "palegreen3") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "orange4") +
  labs(x = "Valores ajustados", y = "Residuales estandarizados") +
  theme_minimal(base_size = 14)

Al observar esta gráfica podemos ver como la nube de puntos no esta completamente constante a lo largo de la línea. Se distribuyen más puntos al final que al principio. Aun así, si es un poco complicado notar un patrón y es por esto que llevamos a cabo la prueba Breusch-Pagan y asegurar claridad.

#Prueba Breusch-Pagan
bptest(modelo2)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo2
## BP = 4.5645, df = 2, p-value = 0.1021

El pvalue en esta prueba es cercano a cero por lo que rechazamos la hiótesis nula y confirmamos que la varianza no es constante. El p-value en esta prueba sugiere heterocedasticidad en la variabilidad de los errores.

Supuesto de Independencia

## Gráfica
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 = "lightpink1") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "orange4") +
  labs(x = "Orden/tiempo", y = "Residuales estandarizados") +
  theme_minimal(base_size = 14)

Al mirar la gráfica, no se logra apreciar un patrón en particular y todos los puntos parecen estar dispersos alrededor de cero. Esto nos indica que los errores son independientes entre sí.

Para confirmarlo, realizaremos nuevamente la prueba Durbin-Watson que analiza la independencia de los residuales estandarizados.

#Prueba Durbin-Watson
dwtest(modelo2)
## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 1.9599, p-value = 0.4015
## alternative hypothesis: true autocorrelation is greater than 0

En esta prueba el pvalue es mayor que cualquier valor razonable de alpha, por ende no se puede rechazar la hipótesis nula y se concluye que los errores no están correlacionados entre sí.

Conclusión

Luego de evaluar ambos modelos, concluimos que ninguno de los dos cumple completamente con los supuestos requeridos para considerar sus resultados como unos válidos. El modelo de Regresión Lineal Simple no cumple con los supuestos de normalidad ni homocedasticidad, mientras que el modelo de Regresión Lineal Múltiple no cumple con el supuesto de homocedasticidad.

Esto es una indicación de que los modelos analizados en este trabajo pueden tener estimaciones poco precisas y mayor sensibilidad a observaciones atípicas. Antes de utilizar cualquiera de estos modelos para la inferencia o predicción de algun suceso, se requiere ajustarlos o transformarlos para que así sea más confiable y seguro el análisis.