n_filas    <- nrow(datos)
n_codigos  <- sum(datos$ingreso %in% c(999999, -1), na.rm = TRUE)
n_blanco   <- sum(is.na(datos$ingreso))
n_ingreso  <- nrow(ing)
n_cero     <- sum(ing$ingreso == 0)

q <- quantile(ing$ingreso, c(0.25, 0.75))
lim_sup    <- q[2] + 1.5 * (q[2] - q[1])
n_atipicos <- sum(ing$ingreso > lim_sup)

summary(ing$ingreso)
# Control: n_filas - n_codigos - n_blanco = n_ingreso

# Media, mediana y CV: sin y con factor de expansión
estad <- tibble(
  medida = c("Media", "Mediana", "Desv. estándar", "CV"),
  sin_factor = c(mean(ing$ingreso), median(ing$ingreso),
                 sd(ing$ingreso), sd(ing$ingreso) / mean(ing$ingreso)),
  con_factor = c(media_pond(ing$ingreso, ing$factor),
                 mediana_pond(ing$ingreso, ing$factor),
                 sd_pond(ing$ingreso, ing$factor),
                 sd_pond(ing$ingreso, ing$factor) /
                   media_pond(ing$ingreso, ing$factor)))

# Correlación de Spearman (el ingreso es muy asimétrico)
cor(ing$instruccion, ing$ingreso, method = "spearman", use = "complete.obs")

# Probabilidad como proporción: ¿gana menos que el salario básico?
SBU <- 470
mean(ing$ingreso < SBU)
prop_pond(ing$ingreso < SBU, ing$factor)

# Intervalo de confianza al 95 % por grupo
ic_sexo <- ing |>
  group_by(sexo) |>
  summarise(n = n(),
            media = mean(ingreso),
            ee = sd(ingreso) / sqrt(n),
            li = media - qt(0.975, n - 1) * ee,
            ls = media + qt(0.975, n - 1) * ee)

# Prueba t de Welch (no usa el factor de expansión)
prueba <- t.test(ingreso ~ sexo, data = ing)
prueba$conf.int
prueba$p.value

set.seed(2026)

seg <- ing |>
  filter(ingreso > 0, !is.na(horas), horas > 0, horas <= 168) |>
  mutate(log_ingreso = log(ingreso))

km_datos <- seg |>
  select(edad, horas, log_ingreso) |>
  scale()

# Método del codo
codo <- tibble(
  k = 1:8,
  sct = map_dbl(1:8, \(k) kmeans(km_datos, k, nstart = 10)$tot.withinss))

ggplot(codo, aes(k, sct)) +
  geom_line() + geom_point() +
  scale_x_continuous(breaks = 1:8) +
  labs(title = "Método del codo",
       x = "Número de grupos (k)",
       y = "Variación dentro de los grupos",
       caption = "Fuente: INEC, ENEMDU anual 2025")

K <- 3   # tu decisión: codo + tamaño + sentido de los grupos
km <- kmeans(km_datos, centers = K, nstart = 25)
seg <- seg |> mutate(grupo = factor(km$cluster))

perfil <- seg |>
  group_by(grupo) |>
  summarise(n = n(),
            pct_pob = 100 * sum(factor) / sum(seg$factor),
            edad = median(edad),
            horas = median(horas),
            ingreso = median(ingreso),
            pct_mujer = 100 * mean(sexo == "Mujer"),
            pct_rural = 100 * mean(area == "Rural"))
perfil

# Estabilidad: ¿otra semilla da la misma partición?
set.seed(99)
km2 <- kmeans(km_datos, centers = K, nstart = 25)
table(km$cluster, km2$cluster)

serie_mes <- ing |>
  group_by(mes = as.numeric(mes)) |>
  summarise(n = n(),
            media = mean(ingreso),
            ee = sd(ingreso) / sqrt(n),
            li = media - qt(0.975, n - 1) * ee,
            ls = media + qt(0.975, n - 1) * ee,
            media_p = media_pond(ingreso, factor))

# Media móvil de 3 meses
serie_mes <- serie_mes |>
  mutate(mm3 = (media_p + lag(media_p) + lag(media_p, 2)) / 3)

# Tendencia lineal y su valor p
modelo <- lm(media_p ~ mes, data = serie_mes)
summary(modelo)

# ¿Qué meses se alejan de la media anual más de lo que explica el ruido?
media_anual <- mean(ing$ingreso)
serie_mes |>
  mutate(fuera = li > media_anual | ls < media_anual) |>
  select(mes, media, li, ls, fuera)

# Proyección de 3 meses con su rango (13 a 15 = ene a mar de 2026)
nuevos <- tibble(mes = 13:15)
predict(modelo, newdata = nuevos, interval = "prediction", level = 0.95)

g_bajo <- which.min(perfil$ingreso)    # perfil de menor ingreso mediano

kpis <- tibble(
  kpi = c("Ingreso mediano (con factor, USD)",
          "% que gana menos que el SBU (con factor)",
          "Brecha de ingreso hombre-mujer (USD)",
          "% de la población en el perfil de menor ingreso",
          "Tendencia mensual del ingreso (USD por mes)"),
  valor = c(mediana_pond(ing$ingreso, ing$factor),
            100 * prop_pond(ing$ingreso < SBU, ing$factor),
            unname(prueba$estimate[1] - prueba$estimate[2]),
            perfil$pct_pob[g_bajo],
            unname(coef(modelo)["mes"])))
print(kpis)