01_preprocessing_data.Rmd

---
title: "Preprocessing Data"
author:
  - name: Brian Peña-Calero
    email: brian.pena@upch.pe
    affiliation: Innovar, UPCH
date: "`r Sys.Date()`"
output: 
  html_document: 
    toc: yes
    toc_depth: 3
    toc_float: yes
    number_sections: yes
    code_folding: show
    code_download: yes
    theme: 
      bootswatch: flatly
    highlight: kate
    highlight_downlit: true
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  dpi = 300,
  fig.align = 'center'
  # message = FALSE,
  # warning = FALSE,
  # error = TRUE
)
```

# Carga de paquetes

```{r}
library(tidyverse)
library(sf)
```

# Importar datos

```{r}
salid1_pubsalid1 <- readr::read_csv("01_data/raw/SALID1_PUBSALID1.csv")
gdb_path <- "01_data/raw/Data Request MS242_09222023/MS242_L1.gdb"
gdb_salurbal <- st_layers(gdb_path)

cities_l1ad <- st_read(
  gdb_path,
  gdb_salurbal$name[[1]]
) %>%
  as_tibble() %>%
  st_as_sf() %>%
  left_join(
    salid1_pubsalid1
  ) %>%
  relocate(SALID1, .after = PUBSALID1)

cities_l1ad_noislands <- st_read(
  gdb_path,
  gdb_salurbal$name[[2]]
) %>%
  as_tibble() %>%
  st_as_sf() %>%
  left_join(
    salid1_pubsalid1
  ) %>%
  relocate(SALID1, .after = PUBSALID1)
```


## Datos Crecimiento Urbano

Los datos de crecimiento urbano tienen columna `ISO2` que identifica al paós, y una codificación propia llamada `SALID1`, con la cual se tiene una data que indica a que ciudades corresponden (`salid1_pubsalid1`).

```{r}
bec_backcast_l1ux <- read_csv("01_data/raw/BEC_BACKCAST_L1UX_20230302.csv")
```

Se adiciona la columna `SALID0` calculado como los 3 primeros dígitos de `SALID1`.

```{r}
bec_backcast_l1ux <- bec_backcast_l1ux %>%
  filter(ISO2 != "NI") %>% 
  mutate(
    SALID0 = str_sub(SALID1, 1, 3),
    .before = SALID1
  ) %>%
  left_join(
    salid1_pubsalid1
  ) %>%
  relocate(PUBSALID1, .before = SALID1)
```

## Datos Poblacionales

De forma similar en los datos poblaciones, se hace el match de `SALID1` para poder identificar a las ciudades que corresponde.

```{r}
population <- st_read(
  gdb_path,
  gdb_salurbal$name[[3]]
) %>%
  as_tibble() %>%
  left_join(
    salid1_pubsalid1
  ) %>%
  relocate(SALID1, .after = PUBSALID1) %>%
  drop_na(PUBSALID1)
```

```{r}
population <- population %>%
  filter(country != "nicaragua") %>% 
  mutate(
    ISO2 = case_match(
      country,
      "argentina" ~ "AR",
      "brazil" ~ "BR",
      "chile" ~ "CL",
      "colombia" ~ "CO",
      "costa rica" ~ "CR",
      "el salvador" ~ "SV",
      "guatemala" ~ "GT",
      "mexico" ~ "MX",
      "nicaragua" ~ "NI",
      "panama" ~ "PA",
      "peru" ~ "PE"
    ),
    .before = country
  ) %>%
  left_join(
    bec_backcast_l1ux %>%
      select(ISO2, SALID0) %>%
      distinct()
  ) %>%
  relocate(ISO2, SALID0) %>%
  arrange(SALID0)
population
```

Se aplica una pequeña corrección a los datos numéricos poblacionales que contienen caracteres como `-` y `" "` (espacios) para representar valores perdidos.

```{r}
population <- population %>%
  mutate(
    across(
      F1945:F2022,
      as.numeric
    )
  )

population
```

> Guardamos la información de las ciudades para futuras referencias

```{r}
city_info <- population %>%
  select(
    ISO2:City, PUBSALID1:SALID1
  )
```


```{r eval=FALSE}
saveRDS(
  city_info,
  file = "01_data/processed/city_info.rds"
)

#city_info <- readRDS("01_data/processed/city_info.rds")
```


# Imputación de datos

Debido a que hay una gran cantidad de valores perdidos a través de los años en distintas ciudades, se realiza un proceso de imputación de datos tanto para los de crecimiento urbano como datos poblacionales.

El análisis de imputación se realiza para cada ciudad, por lo que se usará el enfoque de ***tidy nested*** para organizar los datos de tal manera que cada fila sea una ciudad; y su columna, sea un conjunto de datos que contenga la información de interés a través de los años. Para el proceso de imputación se usará el paquete `imputeTS` wrappeado en la función `safely` (paquete `purrr`) para evitar la detención de errores debido a problemas de convergencia o insuficiencia de datos en algunas ciudades.

```{r}
library(imputeTS)
na_kalman_safe <- safely(.f = na_kalman)
```

## Crecimiento Urbano

Los datos de crecimiento urbano proporcionado tienen información cada 5 años, por lo que se usará la misma técnica de imputación para tener información anual. Para ello primero se completará los años faltantes en el data.frame:

```{r}
bectuareal_tidy <- bec_backcast_l1ux %>%
  distinct(PUBSALID1) %>%
  crossing(YEAR = 1985:2015) %>%
  left_join(
    bec_backcast_l1ux %>%
      select(PUBSALID1, YEAR, BECTUAREAL1UX)
  ) %>%
  janitor::clean_names()
```

Ahora se procede a crear el formato nest:

```{r}
bectuareal_imp <- bectuareal_tidy %>%
  group_nest(pubsalid1) %>%
  mutate(
    ts_data = map(
      data,
      ~ ts(.x$bectuareal1ux, 1985, 2015)
    ),
    kalman_data = map(
      ts_data,
      na_kalman_safe
    )
  )
```

En este caso también se tiene algunos errores (1)

```{r}
bectuareal_imp_errors <- bectuareal_imp %>%
  filter(map_lgl(.$kalman_data, ~ !is.null(.x$error)))
```

Aunque en la mayoría de ciudades se logra imputar (370 ciudades).

```{r}
bectuareal_imp_clean <- bectuareal_imp %>% 
  filter(map_lgl(.$kalman_data, ~ is.null(.x$error))) %>% 
  mutate(kalman_data = map(kalman_data, 'result'))
```

Finalmente, se añadirá la información sobre la cantidad original de crecimiento urbano existente para fines de futuras validaciones:

```{r}
bectuareal_imp_merge <- bectuareal_imp_clean %>%
  mutate(
    kalman_data = map(
      kalman_data,
      ~ as_tibble(.x) %>%
        mutate(
          year = 1985:2015
        )
    ),
    data_imputation = map2(
      data, kalman_data,
      ~ .x %>%
        left_join(.y,
          by = join_by(year)
        )
    )
  )

bectuareal_imp_merge <- bectuareal_imp_merge %>% 
  select(pubsalid1, data_imputation) %>% 
  unnest(cols = c(data_imputation)) %>% 
  rename(bectuareal1ux_imp = x)
```

## Datos poblacionales


Necesitamos ordenar los datos y darle un formato **nested**. Observamos que los datos poblacionales con los que se trabaja tiene 648 ciudades.

```{r}
population_nested <- population %>%
  janitor::clean_names() %>%
  rename_with(
    ~ str_replace(., "f", "x"),
    f1945:f2015
  ) %>%
  select(pubsalid1, x1945:x2015) %>%
  pivot_longer(
    cols = x1945:x2015,
    names_to = "year",
    values_to = "population"
  ) %>%
  group_nest(pubsalid1)
population_nested
```

Se tiene una función que verifica una desviación máxima de valores poblacionales en una ciudad en particular a partir de su mediana a lo largo de los años. Por defecto, se usa el criterio de 20% máximo de variación. En caso algún valor poblacional supere este valor, se ingresará un `NA` en lugar del valor reportado. Así, la imputación podrá también estimar estos valores *anómalos*.

```{r}
fix_population <- function(data, desviacion_max = 0.2) {
  data_wo_na <- na.omit(data)
  
  if (length(data_wo_na) < 2) {
    return(data)
  }
  
  median_population <- median(data_wo_na)
  
  upper_limit <- median_population * (1 + desviacion_max)
  lower_limit <- median_population * (1 - desviacion_max)
  
  data_fix <- ifelse(
    !is.na(data) & (data > upper_limit | data < lower_limit), 
    NA, 
    data
  )
  return(data_fix)
}
```

```{r}
population_imp <- population_nested %>%
  mutate(
    # Aplicación de la corrección
    pop_fix = map(
      data, 
      ~ fix_population(.x$population) %>%
        enframe(
          name = "year",
          value = "population"
        ) %>%
        mutate(
          year = paste0("x", 1945:2015)
        )
    ),
    # Transformar a objeto `time series`
    ts_pop_fix = map(
      pop_fix,
      ~ ts(.x$population, 1945, 2015) %>%
        log()
    ),
    # Aplicación de la imputación
    kalman_pop_fix = map(
      ts_pop_fix, 
      na_kalman_safe
    )
  )
```

Debido a que algunas ciudades tienen valores perdidos en todos los años de análisis, no en todos los casos la imputación es factible. En este caso hubo un total de *208* ciudades con problemas de imputación.

```{r}
population_imp_errors <- population_imp %>%
  filter(map_lgl(.$kalman_pop_fix, ~ !is.null(.x$error)))
```

En tanto, la cantidad de ciudades con las que se podrá trabajar con información poblacional será de *440*.

```{r}
population_imp_clean <- population_imp %>%
  filter(map_lgl(.$kalman_pop_fix, ~ is.null(.x$error))) %>%
  mutate(kalman_pop_fix = map(kalman_pop_fix, "result"))
```

Finalmente, se añadirá la información sobre la cantidad original de dato poblacional existente para fines de futuras validaciones:

```{r}
population_imp_merge <- population_imp_clean %>%
  mutate(
    kalman_pop_fix = map(
      kalman_pop_fix,
      ~ as_tibble(.x) %>%
        mutate(
          x = exp(x), 
          year = paste0("x", 1945:2015)
        )
    ),
    data_imputation = map2(
      pop_fix, kalman_pop_fix,
      ~ .x %>%
        left_join(.y,
          by = join_by(year)
        )
    )
  )

population_imp_merge <- population_imp_merge %>%
  select(pubsalid1, data_imputation) %>%
  unnest(cols = c(data_imputation)) %>%
  rename(population_imp = x) %>%
  mutate(
    year = str_remove(year, "x"),
    year = as.numeric(year)
  )

population_imp_merge
```


# Unificación de datos

En este punto, ya se tiene los datos poblacionales y de crecimiento urbano limpios y ordenados. Se realiza la unificación en base al código `salid1` que deberían tener ambos objetos. 

Se usa la función `full_join` para tener la posibilidad de analizar los datos poblacionales y de crecimiento juntos y separados.

```{r}
pop_bectuareal <- full_join(
  population_imp_merge,
  bectuareal_imp_merge,
  by = join_by(pubsalid1, year)
)
```

# Export data

```{r}
saveRDS(pop_bectuareal,  
        file = "01_data/processed/pop_bectuareal1ux.rds")
```