Mapas Geo-Referenciados

Para este análisis se utilizará la base de datos de la Asignación a la Excelencia Pedagógica (AEP) del año 2014, un programa del Ministerio de Educación (MINEDUC) que consiste en un proceso de evaluación voluntario dirigido a docentes que desean mejorar su remuneración. Este proceso evalúa el desempeño disciplinario y pedagógico de los docentes mediante diversas pruebas y la elaboración de un portafolio. Con base en estos resultados, el docente podría acreditar un Tramo (I, II o III), lo cual se traduciría en un beneficio económico otorgado directamente por el MINEDUC.

La base de datos está compuesta por 2407 observaciones de 24 variables, de donde es de interés estudiar los promedios de los Puntajes en Prueba y Puntaje en Portafolio por comuna y por región en Chile, y además poder visualizar qué región o comuna Acredita más.

Tratamiento de la base de datos

#Librerias
library(readxl)
library(dplyr)
library(ggplot2)
library(leaflet)
library(tidyverse)
library(sf)
library(tidygeocoder)
library(mapview)
library(biscale)
library(cowplot)
library(rmarkdown)

#Base completa
base_aep <- read_delim("20150810_AEP_2014_20150709_PUBL.csv", delim = ";", escape_double = FALSE, trim_ws = TRUE)

Analizamos si posee datos NA visualmente:

library(naniar)
vis_miss(base_aep)

No hay datos NA, por lo que continuamos con esta base de datos.

Ahora, con el fin de obtener los nombres de las comunas, regiones y provincias adecuado junto a su respectiva codificación, para nuestro análisis se hará uso de la base de datos del CENSO 2017, de donde extraemos las columnas de nombres y codificaciones:

#library(censo2017)
#censo_descargar()

# Extraemos columnas de nombres y codificaciones
#ncomunas<-censo_tabla("comunas");ncomunas
#nprov<-censo_tabla("provincias");nprov<-data.frame(nprov[1],nprov[2], nprov[3],nprov[5])
#nreg<-censo_tabla("regiones");nreg<-data.frame(nreg[1],nreg[5])

Ya que descargar y cargar la base de datos completa del CENSO del 2017 toma mucho tiempo, guardamos las tablas de interés en formato “csv”, para ser llamadas de manera más eficiente en el actual análisis.

#write.csv(ncomunas, "comunas.csv", row.names = FALSE)
#write.csv(nprovincias, "provincias.csv", row.names = FALSE)
#write.csv(nregiones, "regiones.csv", row.names = FALSE)

library(readr)
ncomunas <- read_csv("ncomunas.csv")
nprov <- read_csv("nprov.csv")
nreg <- read_csv("nreg.csv")

Combinamos las provincias, comunas y regiones utilizando la respectiva columna de enlace (provincia_ref_id):

comunas<-merge(nprov,ncomunas,by.x="provincia_ref_id",by.y="provincia_ref_id")
comunas<-merge(nreg,comunas,by.x="region_ref_id",by.y="region_ref_id")
names(comunas)[8]=c("codigo_comuna")

Ahora, utilizando la libreria “chilemapas” obtenemos las coordenadas geográficas de cada comuna:

library(chilemapas)
mapacomunas = mapa_comunas

Procedemos a guardar las posibles variables de interés

X <- data.frame(base_aep[12],base_aep[13],base_aep[11], base_aep[10], base_aep[20],base_aep[21],base_aep[24],base_aep[8]) #12 es COD_COM_RBD
names(X)[1]=c("codigo_comuna")

Los codigos comuna, región y provincia están codificados como: “01”, “02”,… , mientras que la codificacion de enlace de nuestra base de datos es “1”, “2”,…, por lo que para que queden iguales se procede a transformarlos en números y luego a character.

mapacomunas$codigo_comuna <- as.numeric(mapacomunas$codigo_comuna)
mapacomunas$codigo_comuna <- as.character(mapacomunas$codigo_comuna)
mapacomunas$codigo_region <- as.numeric(mapacomunas$codigo_region)
mapacomunas$codigo_region <- as.character(mapacomunas$codigo_region)
mapacomunas$codigo_provincia <- as.numeric(mapacomunas$codigo_provincia)
mapacomunas$codigo_provincia <- as.character(mapacomunas$codigo_provincia)
X$codigo_comuna <- as.character(X$codigo_comuna)

Luego combinamos las posibles variables de interés con sus respectivas coordenadas geográficas de su comuna, usando el “codigo_comuna” como enlace.

base_combinada2 <- left_join(X, mapacomunas, by = "codigo_comuna")
table(base_combinada2$NOM_REG_RBD)

## 
##                               DE ANTOFAGASTA 
##                                           24 
##                        DE ARICA Y PARINACOTA 
##                                           40 
##                                   DE ATACAMA 
##                                           24 
## DE AYSÉN DEL GENERAL CARLOS IBAÑEZ DEL CAMPO 
##                                            9 
##                                  DE COQUIMBO 
##                                          120 
##                              DE LA ARAUCANÍA 
##                                          190 
##                                 DE LOS LAGOS 
##                                          131 
##                                  DE LOS RÍOS 
##                                           49 
##      DE MAGALLANES Y DE LA ANTÁRTICA CHILENA 
##                                           30 
##                                  DE TARAPACÁ 
##                                           46 
##                                DE VALPARAÍSO 
##                                          247 
##                                   DEL BIOBÍO 
##                                          398 
##             DEL LIBERTADOR BERNARDO OHIGGINS 
##                                          132 
##                                    DEL MAULE 
##                                          175 
##                    METROPOLITANA DE SANTIAGO 
##                                          792

Analizamos si existen datos NA

vis_miss(base_combinada2)

datos_con_na <- base_combinada2 %>% filter(is.na(codigo_region))
table(datos_con_na$NOM_REG_RBD) #los datos NA corresponden a un 30,4% de la region del biobio

## 
## DEL BIOBÍO 
##        121

121/398

## [1] 0.3040201

Podemos observar que un 30,4% de los datos de la región del Bio Bio (121 de 398) se vuelven NA al combinar la base. Esto se debe a que la librería chile_mapas no contiene las coordenadas geográficas de algunas de las comunas del Bio Bio. Con el fin de continuar con el estudio se utilizará la base sin estos datos.

base_combinada3 <- na.omit(base_combinada2)

Ahora procederemos a agregar su respectiva región a cada comuna de la base:

names(nreg)[1]=c("codigo_region")
nreg$codigo_region <- as.character(nreg$codigo_region)

base_combinada4 <- left_join(base_combinada3, nreg, by = "codigo_region")
vis_miss(base_combinada4)

names(base_combinada4)[2]=c("nombre_comuna")
base_combinada4[4] = base_combinada4[12]
base_combinada4=base_combinada4[1:11]
names(base_combinada4)[4]=c("nombre_region")

base_combinada4 será la base a utilizar para todos los mapas que realizaremos

Promedios en Puntajes por Comuna

Ahora, calcularemos los promedios de Puntajes por comunas

base_regiones <- base_combinada4 %>%
  group_by(codigo_region, nombre_region, codigo_comuna, nombre_comuna, geometry) %>%
  summarize(
    promedio_PTJE_PF = mean(PTJE_PF, na.rm = TRUE),
    promedio_PTJE_PRUEBA = mean(PTJE_PRUEBA, na.rm = TRUE),
    n_observaciones = n()
  ) %>%
  ungroup()

base_regiones <- st_as_sf(base_regiones)
base_regiones <- st_transform(base_regiones, crs = 4326)

Mapa geo-referenciado de los promedios de Puntaje en Portafolio por Comuna

labels4 <- sprintf("<strong>%s</strong> <br/> <strong>Comuna:</strong> %s <br/> <strong>Promedio Puntaje Portafolio:</strong> %g <br/><strong>Número de Observaciones:</strong> %g<br/>"
                  ,base_regiones$nombre_region,base_regiones$nombre_comuna, round(base_regiones$promedio_PTJE_PF,2), base_regiones$n_observaciones) %>% lapply(htmltools::HTML)

leaflet(data = base_regiones) %>%
  addProviderTiles(provider = providers$OpenStreetMap.Mapnik)%>%
  addTiles() %>%
  addPolygons(
    fillColor = ~colorNumeric("RdYlGn", promedio_PTJE_PF)(promedio_PTJE_PF),
    color = "black",
    weight = 1,
    fillOpacity = 0.7,
    label = labels4,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "15px",
      direction = "auto"
    )
  ) %>%
  addLegend(
    pal = colorNumeric("RdYlGn", domain = base_regiones$promedio_PTJE_PF),
    values = ~promedio_PTJE_PF,
    opacity = 0.7,
    title = "Promedio Puntaje Portafolio<br/> (Por Comuna)",
    position = "bottomright"
  ) %>%
  setView(lng = -70.6483, lat = -33.4569, zoom = 10)

Mapa geo-referenciado de los promedios de Puntaje en Prueba por Comuna

labels5 <- sprintf("<strong> %s </strong> <br/> <strong>Comuna:</strong> %s <br/> <strong>Promedio Puntaje Prueba:</strong> %g <br/><strong>Número de Observaciones:</strong> %g<br/>"
                   ,base_regiones$nombre_region,base_regiones$nombre_comuna, round(base_regiones$promedio_PTJE_PRUEBA,2), base_regiones$n_observaciones) %>% lapply(htmltools::HTML)

leaflet(data = base_regiones) %>%
  addProviderTiles(provider = providers$OpenStreetMap.Mapnik)%>%
  addTiles() %>%
  addPolygons(
    fillColor = ~colorNumeric("RdYlGn", promedio_PTJE_PRUEBA)(promedio_PTJE_PRUEBA),
    color = "black",
    weight = 1,
    fillOpacity = 0.7,
    label = labels5,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "15px",
      direction = "auto"
    )
  ) %>%
  addLegend(
    pal = colorNumeric("RdYlGn", domain = base_regiones$promedio_PTJE_PRUEBA),
    values = ~promedio_PTJE_PRUEBA,
    opacity = 0.7,
    title = "Promedio Puntaje Prueba<br/> (Por Comuna)",
    position = "bottomright"
  ) %>%
  setView(lng = -70.6483, lat = -33.4569, zoom = 10)

Promedios en Puntaje por Región

Ahora, con el fin de obtener los Promedios Regionales juntaremos las geometrías (coordenadas geográficas) de cada comuna con su respectiva Región

############### POR REGION ##################
##### Juntando las comunas para que quede una region:
geometria_por_region <- base_combinada4 %>%
  group_by(nombre_region) %>%
  summarize(geometry = st_union(geometry), .groups = 'drop')

Se procede a calcular los Promedios Regionales

promedios_por_region <- base_combinada4 %>%
  group_by(codigo_region, nombre_region) %>%
  summarize(
    promedio_PTJE_PF = mean(PTJE_PF, na.rm = TRUE),
    promedio_PTJE_PRUEBA = mean(PTJE_PRUEBA, na.rm = TRUE),
    n_observaciones = n()
  ) %>%
  ungroup()

Uniendo las geometrías de las regiones con los promedios calculados

region_con_promedios <- geometria_por_region %>%
  left_join(promedios_por_region, by = "nombre_region")

region_con_promedios <- st_as_sf(region_con_promedios)
region_con_promedios <- st_transform(region_con_promedios, crs = 4326)

Mapa geo-referenciado de los promedios de Puntaje en Portafolio por Región

labels7<- sprintf("<strong> %s </strong> <br/> <strong>Promedio Puntaje Prueba:</strong> %g <br/><strong>Número de Observaciones:</strong> %g<br/>"
                  ,region_con_promedios$nombre_region,round(region_con_promedios$promedio_PTJE_PF,2), region_con_promedios$n_observaciones) %>% lapply(htmltools::HTML)

leaflet(data = region_con_promedios) %>%
  addProviderTiles(provider = providers$OpenStreetMap.Mapnik)%>%
  addTiles() %>%
  addPolygons(
    fillColor = ~colorNumeric("RdYlGn", region_con_promedios$promedio_PTJE_PF)(region_con_promedios$promedio_PTJE_PF),
    color = "black",
    weight = 1,
    fillOpacity = 0.7,
    label = labels7,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "15px",
      direction = "auto"
    )
  ) %>%
  addLegend(
    pal = colorNumeric("RdYlGn", domain = region_con_promedios$promedio_PTJE_PF),
    values = ~promedio_PTJE_PF,
    opacity = 0.7,
    title = "Promedio Puntaje Portafolio <br/> (Por Region)",
    position = "bottomright"
  ) %>%
  setView(lng = -70.6483, lat = -33.4569, zoom = 10)

Mapa geo-referenciado de los promedios de Puntaje en Prueba por Región

labels6<- sprintf("<strong> %s </strong> <br/> <strong>Promedio Puntaje Prueba:</strong> %g <br/><strong>Número de Observaciones:</strong> %g<br/>"
                   ,region_con_promedios$nombre_region,round(region_con_promedios$promedio_PTJE_PRUEBA,2), region_con_promedios$n_observaciones) %>% lapply(htmltools::HTML)

leaflet(data = region_con_promedios) %>%
  addProviderTiles(provider = providers$OpenStreetMap.Mapnik)%>%
  addTiles() %>%
  addPolygons(
    fillColor = ~colorNumeric("RdYlGn", region_con_promedios$promedio_PTJE_PRUEBA)(region_con_promedios$promedio_PTJE_PRUEBA),
    color = "black",
    weight = 1,
    fillOpacity = 0.7,
    label = labels6,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "15px",
      direction = "auto"
    )
  ) %>%
  addLegend(
    pal = colorNumeric("RdYlGn", domain = region_con_promedios$promedio_PTJE_PRUEBA),
    values = ~promedio_PTJE_PRUEBA,
    opacity = 0.7,
    title = "Promedio Puntaje Prueba <br/> (Por Region)",
    position = "bottomright"
  ) %>%
  setView(lng = -70.6483, lat = -33.4569, zoom = 10)

Acreditan

Con el objetivo de observar el desempeño en términos de acreditación en general y simplificar la visualización en el mapa interactivo se creará una nueva variable de agrupación con dos categorías, “Acreditan” (que incluye tramos I, II, y III) y “No Acreditan” (que agrupa “No Acredita” y “No Acredita por copia”), lo cual es en base a la variable “TRAMO_ACRE”. Esto permitirá resaltar de manera más efectiva la diferencia entre quienes logran acreditar y quienes no.

table(base_aep$TRAMO_ACRE)

## 
##      ACREDITA TRAMO I     ACREDITA TRAMO II    ACREDITA TRAMO III 
##                   237                   751                   735 
##           NO ACREDITA NO ACREDITA POR COPIA 
##                   572                   112

Creación de la nueva variable de acreditación

base_combinada5 <- st_as_sf(base_combinada4)
base_combinada5 <- st_transform(base_combinada5, crs = 4326)

base_combinada5 <- base_combinada5 %>%
  mutate(acreditacion = case_when(
    TRAMO_ACRE %in% c("ACREDITA TRAMO I", "ACREDITA TRAMO II", "ACREDITA TRAMO III") ~ "Acreditan",
    TRAMO_ACRE %in% c("NO ACREDITA", "NO ACREDITA POR COPIA") ~ "No Acreditan",
    TRUE ~ NA_character_  # Para manejar otros casos que podrían aparecer
  ))

Calculando el porcentaje de acreditación por comuna y region

base_combinada_porcentaje <- base_combinada5 %>%
  group_by(codigo_comuna, nombre_comuna, nombre_region) %>%
  summarise(
    total = n(),
    acreditan = sum(acreditacion == "Acreditan", na.rm = TRUE),
    no_acreditan = sum(acreditacion == "No Acreditan", na.rm = TRUE),
    porcentaje_acreditacion = (acreditan / total) * 100
  )

Mapa geo-referenciado del porcentaje de Acreditación por Comuna

leaflet(base_combinada_porcentaje) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor = ~colorNumeric(palette = "YlGn", domain = base_combinada_porcentaje$porcentaje_acreditacion)(porcentaje_acreditacion),
    fillOpacity = 0.7,
    color = "black",
    weight = 1,
    label =  sprintf("<strong> %s </strong> <br/> <strong> Comuna: </strong> %s <br/> <strong> Porcentaje de Acreditación:</strong> %g <br/> <strong>Docentes que Acreditan:</strong> %g<br/> <strong>Docentes que No Acreditan:</strong> %g <br/>"
                     ,base_combinada_porcentaje$nombre_region,base_combinada_porcentaje$nombre_comuna,base_combinada_porcentaje$porcentaje_acreditacion, base_combinada_porcentaje$acreditan, base_combinada_porcentaje$no_acreditan) %>% lapply(htmltools::HTML),
    labelOptions = labelOptions(direction = "auto"),
    highlightOptions = highlightOptions(weight = 3, color = "blue", bringToFront = TRUE)
  ) %>%
  addLegend(
    pal = colorNumeric(palette = "YlGn", domain = base_combinada_porcentaje$porcentaje_acreditacion),
    values = ~porcentaje_acreditacion,
    title = "Porcentaje de Acreditación",
    position = "bottomright"
  )

Se calcula el porcentaje de acreditación por Región

base_combinada_porcentaje2 <- base_combinada5 %>%
  group_by(codigo_region, nombre_region) %>%
  summarise(
    total = n(),
    acreditan = sum(acreditacion == "Acreditan", na.rm = TRUE),
    no_acreditan = sum(acreditacion == "No Acreditan", na.rm = TRUE),
    porcentaje_acreditacion = (acreditan / total) * 100
  )  %>%
  ungroup()