Sesión 3 y 4

Analítica de datos.

Sesión 3 Modelo de regresión lineal

En esta sesión aprenderemos los aspectos claves del modelo de regresión lineal simple. Repasaremos el análisis de correlación, ejecutaremos correlogramas y estimaremos un modelo de regresión lineal simple.

Veremos a continuación los comandos:

Primero, instalamos las librerías necesarias:

install.packages("tidyverse") 
install.packages("boot") 
install.packages("car") 
install.packages("QuantPsyc") 
install.packages("corrplot")

Ahora, las cargamos:

library(haven)
library(pacman)
library(dplyr)
library(ggplot2)  
library(ggpubr)
library(boot) 
library(car) 
library(corrplot) 
library(QuantPsyc)
library(corrplot) 
library(ggcorrplot) 
p_load(tidyverse, haven)

Importamos la base de datos

Creamos un subconjunto de datos con las variables de interés

new_dataset <- datos_IDERE %>% 
  dplyr::select(nombre, pais, homicidios, alfabetismo, est_educ_sup, 
         matr_educ_ini, anos_educ, mortalidad_inf, cobertura_salud, 
         esperanza_vida, pobreza, indigencia, trab_informal, desempleo_juv,
         gini, pib_pc, desempleo, partic_elecciones, corrupcion, ing_hog_persona)

Convertimos las variables a “númericas” para seguir los cálculos:

class(new_dataset$est_educ_sup) #para mirar que tipo de variable

Convertimos todo a numérico EXCEPTO nombre y pais

new_dataset <- new_dataset %>%
  mutate(across(-c(pais, nombre), as.numeric))

glimpse(new_dataset) #para confirmar el cambio

Identificamos valores nulos NA, (missing values), a veces aparecen como “NA” o “.”**

Comando para contar cuántos datos nulos tenemos por columna

sapply(new_dataset, function(x) 
  sum(is.na(x))) 

Borramos las filas que contienen datos nulos. Genera una base de datos “limpia”

new_dataset2<- na.omit(new_dataset) 

Comando para saber cuántos casos completos o observaciones tenemos:

sum(complete.cases(new_dataset2)) 

Analizamos las estadísticas descriptivas de los datos:

summary(new_dataset2) 

Ahora si los pasos para estimar la matriz de correlaciones y el modelo de regresión lineal

1. Calculamos la matriz de correlaciones para seleccionar las variables que podrían explicar nuestra variable dependiente “Y”. Queremos entender qué variables están correlacionadas con la variable “homicidio”

# Excluimos las dos variables de texto por su nombre
variables.cor <- cor(new_dataset2[, !(names(new_dataset) %in% c("pais", "nombre"))], 
               method = "pearson")
round(variables.cor, digits = 2)

Representación gráfica de la matriz de correlación

opción 1

corrplot(variables.cor)

opción 2

ggcorrplot(variables.cor, 
           method = 'circle', 
           type = 'lower', 
           lab = FALSE) + # lab = TRUE si quieres ver los números dentro
  ggtitle("Correlograma de variables con tasa de homicidios") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

Identificamos las variables que tienen alta correlación con la variable. Por ejemplo, seleccionamos aquellas variables con coef de corr > .4 : alfabetismo, est_educ_sup, cobertura_salud, esperanza_vida, corrupción

Pasos para construir un scatter con el coef. correlación Pearson. Calculamos la correlación entre las variables

ggplot(new_dataset2, aes(x = homicidios, y = alfabetismo)) + 
  geom_point(alpha = 0.5, color = "darkblue") + 
  geom_smooth(method = "lm", color = "red") + 
  stat_cor(method = "pearson", 
           label.x = 3,    # Coordenada X donde quieres el texto
           label.y = 100,  # Coordenada Y donde quieres el texto
           p.accuracy = 0.001, 
           r.accuracy = 0.01) + 
  theme_minimal()

Estimación de modelo de regresión

3. Modelo de regresión lineal simple

modelo_simple=lm(homicidios ~ alfabetismo , data= new_dataset2, na.action = na.exclude)  
summary(modelo_simple)

Interpretación del coeficiente B de la variable estatura:

En este modelo de regresión lineal simple, el coeficiente \(b\) (la pendiente) nos indica cuánto cambia la variable dependiente (Y) por cada unidad que aumenta la variable independiente (x).

La interpretación directa de \(b\) = -2.73, es: por cada punto porcentual adicional que aumente la tasa de alfabetismo (\(x\)), se espera que la tasa de homicidios (\(y\)) disminuya en 2.73 casos por cada 100,000 habitantes, manteniendo constantes los demás factores.

Es importante tener en cuenta las Unidades de Medida:

\(x\) : Puntos porcentuales de alfabetismo.

\(y\) : Homicidios por cada 100,000 personas.

Recordar que correlación no es Causalidad:

Aunque el coeficiente sea de -2.73, esto no garantiza que “enseñar a leer a alguien cause que no cometa un crimen”. El alfabetismo suele estar ligado a mejores ingresos, más oportunidades y mayor presencia del Estado, que son los factores que realmente reducen la violencia.

Guardamos la base de datos para continuar en la próxima sesión y poder compartirla

saveRDS(new_dataset2, "new_dataset2.rds")

Para ver dónde quedó guardado

getwd()

---

Sesión 4

Post-estimación: Verificación de cumplimiento de los supuestos del modelo de regresión lineal MCO

Retomamos la base de datos de la sesión anterior. Disponible en BN/sesión 4

nombre_de_objeto <- readRDS("ubicación del archivo/new_dataset2.rds")

1. Estimamos dos modelos para compararlos: modelo 1 y modelo 2

Modelo 1

modelo_1=lm(homicidios ~ alfabetismo + esperanza_vida + corrupcion + trab_informal,             data= new_dataset2, na.action = na.exclude) 
summary(modelo_1)

Modelo 2

modelo_2=lm(homicidios ~ alfabetismo + cobertura_salud, 
            data= new_dataset2, na.action = na.exclude) 
summary(modelo_2)

1) Gráfico para Diagnóstico de Linealidad

par(mfrow=c(1,2))  
plot(modelo_1, 1, caption = "Modelo 1")  
plot(modelo_2, 1, caption = "Modelo 2")

Este gráfico sugiere elegir el modelo 2 porque su linea roja es menos “errática”.

Para confirmar la normalidad de los errores podemos aplicar la prueba Shapiro-Wilks. En esta prueba se tienen las siguientes hipótesis.

H0 : los residuales tienen distribución normal.

H1: los residuales NO tienen distribución normal.

ei <- residuals(modelo_2) 
shapiro.test(ei)

Interpretación: el p-value es mucho menor a 0.05 (5% de significancia). Se rechaza la hipótesis nula, es decir, los residuos no son normales.

Una opción para reparar es la transformación logaritmica de la variable Y. Como estamos analizando la tasa de homicidios, aplicaremos el truco experto de sumar 1 (+ 1) para evitar errores si algún valor es cero (ya que el logaritmo de cero es indefinido). Corremos el siguiente modelo:

modelo_log=lm(log(homicidios+1) ~ alfabetismo + esperanza_vida + corrupcion + trab_informal,             data= new_dataset2, na.action = na.exclude) 
summary(modelo_log)

Aplicar la prueba Shapiro-Wilks al modelo transformado

ei <- residuals(modelo_log) 
shapiro.test(ei)

En estadística aplicada y ciencias sociales, un valor de 0.018 a menudo se considera “manejable”, especialmente si la muestra es grande (\(N > 30\) o \(50\)).

2. Gráficos para Diagnóstico de Normalidad de los residuos

par(mfrow=c(2,2))  
plot(modelo_1, 2, caption = "Modelo 1")  
plot(modelo_log, 2, caption = "Modelo Transformado log")  
plot(density(modelo_1$residuals))  
plot(density(modelo_log$residuals))

Se eligirá el modelo que se distribuya “parecido” a una distribución normal.

3) Diagnóstico de homocedasticidad

par(mfrow=c(1,2))  
plot(modelo_1, 3, caption = "Modelo 1")  
plot(modelo_log, 3, caption = "Modelo transformado")

Corrección de heterocedasticidad con errores robustos

install.packages("lmtest")
install.packages("sandwich")

library(lmtest)
library(sandwich)

Crear el modelo y aplicar la corrección

# 1. Ajustar el modelo de regresión lineal estándar (OLS)
modelo_robust=lm(homicidios ~ alfabetismo + esperanza_vida + corrupcion + trab_informal,             data= new_dataset2, na.action = na.exclude)

# 2. Calcular los errores estándar robustos (Tipo HC3 es el más recomendado)
# HC3 es robusto incluso con muestras pequeñas y presencia de valores atípicos.
resultados_robustos <- coeftest(modelo_robust, vcov = vcovHC(modelo_robust, type = "HC3"))

# 3. Ver los resultados
print(modelo_robust)

Comparación Visual (Opcional pero recomendado)

Es muy útil comparar los errores estándar originales vs. los robustos para ver cuánto “ruido” había en tus datos.

# Comparar Errores Estándar Normales vs Robustos
summary(modelo_robust)      # Mira la columna 'Std. Error'
resultados_robustos     # Compara con la nueva columna 'Std. Error'

C.4) Diagnóstico de Observaciones Influyentes: “Leverage” (palanca) y Valores Atípicos

par(mfrow=c(1,2))  plot(modelo1, 5, caption = "Modelo 1")  
plot(modelo2, 5, caption = "Modelo 2")

C.5) Predicción con el modelo de regresión

predecir_homicidios <- predict(modelo_log, 
data= body.dat[-c(1, 2, 3, 5, 6, 7, 9, 20, 21)]) 
summary(modelo2)

C.6) Para ver los coeficientes B estimados

modelo2$coefficients 
sqrt(mean(predecir_peso - body.dat$peso)^2)

C.7) Graficar los datos observados con los predecidos

par(mfrow=c(2, 2))  plot(modelo2)

D) ¿Qué hacer cuando nuestra variable dependiente Y no presenta una distribución normal?

hist(x = body.dat$peso, xlab = "Peso",       
     ylab = "Frecuencias",       main = "Histograma de Peso",       
     col = "blue", labels = T)

D.1) Tomamos logaritmo de nuestra variable Y

log_peso <- log(body.dat$peso)  
summary(log_peso)

hist(x = log_peso, xlab = "log Peso",       
     ylab = "Frecuencias",       main = "Histograma de Peso",       
     col = "blue", labels = T)

D.2) El siguiente código presenta en una tabla los modelos estimados como stata

install.packages("stargazer") 
library("stargazer") 

stargazer(body.dat, type = "text",           
title= "Table 1: Summary Statistics", out= "table1.txt")

stargazer(body.dat[c("peso", "estatura", "sexo")], type = "text",        
          title= "Table 1: Summary Statistics", out= "table1.txt")

stargazer(modelo1, modelo2, type="text", df= FALSE,           
title= "Table 2: Regresion results", out= "table1.txt") 

Trash

2) Gráfico boxplot para ver la distribución de los residuos

boxplot(modelo_1$residuals)   
boxplot(modelo_2$residuals)

Including Plots