library("WDI")
library("dplyr")
library("tidyr")
library("readr")
library("ggplot2")

# 1.2 Definir el indicador que queremos
indicadores <- c(
  pop65   = "SP.POP.65UP.TO.ZS",  # % población ≥65 (envejecimiento)
  ghpub   = "SH.XPD.GHED.CH.ZS",  # gasto público en salud (% del gasto en salud)
  hexp    = "SH.XPD.CHEX.GD.ZS",  # gasto en salud como % del PIB
  lifeexp = "SP.DYN.LE00.IN",     # esperanza de vida al nacer
  mortinf = "SH.DYN.MORT"         # mortalidad infantil menores de 5 años (por 1.000 nacidos vivos)
  )
  
#1.3 Definimos los países usando su código ISO de 2 letras
OCEANIA <- c("AU","NZ","FJ","PG")
AFRICA  <- c("ZA","NG","EG","KE","MA","ET","GH","TZ","SN","UG")
ASIA    <- c("JP","KR","CN","IN","TH","ID","PH","VN","MY","SG")
EUROPE  <- c("ES","FR","DE","IT","SE","NL","PL","PT","GB","NO")
LATAM   <- c("CL","AR","BR","MX","CO","PE","UY","VE","CR","EC")
paises <- c(OCEANIA, AFRICA, ASIA, EUROPE, LATAM)

# 1.4 Descargamos los datos desde 2005 hasta 2024
datos_bm <- WDI(
  indicator = indicadores,
  country = paises,
  start = 2005,
  end = 2024,
  extra     = TRUE    # agrega región, nivel de ingreso, etc.
)

#1.5 Utilizamos head(datos_bm) explorar y verificar que los datos estén en el formato data.frame.
head(datos_bm)

# 2.1) Quitar agregados y NAs de región
datos_bm2 <- datos_bm[!is.na(datos_bm$region) & datos_bm$region != "Aggregates", ]

# 2.2) Crear variable 'grupo' (continente) con un bucle fila a fila
grupo_vec <- character(nrow(datos_bm2))  # vector vacío para ir llenando
for (i in seq_len(nrow(datos_bm2))) {
  code <- datos_bm2$iso2c[i]
  if (code %in% OCEANIA) {
    grupo_vec[i] <- "OCEANIA"
  } else if (code %in% AFRICA) {
    grupo_vec[i] <- "AFRICA"
  } else if (code %in% ASIA) {
    grupo_vec[i] <- "ASIA"
  } else if (code %in% EUROPE) {
    grupo_vec[i] <- "EUROPE"
  } else if (code %in% LATAM) {
    grupo_vec[i] <- "LATAM"
  } else {
    grupo_vec[i] <- "OTRO"
  }
}
datos_bm2$grupo <- grupo_vec

# 2.3) Ordenar niveles de ingreso (para que los resúmenes salgan en orden lógico)
niveles_ing <- c("Low income","Lower middle income","Upper middle income","High income")
datos_bm2$income <- factor(datos_bm2$income, levels = niveles_ing)

# 2.4) Seleccionar columnas finales
df_wide <- datos_bm2 %>%
  select(country, iso2c, iso3c, year, grupo, region, income,
         pop65, ghpub, hexp, lifeexp, mortinf) %>%
  arrange(grupo, country, year)


# 2.5) Chequeo básico del panel
cat("Panel ancho -> filas:", nrow(df_wide),
    "| países:", length(unique(df_wide$country)),
    "| años:", min(df_wide$year, na.rm=TRUE), "-", max(df_wide$year, na.rm=TRUE), "\n")


# 2.6) Variables a apilar / medir
vars <- c("pop65","ghpub","hexp","lifeexp","mortinf")

# pasar a formato largo 
df_long_list <- list()
for (v in vars) {
  tmp <- df_wide[, c("country","iso2c","iso3c","year","grupo","income", v)]
  names(tmp)[7] <- "value"
  tmp$indicator <- v
  # tmp <- tmp[!is.na(tmp$value), ]   # <- BORRAR esta línea (no eliminar NA aquí)
  df_long_list[[v]] <- tmp
}
df_long <- do.call(rbind, df_long_list)
cat("\nPanel largo -> filas totales:", nrow(df_long),
    "| con dato:", sum(!is.na(df_long$value)), "\n")

# 2.7) Obtener % de NA por indicador
na_pct <- numeric(length(vars)); names(na_pct) <- vars
for (j in seq_along(vars)) {
  v <- vars[j]
  na_pct[j] <- mean(is.na(df_wide[[v]])) * 100
}
cat("\n% de NA por indicador (2005–2024):\n")
print(round(na_pct, 2))

# 2.8) Conteo de datos válidos por continente e indicador
grupos <- sort(unique(df_wide$grupo)); grupos <- grupos[!is.na(grupos)]
conteo <- data.frame(grupo = grupos, stringsAsFactors = FALSE)
for (v in vars) conteo[[v]] <- 0L

for (g in seq_along(grupos)) {
  gg <- grupos[g]
  sub <- df_wide[df_wide$grupo == gg, ]
  for (v in vars) {
    conteo[g, v] <- sum(!is.na(sub[[v]]))
  }
}
cat("\nCuenta de datos válidos por continente e indicador:\n")
print(conteo)

#3.1 Media y desviación estándar por continente de población, esperanza de vida y mortalidad infantil
desc_disp <- df_wide %>%
  group_by(grupo) %>%
  summarise(
    mean_pop65 = mean(pop65, na.rm=TRUE),
    sd_pop65   = sd(pop65, na.rm=TRUE),
    mean_life  = mean(lifeexp, na.rm=TRUE),
    sd_life    = sd(lifeexp, na.rm=TRUE),
    mean_mort  = mean(mortinf, na.rm=TRUE),
    sd_mort    = sd(mortinf, na.rm=TRUE),
    .groups="drop"
  )

cat("\n[3.2] Media y desviación estándar por continente de población, esperanza de vida y mortalidad infantil:\n")
print(desc_disp)

#3.2 Gasto público y gasto del PIB por continente
gasto_cont <- df_wide %>%
  group_by(grupo) %>%
  summarise(
    mean_hexp  = mean(hexp, na.rm=TRUE),   # gasto total en salud (% PIB)
    sd_hexp    = sd(hexp, na.rm=TRUE),
    mean_ghpub = mean(ghpub, na.rm=TRUE),  # % del gasto financiado por el Estado
    sd_ghpub   = sd(ghpub, na.rm=TRUE),
    .groups="drop"
  )

cat("\n[Descriptivo de gasto en salud por continente]\n")
print(gasto_cont)

# 1) Boxplot con el envejecimiento por continente en todos los años (2005–2024)
ggplot(df_wide, aes(x = grupo, y = pop65, fill = grupo)) +
  geom_boxplot(alpha = 0.7, outlier.alpha = 0.4) +
  labs(
    title = "Distribución del envejecimiento por continente (2005–2024)",
    x = "Continente",
    y = "% población ≥65"
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))
        
#2) Esperanza de vida y mortalidad infantil
# Esperanza de vida (media ± sd)
ggplot(desc_disp, aes(x = grupo, y = mean_life, fill = grupo)) +
  geom_col(alpha = 0.7) +
  geom_errorbar(aes(ymin = mean_life - sd_life, ymax = mean_life + sd_life),
                width = 0.2) +
  labs(
    title = "Esperanza de vida promedio por continente (2005–2024)",
    x = "Continente",
    y = "Esperanza de vida (años)"
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))

# Mortalidad infantil (media ± sd)
ggplot(desc_disp, aes(x = grupo, y = mean_mort, fill = grupo)) +
  geom_col(alpha = 0.7) +
  geom_errorbar(aes(ymin = mean_mort - sd_mort, ymax = mean_mort + sd_mort),
                width = 0.2) +
  labs(
    title = "Mortalidad infantil promedio por continente (2005–2024)",
    x = "Continente",
    y = "Muertes por 1.000 nacidos vivos"
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))

# 3) Resumen transversal de gasto (todos los años)
gasto_allyears <- df_wide %>%
  group_by(grupo) %>%
  summarise(
    mean_hexp  = mean(hexp,  na.rm=TRUE), sd_hexp  = sd(hexp,  na.rm=TRUE),  # % del PIB
    mean_ghpub = mean(ghpub, na.rm=TRUE), sd_ghpub = sd(ghpub, na.rm=TRUE), # % del gasto en salud
    .groups = "drop"
  )
# Largo + etiquetas para unidades
plot_gasto <- gasto_allyears %>%
  tidyr::pivot_longer(
    cols = c(mean_hexp, sd_hexp, mean_ghpub, sd_ghpub),
    names_to = c(".value","indicador"),
    names_pattern = "(mean|sd)_(.*)"
  ) %>%
  mutate(
    unidad = dplyr::case_when(
      indicador == "hexp"  ~ "% del PIB",
      indicador == "ghpub" ~ "% del gasto en salud"
    ),
    etiqueta = dplyr::case_when(
      indicador == "hexp"  ~ "Gasto en salud (% del PIB)",
      indicador == "ghpub" ~ "Financiamiento estatal (% del gasto en salud)"
    )
  )

# Barras + error bars (media ± sd) con TODA la serie
ggplot(plot_gasto, aes(x = grupo, y = mean, fill = grupo)) +
  geom_col(alpha = 0.7) +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2) +
  facet_wrap(~ etiqueta, scales = "free_y") +
  labs(
    title = "Gasto y financiamiento en salud (transversal 2005–2024)",
    subtitle = "Media y dispersión (sd) usando todas las observaciones país-año",
    x = "Continente", y = NULL
  ) +
  theme_minimal() +
  theme(legend.position = "none")