Sesión 3

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 %>% 
  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:

El coeficiente alfabetismo de 2.73 puntos indica que, manteniendo constantes las demás variables del modelo, por cada incremento de 1 centímetro en la estatura, se espera que el peso aumente en promedio 0.3785 kilogramos.

Por cada 1% que aumente el alfabetismo, la tasa de homicidios disminuye en 2.73 puntos.

Interpretación del coeficiente de sexo (-8.58). Recordar que la variable es dicotómica: 1 si hombre, 0 si mujer. El coeficiente de regresión B de sexo indica que, manteniendo constantes las demás variables del modelo, ….

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

boxplot(modelo2$residuals)  
boxplot(modelo1$residuals)

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

Compararemos los dos modelos estimados: modelo 1 y modelo 2

C.1) Gráfico para Diagnóstico de Linealidad

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

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(modelo1) shapiro.test(ei)

C.2) Gráficos para Diagnóstico de Normalidad

par(mfrow=c(2,2))  plot(modelo1, 2, caption = "Modelo 1")  
plot(modelo2, 2, caption = "Modelo 2")  
plot(density(modelo1$residuals))  
plot(density(modelo2$residuals))

C.3) Diagnóstico de homocedasticidad

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

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_peso <- predict(modelo2, 
                         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

Including Plots