# =========================================================
# Sesión 8: Modelo de Regresión lineal
# =========================================================

# ==============================
# 1. PAQUETES
# ==============================

install.packages(c("tidyverse","haven","janitor","skimr","moments","lmtest","sandwich"))

library(tidyverse)
library(haven)
library(janitor)
library(skimr)
library(moments)
library(lmtest)
library(sandwich)

# ==============================
# 2. CARGA DE DATOS
# ==============================

base_enemdu <- read_excel("G:/Mi unidad/3. Trabajos/8. Senamhi  Análisis de datos en R/4. Curso R/Módulo 2/Ejercicios 6/base_simulada_enemdu.xlsx")

# ==============================
# 3. ESTADÍSTICOS DESCRIPTIVOS
# ==============================

moda <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

# ==============================
#  3.1 TENDENCIA CENTRAL
#==============================
  
  mean(base_enemdu$ingreso) # Media
  median(base_enemdu$ingreso, na.rm = TRUE) # Mediana
  moda(base_enemdu$ingreso) # Moda
  
  
  table(base_enemdu$sexo)
  
  base_enemdu$sexo <- factor(base_enemdu$sexo,
                        levels = c(1, 2),
                        labels = c("Hombre", "Mujer"))
  
#==============================
#  3.2 DISPERSIÓN
#==============================
  
    sd(base_enemdu$ingreso) # Desviación estándar
    var(base_enemdu$ingreso) # Varianza
    range(base_enemdu$ingreso) # Rango
    IQR(base_enemdu$ingreso) # IQR

#==============================
#  3.3 FORMA
#==============================
    
    skewness(base_enemdu$ingreso) # Asimetría
    kurtosis(base_enemdu$ingreso) # Curtosis

    #==============================
    #  3.4 Resúmenes estadísticos
    #==============================
    
    summary(base_enemdu$ingreso)
    skim(base_enemdu$ingreso)

#==============================
#  4. ANALISIS POR GRUPOS
#==============================
  
    # Ingreso por sexo (ponderado)
    base_enemdu %>%
    group_by(sexo) %>%
    summarise(
      media = weighted.mean(ingreso, fexp, na.rm = TRUE),
      mediana = median(ingreso, na.rm = TRUE),
      sd = sd(ingreso, na.rm = TRUE)
    )

   # Ingreso por educación
  
      base_enemdu %>%
    group_by(educacion) %>%
    summarise(
      media = weighted.mean(ingreso, fexp, na.rm = TRUE),
      mediana = median(ingreso, na.rm=TRUE),
      sd = sd(ingreso, na.rm = TRUE)
    )
    
# ==============================
#  6. MODELO MCO (BRECHA SALARIAL)
# ==============================
  
  base_enemdu <- base_enemdu %>%
  mutate(ln_ingreso = log(ingreso))

  # Modelo 1 - Brecha Salarial
  modelo1 <- lm(ln_ingreso ~ sexo, data = base_enemdu)
  summary(modelo1)
  
  modelo2 <- lm(ln_ingreso ~ sexo + horas, data = base_enemdu)
  summary(modelo2)
  
  
  modelo3 <- lm(ln_ingreso ~ sexo + horas + experiencia, data = base_enemdu)
  summary(modelo3)
  
  modelo4 <- lm(ln_ingreso ~ sexo + horas + experiencia + edad, data = base_enemdu)
  summary(modelo4)
  
  
  
  cor(base_enemdu$experiencia, base_enemdu$edad)
  
  ggplot(base_enemdu, aes(x = experiencia, y = edad)) +
    geom_point() +
    labs(title = "Relación entre experiencia y edad",
         x = "Experiencia (años)",
         y = "Edad") +
    theme_minimal()