Introducción

La regresión lineal es una técnica estadística ampliamente utilizada que permite modelar y analizar la relación entre una variable dependiente y una o más variables independientes. En su forma más simple, asume una relación lineal entre estas variables y busca encontrar la mejor línea recta que se ajuste a los datos observados.

Paquetes y dataset

library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(dslabs)
library(writexl)
library(lmtest)
## Cargando paquete requerido: zoo
## 
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(plotly)
## 
## Adjuntando el paquete: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(corrplot)
## corrplot 0.92 loaded
data("murders")
murders <- mutate(murders, rate = total / population * 100000)

Coeficiente de correlación de Pearson

El coeficiente de correlación de Pearson indica la fuerza y dirección de la relación lineal entre las dos variables. Un valor cercano a 1 o -1 indica una fuerte relación, mientras que un valor cercano a 0 indica una relación débil o nula.

#Matriz de correlación
# Selecciona solo las variables numéricas
numeric_vars <- murders[, sapply(murders, is.numeric)]

# Calcula la matriz de correlación
correlation_matrix <- cor(numeric_vars)
#correlation_matrix
corrplot(correlation_matrix, method = "number", tl.col = "black")

Calculamos el coeficiente de correlación de Pearson para cada región, lo que nos permitirá entender la relación lineal entre la población y el número de homicidios.

murders %>% group_by(region) %>%
  summarise(cor(log10(population/10^6),log10(total),
                          method = "pearson"))
## # A tibble: 4 × 2
##   region        `cor(log10(population/10^6), log10(total), method = "pearson")`
##   <fct>                                                                   <dbl>
## 1 Northeast                                                               0.970
## 2 South                                                                   0.852
## 3 North Central                                                           0.949
## 4 West                                                                    0.915

Gráfico de Dispersión

A continuación, creamos un gráfico de dispersión para visualizar la relación entre la población (en escala logarítmica) y el total de homicidios (también en escala logarítmica), segmentando por región.

#Gráfico de Dispersión
murders %>% ggplot(aes(x=log10(population/10^6),y=log10(total),
                       color=region,shape=region))+
  geom_point(show.legend = FALSE)+facet_wrap(~region,nrow=2)

Filtrado del Dataset

Filtramos los datos para enfocarnos en la región del Noreste y preparamos los datos para el análisis de regresión.

northeast <- murders %>% filter(region=="Northeast") %>%
  mutate(log10pop=log10(population/10^6),log10tot=log10(total)) %>%
  select(population,log10pop,total,log10tot)
class(northeast)
## [1] "data.frame"

Modelo de Regresión Lineal

Ahora ajustamos un modelo de regresión lineal para predecir el número total de homicidios a partir del tamaño de la población.

#Modelo de regresión lineal
modelo<- northeast %>% select(log10pop,log10tot) %>%
  lm(log10tot~log10pop,data=.)

summary<-summary(modelo)

intercept <- modelo[["coefficients"]][["(Intercept)"]]
slope <- modelo[["coefficients"]][["log10pop"]]

Los coeficientes de la ecuación de regresión nos indican la relación entre las variables independientes y la variable dependiente. Por ejemplo, un coeficiente positivo indica una relación directa, mientras que un coeficiente negativo indica una relación inversa. Además, podemos evaluar la significancia estadística de los coeficientes utilizando pruebas de hipótesis y determinar si el modelo en su conjunto es significativo.

#Gráfico del modelo
graf <- northeast %>% ggplot(aes(x=log10pop,y=log10tot)) +
  geom_point() + ggtitle("Regresion linear model") +
  geom_smooth(method = "lm",color="blue")+
  annotate("text", x = 0.1, y = 3, 
           label = paste("y = ", round(slope, 2), "x +", round(intercept, 2)), 
           color = "blue")
ggplotly(graf)
## `geom_smooth()` using formula = 'y ~ x'

Este gráfico muestra el ajuste del modelo de regresión a los datos del Noreste, junto con la ecuación de la recta que representa la relación estimada entre la población y el total de homicidios.

Análisis de los Residuales

Los valores residuales de un análisis de regresión son las diferencias entre los valores observados del dataset y los valores estimados calculados con la ecuación de regresión.

ei <- residuals(modelo);ei
##           1           2           3           4           5           6 
##  0.27453592  0.01986875 -0.06282747 -0.31627887  0.05054933 -0.17841806 
##           7           8           9 
##  0.06274895  0.34499620 -0.19517475
# Calcular la media de los residuales
media_residuales <- mean(ei)

# Crear el gráfico de residuales
residuals_plot <- ggplot(data = data.frame(residuals = ei, fitted = modelo$fitted.values), 
                         aes(x = fitted, y = residuals)) +
  geom_point(color = "blue") +
  geom_hline(yintercept = media_residuales, color = "red", linetype = "dashed") +
  ggtitle("Gráfico de Residuales con Línea de la Media") +
  xlab("Valores Predichos") +
  ylab("Residuales") +
  coord_cartesian(ylim = c(-1, 1)) +
  theme_minimal()

# Mostrar el gráfico
residuals_plot

pred<-sort(fitted(modelo));pred
##         9         8         4         2         1         3         5         7 
## 0.4962047 0.8591238 1.0152489 1.0215239 1.7122358 2.1347095 2.3403858 2.5971672 
##         6 
## 2.8919086

Detección de Puntos Influyentes

Identificamos puntos influyentes mediante la distancia de Cook.

#Detección de puntos influyentes
#plot(cooks.distance(modelo))

cooks_d <- cooks.distance(modelo)
cooks_df <- data.frame(
  Index = seq_along(cooks_d),
  CooksDistance = cooks_d
)

cooks_plot <- ggplot(cooks_df, aes(x = Index, y = CooksDistance)) +
  geom_point(color = "blue") +  # Puntos de Cook's distance
  geom_hline(yintercept = 4 / length(cooks_d), color = "red", linetype = "dashed") +  # Línea de referencia
  ggtitle("Gráfico de Cook's Distance") +
  xlab("Índice del Observación") +
  ylab("Cook's Distance") +
  coord_cartesian(ylim = c(0, 0.5)) +
  theme_minimal()

# Mostrar el gráfico
cooks_plot

El umbral para la distancia de Cook es una forma de identificar observaciones influyentes en un modelo de regresión. Aunque no hay un valor universalmente aceptado para todos los contextos, el umbral comúnmente usado para detectar puntos influyentes es 4/n. Donde n es el número total de observaciones en el conjunto de datos. Este umbral ayuda a identificar observaciones cuyo Cook’s distance es suficientemente grande como para sugerir que pueden tener una influencia desproporcionada en el ajuste del modelo.

Validación del modelo

La regresión lineal se basa en varios supuestos fundamentales. Estos nos permiten obtener estimaciones precisas y confiables a partir de los datos.

Prueba de la Media de los Errores (Prueba t de Student)

Propósito: Verificar si los errores (residuales) tienen una media de cero, lo cual es un supuesto básico en la regresión lineal.

Hipótesis:

Hipótesis nula (Ho): La media de los errores es igual a cero.

Hipótesis alternativa (Ha): La media de los errores es diferente de cero.

t.test(ei)
## 
##  One Sample t-test
## 
## data:  ei
## t = -1.3876e-16, df = 8, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1664627  0.1664627
## sample estimates:
##     mean of x 
## -1.001682e-17

Si el p-valor de la prueba es mayor que un nivel de significancia de 0.05, no se rechaza la hipótesis nula, lo que sugiere que la media de los errores es cero.

Si el p-valor es menor que 0.05, se rechaza la hipótesis nula, lo que indica que los errores podrían tener una media distinta de cero, sugiriendo un posible sesgo en el modelo.

Prueba de Normalidad de los Errores (Prueba de Shapiro-Wilk n<30)

Propósito: Verificar si los errores siguen una distribución normal, lo cual es otro supuesto importante en la regresión lineal.

El Q-Q plot (quantile-quantile plot) es una herramienta útil para evaluar la normalidad de los residuales de un modelo de regresión. Donde:

Puntos cerca de la línea: Indica que los residuales siguen una distribución normal.

Puntos alejados de la línea: Sugiere que los residuales pueden no seguir una distribución normal, lo que podría afectar la validez de las inferencias del modelo.

# Q-Q Plot de los residuales
qq_plot <- ggplot(data = data.frame(residuals = ei), aes(sample = residuals)) +
  stat_qq() +
  stat_qq_line() +
  ggtitle("Q-Q Plot de los Residuales") +
  coord_cartesian(ylim = c(-0.5, 0.5)) +
  theme_minimal()

# Mostrar el gráfico
qq_plot

Hipótesis:

Hipótesis nula (Ho): Los errores siguen una distribución normal.

Hipótesis alternativa (Ha): Los errores no siguen una distribución normal.

shapiro.test(ei)
## 
##  Shapiro-Wilk normality test
## 
## data:  ei
## W = 0.96004, p-value = 0.7988

Si el p-valor es mayor que 0.05, no se rechaza la hipótesis nula, lo que indica que los errores pueden ser normalmente distribuidos.

Si el p-valor es menor que 0.05, se rechaza la hipótesis nula, sugiriendo que los errores no siguen una distribución normal, lo que podría afectar la validez de las inferencias realizadas a partir del modelo.

Prueba de Homocedasticidad (Prueba de Breusch-Pagan)

Propósito: Verificar si la varianza de los errores es constante a lo largo de todas las observaciones. Este supuesto se llama homocedasticidad.

Hipótesis:

Hipótesis nula (Ho): La varianza de los errores es constante (homocedasticidad).

Hipótesis alternativa (Ha): La varianza de los errores no es constante (heterocedasticidad).

bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 2.3144, df = 1, p-value = 0.1282

Si el p-valor es mayor que 0.05, no se rechaza la hipótesis nula, sugiriendo que hay homocedasticidad, es decir, la varianza de los errores es constante.

Si el p-valor es menor que 0.05, se rechaza la hipótesis nula, lo que indica la presencia de heterocedasticidad, lo que puede afectar la eficiencia de los estimadores del modelo.

Prueba de Independencia de los Errores (Prueba de Durbin-Watson)

Propósito: Evaluar si los errores son independientes entre sí, es decir, si no existe correlación entre los errores de observaciones consecutivas.

Hipótesis:

Hipótesis nula (Ho): Los errores son independientes.

Hipótesis alternativa (Ha): Los errores no son independientes.

dwtest(modelo,alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 2.0058, p-value = 0.864
## alternative hypothesis: true autocorrelation is not 0

La estadística de Durbin-Watson toma valores entre 0 y 4. Un valor cercano a 2 sugiere que no hay correlación entre los errores. Valores cercanos a 0 indican autocorrelación positiva, mientras que valores cercanos a 4 indican autocorrelación negativa.

Si el p-valor asociado es mayor que 0.05, no se rechaza la hipótesis nula, sugiriendo que los errores son independientes.

Si el p-valor es menor que 0.05, se rechaza la hipótesis nula, indicando que los errores pueden estar correlacionados, lo que puede invalidar las inferencias del modelo.

Resumen de la Importancia de las Pruebas:

Media de los Errores: Asegura que no hay sesgo sistemático en las predicciones.

Normalidad de los Errores: Garantiza la validez de los intervalos de confianza y pruebas de significancia.

Homocedasticidad: Asegura que las estimaciones de los coeficientes sean eficientes y que las inferencias sean válidas.

Independencia de los Errores: Previene la autocorrelación, que puede llevar a conclusiones engañosas.

Cumplir con estos supuestos es esencial para que el modelo de regresión sea robusto y confiable. Si alguno de estos supuestos no se cumple, puede ser necesario reconsiderar el modelo, transformar las variables, o aplicar técnicas alternativas.

Predicciones

Finalmente, realizamos una predicción del número total de homicidios para una población específica.

round(log10(20000000/10^6),2)
## [1] 1.3
pred1 <- predict(object = modelo, newdata = data.frame(
  log10pop=log10(20000000/10^6)));pred1
##        1 
## 2.913953
total=10^(pred1);total
##        1 
## 820.2624

Hemos predicho el número total de homicidios para una población de 20 millones utilizando el modelo ajustado.

Conclusión

Este análisis nos ha permitido explorar la relación entre la población y el número total de homicidios en la región del Noreste de Estados Unidos. A través del modelo de regresión lineal y las pruebas de validación, hemos podido confirmar la robustez del modelo y realizar predicciones informadas.