1. INTRODUCCIÓN

La educación es uno de los pilares fundamentales del desarrollo humano y social, al constituirse en un factor determinante para la reducción de la pobreza, la movilidad social y el crecimiento económico sostenible. Sin embargo, su acceso y calidad no se distribuyen de manera equitativa en el territorio, reflejando desigualdades estructurales que responden a dinámicas sociales, económicas y espaciales. En ciudades latinoamericanas como Guayaquil, caracterizadas por un acelerado crecimiento urbano y fuertes contrastes socioeconómicos, estas desigualdades educativas adquieren una dimensión territorial evidente.

El análisis espacial de la escolaridad permite comprender cómo se configuran los patrones de concentración o dispersión del nivel educativo en el espacio urbano, y cómo estos se asocian con otros determinantes como la infraestructura, la accesibilidad o las condiciones socioeconómicas de la población. La estadística espacial ofrece herramientas para identificar la existencia de autocorrelación o dependencia espacial en los datos, es decir, si los sectores con altos o bajos niveles educativos tienden a agruparse geográficamente.

En este contexto, el presente estudio aplica un enfoque de estadística espacial para examinar la desigualdad educativa en la ciudad de Guayaquil, a partir de los datos del Censo de Población y Vivienda 2022 del INEC. Se busca identificar patrones de asociación espacial en la escolaridad por sectores censales. Este enfoque integrador contribuye a comprender la dimensión territorial de la desigualdad educativa y a generar evidencia útil para la planificación y gestión de políticas públicas orientadas hacia una educación más equitativa e inclusiva.

2. OBJETIVOS

Analizar la distribución espacial de la escolaridad en Guayaquil y su relación con variables socioeconómicas, identificando patrones de desigualdad mediante estadística espacial.

2.1. OBJETIVOS ESPECÍFICOS

  • Realizar un análisis exploratorio de datos censales para conocer la distribución de la escolaridad e indicadores claves que describan el fenómeno estudiado.
  • Evaluar la existencia de autocorrelación espacial en los niveles de escolaridad por sectores censales, mediante medidas globales y locales.
  • Interpretar los patrones espaciales obtenidos para aportar evidencia sobre la desigualdad educativa y sus implicaciones territoriales en la ciudad.

3. HIPÓTESIS

La distribución espacial de la escolaridad en Guayaquil no es aleatoria. Se plantea que existe un agrupamiento espacial significativo, donde los niveles educativos altos tienden a concentrarse en sectores con mejores condiciones socioeconómicas, mientras que los niveles educativos bajos se agrupan en zonas periféricas con menor infraestructura educativa y menor calidad de vida.

4. DESARROLLO Y RESULTADOS

El presente estudio se estructura en diferentes estapas de desarrollo utilizando un analisis exploratorio de datos de tipo estadístico y espacial. Luego se trabajó con asociaciones espaciales de tipo global y local. Por último, se estudiaron las variables desde el enfoque de regresión geograficamente ponderada tomando en consideración los supuestos de una regresión. Este enfoque fue utilizado para conocer la distribución de escolaridad en Guayaquil-Ecuador conociendo que su espacialización no es aleatoria sino que se ve influenciada por variables socioeconómicas.

4.1. Preparación del ambiente de trabajo

Antes de iniciar se cargaron los paquetes que permitieron la manipulación y transformación de los datos; así como funcionalidades de estadística espacial requerido para abordar el objetivo general, específico y testear la hipótesis planteada.

#paquetes importantes iniciales
library(readr) #Importar dataset de tipo csv.
library(sf) # gestionar datos espaciales
library(tidyverse) # gestioner datos tabulares
library(spdep) # asociaciones espaciales
library(patchwork) # Mapeo múltiple
library(spgwr) # para GWR.
library(nortest) # Normalidad.
library(bestNormalize) # #libreria para normalizar aplicando la mejor transformación
library(car) #para seleccionar variables VIF y modelo
library(classInt) #crear intervalos

4.2. Fuentes de datos para el estudio

  • Dataset del VIII Censo de Poblacion y VII de Vivienda: CSV con variables del censo, escolaridad, empleo, tipo de vivienda, NBI (índice de pobreza), entre otras variables.
  • Dataset de Sectores Censales: Son unidades censales en fomato shapefile de tipo polígono para todo el Ecuador año 2022.
  • Dataset de Parroquias urbanas: Shapefile con geometría de tipo polígono que contiene los límites administrativos dentro del cantón Guayaquil.

4.2.1. Dataset del VIII Censo de Población y VII de Vivienda

El Instituto Nacional de Estadística y Censos presentó el dataset del censo en cinco versiones diferentes que van en función de la unidad de análisis respectiva. Para el estudio se utilizó el dataset CPV_Población_2022_Nacional_Sector_Editado: que contiene el conteo poblacional, estructura y características sociodemográficas de toda la población residente en el país en el periodo censal. Las características sociodemográficas abarcan datos relacionados con el nivel educativo de la población, su grado de escolaridad, las tasas de asistencia a la educación formal, el uso de tecnologías de la información y comunicación (TIC), y dificultades funcionales.

El dataset CPV_Población_2022_Nacional_Sector_Editado es bastante extenso. Se detallan las variables relevantes para el estudio:

  • I01: Representa el identificador de la provincia. Variable de tipo texto.

  • I02: Representa el identificador del Cantón. Variable de tipo texto.

  • I03: Representa el identificador de la Parroquia. Variable de tipo texto.

  • I04: Representa el identificador de la Zona. Variable de tipo texto.

  • I05: Representa el identidicador del Sector. Variable de tipo texto.

  • I10: Representa el identificador del número de viviendas. Variable de tipo texto.

  • P00: Representa el identificador con el número de personas. Variable de tipo texto.

  • P02: Representa el identificador del sexo al nacer. Variable de tipo categórica.

         1. Hombre
         2. Mujer
  • P17R: Representa el identificador con el nivel de instrucción más alto. Variable de tipo categórica.

         1. Ninguno
         2. Centro de Desarrollo Infantil
         3. Educación Inicial
         4. Alfabetización
         5. Educación General Básica
         6. Bachillerato
         7. Ciclo Postbachillerato (No superior)
         8. Educación Técnica o Tecnológica Superior
         9. Educación Superior (Universidades)
         10. Maestría/Especialización
         11. PHD/Doctorado
  • P29: Representa el identificador sobre la categoría de ocupción. Variable de tipo categórica.

         1. Empleada/o u obrera/o privado?
         2. Empleada/o u obrera/o del Estado, Gobierno, Municipio, Consejo Provincial, Junta Parroquial?
         3. Jornalera/o o peón?
         4. Empleada/o doméstica/o?
         5. Patrona/o?
         6. Cuenta propia?
         7. Socia/o?
         8. Trabajadora/or familiar no remunerada/o?
         9. Se ignora
  • CANTON: Representa el nombre del cantón. Variable de tipo texto.

  • PARROQ: Representa el nombre de la parroquia. Variable de tipo texto.

  • ETAEDAD: Representa el grupo de edad por etapas de vida. Variable de tipo categórica.

         1. Niñas/os de 0 a 11 años 
         2. Adolescentes de 12 a 17 años 
         3. Jóvenes de 18 a 29 años 
         4. Adultas/os de 30 a 64 años 
         5. Adultas/os mayores de 65 años o más 
  • ESCOLA: Representa los años de escolaridad. Variable de tipo numérico.

4.2.2. Dataset de Sectores Censales

La base CPV_Población_2022_Nacional_Sector_Editado al no estár georeferenciada, se necesitan archivos auxiliares para poder representar los datos a nivel territorial. Para ello, como dataset auxiliar, usamos los sectores censales que consiste en una división realizada para investigaciones estadísticas. Estas superficies son perfectamente delimitadas y geográficamente contiguas.

Se detallan las variables relevantes para el estudio:

  • sec_anm: Representa los sectores amanzanados. Variable de tipo texto.

  • nom_par: Representa el nombre de la parroquia. Variable de tipo texto.

  • nom_can: Representa el nombre del cantón. Variable de tipo texto.

  • nom_pro: Representa el nombre de la provincia. Variable de tipo texto.

4.2.3. Dataset de Parroquias Urbanas

Se detallan las variables relevantes para el estudio:

  • Nam: Representa el nombre de la parroquia urbana de Guayaquil. Variable de tipo texto.

  • has: Representa las hectáreas de cada parroquia urbana. Variable de tipo numérico.

4.3. Análisis exploratorio descriptivo

4.3.1. Selección de variables de interés

En esta sección exploramos las variables que tiene el dataset CPV_Población_2022_Nacional_Sector_Editado para comprender mejor su contenido y realizar posteriormente una selección de variables de interés. Para ello, cargados el dataset que tiene formato .csv.

#Cargaos el dataset completo
BDD_POB_CPV2022_SECT <- read_delim("Base de datos/BDD_POB_CPV2022_SECT.csv", 
                                   delim = ";", escape_double = FALSE, trim_ws = TRUE)

El dataset contiene 16’938.986 observaciones con 88 variables. Como se podría intuir, contiene los datos sociodemográficos de toda la población del Ecuador. Estas 88 variables se encuentran codificadas, de tal forma, que se requirió utilizar un diccionario desarrollado por el INEC para entender su signficado.

Revisamos su estructura usando la función glimpse()

#revisión de estructura 
glimpse(BDD_POB_CPV2022_SECT)
## Rows: 16,938,986
## Columns: 88
## $ I01       <chr> "01", "01", "01", "01", "01", "01", "01", "01", "01", "01", …
## $ I02       <chr> "01", "01", "01", "01", "01", "01", "01", "01", "01", "01", …
## $ I03       <dbl> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, …
## $ I04       <chr> "001", "001", "001", "001", "001", "001", "001", "001", "001…
## $ I05       <chr> "001", "001", "001", "001", "001", "001", "001", "001", "001…
## $ I10       <chr> "0001", "0001", "0001", "0001", "0001", "0002", "0002", "000…
## $ INH       <chr> "01", "01", "01", "01", "01", "01", "01", "01", "01", "01", …
## $ P00       <chr> "0001", "0002", "0003", "0004", "0005", "0001", "0002", "000…
## $ P01       <dbl> 1, 2, 3, 3, 3, 1, 2, 6, 1, 1, 2, 1, 2, 3, 3, 1, 2, 3, 3, 1, …
## $ P02       <dbl> 1, 2, 2, 2, 2, 1, 2, 2, 1, 1, 2, 2, 1, 2, 1, 1, 2, 1, 2, 2, …
## $ P03       <dbl> 54, 52, 27, 23, 14, 80, 77, 19, 30, 32, 40, 35, 47, 11, 3, 3…
## $ P05       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ P0601     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ P0602     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P0701     <dbl> 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, …
## $ P0702     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, …
## $ P0703     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, …
## $ P0704     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, …
## $ P0705     <dbl> 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1, 2, …
## $ P0706     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, …
## $ P08       <dbl> 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 2, 1, 1, 1, 1, 2, 2, 2, 1, 1, …
## $ P08P      <chr> "01", "01", "01", "01", "01", "03", "03", "10", "01", "01", …
## $ P08C      <chr> "0101", "0101", "0101", "0101", "0101", "0301", "0301", "100…
## $ P08Q      <chr> "010150", "010150", "010150", "010150", "010150", "030150", …
## $ P0803A    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P09       <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 4, 1, 1, 1, 1, 1, …
## $ P09P      <chr> "01", "01", "01", "01", "01", "01", "01", "10", "01", "01", …
## $ P09C      <chr> "0101", "0101", "0101", "0101", "0101", "0101", "0101", "100…
## $ P09Q      <chr> "010150", "010150", "010150", "010150", "010150", "010150", …
## $ P1001     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P1002     <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ P1003     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3, NA, NA, N…
## $ P1004     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P1005     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P1001I    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P11R      <dbl> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, …
## $ P1108     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P12       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P13       <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ P15       <dbl> 2, 2, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 1, 2, …
## $ P16       <dbl> NA, NA, NA, NA, 2, NA, NA, 2, NA, NA, NA, NA, NA, 2, NA, NA,…
## $ P17R      <dbl> 5, 9, 9, 6, 5, 10, 9, 9, 9, 6, 6, 9, 6, 5, 1, 6, 6, 6, 5, 5,…
## $ P18R      <dbl> 10, 8, 8, 3, 9, 2, 2, 0, 3, 3, 3, 5, 3, 6, NA, 3, 3, 1, 9, 6…
## $ P19       <dbl> 1, NA, NA, NA, 1, NA, NA, NA, NA, NA, NA, NA, NA, 1, NA, NA,…
## $ P20       <dbl> NA, 1, 1, NA, NA, 1, 2, 2, 2, NA, NA, 1, NA, NA, NA, NA, NA,…
## $ P2101     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, NA, 1, 1, 1, 1, 2,…
## $ P2102     <dbl> 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, NA, 1, 1, 1, 1, 2,…
## $ P2103     <dbl> 2, 1, 2, 2, 2, 1, 2, 1, 2, 1, 1, 2, 2, 1, NA, 1, 1, 1, 1, 2,…
## $ P2104     <dbl> 2, 2, 2, 2, 2, 1, 2, 1, 2, 1, 1, 2, 2, 2, NA, 1, 2, 2, 1, 2,…
## $ P22       <dbl> 1, 1, 7, 1, 7, 7, 7, 1, 1, 1, 1, 7, 1, 7, NA, 1, 7, 7, 7, 7,…
## $ P23       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P24       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ P25       <dbl> NA, NA, 1, NA, 2, 2, 2, NA, NA, NA, NA, 2, NA, 2, NA, NA, 2,…
## $ P26       <dbl> NA, NA, NA, NA, 3, 2, 2, NA, NA, NA, NA, 4, NA, 3, NA, NA, 4…
## $ P27       <chr> "5221", "2341", NA, "4311", NA, NA, NA, "9999", "5120", "711…
## $ P28       <chr> "8411", "8411", NA, "4690", NA, NA, NA, "9999", "5610", "410…
## $ P29       <dbl> 2, 2, NA, 1, NA, NA, NA, 1, 1, 1, 4, NA, 2, NA, NA, 1, NA, N…
## $ P30       <dbl> 1, 1, 7, 1, 7, 6, 6, 7, 1, 7, 1, 7, 1, NA, NA, 1, 7, 7, 7, 6…
## $ P31       <dbl> 5, 5, 6, 6, 6, 5, 5, 6, 6, 1, 1, 5, 5, NA, NA, 1, 1, 6, 6, 4…
## $ P3201     <dbl> NA, 3, 0, 0, 0, NA, 2, 0, NA, NA, 0, 1, NA, NA, NA, NA, 1, N…
## $ P3202     <dbl> NA, 0, 0, 0, 0, NA, 3, 0, NA, NA, 1, 1, NA, NA, NA, NA, 1, N…
## $ P3203     <dbl> NA, 3, 0, 0, 0, NA, 5, 0, NA, NA, 1, 2, NA, NA, NA, NA, 2, N…
## $ P3301     <dbl> NA, 3, NA, NA, NA, NA, 2, NA, NA, NA, 0, 1, NA, NA, NA, NA, …
## $ P3302     <dbl> NA, 0, NA, NA, NA, NA, 3, NA, NA, NA, 1, 1, NA, NA, NA, NA, …
## $ P3303     <dbl> NA, 3, NA, NA, NA, NA, 5, NA, NA, NA, 1, 2, NA, NA, NA, NA, …
## $ P34       <dbl> NA, 23, NA, NA, NA, NA, 24, NA, NA, NA, 18, 22, NA, NA, NA, …
## $ P3501     <dbl> NA, 5, NA, NA, NA, NA, 25, NA, NA, NA, 5, 6, NA, NA, NA, NA,…
## $ P3502     <dbl> NA, 9, NA, NA, NA, NA, 8, NA, NA, NA, 4, 3, NA, NA, NA, NA, …
## $ P3503     <dbl> NA, 2008, NA, NA, NA, NA, 1983, NA, NA, NA, 2001, 2019, NA, …
## $ AUR       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ CANTON    <chr> "0101", "0101", "0101", "0101", "0101", "0101", "0101", "010…
## $ PARROQ    <chr> "010150", "010150", "010150", "010150", "010150", "010150", …
## $ ID_VIV    <chr> "0101500010010001", "0101500010010001", "0101500010010001", …
## $ ID_HOG    <chr> "010150001001000101", "010150001001000101", "010150001001000…
## $ ID_PER    <chr> "0101500010010001010001", "0101500010010001010002", "0101500…
## $ GEDAD     <dbl> 12, 12, 7, 6, 4, 18, 17, 5, 8, 8, 10, 9, 11, 4, 2, 9, 9, 5, …
## $ GRANEDAD  <dbl> 2, 2, 2, 2, 1, 3, 3, 2, 2, 2, 2, 2, 2, 1, 1, 2, 2, 2, 1, 3, …
## $ ETAEDAD   <dbl> 4, 4, 3, 3, 2, 5, 5, 3, 4, 4, 4, 4, 4, 1, 1, 4, 4, 2, 2, 5, …
## $ DFUNC     <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 1, …
## $ P10R      <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 7, 2, 2, 2, 2, 2, 2, 2, …
## $ ESCOLA    <dbl> 10, 21, 21, 13, 9, 20, 15, 13, 16, 13, 13, 18, 13, 6, 0, 13,…
## $ ANALF     <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, NA, 2, 2, 2, 2, 2,…
## $ ANALF_DIG <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, NA, 2, 2, 2, 2, 1,…
## $ CONDACT   <dbl> 2, 2, 2, 2, 3, 3, 3, 2, 2, 2, 2, 3, 2, 3, 1, 2, 3, 3, 3, 3, …
## $ CONDACT1  <dbl> 2, 2, 3, 2, 4, 4, 4, 2, 2, 2, 2, 4, 2, 4, 1, 2, 4, 4, 4, 4, …
## $ GRUPO1    <dbl> 5, 2, NA, 4, NA, NA, NA, 99, 5, 7, 5, NA, 8, NA, NA, 5, NA, …
## $ RAMA1     <dbl> 15, 15, NA, 7, NA, NA, NA, 99, 9, 6, 20, NA, 8, NA, NA, 14, …
## $ IMP_VOPA  <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …

La variables I1 - I5 contienen los códigos del INEC correspondiente a la Provincia, Cantón, Parroquia, Zona, Sector los cuales debemos unificar para poder crear el ID_GEOGRAFICO que nos servirá para vincular con sectores censales. Adicionalmente, segregamos el dataset para que contenga solamente la información correspondiente a la parroquia Guayaquil ya sea urbana y rural utilizando el código 090150 filtrato a partir de la variable PARROQ.

#segregamos la base 
BDD_POB_GYE <- BDD_POB_CPV2022_SECT |> 
  select(I01, I02, I03, I04, I05, I10, P00, P02, P17R, P29, CANTON, PARROQ, ETAEDAD, ESCOLA) |> 
  filter(PARROQ == "090150") |> 
  mutate(ID_GEOGRAFICO = paste0(I01, I02, I03, I04, I05))

Del total de observaciones, nos quedamos con 2’665.392 observaciones y 14 variables que corresponden a la parte urbana y rural de la parroquia Guayaquil. El registro censal se dio a nivel de viviendas y personas censadas. Por lo cual, tenemos varios registros por cada vivienda dentro de cada sector censal. A continuación, procedemos con un análisis exploratorio descriptivo de variables de interés para conocer su distribución previo al análisis de estadística espacial.

4.3.2. Análisis estadístico exploratorio descriptivo

Como análisis exploratorio inicial, aplicamos la función summary para conocer la distribuición de las variables de estudio.

summary(BDD_POB_GYE)
##      I01                I02                 I03         I04           
##  Length:2665392     Length:2665392     Min.   :50   Length:2665392    
##  Class :character   Class :character   1st Qu.:50   Class :character  
##  Mode  :character   Mode  :character   Median :50   Mode  :character  
##                                        Mean   :50                     
##                                        3rd Qu.:50                     
##                                        Max.   :50                     
##                                                                       
##      I05                I10                P00                 P02       
##  Length:2665392     Length:2665392     Length:2665392     Min.   :1.000  
##  Class :character   Class :character   Class :character   1st Qu.:1.000  
##  Mode  :character   Mode  :character   Mode  :character   Median :2.000  
##                                                           Mean   :1.511  
##                                                           3rd Qu.:2.000  
##                                                           Max.   :2.000  
##                                                                          
##       P17R             P29             CANTON             PARROQ         
##  Min.   : 1.000   Min.   :1.00      Length:2665392     Length:2665392    
##  1st Qu.: 5.000   1st Qu.:1.00      Class :character   Class :character  
##  Median : 6.000   Median :2.00      Mode  :character   Mode  :character  
##  Mean   : 5.909   Mean   :3.65                                           
##  3rd Qu.: 6.000   3rd Qu.:6.00                                           
##  Max.   :11.000   Max.   :9.00                                           
##                   NA's   :1613419                                        
##     ETAEDAD          ESCOLA      ID_GEOGRAFICO     
##  Min.   :1.000   Min.   : 0.00   Length:2665392    
##  1st Qu.:2.000   1st Qu.: 7.00   Class :character  
##  Median :3.000   Median :12.00   Mode  :character  
##  Mean   :3.091   Mean   :10.08                     
##  3rd Qu.:4.000   3rd Qu.:13.00                     
##  Max.   :5.000   Max.   :25.00                     
## 

El nivel educativo promedio de la población censada en Guayaquil se ubica en el grupo correspondiente a Bachillerato, con una distribución aproximadamente normal (media y mediana = 6). Los valores mínimos registrados corresponden a 1, que representa “ningún nivel educativo”, mientras que el valor máximo es 11, que corresponde a “PhD o Doctorado”.

Respecto a la ocupación, la media de la población se ubica en 3.65, lo que corresponde a los grupos “jornalero o peón” y “empleado/a doméstico/a”. Los valores mínimos (1) representan “empleado/a u obrero/a privado”, y los máximos indican categorías “no especificadas”. En esta variable se registraron 1.613.419 valores faltantes (NA).

Finalmente, al analizar la edad agrupada por etapas de vida, la población promedio pertenece al grupo 3, correspondiente a jóvenes de 18 a 29 años, con un nivel de escolaridad promedio de 10 años.

Para conocer la distribución categóricas bivariadas que pueden explicar el comportamiento de la escolaridad en Guayaquil, iniciamos por la distribución de los encuestados por sexo como se muestra en la figura 1.

#Distribución de la población por sexo
ggplot(BDD_POB_GYE, aes(x = factor(P02, levels = c(1, 2), labels = c("Hombre", "Mujer")))) +
  geom_bar(fill = "#3182bd", color = "black") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.3) +  # agrega los valores encima de cada barra
  labs(title = "Fig.1.Distribución de las personas encuestadas por sexo", 
       subtitle = " Población censada en Guayaquil",
       x = "Sexo", 
       y = "Frecuencia de personas", 
       caption = "Fuente: Censo de Población y Vivienda 2022 - INEC") +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold"))

La figura 1 muestra la frecuencia de personas (el número total de personas) encuestadas o censadas en Guayaquil, diferenciadas por genero: hombre y mujer.

El resultado del Censo de Población y Vivienda 2022 del INEC, muestra que la población de mujeres censadas en Guayaquil es mayor que la de hombres con un diferencia porcentual de 2,22%.

Para continuar con el análisis de la población objeto de estudio en Guayaquil - Ecuador, mostramos en la figura 2 una gráfica de barras bivariada crucial que muestra la distribución del grupo de edad por etapas de vida de la Población censada en Guayaquil, desagregada por sexo.

#Distribución de grupo etario y sexo
ggplot(BDD_POB_GYE, aes(x = factor(ETAEDAD, 
                                   levels = c(1,2,3,4,5),
                                   labels = c("Niñas/os de 0 a 11 años", "Adolescentes de 12 a 17 años ", "Jóvenes de 18 a 29 años",
                                              "Adultas/os de 30 a 64 años", "Adultas/os mayores de 65 años o más")))) +
  geom_bar(fill = "#3182bd", color = "black") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.3) +
  labs(title = "Fig.2.Distribución de grupo de edad por estapas de vida",
       subtitle = "Población censada en Guayaquil",
       x = "Etapas de vida",
       y = "Frecuencia de personas",
       caption = "Fuente: Censo de Población y Vivienda 2022 - INEC") +
  facet_wrap(~ factor(P02, levels = c(1,2), labels = c("Hombres", "Mujeres")), 
             scales = "free_y") +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold"))

La estructura poblacional de Guayaquil está dominada por adultos en edad productiva (30-64 años) y muestra un claro patrón donde los hombres son más frecuentes en la niñez, mientras que las mujeres predominan en la edad adulta y la vejez.

El Grupo Demográfico Dominante

Distribución en los Extremos

Diferencia por Sexo

La figura 3 muestra un gáfico de barras con la distribución del nivel de escolaridad por sexo de la población censada en Guayaquil.

#Distribución del nivel de escolaridad y sexo
ggplot(BDD_POB_GYE, aes(x = factor(P17R, 
                                   levels = c(1,2,3,4,5,6,7,8,9,10,11),
                                   labels = c("Ninguno", "Centro de Desarrollo Infantil", "Educación Inicial",
                                              "Alfabetización", "Educación Básica", "Bachillerato", "Ciclo Postbachillerato",
                                              "Educación Técnica", "Educación Superior", "Maestría", "PHD")))) +
  geom_bar(fill = "#3182bd", color = "black") +
  labs(title = "Fig.3.Distribución de nivel de escolaridad por sexo",
       subtitle = "Población censada en Guayaquil",
       x = "Nivel de escolaridad",
       y = "Frecuencia de personas",
       caption = "Fuente: Censo de Población y Vivienda 2022 - INEC") +
  facet_wrap(~ factor(P02, levels = c(1,2), labels = c("Hombres", "Mujeres")), 
             scales = "free_y") +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold")
  )

Hombres

El nivel de escolaridad más frecuente entre los hombres es Educación Básica, con aproximadamente 466.349 personas, seguido muy de cerca por Bachillerato, con cerca de 431.245 personas. El tercer nivel más común es Educación Superior (Universitaria), con aproximadamente 228.120 personas. Los niveles menos frecuentes son PHD (1.525), Alfabetización (604) y Posgrado/Maestría (18.119).

Mujeres

Al igual que en los hombres, el nivel de escolaridad más frecuente entre las mujeres es Educación Básica, con aproximadamente 474.582 personas. El segundo nivel más común es Bachillerato, con cerca de 422.595 personas. El tercer nivel más común es Educación superior (Universitario), con aproximadamente 278.151 personas, lo que sugiere una mayor proporción de mujeres que han alcanzado este nivel académico en comparación con los hombres. Los niveles menos frecuentes son PHD (1.381), Alfabetización (925) y Educación Técnica (7.219).

La figura 4 muestra un gáfico de barras con la distribución de los años de escolaridad agrupados por sexo de la población.

# Distribución de la escolaridad por sexo
ggplot(BDD_POB_GYE, aes(x = factor(ESCOLA))) +
  geom_bar(fill = "#3182bd", color = "white") +
  labs(title = "Fig.4.Distribución de los años de escolaridad por sexo", 
       subtitle = " Población censada en Guayaquil",
       x = "Años de escolaridad", 
       y = "Frecuencia de personas", 
       caption = "Fuente: Censo de Población y Vivienda 2022 - INEC") +
  facet_wrap(~ P02, scales = "free_y", labeller = as_labeller(c("1" = "Hombres", "2" = "Mujeres"))) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold"))

En ambos grupos se observa una tendencia central marcada hacia los 13 años de estudio, correspondiente al nivel de Bachillerato, lo que refleja que la mayoría de la población ha alcanzado la educación secundaria completa.

Asimismo, se identifican picos secundarios en los valores bajos (0 y 7 años), asociados a personas sin instrucción o con niveles educativos incompletos, lo que evidencia la persistencia de segmentos poblacionales con acceso limitado a la educación formal. A medida que aumenta el nivel educativo, la frecuencia disminuye de manera progresiva, siendo muy reducida la proporción de personas con formación universitaria o de posgrado (más de 16 años de escolaridad).

Al comparar por sexo, las mujeres presentan una distribución similar a la de los hombres, aunque con ligeramente mayor concentración en los niveles medios (entre 14 y 18 años), lo que sugiere una participación educativa relativamente equilibrada entre ambos grupos. Sin embargo, la magnitud de las diferencias entre los extremos (sin instrucción y estudios superiores) podría indicar brechas persistentes en el acceso y continuidad educativa, posiblemente vinculadas a factores socioeconómicos o territoriales.

En conjunto, esta distribución evidencia que, aunque Guayaquil mantiene una base educativa sólida en los niveles medios, existen desigualdades estructurales en los extremos del sistema educativo, que serán posteriormente analizadas mediante técnicas de asociación espacial para identificar patrones locales de desigualdad educativa.

Otra perspectiva de análisis fue conocer como de distribuye los años de escolaridad en función de las etapas de vida como se muestra en la figura 5.

# Distribución de la escolaridad por estapas de vida
ggplot(BDD_POB_GYE, aes(x = factor(ESCOLA))) +
  geom_bar(fill = "#3182bd", color = "white") +
  labs(title = "Fig.5.Distribución de los años de escolaridad por etapas de vida", 
       subtitle = " Población censada en Guayaquil",
       x = "Años de escolaridad", 
       y = "Frecuencia de personas", 
       caption = "Fuente: Censo de Población y Vivienda 2022 - INEC") +
  facet_wrap(~ ETAEDAD, scales = "free_y", labeller = as_labeller(c("1" = "Niñas/os de 0 a 11 años ", "2" = "Adolescentes de 12 a 17 años ", 
                                                                    "3" = "Jóvenes de 18 a 29 años ", "4" = "Adultas/os de 30 a 64 años",
                                                                    "5" = "Adultas/os mayores de 65 años o más "))) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    axis.title.x = element_text(face = "bold"))

En el grupo infantil (0–11 años) se observa una alta concentración en los valores más bajos de escolaridad (0 años), correspondiente a niños y niñas que aún no han iniciado la educación básica. Otro subconjunto se distribuye de manera más uniforme entre 1 y 7 años, reflejando la progresión en los primeros niveles de formación.

Entre los adolescentes (12–17 años), la distribución se concentra principalmente entre 7 y 11 años de escolaridad, lo que indica que la mayoría se encuentra cursando niveles intermedios o finalizando la educación básica. No obstante, un pequeño grupo mantiene niveles bajos de escolaridad, evidenciando que algunos aún no han iniciado formalmente su educación.

En la juventud (18–29 años) se aprecia un pico marcado en los 13 años de escolaridad, equivalente a la culminación del Bachillerato, señalando que este grupo ha completado la educación secundaria. La baja frecuencia de valores superiores a 16 años indica que solo una proporción limitada accede a estudios universitarios o de posgrado.

En los adultos (30–64 años), la distribución es más dispersa, aunque se mantiene un patrón central alrededor de los 13 años de escolaridad. Se observa, además, un grupo con niveles bajos (0–7 años), reflejando brechas educativas respecto a generaciones más jóvenes, así como un pequeño grupo con estudios universitarios completos, centrados en los 18 años de escolaridad.

Finalmente, en el grupo de adultos mayores (65 años o más), se destaca una mayor proporción de personas con escolaridad de 7 años (educación básica), otro grupo con mayor concentración en los 13 años (Bachillerato) y un pequeño grupo con estudios universitarios de grado. Esto refleja las limitaciones históricas en el acceso a la educación formal en décadas anteriores.

En conjunto, los resultados evidencian un progreso gradual en los niveles educativos de las generaciones más jóvenes, aunque persisten desigualdades significativas entre grupos etarios, asociadas al tiempo, al contexto socioeconómico y, posiblemente, al territorio, aspectos que serán analizados en detalle en la sección de análisis espacial.

Por último, fue necesario conocer como de distribuye la ocupación de los encuestados categorizado por sexo como se detalla en la figura 6.

# Distribución de la escolaridad
ggplot(BDD_POB_GYE, aes(x = factor(P29, 
                                   levels = c(1,2,3,4,5,6,7,8,9),
                                   labels = c("Empleada/o privado", "Empleada/o del Estado", "Jornalera/o",
                                              "Empleada/o doméstica/o", "Patrona/o", "Cuenta propia", "Socia/o",
                                              "Trabajadora/o no remunerado", "Se ignora")))) +
  geom_bar(fill = "#3182bd", color = "black") +
  labs(title = "Fig.6.Distribución de la ocupación por sexo",
       subtitle = "Población censada en Guayaquil",
       x = "Categoría de ocupación",
       y = "Frecuencia de personas",
       caption = "Fuente: Censo de Población y Vivienda 2022 - INEC") +
  facet_wrap(~ factor(P02, levels = c(1,2), labels = c("Hombres", "Mujeres")), 
             scales = "free_y") +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold")
  )

Hombres

Mujeres

La estructura de ocupación en Guayaquil muestra una mayor inserción de los hombres en el mercado laboral remunerado (tanto formal como por cuenta propia) y una mayor concentración de las mujeres en roles no remunerados (reflejado en la categoría NA).

4.4. Análisis exploratorio espacial descriptivo.

Cargamos el dataset espacial con los sectores censales del Ecuador

#cargamos sectores censales
sectores <-  st_read("SHP/sectores_censales_GYE.shp")
## Reading layer `sectores_censales_GYE' from data source 
##   `C:\R_directory\4. Geoestadística\Trabajo final - INEC\SHP\sectores_censales_GYE.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6357 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -80.09211 ymin: -2.28721 xmax: -79.85996 ymax: -1.962213
## Geodetic CRS:  WGS 84

Como primera aproximación, vemos que el dataset cargado es de tipo multipolígono y contiene 6.357 geometrias con 5 columnas despcriptivas que corresponden a los sectores censales de Guayquil.

Vemos la estructura completa del dataset cargado.

#estructura del dataset
glimpse (sectores)
## Rows: 6,357
## Columns: 6
## $ fid      <dbl> 12, 21, 23, 25, 28, 31, 33, 36, 37, 40, 51, 96, 97, 98, 99, 1…
## $ sec_anm  <chr> "090150318010", "090150318006", "090150318007", "090150318004…
## $ nom_par  <chr> "GUAYAQUIL", "GUAYAQUIL", "GUAYAQUIL", "GUAYAQUIL", "GUAYAQUI…
## $ nom_can  <chr> "GUAYAQUIL", "GUAYAQUIL", "GUAYAQUIL", "GUAYAQUIL", "GUAYAQUI…
## $ nom_pro  <chr> "GUAYAS", "GUAYAS", "GUAYAS", "GUAYAS", "GUAYAS", "GUAYAS", "…
## $ geometry <POLYGON [°]> POLYGON ((-79.87609 -2.1823..., POLYGON ((-79.87529 -…

Revisada la estructura, tenemos la variable sec_anm que contiene el ID_GEOGRAFICO que necesitamos para unirlo con el dataset del censo pero primero, debemos renombrar la variable sec_anm para luego utilizarlo en la vinculación con la base alfanumérica. Este renombre debe ser igual al del dataset del censo. Esto es ID_GEOGRAFICO.

#filtrado los sectores censales
sectores <-  sectores |> 
  rename(ID_GEOGRAFICO = sec_anm) #cambiamos el nombre de la variable para el join

Tenemos en total 6.415 sectores censales con los cuales vamos a trabajar. Como siguiente punto, es determinar indicadores por sectores censales para continuar con la unión espacial. Estos indicadores se estándarizaron a nivel de nomenclatura de las variables para evitar errores.

Los indicadores Demográficos y Socioeconómicos claves fueron:

  • Indicadores Demográficos y de Género: índice de Masculinidad, índice de Feminidad, Índice de envejecimiento, Índice de Juventud, Índice de adultez, Índice de adultos mayores.
  • Indicadores de Escolaridad (Logro Educativo): Tasa de Analfabetismo, Tasa de Logro Básico, Tasa de Logro Bachillerato, Tasa de Logro Superior.
  • Indicadores del Mercado Laboral (ocupación): Tasa de empleo formal (Privado), Tasa de empleo propio (Informalidad), Tasa de Trabajo No Remunerado, Tasa de ocupación por sexo.
#indicadores censales
indicadores_censales <- BDD_POB_GYE |> 
  group_by(ID_GEOGRAFICO) |> 
  summarise(I_MASCULINIDAD = (sum(P02 == 1) / sum(P02 == 2)) * 100,
            I_FEMINIDAD = (sum(P02 == 2) / sum(P02 == 1)) * 100,
            I_ENVEJECIMIENTO  = (sum(ETAEDAD == 05) / sum(ETAEDAD == 01)) * 100,
            I_JUVENTUD        = (sum(ETAEDAD == 03) / n()) * 100,
            I_ADULTEZ         = (sum(ETAEDAD == 04) / n()) * 100,
            I_ADULTEZ_MAYOR   = (sum(ETAEDAD == 05) / n()) * 100,
            T_ANALFABETISMO = (sum(P17R == 01)/ sum (ETAEDAD %in% c(03, 04, 05))) *100,
            T_LOGRO_BACH = (sum(P17R == 06)/ sum (ETAEDAD %in% c(03, 04, 05))) *100,
            T_LOGRO_SUP = (sum(P17R == 09)/ sum (ETAEDAD %in% c(04, 05))) *100,
            POB_OCUP = sum(P29 %in% c(01, 02, 03, 04, 05, 06, 07, 08), na.rm = TRUE),
            T_EMPLEO_FORMAL = (sum(P29 %in% c(01, 02), na.rm = TRUE) /  POB_OCUP *100),
            C_PROPIA = (sum(P29 == 06, na.rm = TRUE) / POB_OCUP) * 100,
            T_TRABAJO_NO_REM = (sum(P29 == 08, na.rm = TRUE) / POB_OCUP) * 100,
            T_OCUPACION_SEXO = (sum(P29 != 9 & P02 == 2, na.rm = TRUE) / sum(P29 != 9 & P02 == 1, na.rm = TRUE))
            )

Calculado los indicadores censales procedemos con la vinculación de los sectores censales.

#Join con sectores censales
sectores_censales_BDD <- sectores |> 
  left_join(indicadores_censales, by = 'ID_GEOGRAFICO')

Primero graficamos los Indicadores Demográficos para entender su distribución en el territorio.

# Tema base sin ejes ni coordenadas
tema_mapa <- theme_void() +
  theme(
    plot.title = element_text(
      face = "bold",         # negrita
      size = 12,             # tamaño moderado
      hjust = 0.5            # centrado
    ),
    legend.position = "bottom",
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 8)
  )

# Mapas individuales
p1DEMO <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = I_ENVEJECIMIENTO), color = NA) +
  scale_fill_viridis_c(option = "inferno", name = "Porcentaje") +
  labs(title = "I_ENVEJECIMIENTO") +
  tema_mapa

p2DEMO <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = I_JUVENTUD), color = NA) +
  scale_fill_viridis_c(option = "inferno", name = "Porcentaje") +
  labs(title = "I_JUVENTUD") +
  tema_mapa

p3DEMO <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = I_ADULTEZ), color = NA) +
  scale_fill_viridis_c(option = "inferno", name = "Porcentaje") +
  labs(title = "I_ADULTEZ") +
  tema_mapa

p4DEMO <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = I_ADULTEZ_MAYOR), color = NA) +
  scale_fill_viridis_c(option = "inferno", name = "Porcentaje") +
  labs(title = "I_ADULTEZ_MAYOR") +
  tema_mapa


# Combinación y formato del título general
(p1DEMO | p2DEMO | p3DEMO | p4DEMO) +
  plot_annotation(
    title = "Indicadores Demográficos para Guayaquil",
    caption = "Fuente: INEC - CPV 2022",
    theme = theme(
      plot.title = element_text(
        face = "bold",
        size = 18,      # más grande
        hjust = 0.5     # centrado
      ),
      plot.caption = element_text(size = 9, hjust = 0.5)
    )
  )

Los mapas confirman una clara segregación demográfica en Guayaquil: la población joven y la población envejecida se concentran en polos opuestos de la ciudad (periferia vs. centro/áreas históricas), mientras que la población adulta está más distribuida.

Luego, graficamos los Indicadores de Escolaridad para entender su distribución en el territorio.

# Tema base sin ejes ni coordenadas
tema_mapa <- theme_void() +
  theme(
    plot.title = element_text(
      face = "bold",         # negrita
      size = 12,             # tamaño moderado
      hjust = 0.5            # centrado
    ),
    legend.position = "bottom",
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 8)
  )

# Mapas individuales
p2ESC <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = T_ANALFABETISMO), color = NA) +
  scale_fill_viridis_c(option = "cividis", name = "Porcentaje") +
  labs(title = "T_Analfabetismo") +
  tema_mapa

p4ESC <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = T_LOGRO_BACH), color = NA) +
  scale_fill_viridis_c(option = "cividis", name = "Porcentaje") +
  labs(title = "T_Bachillerato") +
  tema_mapa

p5ESC <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = T_LOGRO_SUP), color = NA) +
  scale_fill_viridis_c(option = "cividis", name = "Porcentaje") +
  labs(title = "T_Superior") +
  tema_mapa

# Combinación y formato del título general
(p2ESC | p4ESC | p5ESC) +
  plot_annotation(
    title = "Indicadores de Escolaridad en Guayaquil",
    caption = "Fuente: INEC - CPV 2022",
    theme = theme(
      plot.title = element_text(
        face = "bold",
        size = 16,      # más grande
        hjust = 0.5     # centrado
      ),
      plot.caption = element_text(size = 9, hjust = 0.5)
    )
  )

Los indicadores de escolaridad evidencian una profunda desigualdad socio-espacial. Los sectores con bajos logros educativos y alto analfabetismo están geográficamente aislados de los sectores con alto capital humano y acceso a la educación superior. En conjunto, los mapas reflejan una desigualdad espacial en la educación, donde los sectores centrales y del norte de Guayaquil presentan mayores niveles de escolaridad y educación superior, mientras que los del sur y suroeste muestran mayor analfabetismo y predominio de educación básica.

Por último, graficamos los Indicadores del Mercado Laboral para entender su distribución en el territorio.

# Tema base sin ejes ni coordenadas
tema_mapa <- theme_void() +
  theme(
    plot.title = element_text(
      face = "bold",         # negrita
      size = 12,             # tamaño moderado
      hjust = 0.5            # centrado
    ),
    legend.position = "bottom",
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 8)
  )

# Mapas individuales
p1MLAB <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = POB_OCUP), color = NA) +
  scale_fill_viridis_c(option = "viridis", name = "Recuento") +
  labs(title = "POB_OCUP") +
  tema_mapa

p2MLAB <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = T_EMPLEO_FORMAL ), color = NA) +
  scale_fill_viridis_c(option = "viridis", name = "Porcentaje") +
  labs(title = "EMPLEO_FORMAL") +
  tema_mapa

p3MLAB <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = C_PROPIA ), color = NA) +
  scale_fill_viridis_c(option = "viridis", name = "Porcentaje") +
  labs(title = "C_PROPIA") +
  tema_mapa

p4MLAB <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = T_TRABAJO_NO_REM), color = NA) +
  scale_fill_viridis_c(option = "viridis", name = "Porcentaje") +
  labs(title = "TRABAJO_NO_REM") +
  tema_mapa

p5MLAB <- ggplot(sectores_censales_BDD) +
  geom_sf(aes(fill = T_OCUPACION_SEXO), color = NA) +
  scale_fill_viridis_c(option = "viridis", name = "Razón") +
  labs(title = "OCUPA_SEXO") +
  tema_mapa

# Combinación y formato del título general
(p1MLAB | p2MLAB | p3MLAB | p4MLAB | p5MLAB) +
  plot_annotation(
    title = "Indicadores del Mercado Laboral",
    caption = "Fuente: INEC - CPV 2022",
    theme = theme(
      plot.title = element_text(
        face = "bold",
        size = 16,      # más grande
        hjust = 0.5     # centrado
      ),
      plot.caption = element_text(size = 9, hjust = 0.5)
    )
  )

Los indicadores evidencian una heterogeneidad espacial del mercado laboral: las zonas centrales muestran mejores condiciones de formalidad y ocupación, mientras que las periferias reflejan vulnerabilidad laboral y predominio de empleos informales o no remunerados. Estos patrones son consistentes con procesos de segregación socioeconómica urbana, donde la localización condiciona las oportunidades laborales.

4.5. Análisis de asociación espacial global con I_MORAN

Luego de haber realizado un análisis exploratorio de los datos del censo de población y vivienda, necesitamos conocer como se asocian estas variables analizadas consigo mismas en diferentes ubicaciones del espacio geográfico. La asociación espacial no sólo da una medida y significancia de la autocorrelación espacial, suno que también informan sobre cómo losvalores de una variable se disponen en el espacio (patrón espacial) a escala global y local.

El proceso de I_Moran necesita tener dos variables que representen el análisis. La primera es el código del sector censal y la variable que se quiere autocorrelacionar que sería T_LOGRO_BACH, T_LOGRO_SUP. Adicionalmente, se requieren definir los vecindarios ya sea por:

  • Contiguidad física o adyacencia topológica, dada por compartir bordes.
  • Distancias definida a partir de los centroides de los polígonos.

Esta definición de vecindarios pueden ser de primer nivel u orden o de segundo nivel u orden tomando en consideración la forma de búsqueda de estos vecindarios que puede ser de tipo queen (lados y vértices) y de tipo rook (sólo lados). Para el estudio utilizamos la contiguidad física.

4.5.1. Armar lista con vecinos

Explicado lo anterior, generamos los vecinos utilizando poly2nb y la búsqueda de tipo queen.

#creamos los vecinos para I Moran
nbCenso <- poly2nb(sectores_censales_BDD, row.names = sectores_censales_BDD$ID_GEOGRAFICO, queen = T)

Como primera aproximación de la lista de vecinos, se grafico un mapa de grafos sencillo para identificar su relación.

#mapa de grafos por VECINDAD
plot(st_geometry(sectores_censales_BDD)) #obtengo las geometrías.
plot(nbCenso, 
     coords =  st_coordinates(st_centroid(sectores_censales_BDD)), #obtengo los centroides de los sectores censales.
     add = TRUE, #grafico las relaciones con líneas
     col="blue") #coloreo de azúl.

4.5.2. Asignar pesos a los vecinos

Los pesos espaciales controlan cómo cada unidad geográfica influye sobre sus vecinos cuando calculas autocorrelación espacial. La asignación de pesos a los vecinos se puede realizar de múltiples formas. Los más utilizados son pesos espaciales binarios y pesos espaciales estandarizados.

Para sectores censales y estudios socioeconómicos, la bibliografía recomienda utilizar pesos estandarizados porque los valores de Moran suelen ser más estables y comparables entre datasets.

#creamos los pesos usando el método estandarizado
lwCenso <-  nb2listw(nbCenso, style = 'W')

#revisamos el resultado
lwCenso
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 6357 
## Number of nonzero links: 39786 
## Percentage nonzero weights: 0.0984523 
## Average number of links: 6.258613 
## 
## Weights style: W 
## Weights constants summary:
##      n       nn   S0       S1       S2
## W 6357 40411449 6357 2124.984 26325.41

El resultado arrojado me indica:

Como conclusión del resultado obtenido nos dice que la estructura espacial está bien definida. El promedio de vecinos indica un buen nivel de conectividad. Asi como el porcentaje de pesosno cero muy bajos es muy típico en este tipo de análisis. Obtener miles de enlaces es adecuado para análisis estadístico robusto lo que denota número de regiones alto y son favorables para modelos de Moran globales y locales.

4.5.3. I de Moran (Global)

Una vez obtenidas las relaciones entre unidades espaciales (sectores cenales) procedemos a calcular el I de Moran tomando en consideración las variables que queremos autocorrelacionar como lo son T_LOGRO_BACH, T_LOGRO_SUP.

Como primer paso, reemplazamos todos los valores NA por 0 para evitar que muestre un error en la ejecución y luego ejecutamos el índice de Moran.

T_LOGRO_BACH

# reemplazamo NA por  de T_LOGRO_BACH
sectores_censales_BDD$T_LOGRO_BACH[!is.finite(sectores_censales_BDD$T_LOGRO_BACH)] = 0

# Retenemos sólo el Índice de Moran
I_MORAN_T_LOGRO_BACH <- moran(sectores_censales_BDD$T_LOGRO_BACH, # Variable a autocorrelacionar
                 listw = lwCenso,                                 # Pesos de los vecinos
                 n = length(nbCenso),                             # Cantidad de zonas
                 S0 = Szero(lwCenso))[[1]]                        # Suma de los pesos

I_MORAN_T_LOGRO_BACH
## [1] 0.6961581

El valor obtenido es de 0.6961. Esto explica de que existe autocorrelación positiva Fuerte de la variable T_LOGRO_BACH.

T_LOGRO_SUP

# reemplazamo NA por  de T_LOGRO_SUP
sectores_censales_BDD$T_LOGRO_SUP[!is.finite(sectores_censales_BDD$T_LOGRO_SUP)] = 0

# Retenemos sólo el Índice de Moran
I_MORAN_T_LOGRO_SUP <- moran(sectores_censales_BDD$T_LOGRO_SUP, # Variable a autocorrelacionar
                 listw = lwCenso,                                 # Pesos de los vecinos
                 n = length(nbCenso),                             # Cantidad de zonas
                 S0 = Szero(lwCenso))[[1]]                        # Suma de los pesos

I_MORAN_T_LOGRO_SUP
## [1] 0.7800632

El valor obtenido es de 0.7800. Esto explica de que existe autocorrelación positiva Fuerte autocorrelación espacial de la variable T_LOGRO_SUP.

Obtenido los índices por cada variable se agruparon en un dataframe para mostrar graficamente su resultado y una línea de tendencia de los mismos.

#Agrupamos en una tabla
moran_tabla <- tibble(
  variable = c("T_LOGRO_BACH", "T_LOGRO_SUP"),
  moran_I  = c(I_MORAN_T_LOGRO_BACH, I_MORAN_T_LOGRO_SUP))

#mostramos el resultado
moran_tabla
## # A tibble: 2 × 2
##   variable     moran_I
##   <chr>          <dbl>
## 1 T_LOGRO_BACH   0.696
## 2 T_LOGRO_SUP    0.780

Graficamos el resultado obtenido.

#ggplot para graficar.
ggplot(moran_tabla, aes(x = variable, y = moran_I, group = 1)) +
  geom_point(size = 3) +
  geom_line() +
  geom_text(aes(label = round(moran_I, 3)), vjust = -0.6, size = 3.5) + 
  labs(
    title = "Índice de Moran",
    x = "Variables",
    y = "Moran's I"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold")
  )

Estos valores nos indican que existen patrones agrupados (clústers) para cada variable analizada. Continuamos con la prueba de significancia para medir el grado de autocorrelación obtenido y si este es significativo desde el punto de vista estadístico.

4.5.4. Testear significatividad del I de Moran Global

Lo dividimos por variable estudiada.

T_LOGRO_BACH

#test de significancia de Moran para Bachillerato
moran.test(sectores_censales_BDD$T_LOGRO_BACH, lwCenso)
## 
##  Moran I test under randomisation
## 
## data:  sectores_censales_BDD$T_LOGRO_BACH  
## weights: lwCenso    
## 
## Moran I statistic standard deviate = 96.076, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     0.69615808391    -0.00015733166     0.00005252719

La variable T_LOGRO_BACH presenta una muy fuerte autocorrelación espacial positiva (I = 0.696), estadísticamente significativa (p_value < 0.001). Esto indica que los niveles de logro educativo a nivel de bachillerato tienden a agruparse espacialmente: zonas con valores altos están rodeadas por zonas también altas, y lo mismo ocurre para valores bajos. El patrón espacial no es aleatorio sino altamente estructurado geográficamente.

T_LOGRO_SUP

#test de significancia de Moran para Educación Superior. 
moran.test(sectores_censales_BDD$T_LOGRO_SUP, lwCenso)
## 
##  Moran I test under randomisation
## 
## data:  sectores_censales_BDD$T_LOGRO_SUP  
## weights: lwCenso    
## 
## Moran I statistic standard deviate = 107.64, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     0.78006321235    -0.00015733166     0.00005253504

El análisis del índice de Moran global para la variable T_LOGRO_SUP revela una autocorrelación espacial alta y positivamente significativa (I = 0.7801; p_value < 0.00000000000000022), lo que demuestra que la distribución del logro educativo superior en Guayaquil no es aleatoria. Por el contrario, existe un patrón espacial claramente definido donde los sectores con mayores proporciones de población con estudios superiores tienden a concentrarse entre sí, formando clústeres bien diferenciados dentro del territorio.

4.5.5. Correlograma de I de Moran Global

T_LOGRO_BACH

#Grafico de correlograma para Bachillerato
I_CORR_T_LOGRO_BACH <- sp.correlogram(neighbours = nbCenso,
                         var = sectores_censales_BDD$T_LOGRO_BACH,
                         order = 5, # niveles de correlación
                         method = "I", #método moran
                         style = "W",  #pesos estandarizados
                         zero.policy = TRUE)
#mostramos resultados.
I_CORR_T_LOGRO_BACH
## Spatial correlogram for sectores_censales_BDD$T_LOGRO_BACH 
## method: Moran's I
##               estimate   expectation      variance standard deviate
## 1 (6357)  0.6961580839 -0.0001573317  0.0000525272           96.076
## 2 (6357)  0.5831416540 -0.0001573317  0.0000218030          124.920
## 3 (6357)  0.4889913699 -0.0001573317  0.0000123165          139.379
## 4 (6357)  0.4164964736 -0.0001573317  0.0000083353          144.316
## 5 (6357)  0.3623520067 -0.0001573317  0.0000063564          143.785
##                Pr(I) two sided    
## 1 (6357) < 0.00000000000000022 ***
## 2 (6357) < 0.00000000000000022 ***
## 3 (6357) < 0.00000000000000022 ***
## 4 (6357) < 0.00000000000000022 ***
## 5 (6357) < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El correlograma espacial de Moran para T_LOGRO_BACH muestra una autocorrelación positiva fuerte y significativa en todos los órdenes de vecindad, lo que evidencia un patrón espacial bien definido. Los sectores con altos niveles de logro educativo de bachillerato se agrupan entre sí en los primeros niveles, y lo mismo ocurre con las zonas de bajos niveles, formando clústeres que se mantienen tanto a escala local como regional. Este comportamiento confirma la presencia de desigualdades educativas con una marcada dimensión territorial en Guayaquil.

T_LOGRO_SUP

#Grafico de correlograma para educación superior
I_CORR_T_LOGRO_SUP <- sp.correlogram(neighbours = nbCenso,
                         var = sectores_censales_BDD$T_LOGRO_SUP,
                         order = 5, # niveles de correlación
                         method = "I", #método moran
                         style = "W",  #pesos estandarizados
                         zero.policy = TRUE)
#mostramos resultados.
I_CORR_T_LOGRO_SUP
## Spatial correlogram for sectores_censales_BDD$T_LOGRO_SUP 
## method: Moran's I
##               estimate   expectation      variance standard deviate
## 1 (6357)  0.7800632123 -0.0001573317  0.0000525350           107.64
## 2 (6357)  0.6787113703 -0.0001573317  0.0000218063           145.38
## 3 (6357)  0.5927478191 -0.0001573317  0.0000123183           168.93
## 4 (6357)  0.5269896944 -0.0001573317  0.0000083365           182.57
## 5 (6357)  0.4708121966 -0.0001573317  0.0000063573           186.79
##                Pr(I) two sided    
## 1 (6357) < 0.00000000000000022 ***
## 2 (6357) < 0.00000000000000022 ***
## 3 (6357) < 0.00000000000000022 ***
## 4 (6357) < 0.00000000000000022 ***
## 5 (6357) < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El correlograma espacial de Moran para T_LOGRO_SUP evidencia una autocorrelación positiva muy fuerte y altamente significativa en todos los órdenes de vecindad analizados. Los valores elevados desde el primer hasta el quinto orden indican que los niveles de logro educativo superior presentan un patrón espacial intensamente estructurado, con clústeres amplios y consistentes tanto a nivel local como regional. Esto confirma que las zonas con mayores niveles de educación superior tienden a agruparse entre sí, mientras que los sectores con bajos niveles forman concentraciones igualmente definidas, reflejando desigualdades educativas con una marcada dimensión territorial en Guayaquil.

4.5.6. Diagrama de dispersión de Moran

T_LOGRO_BACH

#Diagrama de dispersión de Moran para Bachillerato
dpT_LOGRO_BACH <- moran.plot(sectores_censales_BDD$T_LOGRO_BACH, 
                 listw = lwCenso, 
                 zero.policy = TRUE,
                 labels = NULL)

La nube de puntos se distribuye mayoritariamente sobre una tendencia positiva, lo que indica que:

  • Sectores con altos valores de T_LOGRO_BACH representan clústeres de alto logro educativo a nivel de Bachillerato. Esto es H-H (Arriba derecha).
  • Sectores con bajo valores de T_LOGRO_BACH representan clústeres de bajo logro educativo a nivel de Bachilletato. Esto es L-L (Abajo izquierda).
  • Sectores con valores altos de Bachillerato rodeados de vecinos con valores bajos de Bachillerato. Esto es H-L (Arriba izquierda).
  • Sectores con valores bajo de Bachillerato rodeados de vecinos con valores altos de Bachilletato. Esto es L-H (Abajo derecha).

En conjunto, el Moran Scatterplot refuerza la conclusión de que la distribución espacial del logro educativo de bachillerato está fuertemente estructurada, con clústeres bien definidos y patrones de desigualdad territorial.

T_LOGRO_SUP

#Diagrama de dispersión de Moran para educación superior
dpT_LOGRO_SUP <- moran.plot(sectores_censales_BDD$T_LOGRO_SUP, 
                 listw = lwCenso, 
                 zero.policy = TRUE,
                 labels = NULL)

La nube de puntos se distribuye mayoritariamente sobre una tendencia positiva, lo que indica que:

  • Representan sectores con altos valores de logro educativo superior rodeados de vecinos también altos. Clústeres H-H.
  • Reflejan sectores donde el logro superior es bajo y los vecinos también presentan bajos niveles. Clústeres L-L.
  • Sectores con alto logro educativo rodeados de áreas de bajo nivel. Aunque pocos, estos puntos representan anomalías espaciales. Outliers H-L.
  • Sectores con bajo nivel de educación superior rodeados de sectores más favorecidos, estos puntos representan anomalías espaciales. Outliers L-H.

La nube de puntos está claramente estructurada en forma de bandas diferenciadas, lo que sugiere gradientes espaciales amplios en el nivel educativo. Además, la aglomeración en valores medios y altos indica que la ciudad presenta sectores con alta concentración de educación superior distribuidos de manera no aleatoria, reforzando la conclusión de que la distribución espacial del logro educativo superior en Guayaquil es altamente desigual y sigue patrones territoriales marcados.

En conjunto, estos resultados muestran que la distribución del logro educativo no es aleatoria, sino que responde a un patrón espacial robusto que refuerza desigualdades socioeducativas territoriales y subraya la importancia de incorporar el enfoque espacial en el análisis y la planificación educativa en la ciudad.

4.6. Análisis de Autocorrelación espacial Local I Moran Local (Ii)

Se calcula el índice para cada una de las variables como se muestra a continuación.

#Valores de moran local para Bachillerato
LOC_MORAN_T_LOGRO_BACH <- localmoran(sectores_censales_BDD$T_LOGRO_BACH,
                                     listw = lwCenso)

#Valores de moran local para educación superior
LOC_MORAN_T_LOGRO_SUP <- localmoran(sectores_censales_BDD$T_LOGRO_SUP,
                                     listw = lwCenso)

Si se requiere analizar el resultado individual (caso que no es aplicado en este momento) cada una de las variables que arroja LISA significa lo siguiente:

  • Ii: Valor del Moran local
  • E.Ii: esperanza del Moran local
  • Var.Ii: varianza del Moran local
  • Z.Ii: desvío estándar del Moran local
  • Pr(…): p-value del Moran local

El resultado de cada índice se lo utiliza para estudiar el comportamiento de la asociación espacial local.

4.6.1. Análisis de asociación espacial local con LISA

T_LOGRO_BACH

#definimos las unidades espaciales respecto al cuadrante LISA
lisaT_LOGRO_BACH <- cbind(sectores_censales_BDD,                   # Base original
                  dpT_LOGRO_BACH[c("x", "wx")],                    # Variable y lag (estandarizado) del diagrama de dispersión
                  LOC_MORAN_T_LOGRO_BACH,                          # Valores de Moran Local
                  attributes(LOC_MORAN_T_LOGRO_BACH)$quadr) |>     # Cuadrantes LISA
  
  # Renombramos la columna "Pr(z != E(Ii))"
  rename(p = Pr.z....E.Ii..) |> 
  
  # Los valores no significativos se diferencian en otra categoría
  mutate(quad = ifelse(p > 0.05, 5, mean),
         quad = factor(quad, levels = 1:5,
                labels = c("Low-Low", "High-Low", 
                           "Low-High", "High-High", 
                           "No Signif")))

#Vemos el resultado de los cuadrantes
table(lisaT_LOGRO_BACH$quad)
## 
##   Low-Low  High-Low  Low-High High-High No Signif 
##      1049        31        39       717      4521

Con esta primera aproximación LISA podemos ver que tenemos 4521 unidades espaciales que no son significativas según LISA que corresponde al 71,118% del dataset.

Como siguiente punto, definimos los colores y realizamos un scaterplot para la variable T_LOGRO_BACH con ello evidenciados los valores de LISA que se encuentran en cuadrantes siendo estos significativos y no significativos.

#definimos los colores
LISA_col <- c("blue","skyblue2", "lightpink", "red", "white")
names(LISA_col) <- levels(lisaT_LOGRO_BACH$quad)

#graficamos scaterplot
lisaT_LOGRO_BACH |> 
  st_drop_geometry() |> 
  ggplot(aes(x = x, y = wx)) + 
  geom_hline(linetype = 2, yintercept = mean(lisaT_LOGRO_BACH$wx)) +
  geom_vline(linetype = 2, xintercept = mean(lisaT_LOGRO_BACH$x)) +
  geom_point(aes(fill = quad), shape = 21) + 
  geom_smooth(method = lm, se = F, linetype = 2, color = "tomato4" ) + 
  labs(x = "Variable", y = "Lag") +
  scale_fill_manual(values = LISA_col, drop = F) +
  theme_minimal()

Como vemos, los valores significativos se encuentran en los extremos tango HH como LL, tendiendo mayor contración los nivels LL. La línea de tendencia positiva indica autocorrelación espacial positiva, es decir: valores similares tienden a agruparse en el espacio. También tenemos valores no significativos ubicados en el centro de los cuadrates. Esto significa que no existe evidencia significativa suficiente de autocorrelación espacial local en esos lugares.

Como último punto del análisis, graficamos un mapa final con los valores de LISA local para la variable de Bachillerato.

# 1. Borde general del área (sin polígonos internos)
borde_general_T_LOGRO_BACH <- st_union(lisaT_LOGRO_BACH)

# 2. Mapa LISA sin líneas internas
ggplot() +
  geom_sf(data = lisaT_LOGRO_BACH, aes(fill = quad), color = NA) +
  geom_sf(data = borde_general_T_LOGRO_BACH, fill = NA, color = "grey30", size = 0.4) +
  scale_fill_manual(values = LISA_col, name = "Asociación local") +
  theme_void() +
  labs(
    title = "Autocorrelación espacial LOCAL",
    subtitle = "Nivel de Bachillerato",
    caption = "Fuente: Censo de Población y Vivienda 2022 - INEC"
  ) +
  theme_void() +
  theme(
    legend.position = "right",
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10))

Se puede observar zonas con valores bajos de bachillerato y zonas con altos valores que son estadísticamente significativos El resto de valores no están representados porque no tienen rigor estadístico para definir un grado de bachillerato.

T_LOGRO_SUP

#definimos las unidades espaciales respecto al cuadrante LISA
lisaT_LOGRO_SUP <- cbind(sectores_censales_BDD,                   # Base original
                  dpT_LOGRO_SUP[c("x", "wx")],                    # Variable y lag (estandarizado) del diagrama de dispersión
                  LOC_MORAN_T_LOGRO_SUP,                          # Valores de Moran Local
                  attributes(LOC_MORAN_T_LOGRO_SUP)$quadr) |>     # Cuadrantes LISA
  
  # Renombramos la columna "Pr(z != E(Ii))"
  rename(p = Pr.z....E.Ii..) |> 
  
  # Los valores no significativos se diferencian en otra categoría
  mutate(quad = ifelse(p > 0.05, 5, mean),
         quad = factor(quad, levels = 1:5,
                labels = c("Low-Low", "High-Low", 
                           "Low-High", "High-High", 
                           "No Signif")))

#Vemos el resultado de los cuadrantes
table(lisaT_LOGRO_SUP$quad)
## 
##   Low-Low  High-Low  Low-High High-High No Signif 
##      1298        26        24      1352      3657

Con esta primera aproximación LISA podemos ver que tenemos 3657 unidades espaciales que no son significativas según LISA que corresponde al 56,953% del dataset. Siendo un valor inferior de aquel con niveles de bachillerato. Esto significa que existe mayor patrones locales que tienen autocorrelación espacial respecto a Bachillerato.

Como siguiente punto, realizamos un scaterplot para la variable T_LOGRO_SUP con ello evidenciados los valores de LISA que se encuentran en cuadrantes siendo estos significativos y no significativos.

#definimos colores para variable
names(LISA_col) <- levels(lisaT_LOGRO_SUP$quad)

#graficamos scaterplot
lisaT_LOGRO_SUP |> 
  st_drop_geometry() |> 
  ggplot(aes(x = x, y = wx)) + 
  geom_hline(linetype = 2, yintercept = mean(lisaT_LOGRO_SUP$wx)) +
  geom_vline(linetype = 2, xintercept = mean(lisaT_LOGRO_SUP$x)) +
  geom_point(aes(fill = quad), shape = 21) + 
  geom_smooth(method = lm, se = F, linetype = 2, color = "tomato4" ) + 
  labs(x = "Variable", y = "Lag") +
  scale_fill_manual(values = LISA_col, drop = F) +
  theme_minimal()

Como vemos, los valores significativos se encuentran en los extremos tango HH como LL pero en mayor proporción. La línea de tendencia positiva indica autocorrelación espacial positiva, es decir: valores similares tienden a agruparse en el espacio. También tenemos valores no significativos ubicados en el centro de los cuadrates. Esto significa que no existe evidencia significativa suficiente de autocorrelación espacial local en esos lugares.

Como último punto del análisis, graficamos un mapa final con los valores de LISA local para la variable educación superior.

# 1. Borde general del área (sin polígonos internos)
borde_general_T_LOGRO_SUP <- st_union(lisaT_LOGRO_SUP)

# 2. Mapa LISA sin líneas internas
ggplot() +
  geom_sf(data = lisaT_LOGRO_SUP, aes(fill = quad), color = NA) +
  geom_sf(data = borde_general_T_LOGRO_SUP, fill = NA, color = "grey30", size = 0.4) +
  scale_fill_manual(values = LISA_col, name = "Asociación local") +
  theme_void() +
  labs(
    title = "Autocorrelación espacial LOCAL",
    subtitle = "Nivel de Educación Superior",
    caption = "Fuente: Censo de Población y Vivienda 2022 - INEC"
  ) +
  theme_void() +
  theme(
    legend.position = "right",
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10))

Se puede observar zonas con valores bajos de educación superior en la parte sur de Guayaquil y noroeste de la urba. También se ve zonas con altos niveles de educación superior en el centro de la ciudad. Estos valores son estadísticamente significativos. El resto de valores no están representados porque no tienen rigor estadístico para definir un grado de educación superior.

Apreciamos también valores atípicos de altos con bajos y bajos con altos que no son muchos.

4.7. GWR de niveles de educación con variables socioeconómicas.

T_LOGRO_BACH

Para medir el grado en que los índices influyen espacialmente en los patrones encontrados con LISA aplicamos una Regresión Geograficamente Ponderada. El flujo de trabajo fue el siguiente:

  • Limpieza de datos removiendo valores NA que posiblemente existan.
  • Test de normalidad con Anderson-Darling.
  • Test de linealidad con Pearson.
  • Test de multicolinedad con VIF y eliminar variables redundantes.
  • Ajuste de un modelo lineal no espacial y evaluar los residuos si tienen autocorrelación espacial.
  • Mapear el R2local para medir explicación del modelo.

4.7.1. Limpieza de datos removiendo valores NA.

Aplicando este algoritmo, pudimos eliminar valores NA. Pasamos de tener un total de 6357 observaciones a 6355.

lisaT_LOGRO_BACH_OK <- lisaT_LOGRO_BACH |>
  select(ID_GEOGRAFICO, T_LOGRO_BACH, I_FEMINIDAD, T_TRABAJO_NO_REM, I_JUVENTUD, T_OCUPACION_SEXO,
         I_ADULTEZ, I_ADULTEZ_MAYOR, POB_OCUP, I_MASCULINIDAD, T_ANALFABETISMO, T_EMPLEO_FORMAL) |> 
  drop_na()

4.7.2. Test de Normalidad con Anderson-Darling.

El resultado del test de Anderson–Darling nos da las siguientes hipótesis:

  • H0: La distribución es normal.
  • H1: La distribución no es normal.
# Test de Anderson-Darling
ad.test(lisaT_LOGRO_BACH_OK$I_FEMINIDAD)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$I_FEMINIDAD
## A = 24.49, p-value < 0.00000000000000022
ad.test(lisaT_LOGRO_BACH_OK$T_TRABAJO_NO_REM)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$T_TRABAJO_NO_REM
## A = 391.11, p-value < 0.00000000000000022
ad.test(lisaT_LOGRO_BACH_OK$I_JUVENTUD)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$I_JUVENTUD
## A = 2.933, p-value = 0.0000002283
ad.test(lisaT_LOGRO_BACH_OK$T_OCUPACION_SEXO)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$T_OCUPACION_SEXO
## A = 67.613, p-value < 0.00000000000000022
ad.test(lisaT_LOGRO_BACH_OK$I_ADULTEZ)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$I_ADULTEZ
## A = 0.7001, p-value = 0.06758
ad.test(lisaT_LOGRO_BACH_OK$I_ADULTEZ_MAYOR)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$I_ADULTEZ_MAYOR
## A = 76.079, p-value < 0.00000000000000022
ad.test(lisaT_LOGRO_BACH_OK$POB_OCUP)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$POB_OCUP
## A = 30.228, p-value < 0.00000000000000022
ad.test(lisaT_LOGRO_BACH_OK$I_MASCULINIDAD)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$I_MASCULINIDAD
## A = 526.83, p-value < 0.00000000000000022
ad.test(lisaT_LOGRO_BACH_OK$T_ANALFABETISMO)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$T_ANALFABETISMO
## A = 36.752, p-value < 0.00000000000000022
ad.test(lisaT_LOGRO_BACH_OK$T_EMPLEO_FORMAL)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_OK$T_EMPLEO_FORMAL
## A = 1.0904, p-value = 0.007372

Solamente la variable I_Adultez cumple con el supuesto de normalidad. El resto de variables muestran un rechazo a la normalidad. Para análisis posteriores, consideramos transformaciones para intentar normalizar las variables sin incluir I_Adultez.

#declaramos las variables que vamos a normalizar
indices <- c("I_FEMINIDAD", "T_TRABAJO_NO_REM", "I_JUVENTUD", "T_OCUPACION_SEXO",
             "I_ADULTEZ_MAYOR", "POB_OCUP", "I_MASCULINIDAD",
             "T_ANALFABETISMO", "T_EMPLEO_FORMAL") 

# Lista donde se guardan los modelos de normalización
set.seed(123)
modelos_norm <- map(indices, ~ bestNormalize(lisaT_LOGRO_BACH_OK[[.x]], standardize = TRUE)) # declaramos también una estandarización de los datos.
#renombrados variables resultante con variables objetivo
names(modelos_norm) <- indices

# Aplicar las transformaciones
lisaT_LOGRO_BACH_TRANS <- lisaT_LOGRO_BACH_OK |> 
  mutate(across(all_of(indices), 
                ~ bestNormalize(.x, standardize = TRUE)$x.t, 
                .names = "{.col}_TRANS"))
#revisamos resultados
modelos_norm
## $I_FEMINIDAD
## Best Normalizing transformation with 6355 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.8174
##  - Box-Cox: 2.8008
##  - Center+scale: 1.8985
##  - Double Reversed Log_b(x+a): 3.7151
##  - Exp(x): 682.4666
##  - Log_b(x+a): 1.8227
##  - orderNorm (ORQ): 1.1582
##  - sqrt(x + a): 1.6402
##  - Yeo-Johnson: 2.715
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6355 nonmissing obs and ties
##  - 4225 unique values 
##  - Original quantiles:
##      0%     25%     50%     75%    100% 
##   7.141  99.059 105.405 112.596 176.923 
## 
## $T_TRABAJO_NO_REM
## Best Normalizing transformation with 6355 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 79.7048
##  - Center+scale: 83.2277
##  - Double Reversed Log_b(x+a): 86.9505
##  - Exp(x): 627.9709
##  - Log_b(x+a): 115.7382
##  - orderNorm (ORQ): 79.1143
##  - sqrt(x + a): 82.7413
##  - Yeo-Johnson: 81.3874
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6355 nonmissing obs and ties
##  - 788 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
##  0.000  0.000  0.714  1.562 21.875 
## 
## $I_JUVENTUD
## Best Normalizing transformation with 6355 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.554
##  - Box-Cox: 1.1164
##  - Center+scale: 1.1167
##  - Double Reversed Log_b(x+a): 1.7244
##  - Exp(x): 570.9429
##  - Log_b(x+a): 1.552
##  - orderNorm (ORQ): 1.1684
##  - sqrt(x + a): 1.175
##  - Yeo-Johnson: 1.123
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## Standardized Box Cox Transformation with 6355 nonmissing obs.:
##  Estimated statistics:
##  - lambda = 0.8972242 
##  - mean (before standardization) = 15.38733 
##  - sd (before standardization) = 2.438967 
## 
## $T_OCUPACION_SEXO
## Best Normalizing transformation with 6355 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 2.6822
##  - Box-Cox: 1.2804
##  - Center+scale: 3.8046
##  - Double Reversed Log_b(x+a): 7.6847
##  - Exp(x): 7.6587
##  - Log_b(x+a): 1.2833
##  - orderNorm (ORQ): 1.1287
##  - sqrt(x + a): 2.1143
##  - Yeo-Johnson: 1.3727
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6355 nonmissing obs and ties
##  - 2736 unique values 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## 0.125 0.549 0.659 0.816 1.632 
## 
## $I_ADULTEZ_MAYOR
## Best Normalizing transformation with 6355 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 2.8086
##  - Center+scale: 5.5623
##  - Double Reversed Log_b(x+a): 7.7189
##  - Exp(x): 662.1039
##  - Log_b(x+a): 3.0703
##  - orderNorm (ORQ): 1.0369
##  - sqrt(x + a): 1.9284
##  - Yeo-Johnson: 1.7908
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6355 nonmissing obs and ties
##  - 4952 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
##  0.000  4.796  8.439 12.080 50.000 
## 
## $POB_OCUP
## Best Normalizing transformation with 6355 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 2.0733
##  - Box-Cox: 1.5224
##  - Center+scale: 2.1267
##  - Double Reversed Log_b(x+a): 5.9099
##  - Log_b(x+a): 2.0775
##  - orderNorm (ORQ): 1.3025
##  - sqrt(x + a): 1.412
##  - Yeo-Johnson: 1.5281
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6355 nonmissing obs and ties
##  - 272 unique values 
##  - Original quantiles:
##   0%  25%  50%  75% 100% 
##   14  120  145  170  888 
## 
## $I_MASCULINIDAD
## Best Normalizing transformation with 6355 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.7933
##  - Box-Cox: 2.7372
##  - Center+scale: 8.9696
##  - Double Reversed Log_b(x+a): 1.307
##  - Log_b(x+a): 1.7925
##  - orderNorm (ORQ): 1.2107
##  - sqrt(x + a): 3.3474
##  - Yeo-Johnson: 2.7132
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6355 nonmissing obs and ties
##  - 4225 unique values 
##  - Original quantiles:
##       0%      25%      50%      75%     100% 
##   56.522   88.813   94.872  100.950 1400.357 
## 
## $T_ANALFABETISMO
## Best Normalizing transformation with 6355 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 2.3611
##  - Box-Cox: 1.0672
##  - Center+scale: 2.4207
##  - Double Reversed Log_b(x+a): 4.0685
##  - Exp(x): 651.2096
##  - Log_b(x+a): 2.5033
##  - orderNorm (ORQ): 1.06
##  - sqrt(x + a): 1.0563
##  - Yeo-Johnson: 1.0525
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## Standardized Yeo-Johnson Transformation with 6355 nonmissing obs.:
##  Estimated statistics:
##  - lambda = 0.3924112 
##  - mean (before standardization) = 3.638938 
##  - sd (before standardization) = 1.016695 
## 
## $T_EMPLEO_FORMAL
## Best Normalizing transformation with 6355 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.7088
##  - Box-Cox: 1.0003
##  - Center+scale: 1.0355
##  - Double Reversed Log_b(x+a): 1.7829
##  - Exp(x): 672.0728
##  - Log_b(x+a): 1.7081
##  - orderNorm (ORQ): 1.1077
##  - sqrt(x + a): 1.2028
##  - Yeo-Johnson: 0.9994
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## Standardized Yeo-Johnson Transformation with 6355 nonmissing obs.:
##  Estimated statistics:
##  - lambda = 1.181104 
##  - mean (before standardization) = 107.461 
##  - sd (before standardization) = 18.21247

Se aplicó el procedimiento bestNormalize() con validación cruzada (10 folds, 5 repeticiones) para seleccionar automáticamente la transformación que mejor aproxima la normalidad en cada variable. El criterio de selección fue el estadístico Pearson P/df, donde valores más bajos indican mejor normalidad.

Como resumen se tiene el siguiente resultado:

  • T_Analfabetismo y $I_JUVENTUD se ajustaron mejor a la transformación Standardized Yeo-Johnson.
  • T_Empleo_Formal se aplicó la transformación Standardized Box Cox.
  • El restante de variables requirieron transformación de tipo OrdenNorm.

Para confirmar la estandarización realizada para aproximar los datos a una distribución normal, aplicamos Anderson-Darling nuevamente sin considerar I_Adultez.

#probamos normalidad con variables transformadas.
ad.test(lisaT_LOGRO_BACH_TRANS$I_FEMINIDAD_TRANS) 
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_TRANS$I_FEMINIDAD_TRANS
## A = 0.013287, p-value = 1
ad.test(lisaT_LOGRO_BACH_TRANS$T_TRABAJO_NO_REM_TRANS)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_TRANS$T_TRABAJO_NO_REM_TRANS
## A = 209.49, p-value < 0.00000000000000022
ad.test(lisaT_LOGRO_BACH_TRANS$I_JUVENTUD_TRANS)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_TRANS$I_JUVENTUD_TRANS
## A = 2.933, p-value = 0.0000002283
ad.test(lisaT_LOGRO_BACH_TRANS$T_OCUPACION_SEXO_TRANS)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_TRANS$T_OCUPACION_SEXO_TRANS
## A = 0.015773, p-value = 1
ad.test(lisaT_LOGRO_BACH_TRANS$I_ADULTEZ_MAYOR_TRANS)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_TRANS$I_ADULTEZ_MAYOR_TRANS
## A = 0.00081602, p-value = 1
ad.test(lisaT_LOGRO_BACH_TRANS$POB_OCUP_TRANS)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_TRANS$POB_OCUP_TRANS
## A = 0.19107, p-value = 0.8976
ad.test(lisaT_LOGRO_BACH_TRANS$I_MASCULINIDAD_TRANS)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_TRANS$I_MASCULINIDAD_TRANS
## A = 0.013287, p-value = 1
ad.test(lisaT_LOGRO_BACH_TRANS$T_ANALFABETISMO_TRANS)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_TRANS$T_ANALFABETISMO_TRANS
## A = 0.97043, p-value = 0.01456
ad.test(lisaT_LOGRO_BACH_TRANS$T_EMPLEO_FORMAL_TRANS)
## 
##  Anderson-Darling normality test
## 
## data:  lisaT_LOGRO_BACH_TRANS$T_EMPLEO_FORMAL_TRANS
## A = 1.0904, p-value = 0.007372

Como resultado, obtuvimos que las variables I_FEMINIDAD, T_OCUPACION_SEXO, I_ADULTEZ_MAYOR, POB_OCUP y I_MASCULINIDAD lograron conseguir una distribución normal. Sin embargo, T_TRABAJO_NO_REM, I_JUVENTUD, T_ANALFABETISMO, T_EMPLEO_FORMAL todavía presetan asimetría tras normalización. Dicho esto, no se van a incluir en el análisis posterior.

Nos quedamos en total con 6 variables que siguen una normalidad.

4.7.3. Test de linealidad con Pearson.

Como siguiente punto del análisis, comprobamos la linealidad de los datos utilizando el método de Pearson como se detalla a continuación

#seleccionamos solo variables normales
lisaT_LOGRO_BACH_TRANS_Final <-  lisaT_LOGRO_BACH_TRANS |> 
  select(T_LOGRO_BACH, I_ADULTEZ, I_FEMINIDAD_TRANS, T_OCUPACION_SEXO_TRANS, I_ADULTEZ_MAYOR_TRANS, POB_OCUP_TRANS, I_MASCULINIDAD_TRANS) |> 
  st_drop_geometry()

#Determinar linealidad entre variables normales
cor_matrix <- cor(lisaT_LOGRO_BACH_TRANS_Final, method = "pearson")

#revisamos resultado
cor_matrix
##                        T_LOGRO_BACH   I_ADULTEZ I_FEMINIDAD_TRANS
## T_LOGRO_BACH             1.00000000 -0.60401086       -0.24094132
## I_ADULTEZ               -0.60401086  1.00000000        0.11322499
## I_FEMINIDAD_TRANS       -0.24094132  0.11322499        1.00000000
## T_OCUPACION_SEXO_TRANS  -0.59622869  0.54497071        0.49933148
## I_ADULTEZ_MAYOR_TRANS   -0.51712312  0.45852694        0.17712915
## POB_OCUP_TRANS           0.02616488  0.04685959        0.05948868
## I_MASCULINIDAD_TRANS     0.24094132 -0.11322499       -1.00000000
##                        T_OCUPACION_SEXO_TRANS I_ADULTEZ_MAYOR_TRANS
## T_LOGRO_BACH                       -0.5962287            -0.5171231
## I_ADULTEZ                           0.5449707             0.4585269
## I_FEMINIDAD_TRANS                   0.4993315             0.1771292
## T_OCUPACION_SEXO_TRANS              1.0000000             0.4876353
## I_ADULTEZ_MAYOR_TRANS               0.4876353             1.0000000
## POB_OCUP_TRANS                      0.1107796            -0.0650856
## I_MASCULINIDAD_TRANS               -0.4993315            -0.1771292
##                        POB_OCUP_TRANS I_MASCULINIDAD_TRANS
## T_LOGRO_BACH               0.02616488           0.24094132
## I_ADULTEZ                  0.04685959          -0.11322499
## I_FEMINIDAD_TRANS          0.05948868          -1.00000000
## T_OCUPACION_SEXO_TRANS     0.11077960          -0.49933148
## I_ADULTEZ_MAYOR_TRANS     -0.06508560          -0.17712915
## POB_OCUP_TRANS             1.00000000          -0.05948868
## I_MASCULINIDAD_TRANS      -0.05948868           1.00000000

Podemos observar que la variable dependiente (T_LOGRO_BACH) tiene una relación negativa con I_ADULTEZ, I_FEMINIDAD_TRANS, T_OCUPACION_SEXO_TRANS, I_ADULTEZ_MAYOR_TRANS. En cambio, tiene una relación positiva con I_MASCULINIDAD_TRANS y no tiene ninguna relación con POB_OCUP_TRANS, lo que sugiere los siguientes supuestos:

  • Sectores con mayor proporción de adultos presentan menor logro educativo.
  • Mayor proporción de gente ocupada menor logro educativo.
  • Envejecimiento poblacional representa menor logro académico.
  • Relación débil negativa de feminidad que representa menor logro académico y viceversa para masculidad.
  • No existe relación con población ocupada.
  • Podemos apreciar multicolinealidad perfecta entre los índices de género (I_FEMINIDAD vs I_MASCULINIDAD). Nos quedamos con I_FEMINIDAD_TRANS.

El logro educativo disminuye con variables ligadas al envejecimiento demográfico y al mercado laboral. Esto sugiere un territorio donde zonas con mayor presión laboral o con menor población joven muestran niveles más bajos de logro académico.

4.7.4. Test de multicolinedad con VIF y eliminar variables redundantes.

Como siguiente paso, integramos criterios VIF (multicolinealidad) para quedarnos con el mejor conjunto de predictores. Eliminamos las variables que no aportan al modelo.

#Aseguramos que solo entren variables numéricas y sin geometría (si fuese sf)
modelo_inical <- lm(T_LOGRO_BACH ~ I_ADULTEZ + I_FEMINIDAD_TRANS + 
                 T_OCUPACION_SEXO_TRANS + I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS,
                    data = lisaT_LOGRO_BACH_TRANS_Final)

#vemos el resultado del modelo
summary(modelo_inical)
## 
## Call:
## lm(formula = T_LOGRO_BACH ~ I_ADULTEZ + I_FEMINIDAD_TRANS + T_OCUPACION_SEXO_TRANS + 
##     I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS, data = lisaT_LOGRO_BACH_TRANS_Final)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -50.171  -4.642   0.563   4.920  26.645 
## 
## Coefficients:
##                        Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)            78.62222    1.08364  72.554 < 0.0000000000000002 ***
## I_ADULTEZ              -0.79936    0.02599 -30.760 < 0.0000000000000002 ***
## I_FEMINIDAD_TRANS      -0.15896    0.10860  -1.464                0.143    
## T_OCUPACION_SEXO_TRANS -3.21349    0.13545 -23.725 < 0.0000000000000002 ***
## I_ADULTEZ_MAYOR_TRANS  -2.08576    0.11007 -18.949 < 0.0000000000000002 ***
## POB_OCUP_TRANS          0.66795    0.09325   7.163    0.000000000000877 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.315 on 6349 degrees of freedom
## Multiple R-squared:  0.5021, Adjusted R-squared:  0.5017 
## F-statistic:  1281 on 5 and 6349 DF,  p-value: < 0.00000000000000022
# Diagnósticos
par(mfrow=c(2,2))
plot(modelo_inical)

El resultado nos indica que todos los predictos son significativos menos I_FEMINIDAD_TRANS. Por lo cual, no seria necesaria incluirla en el estudio. Continuamos con la revisión VIF para identificar si las variables restantes están dentro del margen indicado para continuar.

#revisamos valores vif 
vif_vals <- vif(modelo_inical)
#vemos resultado
vif_vals
##              I_ADULTEZ      I_FEMINIDAD_TRANS T_OCUPACION_SEXO_TRANS 
##               1.595174               1.400531               2.178451 
##  I_ADULTEZ_MAYOR_TRANS         POB_OCUP_TRANS 
##               1.438553               1.032327

Todas las variables tienen un valor menor a 5. Descartamos I_FEMINIDAD_TRANS por no ser significativa y nos queda:

  • I_ADULTEZ
  • T_OCUPACION_SEXO_TRANS
  • I_ADULTEZ_MAYOR_TRANS
  • POB_OCUP_TRANS

4.7.5. Ajuste de modelo no lineal y evaluación de residuos si tienen autocorrelación espacial.

Luego de aplicar un modelo lineal inicial para ajustar multicolinealidad. Ejecutamos el modelo final y evaluamos si existe dependencia espacial en los residuos. Si es favorable, seguimos con el modelo GWR.

#ejecución del modelo final
modelo_final <- lm(T_LOGRO_BACH ~ I_ADULTEZ + T_OCUPACION_SEXO_TRANS + I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS,
                    data = lisaT_LOGRO_BACH_TRANS_Final)

#autocorrelación espacial del error
coords <- st_centroid(st_geometry(lisaT_LOGRO_BACH_OK)) |> st_coordinates()
vecinos <- knearneigh(coords, k=5) |> knn2nb() #usamos k-nearest
lw <- nb2listw(vecinos)

#aplicamos el test de moran a los residuos
moran.test(residuals(modelo_final), lw)
## 
##  Moran I test under randomisation
## 
## data:  residuals(modelo_final)  
## weights: lw    
## 
## Moran I statistic standard deviate = 55.693, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     0.42290392498    -0.00015738118     0.00005770454

El resultado que muestra el test de Moran con los residuos es que el P-valor es extremadamente bajo lo que implica que la autocorrelación espacial no es aleatoria. El valor de Moran’s I es aproximadamente 0.4229, positivo y considerable. Esto se traduce que mis residuos están agrupados espacialmente y esto viola es supuesto de independencia de los residuos.

En otras palabras, mi variable dependiente tiene una estructura espacial que el modelo lineal no está capturando. Por lo tanto, necesitamos avanzar hacia un modelo regresión espacial.

Primero seleccionamos el mejor kernel para ajustar el modelo

set.seed(123)
# 1.- Proyectamos a UTM 17S
LisaT_LOGRO_BACH_sp <- st_transform(lisaT_LOGRO_BACH_TRANS, 32717) |>  as("Spatial")

#2.- Seleccionamos bandwidth óptimo con spgwr
bw_gwr_LOGRO_BACH <- gwr.sel(
  T_LOGRO_BACH ~ I_ADULTEZ + T_OCUPACION_SEXO_TRANS + I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS, 
  data = LisaT_LOGRO_BACH_sp,
  coords = st_coordinate(LisaT_LOGRO_BACH_sp),
  adapt = T,
  gweight = gwr.Gauss,  # kernel gaussiano
  verbose = T, 
  show.error.messages = T)
## Adaptive q: 0.381966 CV score: 310290.4 
## Adaptive q: 0.618034 CV score: 324535.1 
## Adaptive q: 0.236068 CV score: 299727.2 
## Adaptive q: 0.145898 CV score: 289859.1 
## Adaptive q: 0.09016994 CV score: 280057.8 
## Adaptive q: 0.05572809 CV score: 270794 
## Adaptive q: 0.03444185 CV score: 261747.9 
## Adaptive q: 0.02128624 CV score: 252154.6 
## Adaptive q: 0.01315562 CV score: 240512.7 
## Adaptive q: 0.008130619 CV score: 225901.1 
## Adaptive q: 0.005024999 CV score: 217415.5 
## Adaptive q: 0.00310562 CV score: 213654.2 
## Adaptive q: 0.001919379 CV score: 212142.3 
## Adaptive q: 0.001186241 CV score: 215132.3 
## Adaptive q: 0.002284003 CV score: 212510.7 
## Adaptive q: 0.001639345 CV score: 212484.1 
## Adaptive q: 0.001960069 CV score: 212136.3 
## Adaptive q: 0.002000759 CV score: 212143.5 
## Adaptive q: 0.001960069 CV score: 212136.3
#revisamos resultado
bw_gwr_LOGRO_BACH
## [1] 0.001960069

Como resultado obtuvimos el bandwidth óptimo adaptativo (q) de 0.001960069 para mi modelo GWR. Primera apreciación al salir un valor muy bajo es que el modelo detecta fuerte heterogeneidad espacial. Las relaciones entre variables socioeducativas cambian bastante entre sectores censales.

Luego ajustamos el modelo GWR con el bandwidth elegido

set.seed(123)
#GWR con el kernel adaptativo.
modelo_gwr_sp_LOGRO_BACH <- gwr(
    T_LOGRO_BACH ~ I_ADULTEZ + T_OCUPACION_SEXO_TRANS + I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS, 
    data = LisaT_LOGRO_BACH_sp,
    coords = st_coordinates(LisaT_LOGRO_BACH_sp),
    adapt = bw_gwr_LOGRO_BACH)

4.7.6. Mapear el R2 local para medir explicación del modelo.

T_LOGRO_BACH

Para mapear el R2 local debemos primero transformar el modelo GWR en un ST y extraemos los coeficientes.

#guardamos coeficientes sobre st para graficar.
SDF <- modelo_gwr_sp_LOGRO_BACH$SDF
SDF <- st_as_sf(modelo_gwr_sp_LOGRO_BACH$SDF)

#Graficamos coeficientes
ggplot(SDF) +
  geom_sf(aes(fill = localR2), color = NA) +
  geom_sf(data = borde_general_T_LOGRO_SUP, fill = NA, color = "grey30", size = 0.4) +
  scale_fill_viridis_c(name="R2 local", direction=-1) +
  labs(
    title = "R2 Local del modelo GWR",
    fill = "R2 local",
    subtitle = "Nivel de Bachillerato",
    caption = "Fuente: Censo de Población y Vivienda 2022 - INEC") +
  theme_void() +
  theme(
    legend.position = "right",
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10))

El mapa de R2 local muestra con claridad la distribución del modeo GWR donde explica mejor la variación de la variable dependiente Logro Bachillerato indicando: - Los valores de R2 local alto (Azúl oscuro-verde oscuro) explica bien la fuerte variación negativa entre T_LOGRO_BACH y los predictores en esa zona. Esta explicación se concentra mejor en el centro de guayaquil y en la parte oeste. - Los valores de R2 local bajo (Verde claro - Amarillo) explica poco la variación entre L_LOGRO_BACH y los predictores en esa zona. Existiendo otros factores socioeconómicos relevantes.

# Calcular centroides para pasar atributos de una entidad a otra.
centroid_sdf <- st_centroid(SDF)

# Join espacial usando centroides
centroid_sdf <- st_transform(centroid_sdf, 32717) #transformamos al mismo sistema de referencia
lisaT_LOGRO_BACH <- st_transform(lisaT_LOGRO_BACH, 32717) #transformamos al mismo sistema de referencia

#realizamos el spatial join
centroid_sdf <- st_join(centroid_sdf, lisaT_LOGRO_BACH["quad"]) 

#Quitamos geometria para el join por atributo
centroid_data <- centroid_sdf |> 
  st_drop_geometry() |> 
  select(localR2, quad)

# Agregar la variable LISA a SDF
SDF <- SDF |> right_join(centroid_data, by = "localR2")

# Clasificar R2 local en terciles
breaks_R2 <- classIntervals(SDF$localR2, n = 3, style = "quantile")$brks
SDF$R2_cat <- cut(SDF$localR2,
                     breaks = breaks_R2,
                     include.lowest = TRUE,
                     labels = c("Bajo", "Medio", "Alto"))

#Crear variable bivariada combinada
SDF <-  SDF |> drop_na()
SDF$bivar <- paste(SDF$R2_cat, SDF$quad, sep = " | ")

#Paleta de colores
bivar_colors <- c(
  "Bajo | Low-Low"   = "#e8e8e8",
  "Medio | Low-Low"  = "#b5c0da",
  "Alto | Low-Low"   = "#6c83b5",
  
  "Bajo | Low-High"  = "#cbb8d7",
  "Medio | Low-High" = "#9f8bbd",
  "Alto | Low-High"  = "#73579d",
  
  "Bajo | High-Low"  = "#e4acac",
  "Medio | High-Low" = "#c85a5a",
  "Alto | High-Low"  = "#b10101",
  
  "Bajo | High-High" = "#d4a3c6",
  "Medio | High-High"= "#a05a8c",
  "Alto | High-High" = "#7a0177",
   
  "Bajo | No Signif"   = "white",
  "Medio | No Signif"  = "white",
  "Alto | No Signif"   = "white"
)

#Plot final del mapa bivariado
ggplot(SDF) +
  geom_sf(aes(fill = bivar), color = NA) +
  geom_sf(data = borde_general_T_LOGRO_SUP, fill = NA, color = "grey30", size = 0.4) +
  scale_fill_manual(values = bivar_colors) +
  
  labs(
    title = "Mapa Bivariado: R2 Local del GWR vs LISA Local",
    subtitle = "NIVEL BACHILLERATO: Modelo explicativo + estructura espacial",
    fill = "R2 Local | LISA"
  ) +
  
  theme_void() +
  theme(
    legend.position = "right",
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10))

Se puede observar que el modelo GWR con valores altos de R2 local puede explicar bien los clústers Low-Low de LISA y pocas zonas High-High. Gran parte del modelo existen zonas de color blanco que fueron categorizadas como no significativas.

T_LOGRO_SUP

Respecto a la variable dependiente Logro Superior, se aplicó cierta parte del procedimiento efectuado con T_LOGRO_BACH y copiamos ciertos resultado. La metodologísse detalla a continuación:

  • Limpieza de datos removiendo valores NA que posiblemente existan.
  • Test de linealidad con Pearson.
  • Test de multicolinedad con VIF y eliminar variables redundantes.
  • Ajuste de un modelo lineal no espacial y evaluar los residuos si tienen autocorrelación espacial.
  • Mapear el R2local para medir explicación del modelo.

4.7.6.1. Limpieza de datos removiendo valores NA.

Aplicando este algoritmo, pudimos eliminar valores NA. Pasamos de tener un total de 6357 observaciones a 6355.

lisaT_LOGRO_SUP_OK <- lisaT_LOGRO_SUP |>
  select(ID_GEOGRAFICO, T_LOGRO_SUP, I_FEMINIDAD, T_TRABAJO_NO_REM, I_JUVENTUD, T_OCUPACION_SEXO,
         I_ADULTEZ, I_ADULTEZ_MAYOR, POB_OCUP, I_MASCULINIDAD, T_ANALFABETISMO, T_EMPLEO_FORMAL) |> 
  drop_na()

Unificamos los valores transformados previamiente eliminando T_LOGRO_BACH y agregando T_LOGRO_SUP que es motivo de análisis.

# Crear objeto final uniendo capas
lisaT_LOGRO_SUP_TRANS <- lisaT_LOGRO_BACH_TRANS |> 
  select(-T_LOGRO_BACH) |>   #remover T_LOGRO_BACH
  left_join(
    lisaT_LOGRO_SUP_OK |> 
      st_drop_geometry() |>     # quitamos geometría del origen
      select(ID_GEOGRAFICO, T_LOGRO_SUP),   # nos quedamos SOLO con esa variable
    by = "ID_GEOGRAFICO") #join por esa variable

4.7.6.2. Test de linealidad con Pearson.

Como siguiente punto del análisis, comprobamos la linealidad de los datos utilizando el método de Pearson como se detalla a continuación

#seleccionamos solo variables normales
lisaT_LOGRO_SUP_TRANS_Final <-  lisaT_LOGRO_SUP_TRANS |> 
  select(T_LOGRO_SUP, I_ADULTEZ, I_FEMINIDAD_TRANS, T_OCUPACION_SEXO_TRANS, I_ADULTEZ_MAYOR_TRANS, POB_OCUP_TRANS, I_MASCULINIDAD_TRANS) |> 
  st_drop_geometry()

#Determinar linealidad entre variables normales
cor_matrix <- cor(lisaT_LOGRO_SUP_TRANS_Final, method = "pearson")

#revisamos resultado
cor_matrix
##                        T_LOGRO_SUP   I_ADULTEZ I_FEMINIDAD_TRANS
## T_LOGRO_SUP             1.00000000  0.64367330        0.28827614
## I_ADULTEZ               0.64367330  1.00000000        0.11322499
## I_FEMINIDAD_TRANS       0.28827614  0.11322499        1.00000000
## T_OCUPACION_SEXO_TRANS  0.70339728  0.54497071        0.49933148
## I_ADULTEZ_MAYOR_TRANS   0.48754810  0.45852694        0.17712915
## POB_OCUP_TRANS          0.03884235  0.04685959        0.05948868
## I_MASCULINIDAD_TRANS   -0.28827614 -0.11322499       -1.00000000
##                        T_OCUPACION_SEXO_TRANS I_ADULTEZ_MAYOR_TRANS
## T_LOGRO_SUP                         0.7033973             0.4875481
## I_ADULTEZ                           0.5449707             0.4585269
## I_FEMINIDAD_TRANS                   0.4993315             0.1771292
## T_OCUPACION_SEXO_TRANS              1.0000000             0.4876353
## I_ADULTEZ_MAYOR_TRANS               0.4876353             1.0000000
## POB_OCUP_TRANS                      0.1107796            -0.0650856
## I_MASCULINIDAD_TRANS               -0.4993315            -0.1771292
##                        POB_OCUP_TRANS I_MASCULINIDAD_TRANS
## T_LOGRO_SUP                0.03884235          -0.28827614
## I_ADULTEZ                  0.04685959          -0.11322499
## I_FEMINIDAD_TRANS          0.05948868          -1.00000000
## T_OCUPACION_SEXO_TRANS     0.11077960          -0.49933148
## I_ADULTEZ_MAYOR_TRANS     -0.06508560          -0.17712915
## POB_OCUP_TRANS             1.00000000          -0.05948868
## I_MASCULINIDAD_TRANS      -0.05948868           1.00000000

Podemos observar que la variable dependiente (T_LOGRO_SUP) tiene una relación positiva FUERTE con T_OCUPACION_SEXO_TRANS, I_ADULTEZ. Relación positiva moderada con I_ADULTEZ_MAYOR_TRANS y una relación débil positiva con I_FEMINIDAD_tRANS. En cambio, tiene una relación negativa débil con I_MASCULINIDAD_TRANS y no tiene ninguna relación con POB_OCUP_TRANS, lo que sugiere los siguientes supuestos:

  • Sectores con mayor ocupación por sexo explica mucho mejor el logro educativo SUPERIOR.
  • Mayor proporción de adultos mayor logro educativo SUPERIOR.
  • Envejecimiento poblacional representa mayor logro SUPERIOR.
  • Relación débil de feminidad que representa mayor logro SUPERIOR y viceversa para masculidad.
  • No existe relación con población ocupada.
  • Podemos apreciar multicolinealidad perfecta entre los índices de género (I_FEMINIDAD vs I_MASCULINIDAD).

El logro educativo SUPERIOR aumenta con variables ligadas a la ocupación y población adulta. Esto sugiere un territorio donde zonas con mayor presión laboral o con menor población joven muestran niveles más altos de logro académico SUPERIOR. Esta relación, es positiva, respecto a Bachillerato que era negativa.

4.7.6.3. Test de multicolinedad con VIF y eliminar variables redundantes.

Como siguiente paso, integramos criterios VIF (multicolinealidad) para quedarnos con el mejor conjunto de predictores. Eliminamos las variables que no aportan al modelo.

#Aseguramos que solo entren variables numéricas y sin geometría (si fuese sf)
modelo_inicalSUP <- lm(T_LOGRO_SUP ~ I_ADULTEZ + I_FEMINIDAD_TRANS + 
                 T_OCUPACION_SEXO_TRANS + I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS,
                    data = lisaT_LOGRO_SUP_TRANS_Final)

#vemos el resultado del modelo
summary(modelo_inicalSUP)
## 
## Call:
## lm(formula = T_LOGRO_SUP ~ I_ADULTEZ + I_FEMINIDAD_TRANS + T_OCUPACION_SEXO_TRANS + 
##     I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS, data = lisaT_LOGRO_SUP_TRANS_Final)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.780  -8.054  -1.216   6.945  89.238 
## 
## Coefficients:
##                         Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)            -19.41778    1.69940 -11.426 < 0.0000000000000002 ***
## I_ADULTEZ                1.38931    0.04075  34.090 < 0.0000000000000002 ***
## I_FEMINIDAD_TRANS       -0.04611    0.17032  -0.271              0.78658    
## T_OCUPACION_SEXO_TRANS   8.55630    0.21242  40.281 < 0.0000000000000002 ***
## I_ADULTEZ_MAYOR_TRANS    1.79923    0.17262  10.423 < 0.0000000000000002 ***
## POB_OCUP_TRANS          -0.41480    0.14623  -2.837              0.00457 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.47 on 6349 degrees of freedom
## Multiple R-squared:  0.5992, Adjusted R-squared:  0.5989 
## F-statistic:  1899 on 5 and 6349 DF,  p-value: < 0.00000000000000022
# Diagnósticos
par(mfrow=c(2,2))
plot(modelo_inicalSUP)

El resultado nos indica que todos los predictos son significativos menos I_FEMINIDAD_TRANS e I_MASCULINIDAD_TRANS. Por lo cual, no seria necesaria incluirla en el estudio. El gráfico de Q-Q Residuals muestra bien la tendencia normalidad y linealidad de las variables; asi como, Residuals vs Fiteed que están regularmente dispersos sobre la media.

Continuamos con la revisión VIF para identificar si las variables restantes están dentro del margen indicado para continuar.

#revisamos valores vif 
vif_valsSUP <- vif(modelo_inicalSUP)
#vemos resultado
vif_valsSUP
##              I_ADULTEZ      I_FEMINIDAD_TRANS T_OCUPACION_SEXO_TRANS 
##               1.595174               1.400531               2.178451 
##  I_ADULTEZ_MAYOR_TRANS         POB_OCUP_TRANS 
##               1.438553               1.032327

Todas las variables tienen un valor menor a 5. Descartamos I_FEMINIDAD_TRANS por no ser significativa y nos queda:

  • I_ADULTEZ
  • T_OCUPACION_SEXO_TRANS
  • I_ADULTEZ_MAYOR_TRANS
  • POB_OCUP_TRANS

4.7.6.4. Ajuste de modelo no lineal y evaluación de residuos si tienen autocorrelación espacial.

Luego de aplicar un modelo lineal inicial para ajustar multicolinealidad. Ejecutamos el modelo final y evaluamos si existe dependencia espacial en los residuos. Si es favorable, seguimos con el modelo GWR.

#ejecución del modelo final
modelo_finalSUP <- lm(T_LOGRO_SUP ~ I_ADULTEZ + T_OCUPACION_SEXO_TRANS + I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS,
                    data = lisaT_LOGRO_SUP_TRANS_Final)

#autocorrelación espacial del error
coordsSUP <- st_centroid(st_geometry(lisaT_LOGRO_SUP_OK)) |> st_coordinates()
vecinosSUP <- knearneigh(coords, k=5) |> knn2nb() #usamos k-nearest
lwSUP <- nb2listw(vecinosSUP)

#aplicamos el test de moran a los residuos
moran.test(residuals(modelo_finalSUP), lwSUP)
## 
##  Moran I test under randomisation
## 
## data:  residuals(modelo_finalSUP)  
## weights: lwSUP    
## 
## Moran I statistic standard deviate = 38.994, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     0.29604209208    -0.00015738118     0.00005770109

El resultado que muestra el test de Moran con los residuos es que el P-valor es extremadamente bajo lo que implica que la autocorrelación espacial no es aleatoria. El valor de Moran’s I es aproximadamente 0.2960, positivo pero menor que Bachillerato. Esto se traduce que mis residuos están agrupados espacialmente y esto viola es supuesto de independencia de los residuos.

En otras palabras, mi variable dependiente tiene una estructura espacial que el modelo lineal no está capturando. Por lo tanto, necesitamos avanzar hacia un modelo regresión espacial.

Primero seleccionamos el mejor kernel para ajustar el modelo

set.seed(123)
# 1.- Proyectamos a UTM 17S
LisaT_LOGRO_SUP_sp <- st_transform(lisaT_LOGRO_SUP_TRANS, 32717) |>  as("Spatial")

#2.- Seleccionamos bandwidth óptimo con spgwr
bw_gwr_LOGRO_SUP <- gwr.sel(
  T_LOGRO_SUP ~ I_ADULTEZ + T_OCUPACION_SEXO_TRANS + I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS, 
  data = LisaT_LOGRO_SUP_sp,
  coords = st_coordinate(LisaT_LOGRO_SUP_sp),
  adapt = T,
  gweight = gwr.Gauss,  # kernel gaussiano
  verbose = T, 
  show.error.messages = T)
## Adaptive q: 0.381966 CV score: 781072.4 
## Adaptive q: 0.618034 CV score: 809685.2 
## Adaptive q: 0.236068 CV score: 759766.7 
## Adaptive q: 0.145898 CV score: 736115.1 
## Adaptive q: 0.09016994 CV score: 708577.8 
## Adaptive q: 0.05572809 CV score: 679295.3 
## Adaptive q: 0.03444185 CV score: 649203.3 
## Adaptive q: 0.02128624 CV score: 623254.8 
## Adaptive q: 0.01315562 CV score: 594229.4 
## Adaptive q: 0.008130619 CV score: 559313.4 
## Adaptive q: 0.005024999 CV score: 531371.9 
## Adaptive q: 0.00310562 CV score: 507091.4 
## Adaptive q: 0.001919379 CV score: 493459.1 
## Adaptive q: 0.001186241 CV score: 493625.9 
## Adaptive q: 0.001571441 CV score: 491447.1 
## Adaptive q: 0.001530751 CV score: 491261.9 
## Adaptive q: 0.00139916 CV score: 491486.4 
## Adaptive q: 0.001488437 CV score: 491179.1 
## Adaptive q: 0.001447747 CV score: 491234.2 
## Adaptive q: 0.001488437 CV score: 491179.1
#revisamos resultado
bw_gwr_LOGRO_SUP
## [1] 0.001488437

Como resultado obtuvimos el bandwidth óptimo adaptativo (q) de 0.001488437 para mi modelo GWR de educación superior. Primera apreciación al salir un valor muy bajo es que el modelo detecta fuerte heterogeneidad espacial, incluso más que Bachillerato.

Luego ajustamos el modelo GWR con el bandwidth elegido

set.seed(123)
#GWR con el kernel adaptativo.
modelo_gwr_sp_LOGRO_SUP <- gwr(
    T_LOGRO_SUP ~ I_ADULTEZ + T_OCUPACION_SEXO_TRANS + I_ADULTEZ_MAYOR_TRANS + POB_OCUP_TRANS, 
    data = LisaT_LOGRO_SUP_sp,
    coords = st_coordinates(LisaT_LOGRO_SUP_sp),
    adapt = bw_gwr_LOGRO_SUP)

4.7.6.5. Mapear el R2 local para medir explicación del modelo.

Para mapear el R2 local debemos primero transformar el modelo GWR en un ST y extraemos los coeficientes.

#guardamos coeficientes sobre st para graficar.
SDFSUP <- modelo_gwr_sp_LOGRO_SUP$SDF
SDFSUP <- st_as_sf(modelo_gwr_sp_LOGRO_SUP$SDF)

#Graficamos coeficientes
ggplot(SDFSUP) +
  geom_sf(aes(fill = localR2), color = NA) +
  geom_sf(data = borde_general_T_LOGRO_SUP, fill = NA, color = "grey30", size = 0.4) +
  scale_fill_viridis_c(name="R2 local", direction=-1) +
  labs(
    title = "R2 Local del modelo GWR",
    fill = "R2 local",
    subtitle = "Nivel Superior",
    caption = "Fuente: Censo de Población y Vivienda 2022 - INEC") +
  theme_void() +
  theme(
    legend.position = "right",
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10))

El mapa de R2 local muestra con claridad la distribución del modeo GWR donde explica mejor la variación de la variable dependiente Logro Superior indicando: - Los valores de R2 local alto (Azúl oscuro-verde oscuro) explica bien la fuerte variación Positiva entre T_LOGRO_SUP y los predictores en esa zona. Estas zonas mejor explicadas se concentra en la parte oeste de Guayaquil y centro norte. - Los valores de R2 local bajo (Verde claro - Amarillo) explica poco la variación entre T_LOGRO_SUP y los predictores en esa zona. Existiendo otros factores socioeconómicos relevantes.

# Calcular centroides para pasar atributos de una entidad a otra.
centroid_sdfsup <- st_centroid(SDFSUP)

# Join espacial usando centroides
centroid_sdfsup <- st_transform(centroid_sdfsup, 32717) #transformamos al mismo sistema de referencia
lisaT_LOGRO_SUP <- st_transform(lisaT_LOGRO_SUP, 32717) #transformamos al mismo sistema de referencia

#realizamos el spatial join
centroid_sdfsup <- st_join(centroid_sdfsup, lisaT_LOGRO_SUP["quad"])


# Agregar la variable LISA a SDFSUP
SDFSUP$quad <- centroid_sdfsup$quad

# Clasificar R2 local en terciles
breaks_R2 <- classIntervals(SDFSUP$localR2, n = 3, style = "quantile")$brks
SDFSUP$R2_cat <- cut(SDFSUP$localR2,
                     breaks = breaks_R2,
                     include.lowest = TRUE,
                     labels = c("Bajo", "Medio", "Alto"))

#Normalizar los nombres de LISA
#SDFSUP$quad <- factor(SDFSUP$quad,
                  #    levels = c("Low-Low","High-Low","Low-High","High-High","Not significant"))

#Crear variable bivariada combinada
SDFSUP <-  SDFSUP |> drop_na()

SDFSUP$bivar <- paste(SDFSUP$R2_cat, SDFSUP$quad, sep = " | ")

#Paleta de colores
bivar_colors <- c(
  "Bajo | Low-Low"   = "#e8e8e5",
  "Medio | Low-Low"  = "#b5c0da",
  "Alto | Low-Low"   = "#6c83b5",
  
  "Bajo | Low-High"  = "#cbb8d7",
  "Medio | Low-High" = "#9f8bbd",
  "Alto | Low-High"  = "#73579d",
  
  "Bajo | High-Low"  = "#e4acac",
  "Medio | High-Low" = "#c85a5a",
  "Alto | High-Low"  = "#b10101",
  
  "Bajo | High-High" = "#d4a3c6",
  "Medio | High-High"= "#a05a8c",
  "Alto | High-High" = "#7a0177",
  
  "Bajo | No Signif"   = "white",
  "Medio | No Signif"  = "white",
  "Alto | No Signif"   = "white"
)

#Plot final del mapa bivariado
ggplot(SDFSUP) +
  geom_sf(aes(fill = bivar), color = NA) +
  geom_sf(data = borde_general_T_LOGRO_SUP, fill = NA, color = "grey30", size = 0.4) +
  scale_fill_manual(values = bivar_colors) +
  labs(
    title = "Mapa Bivariado: R2 Local del GWR vs LISA Local",
    subtitle = "EDUCACIÓN SUPERIOR: Modelo explicativo + estructura espacial",
    fill = "R2 Local | LISA"
  ) +

  theme_void() +
  theme(
    legend.position = "right",
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10))

Se puede observar que el modelo GWR con valores altos de R2 local puede explicar los clústers High-High que se encuentran con mayor proporción en el territorio y también aquellos con valores Low-Low de LISA. Existen zonas que no son significativas donde el modelo explica poco las cuales están en color blanco.

5.0 Discusión

Aunque se realizaron las transformaciones pertinentes para normalizar variables algunas no se pudieron utilizar para generar una regresión geográficamente ponderada por no tener una distribución normal, las cuales fueron descartadas.

Las relaciones encontradas entre T_LOGRO_BACH y T_LOGRO_SUP son inversamente proporcionales. Es decir, que las variables I_ADULTEZ, T_OCUPACION_SEXO_TRANS, I_ADULTEZ_MAYOR_TRANS, POB_OCUP_TRANS tienen relación positiva con T_LOGRO_SUP y relación negativa con T_LOGRO_BACH. Esta relación a nivel geográfico también se encuentra zonificada para cada variable dependiente. El modelo GWR entrenado para T_LOGRO_BACH explica mejor la región centro y oeste de Guayaquil donde tenemos clústeres Low-Low principalmente.

Por el contrario, El modelo GWR entrenado para T_LOGRO_SUP explica mejor la región norte y oeste de Guayaquil donde tenemos clústeres High-High principalmente. Esto puede deberse pues a las relaciones lineales encontradas y la distribución de los predictores.

Esto nos indica que la relación entre los niveles de educación es inversamente proporcional, así como la explicación del modelo R2 local en cada caso que varía según el fenómeno estudiado y la localización. Esta primera aproximación nos permite conocer la escolaridad y su influencia en las variables predictoras estudiadas y nos confirma que existen clusteres asociados por sectores censales.

Se considera oportuno tratar de añadir variables socioeconómicas con distribución normal para inferir en regiones donde el R2 local fue bajo.

6.0 Referencias mínimas.