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.
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.
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.
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.
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 intervalosEl 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. MujerP17R: 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/DoctoradoP29: 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 ignoraCANTON: 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.
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.
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.
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()
## 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.
Como análisis exploratorio inicial, aplicamos la función summary para conocer la distribuición de las variables de estudio.
## 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
Población Joven (0 a 29 años): La suma de Niños/as, Adolescentes y Jóvenes es significativa. En el grupo de Niños/as (0 a 11 años), hay más hombres (261,015) que mujeres (248,221), lo cual es biológicamente común. En el grupo de Jóvenes (18 a 29 años), hay una ligera mayoría de mujeres (275,887 vs 269,730).
Población Mayor (65 años o más): Hay una diferencia notable, con más mujeres (125,629) que hombres (100,756). Esto refleja la tendencia general de las mujeres a tener una mayor esperanza de vida en Guayaquil.
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
Principal Categoría (NA): La barra dominante, con 672,970 personas, está marcada como NA (No Aplica / No Especifica). En el contexto de los censos, esto a menudo incluye a personas fuera de la Población Económicamente Activa (PEA), como estudiantes, jubilados, personas dedicadas al trabajo doméstico no remunerado, o personas que simplemente no declararon su ocupación.
Empleado/a privado/a: Es la categoría de empleo formal más grande, con 301,799 hombres.
Cuenta propia: La segunda categoría más grande, con 182,823 hombres.
Ocupaciones Minoritarias: Las categorías de Jornalero/a (14,179), Socio/a (2,557), Trabajador/a no remunerado/a (4,364) y Empleador/a del Estado (42,237) son considerablemente menores.
Mujeres
Principal Categoría (NA): Al igual que en los hombres, esta es la categoría más grande, pero con un número significativamente mayor: 940,449 mujeres. Esta diferencia acentuada probablemente incluye a la gran mayoría de mujeres dedicadas a las labores domésticas y de cuidado no remuneradas, que son demográficamente más comunes entre las mujeres.
Empleado/a privado/a: Es la categoría de empleo formal más grande, con 177,937 mujeres. Es notablemente inferior al número de hombres en esta misma categoría.
Cuenta propia: La segunda categoría, con 114,459 mujeres, también inferior a la de los hombres.
Ocupaciones Minoritarias: Similarmente, las otras categorías de empleo tienen frecuencias bajas.
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).
Cargamos el dataset espacial con los sectores censales del Ecuador
## 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.
## 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 joinTenemos 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 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.
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:
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.
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.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.
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.
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.
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.
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:
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:
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.
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:
El resultado de cada índice se lo utiliza para estudiar el comportamiento de la asociación espacial local.
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.
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:
Aplicando este algoritmo, pudimos eliminar valores NA. Pasamos de tener un total de 6357 observaciones a 6355.
El resultado del test de Anderson–Darling nos da las siguientes hipótesis:
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_OK$I_FEMINIDAD
## A = 24.49, p-value < 0.00000000000000022
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_OK$T_TRABAJO_NO_REM
## A = 391.11, p-value < 0.00000000000000022
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_OK$I_JUVENTUD
## A = 2.933, p-value = 0.0000002283
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_OK$T_OCUPACION_SEXO
## A = 67.613, p-value < 0.00000000000000022
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_OK$I_ADULTEZ
## A = 0.7001, p-value = 0.06758
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_OK$I_ADULTEZ_MAYOR
## A = 76.079, p-value < 0.00000000000000022
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_OK$POB_OCUP
## A = 30.228, p-value < 0.00000000000000022
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_OK$I_MASCULINIDAD
## A = 526.83, p-value < 0.00000000000000022
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_OK$T_ANALFABETISMO
## A = 36.752, p-value < 0.00000000000000022
##
## 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"))## $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:
Para confirmar la estandarización realizada para aproximar los datos a una distribución normal, aplicamos Anderson-Darling nuevamente sin considerar I_Adultez.
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_TRANS$I_FEMINIDAD_TRANS
## A = 0.013287, p-value = 1
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_TRANS$T_TRABAJO_NO_REM_TRANS
## A = 209.49, p-value < 0.00000000000000022
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_TRANS$I_JUVENTUD_TRANS
## A = 2.933, p-value = 0.0000002283
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_TRANS$T_OCUPACION_SEXO_TRANS
## A = 0.015773, p-value = 1
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_TRANS$I_ADULTEZ_MAYOR_TRANS
## A = 0.00081602, p-value = 1
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_TRANS$POB_OCUP_TRANS
## A = 0.19107, p-value = 0.8976
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_TRANS$I_MASCULINIDAD_TRANS
## A = 0.013287, p-value = 1
##
## Anderson-Darling normality test
##
## data: lisaT_LOGRO_BACH_TRANS$T_ANALFABETISMO_TRANS
## A = 0.97043, p-value = 0.01456
##
## 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.
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:
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.
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
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.
## 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:
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
## [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
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:
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 variableComo 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:
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.
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
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.
## 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:
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
## [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
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.
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.