Evaluación de impacto para la atención de LLA en la población infantil de la CDMX (1998-2024)
JHONATAN GAONA
2026-03-09
# Configuración para evitar problemas
knitr::opts_chunk$set(
echo = TRUE,
warning = FALSE,
message = FALSE,
cache = FALSE # ← Importante: deshabilita caché
)BLOQUE 0:DATASETS
- Encuesta de defunciones registradas (EDR)
Se importa desde la subcarpeta localizada en el ordenador la Encuesta de Defunciones Registradas publicadas al 2024 por parte del INEGI
library(plyr)
library(foreign)
# Read files directly without changing working directory
file_path <- "/Users/jhonatangaona/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/EDR/DEFUN"
# Get list of DBF files
dbf_files <- list.files(path = file_path,
pattern = "\\.dbf$",
full.names = TRUE)
# Read all DBF files
dataset <- ldply(dbf_files, read.dbf)Los atributos que permanecen a lo largo de los ciclos publicados, con su correspondiente modificación en la codificación son:
- Escolaridad: 1990-2003, 2004-2011, 2012-2023
- Lugar de ocurrencia: 1990-2002,2003, 2004-2011, 2012-2023,
- Derechohabiente: 1990-2003, 2004-2011, 2012-2017, 2018-2021, 2022-2023
- sitio de ocurrencia: 1990-2003, 2004-2023
- causa de defunción: 1990-1997 (2040), 1998-2023 (c910)
colnames(dataset)## [1] "ENT_REGIS" "MUN_REGIS" "ENT_RESID" "MUN_RESID" "TLOC_RESID"
## [6] "ENT_OCURR" "MUN_OCURR" "CAUSA_DEF" "LISTA_MEX" "SEXO"
## [11] "EDAD" "DIA_OCURR" "MES_OCURR" "ANIO_OCUR" "MES_REGIS"
## [16] "ANIO_REGIS" "DIA_NACIM" "MES_NACIM" "ANIO_NACIM" "OCUPACION"
## [21] "ESCOLARIDA" "EDO_CIVIL" "PRESUNTO" "OCURR_TRAB" "LUGAR_OCUR"
## [26] "NECROPSIA" "ASIST_MEDI" "SITIO_OCUR" "COND_CERT" "CERT_NOMED"
## [31] "NACIONALID" "DERECHOHAB" "EMBARAZO" "REL_EMBA" "HORAS"
## [36] "MINUTOS" "CAPITULO" "GRUPO" "LISTA1" "GR_LISMEX"
## [41] "VIO_FAMI" "EDAD_AGRU" "MATERNAS" "DIS_RE_OAX" "LOC_RESID"
## [46] "TLOC_OCURR" "LOC_OCURR" "DIA_REGIS" "AREA_UR" "COMPLICARO"
## [51] "DIA_CERT" "MES_CERT" "ANIO_CERT" "PESO" "LENGUA"
## [56] "COND_ACT" "PAR_AGRE" "ENT_OCULES" "MUN_OCULES" "LOC_OCULES"
## [61] "RAZON_M" "TLOC_REGIS" "LOC_REGIS" "COD_ADICIO" "ENT_NAC"
## [66] "AFROMEX" "CONINDIG" "CVE_LENGUA" "NACESP_CVE" "SEM_GEST"
## [71] "GRAMOS" "TIPO_DEFUN" "CIRUGIA" "NATVIOLE" "USONECROPS"
## [76] "ENCEFALICA" "DONADOR" "ID_INEGI" "LISTA_BAS"Seleccion de atributos
# seleccion de atributos por posición
df<-dataset[,c(3,4,5,8,10,12,13,14,21,27,28,32,42,49)]
colnames(df)## [1] "ENT_RESID" "MUN_RESID" "TLOC_RESID" "CAUSA_DEF" "SEXO"
## [6] "DIA_OCURR" "MES_OCURR" "ANIO_OCUR" "ESCOLARIDA" "ASIST_MEDI"
## [11] "SITIO_OCUR" "DERECHOHAB" "EDAD_AGRU" "AREA_UR"# changue data type
str(df)## 'data.frame': 20402630 obs. of 14 variables:
## $ ENT_RESID : int 1 1 1 1 1 1 1 1 1 1 ...
## $ MUN_RESID : int 1 1 1 1 1 1 1 1 1 1 ...
## $ TLOC_RESID: int 14 14 14 14 14 14 14 14 14 14 ...
## $ CAUSA_DEF : Factor w/ 11209 levels "A010","A014",..: 73 2243 3346 1186 797 1510 1557 2409 401 279 ...
## $ SEXO : int 2 2 2 2 2 1 1 1 2 1 ...
## $ DIA_OCURR : int 5 99 1 1 2 29 1 1 19 11 ...
## $ MES_OCURR : int 1 1 1 1 1 1 1 1 1 1 ...
## $ ANIO_OCUR : int 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 ...
## $ ESCOLARIDA: int 3 9 1 9 4 7 4 9 6 2 ...
## $ ASIST_MEDI: int 1 1 2 1 1 1 2 1 1 1 ...
## $ SITIO_OCUR: int 3 1 4 3 1 2 3 1 3 3 ...
## $ DERECHOHAB: int 1 2 1 99 2 3 2 2 1 1 ...
## $ EDAD_AGRU : Factor w/ 30 levels "01","02","03",..: 23 1 30 30 20 15 20 1 21 24 ...
## $ AREA_UR : int NA NA NA NA NA NA NA NA NA NA ...- Catálogo Único de Claves de Áreas Geoestadísticas Estatales, Municipales y Localidades.
Importamos un catalogo de entidades federativas proporcionado por el INEGI para poder realizar una relación de códigos con sus respectivos nombres.
library(readxl)
# catalogo de claves, entidades y abreviaturas
AGEEML_INEGI <- read_excel("~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/Bases/INEGI/AGEEML_INEGI.xlsx",
skip = 3)
# transformar la clave de entidad
AGEEML_INEGI$CVE_ENT<-as.numeric(AGEEML_INEGI$CVE_ENT)
# seleccion por posicion de columnas
AGEEML_INEGI<-AGEEML_INEGI[,c(2,4)]
# eliminar duplicados en todo el conjunto
AGEEML_INEGI<-AGEEML_INEGI[!duplicated(AGEEML_INEGI), ]
# renombrar columnas
names(AGEEML_INEGI)[1] <- "Entidad"
head(AGEEML_INEGI)## # A tibble: 6 × 2
## Entidad NOM_ABR
## <dbl> <chr>
## 1 1 Ags.
## 2 2 BC
## 3 3 BCS
## 4 4 Camp.
## 5 5 Coah.
## 6 6 Col.- Catálogo Único de Claves de Áreas Geoestadísticas Estatales, Municipales y Localidades.
Importamos un catalogo de municipios proporcionado por el INEGI para poder realizar una relación de códigos con sus respectivos nombres.
# catalogo de claves, entidades, municipios y abreviaturas
AGEEMLM_INEGI <- read_excel("~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/Bases/INEGI/AGEEML_M_INEGI.xlsx",
skip = 3)
# cambio del tipo de datos
AGEEMLM_INEGI$CVE_ENT<-as.numeric(AGEEMLM_INEGI$CVE_ENT)
AGEEMLM_INEGI$CVE_MUN<-as.numeric(AGEEMLM_INEGI$CVE_MUN)
# solo CDMX
AGEEMLM_INEGI<-AGEEMLM_INEGI[AGEEMLM_INEGI$CVE_ENT==9,]
# seleccion de columnas
AGEEMLM_INEGI<-AGEEMLM_INEGI[,c(5,6)]
# renombrar columnas
names(AGEEMLM_INEGI)[1] <- "MUN_RESID"
head(AGEEMLM_INEGI)## # A tibble: 6 × 2
## MUN_RESID NOM_MUN
## <dbl> <chr>
## 1 2 Azcapotzalco
## 2 3 Coyoacán
## 3 4 Cuajimalpa de Morelos
## 4 5 Gustavo A. Madero
## 5 6 Iztacalco
## 6 7 Iztapalapa- Transparencia presupuestaria: Fideicomisos al 4T del 2025.
library(data.table)
# Leer sin mensajes
fideicomiso <- fread("~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/Bases/SHCP/fideicomisos_4t2025.csv")
# columnas de interés
fideicomiso<- fideicomiso[,.(CICLO, TRIMESTRE,ID_RAMO,DESC_RAMO,
ID_UR,DESC_GRUPO_TEMATICO,
MONTO_INGRESOS,
MONTO_EGRESOS,MONTO_DISPONIBILIDAD
)]
# filtrado de atributos
fideicomiso<-fideicomiso[ID_UR== "M7B"]
# Convertir a numérico y en miles de millones
fideicomiso[, `:=`(
MONTO_INGRESOS = round(as.numeric(MONTO_INGRESOS) / 1e9,0),
MONTO_EGRESOS = round(as.numeric(MONTO_EGRESOS)/ 1e9,0),
MONTO_DISPONIBILIDAD = round(as.numeric(MONTO_DISPONIBILIDAD) / 1e9,0)
)]
# Ver estructura
head(fideicomiso)## CICLO TRIMESTRE ID_RAMO DESC_RAMO ID_UR DESC_GRUPO_TEMATICO MONTO_INGRESOS
## <int> <int> <int> <char> <char> <char> <num>
## 1: 2005 1 12 SALUD M7B APOYOS FINANCIEROS 1
## 2: 2005 2 12 SALUD M7B APOYOS FINANCIEROS 2
## 3: 2005 3 12 SALUD M7B APOYOS FINANCIEROS 3
## 4: 2005 4 12 SALUD M7B APOYOS FINANCIEROS 3
## 5: 2006 1 12 SALUD M7B APOYOS FINANCIEROS 1
## 6: 2006 2 12 SALUD M7B APOYOS FINANCIEROS 2
## MONTO_EGRESOS MONTO_DISPONIBILIDAD
## <num> <num>
## 1: 0 2
## 2: 0 2
## 3: 1 3
## 4: 1 3
## 5: 0 4
## 6: 1 5- INPC
# DEFINIR la ruta del archivo primero
archivo_excel <- "~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/Bases/INEGI/ca61_2018a.xlsx"
# Ahora SÍ puedes leer
inpc_cdmx <- read_excel(archivo_excel,
sheet = "cdmx")
setDT(inpc_cdmx)
inpc_cdmx## fecha_mes INPC_CDMX mes año
## <char> <num> <char> <num>
## 1: Ene 1998 33.36190 Ene 1998
## 2: Feb 1998 33.94154 Feb 1998
## 3: Mar 1998 34.33677 Mar 1998
## 4: Abr 1998 34.66434 Abr 1998
## 5: May 1998 34.97156 May 1998
## ---
## 320: Ago 2024 133.13900 Ago 2024
## 321: Sep 2024 133.14300 Sep 2024
## 322: Oct 2024 133.51800 Oct 2024
## 323: Nov 2024 133.51100 Nov 2024
## 324: Dic 2024 134.10700 Dic 2024- PIBE-CDMX
# Ruta del archivo Excel
archivo_excel <- "~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/Bases/INEGI/pibe_r_17.xlsx"
pibe_cdmx<- read_excel(archivo_excel, sheet = "Tabulado", skip = 4)
setDT(pibe_cdmx)
# PASO 2: Identificar y filtrar la fila de PIB a precios de mercado
# Buscar la fila que contiene el texto específico
fila_pib <- pibe_cdmx[grepl("Producto interno bruto.*precios de mercado",
Concepto, ignore.case = TRUE)]
if(nrow(fila_pib) > 0) {
# Si hay múltiples coincidencias, tomar la primera
fila_pib <- fila_pib[1]
cat(" ✓ Fila encontrada:", fila_pib$Concepto, "\n\n")
} else {
stop("no se encontró la fila de 'Producto interno bruto, a precios de mercado'")
}## ✓ Fila encontrada: ____aProducto interno bruto, a precios de mercado# PASO 3: Seleccionar solo columnas de años (1998-2024)
# Identificar columnas que son años
cols_años <- names(pibe_cdmx)[grepl("^19[0-9]{2}$|^20[0-9]{2}$", names(pibe_cdmx))]
# Filtrar solo años entre 1998 y 2024
años_validos <- as.numeric(cols_años)
cols_años_filtrados <- cols_años[años_validos >= 1998 & años_validos <= 2024]
# Seleccionar solo esas columnas
pib_ancho <- fila_pib[, c("Concepto", cols_años_filtrados), with = FALSE]
# PASO 4: Convertir a formato largo
pib_nominal_cdmx <- melt(
pib_ancho,
id.vars = "Concepto",
variable.name = "año",
value.name = "PIB_nominal"
)
# Eliminar la columna Concepto (ya no la necesitamos)
pib_nominal_cdmx[, Concepto := NULL]
# PASO 5: Convertir año a numérico y PIB_nominal a numérico
pib_nominal_cdmx[, `:=`(
año = as.integer(as.character(año)),
PIB_nominal = as.numeric(PIB_nominal)
)]
# PASO 6: Ordenar por año
setorder(pib_nominal_cdmx, año)
pib_nominal_cdmx## año PIB_nominal
## <int> <num>
## 1: 1998 2657943
## 2: 1999 2739399
## 3: 2000 2903709
## 4: 2001 2870934
## 5: 2002 2893524
## 6: 2003 2842789
## 7: 2004 2962342
## 8: 2005 2991075
## 9: 2006 3134225
## 10: 2007 3142464
## 11: 2008 3183429
## 12: 2009 3008812
## 13: 2010 3125717
## 14: 2011 3231440
## 15: 2012 3332209
## 16: 2013 3363394
## 17: 2014 3412110
## 18: 2015 3503520
## 19: 2016 3571368
## 20: 2017 3646318
## 21: 2018 3694575
## 22: 2019 3679917
## 23: 2020 3293855
## 24: 2021 3496056
## 25: 2022 3650378
## 26: 2023 3805614
## 27: 2024 3885957
## año PIB_nominal
## <int> <num>7.PIBE + INPC
Con el objetio de transformar los precios nominales/corrientes a precios constantes/reales, se hace un merge de del PIB nominal de la CDMX con el INPC a nivel CDMX.
# INPC a nivel año
inpc_anual <- inpc_cdmx[, .(
INPC_CDMX_promedio = mean(INPC_CDMX, na.rm = TRUE)), by = año]
# union de datasets
economia_cdmx <- merge(pib_nominal_cdmx,
inpc_anual,
by = "año",
all.x = TRUE)
# PIB real calculado
economia_cdmx[,PIB_real := (PIB_nominal/INPC_CDMX_promedio)*100]
economia_cdmx## Key: <año>
## año PIB_nominal INPC_CDMX_promedio PIB_real
## <int> <num> <num> <num>
## 1: 1998 2657943 35.71249 7442614
## 2: 1999 2739399 41.59066 6586573
## 3: 2000 2903709 45.36122 6401302
## 4: 2001 2870934 47.93227 5989564
## 5: 2002 2893524 50.44800 5735657
## 6: 2003 2842789 52.81829 5382205
## 7: 2004 2962342 55.32851 5354098
## 8: 2005 2991075 57.45369 5206063
## 9: 2006 3134225 59.51390 5266375
## 10: 2007 3142464 61.91556 5075402
## 11: 2008 3183429 65.22721 4880523
## 12: 2009 3008812 68.66254 4382028
## 13: 2010 3125717 71.74483 4356713
## 14: 2011 3231440 74.42420 4341921
## 15: 2012 3332209 77.48259 4300590
## 16: 2013 3363394 80.76191 4164580
## 17: 2014 3412110 84.35400 4044989
## 18: 2015 3503520 86.98035 4027944
## 19: 2016 3571368 89.54426 3988383
## 20: 2017 3646318 94.89375 3842527
## 21: 2018 3694575 99.90812 3697973
## 22: 2019 3679917 103.79942 3545220
## 23: 2020 3293855 107.31267 3069400
## 24: 2021 3496056 112.67458 3102790
## 25: 2022 3650378 120.33992 3033389
## 26: 2023 3805614 126.43025 3010051
## 27: 2024 3885957 132.14383 2940702
## año PIB_nominal INPC_CDMX_promedio PIB_real
## <int> <num> <num> <num>8.PIB per-capita
Calculo de edad promedio de vida para la CDMX desde 1998-2025
archivo_excel <- "~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/Bases/CONAPO/6_Esperanza_Vida_Nacer_1950_2070.xlsx"
ev_nacional <- read_excel(archivo_excel,
sheet = "EV_completa")
setDT(ev_nacional)
# Filtrado
df_cdmx <- ev_nacional[CVE_GEO == 9]
df_cdmx <- ev_nacional[AÑO >= 1998 & AÑO <= 2025]
df_cdmx <- ev_nacional[SEXO %in% c("Hombres", "Mujeres")]
# Calculos agregados
ev_cdmx <- df_cdmx[, .(esperanza_vida = mean(EV, na.rm = TRUE)), by = AÑO]
# Renombrar AÑO a año
setnames(ev_cdmx, "AÑO", "año")
# Ordenar por año
setorder(ev_cdmx, año)
ev_cdmx## año esperanza_vida
## <num> <num>
## 1: 1950 46.30306
## 2: 1951 45.72893
## 3: 1952 48.09234
## 4: 1953 47.74312
## 5: 1954 51.47421
## ---
## 117: 2066 82.53377
## 118: 2067 82.68164
## 119: 2068 82.82903
## 120: 2069 82.97596
## 121: 2070 83.122459.Población a mitad de año
Se construye la tabla de contingencia a partir de las estimación poblacional a mitad de año para la CDMX publicado por la CONAPO en 2024.
PMA = (Población Inicial + Nacimientos - Muertes + Migración Neta)
archivo_csv <- "~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/Bases/CONAPO/Pob_Mitad_1950_2070.csv"
# Leer archivo CSV
pob <- fread(archivo_csv)
# Seleccionar columnas y filtrar
pob_total_cdmx <- pob[CVE_GEO == 9 &
AÑO >= 1998 & AÑO <= 2025 &
EDAD <=150,
.(poblacion_cdmx = sum(POBLACION, na.rm = TRUE)),
by = AÑO]
# Renombrar AÑO a año
setnames(pob_total_cdmx, "AÑO", "año")
# Ordenar por año
setorder(pob_total_cdmx, año)
pob_total_cdmx## año poblacion_cdmx
## <int> <int>
## 1: 1998 8674660
## 2: 1999 8682629
## 3: 2000 8706087
## 4: 2001 8758667
## 5: 2002 8819709
## 6: 2003 8878329
## 7: 2004 8931633
## 8: 2005 8977190
## 9: 2006 9001741
## 10: 2007 9011907
## 11: 2008 9022695
## 12: 2009 9029228
## 13: 2010 9044151
## 14: 2011 9074286
## 15: 2012 9102171
## 16: 2013 9121489
## 17: 2014 9131261
## 18: 2015 9152227
## 19: 2016 9196836
## 20: 2017 9246758
## 21: 2018 9296268
## 22: 2019 9343071
## 23: 2020 9332593
## 24: 2021 9272649
## 25: 2022 9237644
## 26: 2023 9221637
## 27: 2024 9204018
## 28: 2025 9183658
## año poblacion_cdmx
## <int> <int>- Macro indicadores CDMX Procederemos a realizar un dataframe donde en su conjunto contenga los principales macro indicadores calculados previamente
lista_datasets <- list(economia_cdmx, ev_cdmx, pob_total_cdmx)
# Unir todos por 'año'
macros_cdmx <- Reduce(function(x, y) merge(x, y,
by = "año",
all = FALSE), lista_datasets)
# PIB_percapita calculado
macros_cdmx[,PIB_per_capita := (PIB_real*1000000/poblacion_cdmx)]
macros_cdmx## Key: <año>
## año PIB_nominal INPC_CDMX_promedio PIB_real esperanza_vida poblacion_cdmx
## <int> <num> <num> <num> <num> <int>
## 1: 1998 2657943 35.71249 7442614 73.08981 8674660
## 2: 1999 2739399 41.59066 6586573 73.48701 8682629
## 3: 2000 2903709 45.36122 6401302 74.06916 8706087
## 4: 2001 2870934 47.93227 5989564 74.38121 8758667
## 5: 2002 2893524 50.44800 5735657 74.42981 8819709
## 6: 2003 2842789 52.81829 5382205 74.49471 8878329
## 7: 2004 2962342 55.32851 5354098 74.86737 8931633
## 8: 2005 2991075 57.45369 5206063 74.76224 8977190
## 9: 2006 3134225 59.51390 5266375 75.14322 9001741
## 10: 2007 3142464 61.91556 5075402 75.05439 9011907
## 11: 2008 3183429 65.22721 4880523 74.83693 9022695
## 12: 2009 3008812 68.66254 4382028 74.63739 9029228
## 13: 2010 3125717 71.74483 4356713 74.39732 9044151
## 14: 2011 3231440 74.42420 4341921 74.85093 9074286
## 15: 2012 3332209 77.48259 4300590 75.06629 9102171
## 16: 2013 3363394 80.76191 4164580 75.29811 9121489
## 17: 2014 3412110 84.35400 4044989 75.20473 9131261
## 18: 2015 3503520 86.98035 4027944 75.16729 9152227
## 19: 2016 3571368 89.54426 3988383 74.91849 9196836
## 20: 2017 3646318 94.89375 3842527 74.91544 9246758
## 21: 2018 3694575 99.90812 3697973 74.98876 9296268
## 22: 2019 3679917 103.79942 3545220 74.93456 9343071
## 23: 2020 3293855 107.31267 3069400 69.51042 9332593
## 24: 2021 3496056 112.67458 3102790 69.25821 9272649
## 25: 2022 3650378 120.33992 3033389 75.06909 9237644
## 26: 2023 3805614 126.43025 3010051 75.23577 9221637
## 27: 2024 3885957 132.14383 2940702 75.43892 9204018
## año PIB_nominal INPC_CDMX_promedio PIB_real esperanza_vida poblacion_cdmx
## <int> <num> <num> <num> <num> <int>
## PIB_per_capita
## <num>
## 1: 857971.9
## 2: 758592.1
## 3: 735267.4
## 4: 683844.3
## 5: 650322.7
## 6: 606218.2
## 7: 599453.4
## 8: 579921.2
## 9: 585039.6
## 10: 563188.5
## 11: 540916.4
## 12: 485315.9
## 13: 481716.1
## 14: 478486.3
## 15: 472479.6
## 16: 456568.0
## 17: 442982.6
## 18: 440105.4
## 19: 433669.1
## 20: 415554.0
## 21: 397791.1
## 22: 379449.1
## 23: 328890.4
## 24: 334617.4
## 25: 328372.6
## 26: 326411.7
## 27: 319502.0
## PIB_per_capita
## <num>BLOQUE 1: QAD
Con el objeto de homologar términos y atributos, se toma la última actualización del diccionario de dato de la EDR-24 publicado por el INEGI. Primeramente, las valores nominales se les asigna con caracteres de letra (genero,); segundo, algunas clasificaciones dentro de los atributos se han agrupado bajo términos de mayor amplitud (sitio de ocurrencia, tipo de seguro medico, escolaridad); y tercero, se traslada la codificación de valores a caracteres alfanuméricos (densidad poblacional,asistencia médica,edad agrupada, área urbana)
Calidad de los datoss
library(data.table)
library(validate)
library(naniar)library(data.table)
library(validate)
# ═══════════════════════════════════════════════════════════════════
# FUNCIONES AUXILIARES - VERSIÓN OPTIMIZADA PARA MEMORIA
# ═══════════════════════════════════════════════════════════════════
print_section <- function(title, char = "-", width = 70) {
cat(rep(char, width), "\n", sep = "")
cat(title, "\n")
cat(rep(char, width), "\n\n", sep = "")
}
# Función OPTIMIZADA para calcular completitud (evita copias de data)
calcular_completitud <- function(dataset, vars_criticas, vars_importantes, vars_opcionales) {
cat("Calculando completitud...\n")
# Calcular estadísticas sin crear matrices grandes
n_total <- nrow(dataset)
vars <- names(dataset)
# Procesar una variable a la vez para no sobrecargar memoria
completitud <- data.table(Variable = vars)
completitud[, Tipo := sapply(vars, function(v) class(dataset[[v]])[1])]
completitud[, Total := n_total]
completitud[, NAs := sapply(vars, function(v) sum(is.na(dataset[[v]])))]
completitud[, Completos := Total - NAs]
completitud[, Pct_Completo := round((Completos / Total) * 100, 2)]
setorder(completitud, -NAs)
completitud[, Criticidad := fcase(
Variable %in% vars_criticas, "CRITICA",
Variable %in% vars_importantes, "IMPORTANTE",
Variable %in% vars_opcionales, "OPCIONAL",
default = "OTRA"
)]
gc() # Limpiar memoria
list(
tabla = completitud,
score_general = mean(completitud$Pct_Completo),
score_criticas = completitud[Variable %in% vars_criticas, mean(Pct_Completo)],
score_importantes = completitud[Variable %in% vars_importantes, mean(Pct_Completo)]
)
}
# Función OPTIMIZADA para validaciones (procesa en lotes)
ejecutar_validaciones <- function(dataset, reglas, nombre_dimension) {
resultado <- tryCatch({
cat("Ejecutando validaciones de", nombre_dimension, "...\n")
# Configurar validate para usar menos memoria
conf_result <- confront(dataset, reglas, raise = "all")
resumen <- as.data.table(summary(conf_result))
score <- mean(resumen$passes / resumen$items, na.rm = TRUE) * 100
cat("RESUMEN DE", toupper(nombre_dimension), ":\n")
cat(" Score:", sprintf("%.2f%%", score), "\n")
cat(" Reglas pasadas:", sum(resumen$passes == resumen$items, na.rm = TRUE),
"de", nrow(resumen), "\n\n")
rm(conf_result) # Liberar memoria
gc()
list(score = score, resumen = resumen)
}, error = function(e) {
cat("Error en", nombre_dimension, ":", e$message, "\n\n")
gc()
list(score = 0, resumen = NULL)
})
return(resultado)
}
# Función OPTIMIZADA para accuracy (evita copias grandes)
analizar_accuracy <- function(dataset) {
resultado <- tryCatch({
resumen <- data.table(
Metrica = character(),
Categoria = character(),
Valor = character(),
Detalle = character()
)
# ===== EDAD - Extraer solo lo necesario =====
cat("Analizando edad...\n")
# Solo extraer edades válidas (no crear vector completo)
edad_anos <- dataset[EDAD >= 4001 & EDAD <= 4998, EDAD - 4000]
if(length(edad_anos) > 0) {
# Calcular estadísticas directamente
media_edad <- mean(edad_anos, na.rm = TRUE)
mediana_edad <- median(edad_anos, na.rm = TRUE)
q1_edad <- quantile(edad_anos, 0.25, na.rm = TRUE)
q3_edad <- quantile(edad_anos, 0.75, na.rm = TRUE)
iqr_edad <- q3_edad - q1_edad
# Outliers
limite_inf <- q1_edad - 1.5 * iqr_edad
limite_sup <- q3_edad + 1.5 * iqr_edad
outliers <- sum(edad_anos < limite_inf | edad_anos > limite_sup, na.rm = TRUE)
pct_outliers <- round(outliers / length(edad_anos) * 100, 2)
cat(sprintf(" Media: %.1f | Mediana: %.0f | Outliers: %.2f%%\n\n",
media_edad, mediana_edad, pct_outliers))
score_edad <- max(0, 100 - pct_outliers)
# Agregar a resumen
resumen <- rbindlist(list(
resumen,
data.table(
Metrica = "Edad",
Categoria = c("Media", "Mediana", "Q1-Q3", "Outliers", "Score"),
Valor = as.character(c(
round(media_edad, 1),
mediana_edad,
paste0(q1_edad, " - ", q3_edad),
format(outliers, big.mark = ","),
sprintf("%.2f%%", score_edad)
)),
Detalle = c(
"Promedio de edad",
"Mediana de edad",
"Rango intercuartilico",
paste0(pct_outliers, "% del total"),
"100 - % outliers"
)
)
))
rm(edad_anos) # Liberar memoria
gc()
} else {
score_edad <- 100
}
# ===== SEXO - Agregación directa =====
cat("Analizando sexo...\n")
dist_sexo <- dataset[SEXO %in% c(1, 2), .N, by = SEXO]
if(nrow(dist_sexo) > 0) {
dist_sexo[, `:=`(
Pct = round(N / sum(N) * 100, 2),
Sexo_Desc = ifelse(SEXO == 1, "Hombre", "Mujer")
)]
cat("DISTRIBUCION DE SEXO:\n")
print(dist_sexo[, .(Sexo_Desc, N, Pct)])
pct_masc <- dist_sexo[SEXO == 1, Pct]
if(length(pct_masc) > 0) {
desviacion <- abs(pct_masc - 52)
cat(sprintf("\n Desviacion: %.2f%%\n\n", desviacion))
score_sexo <- max(0, 100 - (desviacion * 5))
# Agregar a resumen
resumen <- rbindlist(list(
resumen,
dist_sexo[, .(
Metrica = "Sexo",
Categoria = Sexo_Desc,
Valor = format(N, big.mark = ","),
Detalle = paste0(Pct, "%")
)],
data.table(
Metrica = "Sexo",
Categoria = c("Desviacion", "Score"),
Valor = as.character(c(sprintf("%.2f%%", desviacion), sprintf("%.2f%%", score_sexo))),
Detalle = c("vs. 52% esperado", "100 - (desviacion x 5)")
)
))
} else {
score_sexo <- 100
}
} else {
score_sexo <- 100
}
# Score total
score_total <- mean(c(score_edad, score_sexo), na.rm = TRUE)
resumen <- rbindlist(list(
resumen,
data.table(
Metrica = "TOTAL",
Categoria = "Score Accuracy",
Valor = sprintf("%.2f%%", score_total),
Detalle = "Promedio edad y sexo"
)
))
cat("SCORE DE ACCURACY:", sprintf("%.2f%%", score_total), "\n\n")
gc() # Limpiar memoria
list(score = score_total, resumen = resumen)
}, error = function(e) {
cat("Error en accuracy:", e$message, "\n\n")
gc()
list(
score = 0,
resumen = data.table(
Metrica = "ERROR",
Categoria = "Error",
Valor = "0%",
Detalle = e$message
)
)
})
return(resultado)
}
# ═══════════════════════════════════════════════════════════════════
# FUNCIÓN PRINCIPAL - OPTIMIZADA
# ═══════════════════════════════════════════════════════════════════
evaluar_calidad_mortalidad_oficial <- function(dataset, ruta_salida = NULL) {
# Asegurar que es data.table (sin copiar si ya lo es)
if(!is.data.table(dataset)) {
setDT(dataset)
}
cat("\n")
print_section("EVALUACION DE CALIDAD - DEFUNCIONES MEXICO", "=")
cat("Dataset:", format(nrow(dataset), big.mark = ","), "registros\n")
cat("Memoria en uso:", format(object.size(dataset), units = "Gb"), "\n\n")
# ═══════════════════════════════════════════════════════════════════
# 1. COMPLETENESS
# ═══════════════════════════════════════════════════════════════════
print_section("DIMENSION 1: COMPLETENESS")
vars_criticas <- c("ENT_REGIS", "MUN_REGIS", "ENT_RESID", "MUN_RESID",
"ENT_OCURR", "MUN_OCURR", "ANIO_OCUR", "MES_OCURR",
"SEXO", "EDAD", "CAUSA_DEF")
vars_importantes <- c("DIA_OCURR", "ESCOLARIDA", "EDO_CIVIL", "ASIST_MEDI",
"SITIO_OCUR", "DERECHOHAB", "NACIONALID")
vars_opcionales <- c("OCUPACION", "LENGUA", "AFROMEX", "CONINDIG",
"NECROPSIA", "PRESUNTO", "OCURR_TRAB")
result_completitud <- calcular_completitud(dataset, vars_criticas, vars_importantes, vars_opcionales)
cat("RESUMEN:\n")
cat(sprintf(" General: %.2f%% | Criticas: %.2f%% | Importantes: %.2f%%\n",
result_completitud$score_general,
result_completitud$score_criticas,
result_completitud$score_importantes))
cat(" 100% completas:", sum(result_completitud$tabla$Pct_Completo == 100),
"de", nrow(result_completitud$tabla), "\n\n")
# ═══════════════════════════════════════════════════════════════════
# 2. INTEGRITY - REDUCIDO (solo reglas críticas)
# ═══════════════════════════════════════════════════════════════════
print_section("DIMENSION 2: INTEGRITY")
# REDUCIR REGLAS para ahorrar memoria
reglas_integridad <- validator(
# Solo geografía crítica
ent_regis_rango = ENT_REGIS >= 1 & ENT_REGIS <= 32,
ent_resid_rango = ENT_RESID %in% c(1:32, 35, 99),
mun_regis_rango = MUN_REGIS >= 1 & MUN_REGIS <= 570,
# Solo fechas críticas
mes_ocurr_rango = MES_OCURR %in% c(1:12, 99) | is.na(MES_OCURR),
anio_ocur_rango = (ANIO_OCUR >= 1900 & ANIO_OCUR <= 2024) | ANIO_OCUR == 9999,
# Demográficas
sexo_rango = SEXO %in% c(1, 2, 9),
edad_formato = (EDAD >= 1001 & EDAD <= 1098) | (EDAD >= 2001 & EDAD <= 2098) |
(EDAD >= 3001 & EDAD <= 3098) | (EDAD >= 4001 & EDAD <= 4998) | is.na(EDAD)
)
result_integridad <- ejecutar_validaciones(dataset, reglas_integridad, "integridad")
# ═══════════════════════════════════════════════════════════════════
# 3. CONSISTENCY - REDUCIDO
# ═══════════════════════════════════════════════════════════════════
print_section("DIMENSION 3: CONSISTENCY")
# Solo reglas más importantes
reglas_consistencia <- validator(
nacim_antes_ocurr = ANIO_NACIM <= ANIO_OCUR | ANIO_NACIM == 9999 | is.na(ANIO_NACIM),
febrero_29dias = !(MES_OCURR == 2 & DIA_OCURR > 29) | is.na(MES_OCURR) | is.na(DIA_OCURR),
embarazo_solo_mujeres = !(SEXO == 1 & EMBARAZO %in% 1:5) | is.na(EMBARAZO)
)
result_consistencia <- ejecutar_validaciones(dataset, reglas_consistencia, "consistencia")
# ═══════════════════════════════════════════════════════════════════
# 4. ACCURACY
# ═══════════════════════════════════════════════════════════════════
print_section("DIMENSION 4: ACCURACY")
result_accuracy <- analizar_accuracy(dataset)
# ═══════════════════════════════════════════════════════════════════
# SCORE TOTAL
# ═══════════════════════════════════════════════════════════════════
score_total <- mean(c(
result_completitud$score_general,
result_integridad$score,
result_consistencia$score,
result_accuracy$score
), na.rm = TRUE)
print_section("SCORE TOTAL DE CALIDAD", "=")
cat(sprintf(" 1. Completitud : %6.2f%%\n", result_completitud$score_general))
cat(sprintf(" 2. Integridad : %6.2f%%\n", result_integridad$score))
cat(sprintf(" 3. Consistencia : %6.2f%%\n", result_consistencia$score))
cat(sprintf(" 4. Accuracy : %6.2f%%\n\n", result_accuracy$score))
cat(sprintf(" TOTAL: %.2f%%\n\n", score_total))
clasificacion <- fcase(
score_total >= 95, "EXCELENTE",
score_total >= 85, "BUENA",
score_total >= 70, "ACEPTABLE",
score_total >= 50, "REGULAR",
default = "REQUIERE MEJORA URGENTE"
)
cat(" CLASIFICACION:", clasificacion, "\n\n")
# ═══════════════════════════════════════════════════════════════════
# EXPORTAR
# ═══════════════════════════════════════════════════════════════════
ruta_reportes <- if(is.null(ruta_salida)) "reportes_calidad_inegi" else ruta_salida
dir.create(ruta_reportes, showWarnings = FALSE, recursive = TRUE)
cat("Exportando archivos...\n")
# Exportar uno por uno y limpiar memoria
fwrite(result_completitud$tabla,
file.path(ruta_reportes, "1_completitud_detalle.csv"),
bom = TRUE)
cat(" 1_completitud_detalle.csv\n")
if(!is.null(result_integridad$resumen)) {
fwrite(result_integridad$resumen,
file.path(ruta_reportes, "2_integridad_detalle.csv"),
bom = TRUE)
cat(" 2_integridad_detalle.csv\n")
}
if(!is.null(result_consistencia$resumen)) {
fwrite(result_consistencia$resumen,
file.path(ruta_reportes, "3_consistencia_detalle.csv"),
bom = TRUE)
cat(" 3_consistencia_detalle.csv\n")
}
if(!is.null(result_accuracy$resumen)) {
fwrite(result_accuracy$resumen,
file.path(ruta_reportes, "4_accuracy_detalle.csv"),
bom = TRUE)
cat(" 4_accuracy_detalle.csv\n")
}
resumen_ejecutivo <- data.table(
Dimension = c("Completitud General", "Completitud Criticas",
"Integridad", "Consistencia", "Accuracy", "TOTAL"),
Score = c(result_completitud$score_general, result_completitud$score_criticas,
result_integridad$score, result_consistencia$score,
result_accuracy$score, score_total),
Clasificacion = c(rep("", 5), clasificacion)
)
fwrite(resumen_ejecutivo,
file.path(ruta_reportes, "0_resumen_ejecutivo.csv"),
bom = TRUE)
cat(" 0_resumen_ejecutivo.csv\n\n")
gc() # Limpiar memoria final
# ═══════════════════════════════════════════════════════════════════
# RETURN
# ═══════════════════════════════════════════════════════════════════
list(
completitud = result_completitud$tabla,
score_completitud = result_completitud$score_general,
score_criticas = result_completitud$score_criticas,
score_integridad = result_integridad$score,
score_consistencia = result_consistencia$score,
score_accuracy = result_accuracy$score,
score_total = score_total,
clasificacion = clasificacion,
resumen_integridad = result_integridad$resumen,
resumen_consistencia = result_consistencia$resumen,
resumen_accuracy = result_accuracy$resumen,
ruta_reportes = ruta_reportes
)
}
# ═══════════════════════════════════════════════════════════════════
# EJECUTAR CON CONFIGURACIÓN DE MEMORIA
# ═══════════════════════════════════════════════════════════════════
# Aumentar límite de memoria si es posible
#memory.limit(size = 32000) # 32GB (solo Windows)
# Ejecutar
resultados <- evaluar_calidad_mortalidad_oficial(
dataset,
ruta_salida = "~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/DQA_ITS"
)##
## ======================================================================
## EVALUACION DE CALIDAD - DEFUNCIONES MEXICO
## ======================================================================
##
## Dataset: 20,402,630 registros
## Memoria en uso: 6.1 Gb
##
## ----------------------------------------------------------------------
## DIMENSION 1: COMPLETENESS
## ----------------------------------------------------------------------
##
## Calculando completitud...
## RESUMEN:
## General: 65.67% | Criticas: 100.00% | Importantes: 97.63%
## 100% completas: 29 de 79
##
## ----------------------------------------------------------------------
## DIMENSION 2: INTEGRITY
## ----------------------------------------------------------------------
##
## Ejecutando validaciones de integridad ...
## RESUMEN DE INTEGRIDAD :
## Score: 97.53%
## Reglas pasadas: 5 de 7
##
## ----------------------------------------------------------------------
## DIMENSION 3: CONSISTENCY
## ----------------------------------------------------------------------
##
## Ejecutando validaciones de consistencia ...
## RESUMEN DE CONSISTENCIA :
## Score: 99.56%
## Reglas pasadas: 1 de 3
##
## ----------------------------------------------------------------------
## DIMENSION 4: ACCURACY
## ----------------------------------------------------------------------
##
## Analizando edad...
## Media: 69.3 | Mediana: 68 | Outliers: 2.33%
##
## Analizando sexo...
## DISTRIBUCION DE SEXO:
## Sexo_Desc N Pct
## <char> <int> <num>
## 1: Mujer 8917688 43.74
## 2: Hombre 11471519 56.26
##
## Desviacion: 4.26%
##
## SCORE DE ACCURACY: 88.19%
##
## ======================================================================
## SCORE TOTAL DE CALIDAD
## ======================================================================
##
## 1. Completitud : 65.67%
## 2. Integridad : 97.53%
## 3. Consistencia : 99.56%
## 4. Accuracy : 88.19%
##
## TOTAL: 87.74%
##
## CLASIFICACION: BUENA
##
## Exportando archivos...
## 1_completitud_detalle.csv
## 2_integridad_detalle.csv
## 3_consistencia_detalle.csv
## 4_accuracy_detalle.csv
## 0_resumen_ejecutivo.csvlibrary(data.table)
# ═══════════════════════════════════════════════════════════════════
# TABLA RESUMEN DE REGLAS DE VALIDACIÓN
# ═══════════════════════════════════════════════════════════════════
ruta_reportes = "~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/DQA_ITS"
crear_tabla_resumen_reglas <- function(ruta_salida = NULL) {
if(is.null(ruta_salida)) {
ruta_salida <- ruta_reportes
}
dir.create(ruta_salida, showWarnings = FALSE, recursive = TRUE)
# ═══════════════════════════════════════════════════════════════
# DIMENSIÓN 1: COMPLETENESS
# ═══════════════════════════════════════════════════════════════
reglas_completitud <- data.table(
Dimension = "1. COMPLETENESS",
Categoria = c(
"Variables Críticas",
"Variables Críticas",
"Variables Críticas",
"Variables Importantes",
"Variables Opcionales"
),
Regla = c(
"0% valores faltantes permitidos",
"Variables geográficas completas",
"Variables demográficas básicas completas",
"Máximo 5-10% valores faltantes",
"Máximo 20% valores faltantes"
),
Variables_Aplicables = c(
"ENT_REGIS, MUN_REGIS, ENT_RESID, MUN_RESID, ENT_OCURR, MUN_OCURR, ANIO_OCUR, MES_OCURR, SEXO, EDAD, CAUSA_DEF",
"Todas las variables geográficas de registro, residencia y ocurrencia",
"SEXO, EDAD, CAUSA_DEF, ANIO_OCUR, MES_OCURR",
"DIA_OCURR, ESCOLARIDA, EDO_CIVIL, ASIST_MEDI, SITIO_OCUR, DERECHOHAB, NACIONALID",
"OCUPACION, LENGUA, AFROMEX, CONINDIG, NECROPSIA, PRESUNTO, OCURR_TRAB"
),
Criterio_Evaluacion = c(
"Porcentaje de completitud = 100%",
"Sin valores NA en variables geográficas",
"Sin valores NA en variables demográficas básicas",
"Porcentaje de completitud ≥ 90%",
"Porcentaje de completitud ≥ 80%"
)
)
# ═══════════════════════════════════════════════════════════════
# DIMENSIÓN 2: INTEGRITY (Integridad de Dominios)
# ═══════════════════════════════════════════════════════════════
reglas_integridad <- data.table(
Dimension = "2. INTEGRITY",
Categoria = c(
"Geografía - Entidades",
"Geografía - Municipios",
"Fechas - Día",
"Fechas - Mes",
"Fechas - Año ocurrencia",
"Fechas - Año nacimiento",
"Demográficas - Sexo",
"Demográficas - Edad",
"Sociales - Escolaridad",
"Sociales - Estado civil",
"Médicas - Asistencia",
"Médicas - Sitio ocurrencia",
"Médicas - Derechohabiencia"
),
Regla = c(
"Entidades en rango válido según catálogo INEGI",
"Municipios en rango válido según catálogo INEGI",
"Día entre 1-31 o 99 (no especificado)",
"Mes entre 1-12 o 99 (no especificado)",
"Año entre 1900-2024 o 9999 (no especificado)",
"Año entre 1890-2024 o 9999 (no especificado)",
"Valores: 1=Hombre, 2=Mujer, 9=No especificado",
"Formato codificado: 1xxx=horas, 2xxx=días, 3xxx=meses, 4xxx=años",
"Valores 1-10 (niveles educativos), 88=No aplica, 99=No especificado",
"Valores 1-6 (estados civiles), 8=No aplica, 9=No especificado",
"Valores: 1=Con asistencia, 2=Sin asistencia, 9=No especificado",
"Valores 1-13 (lugares de defunción), 99=No especificado",
"Valores 1-10, 12 (instituciones de salud), 99=No especificado"
),
Variables_Aplicables = c(
"ENT_REGIS, ENT_RESID, ENT_OCURR",
"MUN_REGIS, MUN_RESID, MUN_OCURR",
"DIA_OCURR, DIA_NACIM",
"MES_OCURR, MES_NACIM",
"ANIO_OCUR",
"ANIO_NACIM",
"SEXO",
"EDAD",
"ESCOLARIDA",
"EDO_CIVIL",
"ASIST_MEDI",
"SITIO_OCUR",
"DERECHOHAB"
),
Criterio_Evaluacion = c(
"ENT_REGIS (1-32), ENT_RESID (1-32,35,99), ENT_OCURR (1-32,99)",
"MUN_REGIS (1-570), MUN_RESID (1-570,999), MUN_OCURR (1-570,999)",
"Rango: 1-31, 99 o NA",
"Rango: 1-12, 99 o NA",
"Rango: 1900-2024, 9999",
"Rango: 1890-2024, 9999 o NA",
"Dominio: {1, 2, 9}",
"Rangos: 1001-1098, 2001-2098, 3001-3098, 4001-4998 o NA",
"Dominio: {1,2,3,4,5,6,7,8,9,10,88,99} o NA",
"Dominio: {1,2,3,4,5,6,8,9} o NA",
"Dominio: {1, 2, 9} o NA",
"Dominio: {1-13, 99} o NA",
"Dominio: {1-10, 12, 99} o NA"
)
)
# ═══════════════════════════════════════════════════════════════
# DIMENSIÓN 3: CONSISTENCY (Coherencia entre variables)
# ═══════════════════════════════════════════════════════════════
reglas_consistencia <- data.table(
Dimension = "3. CONSISTENCY",
Categoria = c(
"Coherencia Temporal",
"Coherencia Temporal",
"Coherencia de Fechas",
"Coherencia de Fechas",
"Coherencia Demográfica",
"Coherencia Demográfica",
"Coherencia Educativa",
"Coherencia Social"
),
Regla = c(
"Año de nacimiento debe ser anterior o igual al año de ocurrencia",
"Año de registro debe ser posterior o igual al año de ocurrencia",
"Febrero no puede tener más de 29 días",
"Meses de 30 días (abril, junio, sept, nov) no pueden tener día 31",
"Embarazo solo puede ocurrir en mujeres (SEXO=2)",
"Embarazo solo en edad fértil (10-54 años)",
"Menores de 3 años no deben tener escolaridad 2-10",
"Menores de 12 años no deben tener estado civil 1-6"
),
Variables_Aplicables = c(
"ANIO_NACIM, ANIO_OCUR",
"ANIO_REGIS, ANIO_OCUR",
"MES_OCURR, DIA_OCURR",
"MES_OCURR, DIA_OCURR",
"SEXO, EMBARAZO",
"EDAD, EMBARAZO",
"EDAD, ESCOLARIDA",
"EDAD, EDO_CIVIL"
),
Criterio_Evaluacion = c(
"ANIO_NACIM ≤ ANIO_OCUR (excepto 9999=no especificado)",
"ANIO_REGIS ≥ ANIO_OCUR",
"Si MES_OCURR=2 entonces DIA_OCURR ≤ 29",
"Si MES_OCURR ∈ {4,6,9,11} entonces DIA_OCURR ≤ 30",
"Si EMBARAZO ∈ {1,2,3,4,5} entonces SEXO=2",
"Si EMBARAZO ∈ {1-5} entonces EDAD entre 4010-4054 (10-54 años)",
"Si EDAD<4003 (0-2 años) entonces ESCOLARIDA ∉ {2-10}",
"Si EDAD<4012 (<12 años) entonces EDO_CIVIL ∉ {1-6}"
)
)
# ═══════════════════════════════════════════════════════════════
# DIMENSIÓN 4: ACCURACY (Precisión y distribuciones esperadas)
# ═══════════════════════════════════════════════════════════════
reglas_accuracy <- data.table(
Dimension = "4. ACCURACY",
Categoria = c(
"Distribución de Edad",
"Distribución de Edad",
"Distribución de Sexo",
"Distribución de Sexo",
"Valores Extremos"
),
Regla = c(
"Detección de outliers mediante método IQR",
"Edad razonable: menores de 120 años",
"Distribución esperada: ~52% hombres, ~48% mujeres",
"Desviación máxima aceptable: ±5% de distribución esperada",
"Años de ocurrencia dentro del periodo esperado (1950-2024)"
),
Variables_Aplicables = c(
"EDAD (años: códigos 4001-4998)",
"EDAD (años: códigos 4001-4998)",
"SEXO",
"SEXO",
"ANIO_OCUR"
),
Criterio_Evaluacion = c(
"Outliers: valores fuera de [Q1-1.5×IQR, Q3+1.5×IQR]. Score = 100 - %outliers",
"Alertas para EDAD > 4110 (>110 años)",
"Comparación con proporción poblacional histórica de México",
"Score = 100 - (|%masculino - 52| × 5). Penalización por desviación",
"Alertas para ANIO_OCUR < 1950 o > 2024"
)
)
# ═══════════════════════════════════════════════════════════════
# COMBINAR TODAS LAS DIMENSIONES
# ═══════════════════════════════════════════════════════════════
tabla_completa <- rbindlist(list(
reglas_completitud,
reglas_integridad,
reglas_consistencia,
reglas_accuracy
), use.names = TRUE, fill = TRUE)
# ═══════════════════════════════════════════════════════════════
# TABLA RESUMEN EJECUTIVO
# ═══════════════════════════════════════════════════════════════
resumen_ejecutivo <- data.table(
Dimension = c("1. COMPLETENESS", "2. INTEGRITY", "3. CONSISTENCY", "4. ACCURACY"),
Objetivo = c(
"Medir la presencia de valores en los datos",
"Validar que los valores estén dentro de rangos permitidos",
"Verificar coherencia lógica entre variables relacionadas",
"Evaluar precisión mediante distribuciones estadísticas esperadas"
),
Total_Reglas = c(
nrow(reglas_completitud),
nrow(reglas_integridad),
nrow(reglas_consistencia),
nrow(reglas_accuracy)
),
Metrica_Principal = c(
"Porcentaje de completitud por variable",
"Porcentaje de valores dentro de dominio válido",
"Porcentaje de registros que cumplen reglas de coherencia",
"Análisis de outliers y desviación de distribuciones esperadas"
),
Fuente_Reglas = c(
"Clasificación por criticidad de variables (INEGI)",
"Guía Oficial INEGI 2024 (páginas 5-24)",
"Lógica de negocio y relaciones entre variables",
"Análisis estadístico y distribuciones poblacionales"
)
)
# ═══════════════════════════════════════════════════════════════
# EXPORTAR TABLAS
# ═══════════════════════════════════════════════════════════════
fwrite(tabla_completa,
file.path(ruta_salida, "REGLAS_VALIDACION_DETALLE.csv"),
bom = TRUE)
fwrite(resumen_ejecutivo,
file.path(ruta_salida, "REGLAS_VALIDACION_RESUMEN.csv"),
bom = TRUE)
# También crear por separado
fwrite(reglas_completitud,
file.path(ruta_salida, "REGLAS_1_COMPLETITUD.csv"),
bom = TRUE)
fwrite(reglas_integridad,
file.path(ruta_salida, "REGLAS_2_INTEGRIDAD.csv"),
bom = TRUE)
fwrite(reglas_consistencia,
file.path(ruta_salida, "REGLAS_3_CONSISTENCIA.csv"),
bom = TRUE)
fwrite(reglas_accuracy,
file.path(ruta_salida, "REGLAS_4_ACCURACY.csv"),
bom = TRUE)
cat("═══════════════════════════════════════════════════════════════════\n")
cat("✅ TABLAS DE REGLAS DE VALIDACIÓN GENERADAS\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
cat("Archivos generados en:", ruta_salida, "\n\n")
cat("📄 Archivos creados:\n")
cat(" 1. REGLAS_VALIDACION_RESUMEN.csv (", nrow(resumen_ejecutivo), "dimensiones)\n")
cat(" 2. REGLAS_VALIDACION_DETALLE.csv (", nrow(tabla_completa), "reglas totales)\n")
cat(" 3. REGLAS_1_COMPLETITUD.csv (", nrow(reglas_completitud), "reglas)\n")
cat(" 4. REGLAS_2_INTEGRIDAD.csv (", nrow(reglas_integridad), "reglas)\n")
cat(" 5. REGLAS_3_CONSISTENCIA.csv (", nrow(reglas_consistencia), "reglas)\n")
cat(" 6. REGLAS_4_ACCURACY.csv (", nrow(reglas_accuracy), "reglas)\n\n")
# Mostrar resumen en consola
cat("═══════════════════════════════════════════════════════════════════\n")
cat("RESUMEN EJECUTIVO DE REGLAS DE VALIDACIÓN\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
print(resumen_ejecutivo)
cat("\n")
cat("═══════════════════════════════════════════════════════════════════\n")
cat("DETALLE POR DIMENSIÓN (primeras 3 reglas de cada una)\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
for(dim in unique(tabla_completa$Dimension)) {
cat("\n", dim, "\n")
cat(strrep("─", 70), "\n")
reglas_dim <- tabla_completa[Dimension == dim]
print(reglas_dim[1:min(3, nrow(reglas_dim)), .(Categoria, Regla)])
cat(" ... (", nrow(reglas_dim), "reglas totales)\n")
}
cat("\n")
return(list(
resumen = resumen_ejecutivo,
detalle_completo = tabla_completa,
completitud = reglas_completitud,
integridad = reglas_integridad,
consistencia = reglas_consistencia,
accuracy = reglas_accuracy
))
}
# ═══════════════════════════════════════════════════════════════════
# EJECUTAR
# ═══════════════════════════════════════════════════════════════════
tablas_reglas <- crear_tabla_resumen_reglas(
ruta_salida = ruta_reportes
)## ═══════════════════════════════════════════════════════════════════
## ✅ TABLAS DE REGLAS DE VALIDACIÓN GENERADAS
## ═══════════════════════════════════════════════════════════════════
##
## Archivos generados en: ~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/DQA_ITS
##
## 📄 Archivos creados:
## 1. REGLAS_VALIDACION_RESUMEN.csv ( 4 dimensiones)
## 2. REGLAS_VALIDACION_DETALLE.csv ( 31 reglas totales)
## 3. REGLAS_1_COMPLETITUD.csv ( 5 reglas)
## 4. REGLAS_2_INTEGRIDAD.csv ( 13 reglas)
## 5. REGLAS_3_CONSISTENCIA.csv ( 8 reglas)
## 6. REGLAS_4_ACCURACY.csv ( 5 reglas)
##
## ═══════════════════════════════════════════════════════════════════
## RESUMEN EJECUTIVO DE REGLAS DE VALIDACIÓN
## ═══════════════════════════════════════════════════════════════════
##
## Dimension
## <char>
## 1: 1. COMPLETENESS
## 2: 2. INTEGRITY
## 3: 3. CONSISTENCY
## 4: 4. ACCURACY
## Objetivo
## <char>
## 1: Medir la presencia de valores en los datos
## 2: Validar que los valores estén dentro de rangos permitidos
## 3: Verificar coherencia lógica entre variables relacionadas
## 4: Evaluar precisión mediante distribuciones estadísticas esperadas
## Total_Reglas Metrica_Principal
## <int> <char>
## 1: 5 Porcentaje de completitud por variable
## 2: 13 Porcentaje de valores dentro de dominio válido
## 3: 8 Porcentaje de registros que cumplen reglas de coherencia
## 4: 5 Análisis de outliers y desviación de distribuciones esperadas
## Fuente_Reglas
## <char>
## 1: Clasificación por criticidad de variables (INEGI)
## 2: Guía Oficial INEGI 2024 (páginas 5-24)
## 3: Lógica de negocio y relaciones entre variables
## 4: Análisis estadístico y distribuciones poblacionales
##
## ═══════════════════════════════════════════════════════════════════
## DETALLE POR DIMENSIÓN (primeras 3 reglas de cada una)
## ═══════════════════════════════════════════════════════════════════
##
##
## 1. COMPLETENESS
## ──────────────────────────────────────────────────────────────────────
## Categoria Regla
## <char> <char>
## 1: Variables Críticas 0% valores faltantes permitidos
## 2: Variables Críticas Variables geográficas completas
## 3: Variables Críticas Variables demográficas básicas completas
## ... ( 5 reglas totales)
##
## 2. INTEGRITY
## ──────────────────────────────────────────────────────────────────────
## Categoria Regla
## <char> <char>
## 1: Geografía - Entidades Entidades en rango válido según catálogo INEGI
## 2: Geografía - Municipios Municipios en rango válido según catálogo INEGI
## 3: Fechas - Día Día entre 1-31 o 99 (no especificado)
## ... ( 13 reglas totales)
##
## 3. CONSISTENCY
## ──────────────────────────────────────────────────────────────────────
## Categoria
## <char>
## 1: Coherencia Temporal
## 2: Coherencia Temporal
## 3: Coherencia de Fechas
## Regla
## <char>
## 1: Año de nacimiento debe ser anterior o igual al año de ocurrencia
## 2: Año de registro debe ser posterior o igual al año de ocurrencia
## 3: Febrero no puede tener más de 29 días
## ... ( 8 reglas totales)
##
## 4. ACCURACY
## ──────────────────────────────────────────────────────────────────────
## Categoria Regla
## <char> <char>
## 1: Distribución de Edad Detección de outliers mediante método IQR
## 2: Distribución de Edad Edad razonable: menores de 120 años
## 3: Distribución de Sexo Distribución esperada: ~52% hombres, ~48% mujeres
## ... ( 5 reglas totales)# Ver en RStudio
View(tablas_reglas$resumen)
View(tablas_reglas$detalle_completo)Codificación del dataset
Para optimizar memoria, se utiliza el setDT, y al finalizar toda la etapa de procesamiento hay que revertirlo mediante setDF
df <- copy(df) # df_backup
# dataset en memoria
setDT(df)Deseleccionar los registros sin fecha de ocurrencia.
df <- df[DIA_OCURR != 99 &
MES_OCURR != 99 &
ANIO_OCUR != 9999]Creación de columna fecha
library(data.table)
# ISOdate es más rápido que paste + as.Date
df[, fecha_mes := as.Date(ISOdate(ANIO_OCUR, MES_OCURR, 1))]
# Comprobación
df[is.na(fecha_mes)]## Empty data.table (0 rows and 15 cols): ENT_RESID,MUN_RESID,TLOC_RESID,CAUSA_DEF,SEXO,DIA_OCURR...Conversión de valores: numéricos y carácter
library(data.table)
# Convertir factores a character
df[, CAUSA_DEF := as.character(CAUSA_DEF)]
df[, EDAD_AGRU := as.character(EDAD_AGRU)]
# Convertir variables categóricas de int a character (antes de recodificar)
df[, SEXO := as.character(SEXO)]
df[, ESCOLARIDA := as.character(ESCOLARIDA)]
df[, SITIO_OCUR := as.character(SITIO_OCUR)]
df[, DERECHOHAB := as.character(DERECHOHAB)]
df[, ASIST_MEDI := as.character(ASIST_MEDI)]
# Verificar cambios
cat("Tipos de datos después de la conversión:\n")## Tipos de datos después de la conversión:sapply(df, class)## ENT_RESID MUN_RESID TLOC_RESID CAUSA_DEF SEXO DIA_OCURR
## "integer" "integer" "integer" "character" "character" "integer"
## MES_OCURR ANIO_OCUR ESCOLARIDA ASIST_MEDI SITIO_OCUR DERECHOHAB
## "integer" "integer" "character" "character" "character" "character"
## EDAD_AGRU AREA_UR fecha_mes
## "character" "integer" "Date"Cambio de causa de defunción
df[CAUSA_DEF == "2040", CAUSA_DEF := "C910"]Género
# Carga dplyr después de plyr
library(plyr)
library(dplyr)unique(df$SEXO)## [1] "2" "1" "9"# Renombrar columna
setnames(df, "SEXO", "GENERO")
# Recodificar valores
df[, GENERO := fcase(
GENERO == 1, "masculino",
GENERO == 2, "femenino",
default = NA_character_
)]
# Verificar
table(df$GENERO)##
## femenino masculino
## 7443095 9539923Densidad poblacional en la localidad
unique(df$TLOC_RESID)## [1] 14 5 1 11 15 7 2 13 99 16 10 9 6 4 3 12 17 8# Crear mapeo
mapeo_tloc <- c(
"1" = "1-999",
"4" = "2.5K < 5k",
"5" = "5k < 10k",
"6" = "10k < 15k",
"7" = "15k < 20k",
"8" = "20k < 30k",
"9" = "30k < 40k",
"13" = "100k < 250k",
"14" = "250k < 500k",
"15" = "500k < 1M",
"16" = "1M < 2M",
"17" = "> 2M",
"99" = "n/esp"
)
# Recodificar (valores no mapeados = NA automáticamente)
df[, TLOC_RESID := mapeo_tloc[as.character(TLOC_RESID)]]
# Verificar
table(df$TLOC_RESID)##
## > 2M 1-999 100k < 250k 10k < 15k 15k < 20k 1M < 2M
## 759469 2365657 1239154 561867 383569 1538878
## 2.5K < 5k 20k < 30k 250k < 500k 30k < 40k 500k < 1M 5k < 10k
## 925042 602825 2082390 419003 2430513 949478
## n/esp
## 320091Asistencia Medica
unique(df$ASIST_MEDI)## [1] "1" "2" "9"df[, ASIST_MEDI := fcase(
ASIST_MEDI == "1", "Si",
ASIST_MEDI == "2", "No",
ASIST_MEDI == "9", "n/esp"
)]
# Verificar
table(df$ASIST_MEDI)##
## n/esp No Si
## 1052393 2435016 13504932Edad agrupado
df[, EDAD_AGRU := as.numeric(EDAD_AGRU)]
unique(df$EDAD_AGRU)## [1] 23 30 20 15 1 21 24 19 17 8 9 25 18 16 12 22 14 10 11 13 6 5 7 2 4
## [26] 3 26 27 28 29mapeo_edad_agrup <- c("1" = "< 1",
"2" = "1",
"3" = "2",
"4" = "3",
"5" = "4",
"6" = "5-9",
"7" = "10-14"
)
# Recodificar: si no está en el mapeo, asignar "Fuera de rango"
df[, EDAD_AGRU := ifelse(
as.character(EDAD_AGRU) %in% names(mapeo_edad_agrup),
mapeo_edad_agrup[as.character(EDAD_AGRU)],
"Fuera de rango"
)]
# Verificar
table(df$EDAD_AGRU, useNA = "always")##
## < 1 1 10-14 2 3
## 763960 69815 94545 37103 25025
## 4 5-9 Fuera de rango <NA>
## 19903 77244 15904746 0Sitio de ocurrencia
Cada año contiene sitios de ocurrencia diferentes, por lo que se recurre realizar una segmentación para una homologación de términos similar al último diccionario de datos de la EDR-24.
Por ejemplo:
- los nombres de los programas gubernamentales han cambiado al paso de los años, y para ello se les define en la entidad/institución que opera el programa.
- todo sitio de ocurrencia registrado en alguna unidad medica publica, se le asigna a la dependencia de la Secretaria de Salud.
unique(df$ANIO_OCUR)## [1] 2000 1999 1994 1997 1900 1992 1996 1920 1998 1990 1985 1993 1970 1995 1965
## [16] 1987 1912 1991 1986 1961 1989 1967 1951 1978 1981 1980 1966 1988 1984 1964
## [31] 1983 1974 1937 1972 1982 1977 1949 1933 1975 1973 1969 1971 1954 1931 1979
## [46] 1952 1976 1962 1959 1944 1940 1953 1955 1968 1963 1922 1957 1960 1921 1909
## [61] 1939 1935 1913 1910 1947 1925 1956 1958 1928 1901 1924 1902 1948 1917 1930
## [76] 1936 1915 1946 1942 1941 1938 1906 1919 2001 1945 1914 1934 1932 1911 1923
## [91] 1907 1929 1943 1950 1918 1926 2002 1927 1916 2003 1903 2004 1908 1904 2005
## [106] 1905 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
## [121] 2020 2021 2022 2023 2024# 1) PRIMERO: Corregir años 91-99 → 1991-1999
df[ANIO_OCUR >= 90 & ANIO_OCUR <=99, ANIO_OCUR := 1900 + ANIO_OCUR]
# filtrar valores desde 1990
df <- df[ANIO_OCUR >= 1990]
# Verificar
sort(unique(df$ANIO_OCUR))## [1] 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004
## [16] 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
## [31] 2020 2021 2022 2023 2024unique(df$SITIO_OCUR)## [1] "3" "4" "1" "2" "9" "11" "8" "10" "99" "12" "6" "7" "5" "13"La combicación de codficación es la siguiente
PERÍODO ≤ 2003 (Hasta 2003) ═══════════════════════════════════════════════════════════════════════ Código → Significado | Descripción 1 → SSA | Unidad Médica Pública 2 → UMPriv | Unidad Médica Privada 3 → Hogar | En el hogar 4 → Otr.lug | Otro lugar 9 → n/esp | No especificado
PERÍODO ≤ 2003 (Hasta 2003) ═══════════════════════════════════════════════════════════════════════ Código → Significado | Descripción 1 → SSA | Unidad Médica Pública 2 → UMPriv | Unidad Médica Privada 3 → Hogar | En el hogar 4 → Otr.lug | Otro lugar 9 → n/esp | No especificado
PERÍODO ≥ 2005 (De 2005 en adelante) ═══════════════════════════════════════════════════════════════════════ Código → Significado | Descripción 1 → SSA | Secretaría de Salud 2 → IMSS | Instituto Mexicano del Seguro Social 3 → IMSS | IMSS (reclasificado desde Hogar) 4 → ISSSTE | Instituto de Seguridad y Servicios Sociales 5 → PEMEX | Petróleos Mexicanos 6 → SEDENA | Secretaría de la Defensa Nacional 7 → SEMAR | Secretaría de Marina 8 → SSA | Otra Unidad Médica 9 → UMPriv | Unidad Médica Privada 10 → Via Pub | Vía Pública 11 → Hogar | En el hogar 12 → Otr.lug | Otro lugar 99 → n/esp | No especificado
library(data.table)
df[, SITIO_OCUR := as.character(SITIO_OCUR)]
df[, SITIO_ORIGINAL := SITIO_OCUR]
# <= 2003
mapeo_pre2004 <- c("1"="SSA",
"2"="UMPriv",
"3"="Hogar",
"4"="Otr.lug",
"9"="n/esp")
df[ANIO_OCUR <= 2003, SITIO_OCUR := mapeo_pre2004[SITIO_ORIGINAL]]
# >= 2004 (incluye 2004 y adelante)
mapeo_post2004 <- c("1"="SSA",
"2"="IMSS",
"3"="IMSS",
"4"="ISSSTE",
"5"="PEMEX",
"6"="SEDENA",
"7"="SEMAR",
"8"="SSA",
"9"="UMPriv",
"10"="Via Pub",
"11"="Hogar",
"12"="Otr.lug",
"99"="n/esp")
df[ANIO_OCUR >= 2004, SITIO_OCUR := mapeo_post2004[SITIO_ORIGINAL]]
df[, SITIO_ORIGINAL := NULL]
table(df$SITIO_OCUR)##
## Hogar IMSS ISSSTE n/esp Otr.lug PEMEX SEDENA SEMAR SSA UMPriv
## 7564215 2757397 502128 576988 766346 57788 55372 19173 3247477 807408
## Via Pub
## 612668Derechohabiente de seguro medico
Cada año contiene clasificaciones diferentes en el tipo de seguro medico, por lo que se recurre realizar una segmentación para una homologación de términos similar al último diccionario de datos de la EDR-24.
Por ejemplo:
- existen casos en donde se registran seguros en cojunto IMSS/ISSTE o ISSTE/IMSS, por lo que se asigna a la institución que aparece primeramente en la relacion.
- el seguro medico de las fuerzas armadas (SEDENA o SEMARNAT) se le asigna al tipo de institución medica que atiende a ambos grupos
- los programas federales auspiciados por alguna institución de salud (IMSS o ISSTE), se les define bajo el concepto de Aseguramiento Público
#df[, DERECHOHAB := as.numeric(DERECHOHAB)]
unique(df$DERECHOHAB)## [1] "1" "99" "2" "3" "4" "6" "5" "23" "25" "26" "35" "24" "36" "34" "46"
## [16] "56" "45" "7" "8" "9" "10" "12"library(data.table)
# Guardar original
df[, DERECHO_ORIG := DERECHOHAB]
# ≤2003
m1 <- c("1"="Ninguna", "2"="IMSS", "3"="ISSSTE", "4"="PEMEX",
"5"="ISSFAM", "6"="IMSS", "36"="ISSSTE", "99"="Ninguna")
df[ANIO_OCUR <= 2003, DERECHOHAB := m1[DERECHO_ORIG]]
# 2004-2011
m2 <- c("1"="Ninguna", "2"="IMSS", "3"="ISSSTE", "4"="PEMEX",
"5"="ISSFAM", "6"="ISSFAM", "7"="ASP", "8"="Otra", "99"="Ninguna")
df[ANIO_OCUR >= 2004 & ANIO_OCUR <= 2011, DERECHOHAB := m2[DERECHO_ORIG]]
# 2012-2017
m3 <- c("1"="Ninguna", "2"="IMSS", "3"="ISSSTE", "4"="PEMEX",
"5"="ISSFAM", "6"="ISSFAM", "7"="ASP", "8"="Otra",
"9"="ASP", "99"="Ninguna")
df[ANIO_OCUR >= 2012 & ANIO_OCUR <= 2017, DERECHOHAB := m3[DERECHO_ORIG]]
# 2018-2021
m4 <- c("1"="Ninguna", "2"="IMSS", "3"="ISSSTE", "4"="PEMEX",
"5"="ISSFAM", "6"="ISSFAM", "7"="ASP", "8"="Otra",
"9"="ASP", "99"="Ninguna")
df[ANIO_OCUR >= 2018 & ANIO_OCUR <= 2021, DERECHOHAB := m4[DERECHO_ORIG]]
# 2022-2023
m5 <- c("1"="Ninguna", "2"="IMSS", "3"="ISSSTE", "4"="PEMEX",
"5"="ISSFAM", "6"="ISSFAM", "7"="ASP", "8"="Otra",
"9"="ASP", "10"="ISSFAM", "99"="Ninguna")
df[ANIO_OCUR >= 2022, DERECHOHAB := m5[DERECHO_ORIG]]
# Verificar
table(df$DERECHOHAB)##
## ASP IMSS ISSFAM ISSSTE Ninguna Otra PEMEX
## 2039548 5804511 127628 1216061 7334277 322275 137652Escolaridad
Cada año contiene clasificaciones diferentes sobre el nivel educativo, por lo que se recurre realizar una segmentación para una homologación de términos similar al último diccionario de datos de la EDR-24.
Por ejemplo:
- escolaridad desconocida, se define como s/escolaridad (pero podria asignarse al nivel de edad)
- primer, segundo o tercer tramo de primaria se define simplemente como nivel primara. Mismo caso aplica al nivel inicial.
- en los casos donde se indica el nivel educativo trunco, se le asigna al nivel educativo sin importar el status de avance.
- la edad menor a 3 años no aplica una valoración en el sistema educativo, por lo que se le asigna al criterio de s/escolaridad.
#df[, ESCOLARIDA := as.numeric(ESCOLARIDA)]
unique(df$ESCOLARIDA)## [1] "3" "1" "9" "4" "7" "6" "2" "5" "8" "99" "88" "10"library(data.table)
setDT(df)
# Convertir a character
df[, ESCOLARIDA := as.character(ESCOLARIDA)]
# Guardar original
df[, ESCOLAR_ORIG := ESCOLARIDA]
# MAPEOS COMPLETOS (todos los códigos en todos los períodos)
# Período ≤2003
mapeo_pre2004 <- c(
"1" = "s/escol",
"2" = "Primaria",
"3" = "Primaria",
"4" = "Primaria",
"5" = "Secundaria",
"6" = "Bachillerato",
"7" = "Bachillerato", # ← AÑADIDO
"8" = "Preescolar",
"9" = "s/escol",
"10" = "Superior",
"88" = "s/escol",
"99" = "s/escol"
)
df[ANIO_OCUR <= 2003, ESCOLARIDA := mapeo_pre2004[ESCOLAR_ORIG]]
# Período 2004-2011
mapeo_2004_2011 <- c(
"1" = "s/escol",
"2" = "Primaria",
"3" = "Primaria",
"4" = "Secundaria",
"5" = "Secundaria",
"6" = "Bachillerato",
"7" = "Bachillerato",
"8" = "Preescolar",
"9" = "s/escol",
"10" = "Superior",
"88" = "s/escol",
"99" = "s/escol"
)
df[ANIO_OCUR >= 2004 & ANIO_OCUR <= 2011, ESCOLARIDA := mapeo_2004_2011[ESCOLAR_ORIG]]
# Período 2012-2023
mapeo_2012_2023 <- c(
"1" = "s/escol",
"2" = "Preescolar",
"3" = "Primaria",
"4" = "Primaria",
"5" = "Secundaria",
"6" = "Secundaria",
"7" = "Bachillerato",
"8" = "Bachillerato",
"9" = "s/escol",
"10" = "Superior",
"88" = "s/escol",
"99" = "s/escol"
)
df[ANIO_OCUR >= 2012, ESCOLARIDA := mapeo_2012_2023[ESCOLAR_ORIG]]
# Limpiar
df[, ESCOLAR_ORIG := NULL]
# Verificar
table(df$ESCOLARIDA, useNA = "always")##
## Bachillerato Preescolar Primaria s/escol Secundaria Superior
## 1477955 400524 7527670 5512931 2011259 53946
## <NA>
## 0BLOQUE 1: ESTADISTICA DESCIPTIVA
Conjunto de datos CDMX-LLA
Se considera las observaciones por atributo de Residencia ya que esto aísla los casos exclusivos de la CDMX (entidad #9), esto porque cada caso se vincula conforme a su domicilio particular, principal o permanente de la persona. Caso contrario, cuando se utiliza el atributo Registro que inscribe el hecho de defunción en el territorio donde sucedió el caso, sin importar el lugar del Residencia de la persona.
Por otro lado, de acuerdo a la legislación mexicana la poblacion de primera infancia se considera solo los menores de 15 años (agrupaciones del 1 al 7)
Y por ultimo, la LLA se clasifica en el criterio de enfermedades con la sigla C910.
cdmx <- df[ENT_RESID == 9 &
CAUSA_DEF == "C910" &
EDAD_AGRU %in% c("< 1", "1-4", "5-9", "10-14") &
ANIO_OCUR %in% c(1998:2024)
]
setDT(cdmx)Estado general de LLA
library(validate)
library(ggplot2) ruta_graficos<-"~/Library/CloudStorage/GoogleDrive-jhonagaona@gmail.com/My Drive/UBA_FCE/TESIS/Graficos_ITS"Como primera aproximación, se realiza un par de visualizaciones con el objeto de dimensionar la problemática de salud pública que guarda la LLA en la población infantil.
- Nacional
A nivel nacional (32 entidades federativas), se observa una ascendencia de casos por LLA a partir del 2017
# seleccion de entidades con los mismos criterios de interes
nacional<-df[df$CAUSA_DEF=="C910" &
EDAD_AGRU %in% c("< 1", "1-4", "5-9", "10-14") &
ANIO_OCUR %in% c(1998:2024)
]
colSums(is.na(nacional))## ENT_RESID MUN_RESID TLOC_RESID CAUSA_DEF GENERO DIA_OCURR
## 0 0 1689 0 3 0
## MES_OCURR ANIO_OCUR ESCOLARIDA ASIST_MEDI SITIO_OCUR DERECHOHAB
## 0 0 0 0 4 1
## EDAD_AGRU AREA_UR fecha_mes DERECHO_ORIG
## 0 1682 0 0# construir tabla resumen por año
nacional <- with(nacional, table(ANIO_OCUR))
# construir el dataframe
nacional<- as.data.frame(nacional)
names(nacional)[names(nacional) == "ANIO_OCUR"] <- "Año"
names(nacional)[names(nacional) == "Freq"] <- "Frecuencia"
# grafico
casos_nacionales<-
ggplot(data=nacional,
aes(x=Año, y=Frecuencia)) +
geom_bar(stat="identity", fill= "#FD8D3C") +
geom_hline(yintercept = mean(nacional$Frecuencia),
color = "#A52A2A",
linetype = "dashed",
linewidth = 1) +
geom_text(aes(label=Frecuencia),
vjust=1.6,
color="black",
size=2.3) +
scale_x_discrete(guide = guide_axis(angle = 90)) +
labs(x = "Periodo",
y = "Defunciones totales",
title = "Defunciones registradas a nivel nacional",
subtitle = "1998-2024") +
theme_classic()
casos_nacionales
# exportar graficos en directorio
ggsave(filename = "casos_nacionales.png",
plot = casos_nacionales,
width = 159,
height = 100,
units = "mm",
dpi = "retina",
path = ruta_graficos)Cambios de puntos observados a nivel nacional en la EDR por defunciones registradas de LLA
library(tidychangepoint)
nacional_ts <- ts(nacional$Frecuencia,
start = c(1998,1),
frequency=12)
x <- segment(nacional_ts,
method = "pelt",
penalty = "BIC")
tidy(x)## # A tibble: 4 × 10
## region num_obs min max mean sd begin end param_mu param_sigma_hatsq
## <chr> <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 [1,13) 12 395 453 426. 16.1 1 13 426. 238.
## 2 [13,15) 2 409 415 412 4.24 13 15 412 9
## 3 [15,20) 5 414 455 436. 14.6 15 20 436. 170.
## 4 [20,28) 8 365 461 397. 32.2 20 28 397. 909.library(psych)
nacional %>%
filter(Año %in% c(2017:2024)) %>%
pull(Frecuencia) %>%
describe() %>%
t() %>%
as.data.frame() %>%
setNames("Valor")## Valor
## vars 1.0000000
## n 8.0000000
## mean 397.2500000
## sd 32.2390624
## median 390.0000000
## trimmed 397.2500000
## mad 31.1346000
## min 365.0000000
## max 461.0000000
## range 96.0000000
## skew 0.7401734
## kurtosis -0.8304920
## se 11.3982298- Entidades federativas
En particular, la CDMX es la quinta entidad con mayor ocurrencia de casos en el periodo de estudio, antes de ella se encuentra Estado de México, Veracruz, Puebla y Jalisco.
# seleccion de entidades con los mismos criterios de interes
entidades<-df[df$CAUSA_DEF=="C910" &
EDAD_AGRU %in% c("< 1", "1-4", "5-9", "10-14")]
# construir tabla resumen por año
entidades <- with(entidades, table(ENT_RESID))
# construir el dataframe
entidades<-as.data.frame(entidades)
names(entidades)[names(entidades) == "ENT_RESID"] <- "Entidad"
names(entidades)[names(entidades) == "Freq"] <- "Frecuencia"
# relacionar las claves de entidad con sus abreviaturs
entidades<-merge(x = entidades,
y = AGEEML_INEGI,
by = "Entidad",
all = TRUE)
colSums(is.na(entidades))## Entidad Frecuencia NOM_ABR
## 0 0 4entidades <- na.omit(entidades)# grafico
casos_entidades<-
ggplot(data=entidades,
aes(x = reorder(NOM_ABR, -Frecuencia), y = Frecuencia)) +
geom_bar(stat="identity", fill= "#FD8D3C") +
geom_hline(yintercept = mean(entidades$Frecuencia),
color = "#A52A2A",
linetype = "dashed",
linewidth = 1) +
geom_text(aes(label = Frecuencia),
color="black",
size=3,
angle = 90) +
scale_x_discrete(guide = guide_axis(angle = 90)) +
labs(x = "Entidad",
y = "Defunciones totales",
title = "Defunciones registradas acumuladas por entidad",
subtitle = "1998-2024") +
theme_classic()
casos_entidades
# exportar graficos en directorio
ggsave(filename = "casos_entidades.png",
plot = casos_entidades,
width = 159,
height = 100,
units = "mm",
dpi = "retina",
path = ruta_graficos)sum(entidades$Frecuencia)## [1] 11274quantile(entidades$Frecuencia)## 0% 25% 50% 75% 100%
## 60.00 167.50 270.00 438.25 1328.00# asignar la etiqueta quantil sobre el dataframe
entidades_qt<-
entidades %>%
mutate(quantile = ntile(Frecuencia, 3)) %>%
arrange(Frecuencia)
entidades_qt## Entidad Frecuencia NOM_ABR quantile
## 1 6 60 Col. 1
## 2 3 83 BCS 1
## 3 18 96 Nay. 1
## 4 10 100 Dgo. 1
## 5 4 111 Camp. 1
## 6 29 132 Tlax. 1
## 7 32 132 Zac. 1
## 8 1 139 Ags. 1
## 9 25 177 Sin. 1
## 10 17 180 Mor. 1
## 11 23 182 Q. Roo 1
## 12 22 199 Qro. 2
## 13 24 216 SLP 2
## 14 31 238 Yuc. 2
## 15 12 260 Gro. 2
## 16 13 264 Hgo. 2
## 17 26 276 Son. 2
## 18 5 280 Coah. 2
## 19 28 284 Tamps. 2
## 20 8 304 Chih. 2
## 21 27 310 Tab. 2
## 22 2 336 BC 2
## 23 16 372 Mich. 3
## 24 20 436 Oax. 3
## 25 19 445 NL 3
## 26 11 569 Gto. 3
## 27 7 666 Chis. 3
## 28 9 706 CDMX 3
## 29 14 762 Jal. 3
## 30 21 792 Pue. 3
## 31 30 839 Ver. 3
## 32 15 1328 Mex. 3grafico_qt <- ggplot(entidades_qt,
aes(x = quantile,
y = NOM_ABR)) +
geom_point(size = 4, color = "steelblue") +
geom_segment(aes(x = 0,
xend = quantile,
y = NOM_ABR,
yend = NOM_ABR),
color = "gray70") +
labs(title = "Cuantiles por Estado",
x = "Cuantil",
y = "Estado") +
theme_minimal() +
theme(
axis.text.y = element_text(size = 9),
panel.grid.major.y = element_blank(), # ← Elimina líneas horizontales
panel.grid.minor.y = element_blank() # ← Elimina líneas menores
)
grafico_qt
# exportar graficos en directorio
ggsave(filename = "cuantiles_estados.png",
plot = grafico_qt,
width = 159,
height = 100,
units = "mm",
dpi = "retina",
path = ruta_graficos)table(entidades_qt$quantile)##
## 1 2 3
## 11 11 10- Ciudad de México (CDMX)
library(ggplot2)
library(dplyr)
# construir un resumen de frecuencias por año
capital <- with(cdmx, table(ANIO_OCUR))
# construir el dataframe
capital<-as.data.frame(capital)
names(capital)[names(capital) == "ANIO_OCUR"] <- "Año"
names(capital)[names(capital) == "Freq"] <- "Frecuencia"
# grafico
casos_cdmx<-
ggplot(data=capital,
aes(x=Año, y=Frecuencia)) +
geom_bar(stat="identity", fill= "#FD8D3C") +
geom_hline(yintercept = mean(capital$Frecuencia),
color = "#A52A2A",
linetype = "dashed",
linewidth = 1) +
geom_text(aes(label=Frecuencia),
vjust=1.6, color="black", size=3) +
scale_x_discrete(guide = guide_axis(angle = 90)) +
labs(x = "Año",
y = "Defunciones totales",
title = "Defunciones registradas en la CDMX",
subtitle = "1998-2024") +
theme_classic()
casos_cdmx
# exportar graficos en directorio
ggsave(filename = "casos_cdmx.png", plot = casos_cdmx,
width = 159, height = 100, units = "mm", dpi = "retina",
path = ruta_graficos)quantile(capital$Frecuencia)## 0% 25% 50% 75% 100%
## 13.0 23.5 26.0 29.5 40.0# asignar la etiqueta quantil sobre el dataframe
capital %>%
mutate(quantile = ntile(Frecuencia, 3)) %>%
arrange(Frecuencia)## Año Frecuencia quantile
## 1 2020 13 1
## 2 2003 18 1
## 3 2018 18 1
## 4 2009 19 1
## 5 2024 19 1
## 6 2022 20 1
## 7 2015 23 1
## 8 1998 24 1
## 9 2008 24 1
## 10 2011 24 2
## 11 2014 25 2
## 12 2006 26 2
## 13 2012 26 2
## 14 2023 26 2
## 15 2001 27 2
## 16 2004 28 2
## 17 2016 28 2
## 18 2013 29 2
## 19 2017 29 3
## 20 2021 29 3
## 21 2002 30 3
## 22 2007 30 3
## 23 2019 30 3
## 24 1999 33 3
## 25 2000 34 3
## 26 2010 34 3
## 27 2005 40 3cdmx_anual_ts <- ts(capital$Frecuencia,
start = c(1998,1),
frequency=12)
x <- segment(cdmx_anual_ts,
method = "pelt",
penalty = "BIC")
tidy(x)## # A tibble: 1 × 10
## region num_obs min max mean sd begin end param_mu param_sigma_hatsq
## <chr> <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 [1,28) 27 13 40 26.1 5.91 1 28 26.1 33.73.1 CDMX-Alcaldias
# libreria para el calculo de escalas
library(scales)# merge dataset + municipios
cdmx<-merge(x = cdmx,
y = AGEEMLM_INEGI,
by = "MUN_RESID",
all = TRUE)# construir tabla resumen por año
alcaldias <- with(cdmx, table(NOM_MUN))
# construir el dataframe
alcaldias<-as.data.frame(alcaldias)
names(alcaldias)[names(alcaldias) == "NOM_MUN"] <- "Alcaldia"
names(alcaldias)[names(alcaldias) == "Freq"] <- "Frecuencia"
# grafico
alcaldias_cdmx<-
ggplot(data=alcaldias,
aes(x = reorder(Alcaldia, -Frecuencia),
y = Frecuencia)) +
geom_bar(stat="identity",
fill= "#FD8D3C") +
geom_hline(yintercept = mean(alcaldias$Frecuencia),
color = "#A52A2A",
linetype = "dashed",
linewidth = 1) +
geom_text(aes(label = Frecuencia),
color="black",
size=3,
angle = 90) +
scale_x_discrete(guide = guide_axis(angle = 90)) +
labs(x = "Alcaldia",
y = "Defunciones totales",
title = "Defunciones registradas acumuladas por alcaldia",
subtitle = "1998-2024") +
theme_classic()
alcaldias_cdmx
# exportar grafico en directorio
ggsave(filename = "alcaldias_cdmx.png", plot = alcaldias_cdmx,
width = 159, height = 100, units = "mm", dpi = "retina",
path = ruta_graficos)quantile(alcaldias$Frecuencia)## 0% 25% 50% 75% 100%
## 6.00 27.75 31.50 51.75 159.00# asignar la etiqueta quantil sobre el dataframe
alcaldias_qt<-
alcaldias %>%
mutate(quantile = ntile(Frecuencia, 3)) %>%
arrange(Frecuencia)
alcaldias_qt## Alcaldia Frecuencia quantile
## 1 Milpa Alta 6 1
## 2 Cuajimalpa de Morelos 12 1
## 3 La Magdalena Contreras 15 1
## 4 Venustiano Carranza 27 1
## 5 Tláhuac 28 1
## 6 Miguel Hidalgo 29 1
## 7 Azcapotzalco 31 2
## 8 Iztacalco 31 2
## 9 Xochimilco 32 2
## 10 Benito Juárez 33 2
## 11 Cuauhtémoc 43 2
## 12 Álvaro Obregón 51 3
## 13 Coyoacán 54 3
## 14 Tlalpan 57 3
## 15 Gustavo A. Madero 98 3
## 16 Iztapalapa 159 3grafico_alcaldias_qt <- ggplot(alcaldias_qt,
aes(x = quantile,
y = Alcaldia)) +
geom_point(size = 4, color = "steelblue") +
geom_segment(aes(x = 0,
xend = quantile,
y = Alcaldia,
yend = Alcaldia),
color = "gray70") +
labs(title = "Cuantiles por Alcaldia",
x = "Cuantil",
y = "Alcaldia") +
theme_minimal() +
theme(
axis.text.y = element_text(size = 9),
panel.grid.major.y = element_blank(), # ← Elimina líneas horizontales
panel.grid.minor.y = element_blank() # ← Elimina líneas menores
)
grafico_alcaldias_qt
# exportar graficos en directorio
ggsave(filename = "cuantiles_alcaldias.png",
plot = grafico_alcaldias_qt,
width = 159,
height = 100,
units = "mm",
dpi = "retina",
path = ruta_graficos)alcaldias_qt |> count(quantile)## quantile n
## 1 1 6
## 2 2 5
## 3 3 5- Asignación de Gasto Federalizado
# Tabla resumen formato ancho
fondo <- fideicomiso[, .(
Ingresos = sum(MONTO_INGRESOS, na.rm = TRUE),
Egresos = sum(MONTO_EGRESOS, na.rm = TRUE),
Disponibilidad = sum(MONTO_DISPONIBILIDAD, na.rm = TRUE)
), by = CICLO]
# Transformar a formato largo (TABLA LARGA)
tabla_larga <- melt(fondo,
id.vars = "CICLO",
variable.name = "Tipo",
value.name = "Monto")
# Gráfico usando tabla_larga
fondo_gasto <- ggplot(tabla_larga,
aes(x = as.factor(CICLO),
y = Monto,
fill = Tipo)) +
geom_bar(stat = "identity") +
facet_wrap(~Tipo,
scales = "free_y",
ncol = 1) +
scale_y_continuous(labels = comma) +
geom_text(aes(label = scales::comma(Monto)),
vjust = 0.5,
hjust = -0.5,
color = "black",
size = 3,
angle = 90,
position = position_stack(vjust = 0.5)) +
scale_fill_manual(values = c(
"Ingresos" = "#2E7D32",
"Egresos" = "#D32F2F",
"Disponibilidad" = "#1976D2"
)) +
labs(x = "Ciclo",
y = "Miles de Millones ($)",
title = "Montos por Año") +
scale_x_discrete(guide = guide_axis(angle = 90)) +
theme_classic() +
theme(legend.position = "none")
print(fondo_gasto)
# Exportar
ggsave(filename = "fondo.png",
plot = fondo_gasto,
width = 159,
height = 100,
units = "mm",
dpi = "retina",
path = ruta_graficos)# 1. Preparar datos
fondo_largo <- melt(fondo,
id.vars = "CICLO",
variable.name = "Tipo",
value.name = "Monto")
# 2. Función CAGR simple
cagr <- function(v_ini, v_fin, años) {
ifelse(v_ini > 0, round(((v_fin/v_ini)^(1/años) - 1) * 100, 2), NA)
}
# 3. Calcular para todos los tipos
resultados <- fondo_largo[, .(
Periodo = c("2005-2021", "2021-2023", "2005-2023"),
Año_Inicial = c(2005, 2021, 2005),
Año_Final = c(2021, 2023, 2023),
CAGR = c(
cagr(Monto[CICLO == 2005], Monto[CICLO == 2021], 16),
cagr(Monto[CICLO == 2021], Monto[CICLO == 2023], 2),
cagr(Monto[CICLO == 2005], Monto[CICLO == 2023], 18)
)
), by = Tipo]
print(resultados)## Tipo Periodo Año_Inicial Año_Final CAGR
## <fctr> <char> <num> <num> <num>
## 1: Ingresos 2005-2021 2005 2021 21.54
## 2: Ingresos 2021-2023 2021 2023 -79.00
## 3: Ingresos 2005-2023 2005 2023 0.00
## 4: Egresos 2005-2021 2005 2021 37.62
## 5: Egresos 2021-2023 2021 2023 -60.75
## 6: Egresos 2005-2023 2005 2023 19.71
## 7: Disponibilidad 2005-2021 2005 2021 23.50
## 8: Disponibilidad 2021-2023 2021 2023 -37.90
## 9: Disponibilidad 2005-2023 2005 2023 14.42Analisis Univariado
A cotinuación se realizará un análisis univariado y multivariado de toda las varaibles para un mejor entendimiento del contexto sobre los casos de defunciones regitrados por LLA en la población infantil de la CDMX. Se realizará en lo máximo posible una desagregación de las variables ya que la mayor parte del conjunto de datos son de orden categorico.
- Frecuencias relativas por genero-periodo
library(dplyr)
library(ggplot2)
library(scales)# Calcular frecuencias y porcentajes
gen_edad <- cdmx[, .N, by = .(ANIO_OCUR, GENERO)]
# Calcular porcentaje por año
gen_edad[, `:=`(
pct = N / sum(N),
lbl = paste0(round(N / sum(N) * 100, 0))
), by = ANIO_OCUR]
# Gráfico
periodo_genero_cdmx <- ggplot(gen_edad,
aes(x = factor(ANIO_OCUR),
y = pct,
fill = factor(GENERO))) +
geom_bar(stat = "identity", position = "fill") +
scale_y_continuous(breaks = seq(0, 1, .5), label = percent) +
geom_text(aes(label = lbl),
size = 3,
position = position_fill(vjust = 0.5),
color = "black",
angle = 90) +
scale_x_discrete(guide = guide_axis(angle = 90)) +
scale_fill_manual(values = c("#FD8D3C", "#d16a0c")) +
labs(
y = "Mortalidad relativa",
fill = "Género",
x = "Periodo",
title = "Defunciones registradas en la CDMX por género y periodo",
subtitle = "1998-2024"
) +
theme_classic()
print(periodo_genero_cdmx)
# Exportar
ggsave(
filename = "genero_cdmx.png",
plot = periodo_genero_cdmx,
width = 159,
height = 100,
units = "mm",
dpi = "retina",
path = ruta_graficos
)- Frecuencia relativa a total periodo por género
# Calcular frecuencias y porcentajes
genero_lla <- cdmx[, .N, by = .(GENERO)]
# Calcular porcentaje
genero_lla[, `:=`(
pct = N / sum(N),
lbl = paste0(round(N / sum(N) * 100, 0), "%") # ← Agregar % aquí
)]
# Gráfico
genero_lla_cdmx <- ggplot(genero_lla, aes(x = reorder(GENERO, -pct), y = pct)) +
geom_bar(stat = "identity", fill = "#FD8D3C", color = "black") +
geom_text(aes(label = lbl), vjust = -0.25) + # ← Cambiar pctlabel por lbl
scale_y_continuous(labels = percent) +
labs(
x = "Género",
y = "Mortalidad relativa",
title = "Defunciones totales registradas en la CDMX por género",
subtitle = "1998-2024"
) +
theme_classic()
print(genero_lla_cdmx)
# Exportar
ggsave(
filename = "genero_lla_cdmx.png",
plot = genero_lla_cdmx,
width = 159,
height = 100,
units = "mm",
dpi = "retina",
path = ruta_graficos
)- Frecuencia relativa a total periodo por grado escolar
# Calcular frecuencias y porcentajes
edu_lla <- cdmx[, .N, by = .(ANIO_OCUR, ESCOLARIDA)]
# Calcular porcentaje por año
edu_lla[, `:=`(
pct = N / sum(N),
lbl = paste0(round(N / sum(N) * 100, 0))
), by = ANIO_OCUR]
# niveles de factores
periodos<- factor(cdmx$ANIO_OCUR)
periodos<-levels(periodos)
escolar<-factor(cdmx$ESCOLARIDA)
escolar<-levels(escolar)
# grafico de las relaciones
edu_lla_cdmx<-
ggplot(edu_lla,aes(x = factor(ANIO_OCUR, levels = periodos),
y = pct,
fill = factor(ESCOLARIDA,levels = escolar))) +
geom_bar(stat = "identity", position = "fill") +
scale_y_continuous(breaks = seq(0, 1,.5), label = percent) +
geom_text(aes(label = lbl),
size = 3,
position = position_fill(vjust =0.5),
color="black",
angle = 90) +
scale_x_discrete(guide = guide_axis(angle = 90)) +
scale_fill_manual(values = c("#FD8D3C", "#c68251","#ffc08c","#d16a0c"))+
labs(y = "Mortalidad relativa",
fill="Escolaridad",
x = "Periodo",
title = "Defunciones registradas en la CDMX por grado escolar",
subtitle = "1998-2024") +
theme_classic()
edu_lla_cdmx
# exportar graficos en directorio
ggsave(filename = "edu_lla_cdmx.png", plot = edu_lla_cdmx,
width = 159, height = 100, units = "mm", dpi = "retina",
path = ruta_graficos)escolaridad_resumen <- cdmx[, .(
n = .N,
pct = .N / nrow(cdmx),
lbl = percent(.N / nrow(cdmx))
), by = ESCOLARIDA]
print(escolaridad_resumen)## ESCOLARIDA n pct lbl
## <char> <int> <num> <char>
## 1: Primaria 411 0.5821530 58%
## 2: Secundaria 117 0.1657224 17%
## 3: Preescolar 79 0.1118980 11%
## 4: s/escol 99 0.1402266 14%- Frecuencia relativa a total periodo por rango de edades
# Calcular frecuencias y porcentajes
edad_lla <- cdmx[, .N, by = .(ANIO_OCUR, EDAD_AGRU)]
# Calcular porcentaje por año
edad_lla[, `:=`(
pct = N / sum(N),
lbl = paste0(round(N / sum(N) * 100, 0))
), by = ANIO_OCUR]
# niveles de factores
periodo<- factor(cdmx$ANIO_OCUR)
periodo<-levels(periodo)
edad<-factor(cdmx$EDAD_AGRU)
edad<-levels(edad)
# grafico de las relaciones
edad_lla_cdmx<-
ggplot(edad_lla,aes(x = factor(ANIO_OCUR, levels = periodo),
y = pct,
fill = factor(EDAD_AGRU,levels = edad))) +
geom_bar(stat = "identity", position = "fill") +
scale_y_continuous(breaks = seq(0, 1,.5), label = percent) +
geom_text(aes(label = lbl),
size = 3,
position = position_fill(vjust =0.5),
color="black",
angle = 90) +
scale_x_discrete(guide = guide_axis(angle = 90)) +
scale_fill_manual(values = c("#FD8D3C", "#c68251","#ffc08c","#d16a0c",
"#FFBF00", "#FA5F55", "#E3735E")) +
labs(y = "Mortalidad relativa",
fill="Edad agrupada",
x = "Periodo",
title = "Defunciones registradas en la CDMX por edad agrupada",
subtitle = "1998-2024") +
theme_classic()
edad_lla_cdmx
# exportar graficos en directorio
ggsave(filename = "edad_lla_cdmx.png", plot = edad_lla_cdmx,
width = 159, height = 100, units = "mm", dpi = "retina",
path = ruta_graficos)edad_resumen <- cdmx[, .(
n = .N,
pct = .N / nrow(cdmx),
lbl = percent(.N / nrow(cdmx))
), by = EDAD_AGRU]
print(edad_resumen)## EDAD_AGRU n pct lbl
## <char> <int> <num> <char>
## 1: 5-9 348 0.49291785 49%
## 2: 10-14 343 0.48583569 49%
## 3: < 1 15 0.02124646 2%5. Frecuencia relativa a total periodo por afiliación médica
# tasa relativa agrupado por segmento
afiliacion_lla <- cdmx %>%
group_by(ANIO_OCUR, DERECHOHAB) %>%
summarize(n = n(),.groups="drop_last") %>%
mutate(pct = n/sum(n),
lbl = scales::percent(pct, accuracy = 1))
# drop %symbol
afiliacion_lla$lbl <- gsub("%", "", afiliacion_lla$lbl)
# niveles de factor
periodo<-factor(cdmx$ANIO_OCUR)
periodo<-levels(periodo)
afiliacion<- factor(cdmx$DERECHOHAB)
afiliacion<-levels(afiliacion)
#grafico de % en orden ascendente
segmed_lla_cdmx<-
ggplot(afiliacion_lla,aes(x = factor(ANIO_OCUR, levels = periodo),
y = pct,
fill = factor(DERECHOHAB,levels = afiliacion))) +
geom_bar(stat = "identity", position = "fill") +
scale_y_continuous(breaks = seq(0, 1,.5), label = percent) +
geom_text(aes(label = lbl),
size = 3,
position = position_fill(vjust =0.5),
color="black",
angle = 90) +
scale_x_discrete(guide = guide_axis(angle = 90)) +
scale_fill_manual(values = c("#ff8000", "#c68251","#ffc08c","#d16a0c",
"#FFBF00", "#C04000", "#E3735E"))+
labs(y = "Mortalidad relativa",
fill="Seguro médico",
x = "Periodo",
title = "Defunciones registradas en la CDMX seguro médico",
subtitle = "1998-2024") +
theme_classic()
segmed_lla_cdmx
# exportar graficos en directorio
ggsave(filename = "segmed_lla_cdmx.png", plot = segmed_lla_cdmx,
width = 159, height = 100, units = "mm", dpi = "retina",
path = ruta_graficos)cdmx %>%
group_by(DERECHOHAB) %>%
summarize(n = n()) %>%
mutate(pct = n/sum(n), lbl = scales::percent(pct)) %>%
arrange(n)## # A tibble: 7 × 4
## DERECHOHAB n pct lbl
## <chr> <int> <dbl> <chr>
## 1 PEMEX 2 0.00283 0.28%
## 2 Otra 14 0.0198 1.98%
## 3 ISSFAM 16 0.0227 2.27%
## 4 ISSSTE 61 0.0864 8.64%
## 5 ASP 94 0.133 13.31%
## 6 Ninguna 216 0.306 30.59%
## 7 IMSS 303 0.429 42.92%- Frecuencia relativa a total periodo por asistencia médica
# tasa relativa agrupado por segmento
asistencia_lla <- cdmx %>%
group_by(ANIO_OCUR, ASIST_MEDI) %>%
summarize(n = n(),.groups="drop_last") %>%
mutate(pct = n/sum(n),
lbl = scales::percent(pct, accuracy =1))
# drop %symbol
asistencia_lla$lbl <- gsub("%", "", asistencia_lla$lbl)
# niveles de factor
periodo<-factor(cdmx$ANIO_OCUR)
periodo<-levels(periodo)
asistmed<- factor(cdmx$ASIST_MEDI)
asistmed<-levels(asistmed)
#grafico de % en orden ascendente
asistencia_lla_cdmx<-
ggplot(asistencia_lla,aes(x = factor(ANIO_OCUR, levels = periodo),
y = pct,
fill = factor(ASIST_MEDI,levels = asistmed))) +
geom_bar(stat = "identity", position = "fill") +
scale_y_continuous(breaks = seq(0, 1,.5), label = percent) +
geom_text(aes(label = lbl),
size = 3,
position = position_fill(vjust =0.5),
color="black",
angle = 90) +
scale_x_discrete(guide = guide_axis(angle = 90)) +
scale_fill_manual(values = c("#FD8D3C", "#C76E00", "#ffc08c"))+
labs(y = "Mortalidad relativa",
fill="Asistencia médica",
x = "periodo",
title = "Defunciones registradas en la CDMX por asistencia médica",
subtitle = "1998-2024") +
theme_classic()
asistencia_lla_cdmx
# exportar graficos en directorio
ggsave(filename = "asistencia_lla_cdmx.png", plot = asistencia_lla_cdmx,
width = 159, height = 100, units = "mm", dpi = "retina",
path = ruta_graficos)cdmx %>%
group_by(ASIST_MEDI) %>%
summarize(n = n()) %>%
mutate(pct = n/sum(n), lbl = scales::percent(pct))## # A tibble: 3 × 4
## ASIST_MEDI n pct lbl
## <chr> <int> <dbl> <chr>
## 1 No 6 0.00850 1%
## 2 Si 694 0.983 98%
## 3 n/esp 6 0.00850 1%- Frecuencia relativa a total periodo por sitio de ocurrencia
# tasa relativa agrupado por segmento
ocurrencia_lla <- cdmx %>%
group_by(SITIO_OCUR, GENERO) %>%
summarize(n = n(),.groups="drop_last") %>%
mutate(pct = n/sum(n),
lbl = scales::percent(pct, accuracy=1))
# drop %symbol
ocurrencia_lla$lbl <- gsub("%", "", ocurrencia_lla$lbl)
# el periodo deber considerarse como factor + nivel
ocurrencia<- factor(cdmx$SITIO_OCUR)
ocurrencia<-levels(ocurrencia)
#grafico de % en orden ascendente
ocurrencia_lla_cdmx<-
ggplot(ocurrencia_lla,aes(x = factor(SITIO_OCUR, levels = ocurrencia),
y = pct,
fill = factor(GENERO,
levels = c("femenino",
"masculino")))) +
geom_bar(stat = "identity", position = "fill") +
scale_y_continuous(breaks = seq(0, 1,.5), label = percent) +
geom_text(aes(label = lbl),
size = 2.5,
position = position_fill(vjust =0.5),
color="black") +
scale_x_discrete(guide = guide_axis(angle = 90)) +
scale_fill_manual(values = c("#FD8D3C", "#C76E00"))+
labs(y = "Mortalidad relativa",
fill="Género",
x = "Sitio de ocurrencia",
title = "Defunciones registradas en la CDMX por género y sitio de ocurrencia",
subtitle = "1998-2024") +
theme_classic()
ocurrencia_lla_cdmx
# exportar graficos en directorio
ggsave(filename = "ocurrencia_cdmx.png", plot = ocurrencia_lla_cdmx,
width = 159, height = 100, units = "mm", dpi = "retina",
path = ruta_graficos)cdmx %>%
group_by(SITIO_OCUR) %>%
summarize(n = n()) %>%
mutate(pct = n/sum(n),lbl = scales::percent(pct)) %>%
arrange(n)## # A tibble: 10 × 4
## SITIO_OCUR n pct lbl
## <chr> <int> <dbl> <chr>
## 1 PEMEX 1 0.00142 0.14%
## 2 SEMAR 2 0.00283 0.28%
## 3 n/esp 3 0.00425 0.42%
## 4 Otr.lug 11 0.0156 1.56%
## 5 SEDENA 12 0.0170 1.70%
## 6 ISSSTE 32 0.0453 4.53%
## 7 UMPriv 53 0.0751 7.51%
## 8 Hogar 100 0.142 14.16%
## 9 IMSS 177 0.251 25.07%
## 10 SSA 315 0.446 44.62%- Frecuencia relativa a total periodo por densidad poblacional
# tasa relativa agrupado por segmento
loc_lla <- cdmx %>%
group_by(TLOC_RESID,GENERO) %>%
summarize(n = n(),.groups= "drop_last") %>%
mutate(pct = n/sum(n),
lbl = scales::percent(pct, accuracy = 1))
# drop %symbol
loc_lla$lbl <- gsub("%", "", loc_lla$lbl)
# el periodo deber considerarse como factor + nivel
localidad<- factor(cdmx$TLOC_RESID)
localidad<-levels(localidad)
#grafico de % en orden ascendente
localidad_lla_cdmx<-
ggplot(loc_lla,aes(x = factor(TLOC_RESID, levels = localidad),
y = pct,
fill = factor(GENERO,levels = c("femenino",
"masculino")))) +
geom_bar(stat = "identity", position = "fill") +
scale_y_continuous(breaks = seq(0, 1,.5), label = percent) +
geom_text(aes(label = lbl),
size = 3,
position = position_fill(vjust =0.5),
angle = 90,
color="black") +
scale_x_discrete(guide = guide_axis(angle = 90)) +
scale_fill_manual(values = c("#FD8D3C", "#C76E00"))+
labs(y = "Mortalidad relativa",
fill="Género",
x = "Tamaño de localidad de residencia",
title = "Defunciones registradas en la CDMX por género y tamaño de localidad",
subtitle = "1998-2024") +
theme_classic()
localidad_lla_cdmx
# exportar graficos en directorio
ggsave(filename = "localidad_lla_cdmx.png", plot = localidad_lla_cdmx,
width = 159, height = 100, units = "mm", dpi = "retina",
path = ruta_graficos)cdmx %>%
group_by(TLOC_RESID) %>%
summarize(n = n()) %>%
mutate(pct = n/sum(n),lbl = scales::percent(pct)) %>%
arrange(n)## # A tibble: 11 × 4
## TLOC_RESID n pct lbl
## <chr> <int> <dbl> <chr>
## 1 2.5K < 5k 1 0.00142 0.14%
## 2 5k < 10k 1 0.00142 0.14%
## 3 15k < 20k 3 0.00425 0.42%
## 4 10k < 15k 4 0.00567 0.57%
## 5 20k < 30k 6 0.00850 0.85%
## 6 1-999 7 0.00992 0.99%
## 7 100k < 250k 25 0.0354 3.54%
## 8 1M < 2M 112 0.159 15.86%
## 9 > 2M 145 0.205 20.54%
## 10 500k < 1M 196 0.278 27.76%
## 11 250k < 500k 206 0.292 29.18%- Frecuencia relativa a total periodo por municipios de la CDMX
# tasa relativa agrupado por segmento
mun_lla <- cdmx %>%
group_by(NOM_MUN, GENERO) %>%
summarize(n = n(),.groups= "drop_last") %>%
mutate(pct = n/sum(n),
lbl = scales::percent(pct,accuracy = 1))
# drop %symbol
mun_lla$lbl <- gsub("%", "", mun_lla$lbl)
# el periodo deber considerarse como factor + nivel
municipio<- factor(cdmx$NOM_MUN)
municipio<-levels(municipio)
#grafico de % en orden ascendente
municipio_lla_cdmx<-
ggplot(mun_lla,aes(x = factor(NOM_MUN, levels = municipio),
y = pct,
fill = factor(GENERO))) +
geom_bar(stat = "identity", position = "fill") +
scale_y_continuous(breaks = seq(0, 1,.5), label = percent) +
geom_text(aes(label = lbl),
size = 3,
position = position_fill(vjust =0.5),
angle = 90,
color="black") +
scale_x_discrete(guide = guide_axis(angle = 90)) +
scale_fill_manual(values = c("#FD8D3C", "#C76E00")) +
labs(y = "Mortalidad relativa",
x = "Municipios",
fill="Género",
title = "Defunciones registradas en la CDMX por género y municipio",
subtitle = "1998-2024") +
theme_classic()
municipio_lla_cdmx
# exportar graficos en directorio
ggsave(filename = "municipio_lla_cdmx.png", plot = municipio_lla_cdmx,
width = 159, height = 100, units = "mm", dpi = "retina",
path = ruta_graficos)Analisis Multivariado
1. Test de independencia
Se procederá a realizarán los test de variables categóricas: independencia y exactitud de fisher. Esto, con el objeto de conocer posibles relaciones entre las variables y, consecuentemente, buscar literatura que atienda dichos resultados. Se considera el “genero” como la variable grupo, y salida (outcome) las otras variables categóricas
El test de independencia, de uso más frecuente para el análisis de datos categóricos, refiere al supuesto de que dos variables son independientes entre si. Transversalmente, la pregunta de investigación es: saber si los resultados son independientes del grupo al que se encontraba referido. Para ello se tienen dos hipotesis:
- HO: las variables son independientes, y son de proporciones similares entre ellas (p-value > 0.05). Es decir, la “[Variable 1] no esta asociada con la [Variable 2]”
- H1: las variables no son independientes, no son de proporciones similares entre ellas (p-value < 0.05). Es decir, la “[Variable 1] esta asociada con la [Variable 2]”
req <- substitute(require(x, character.only = TRUE))
libs<-c("sjPlot")
sapply(libs, function(x) eval(req) || {install.packages(x); eval(req)})## sjPlot
## TRUELas salidad que proporciona el summary son: χ² = Chi-cuadrado df= grados de libertad (filas x columnas) Fisher’s p = P-valor < 0.05 Cramer’s V = Fuerza de asociación | 0.10 > 0.50 |
library(sjPlot)
# Definir ruta
# Lista de combinaciones de variables
tablas <- list(
#1
list(var1 = "DERECHOHAB", var2 = "GENERO",
nombre = "Derecho habiente-Género"),
#2
list(var1 = "DERECHOHAB", var2 = "TLOC_RESID",
nombre = "Tamaño de localidad-Derecho habiente"),
#3
list(var1 = "DERECHOHAB", var2 = "ESCOLARIDA",
nombre = "Derecho habiente-Escolaridad"),
#4
list(var1 = "DERECHOHAB", var2 = "EDAD_AGRU",
nombre = "Derecho habiente-Edad"),
#5
list(var1 = "SITIO_OCUR", var2 = "GENERO",
nombre = "Sitio de ocurrencia-Género"),
#6
list(var1 = "NOM_MUN", var2 = "GENERO",
nombre = "Alcaldīa-Género")
)
# Crear todas las tablas
ver_tabla <- function(n) {
tab <- tablas[[n]]
tab_xtab(
var.row = cdmx[[tab$var1]],
var.col = cdmx[[tab$var2]],
title = paste(tab$var1, "×", tab$var2),
show.row.prc = TRUE,
show.cell.prc = TRUE,
show.col.prc = TRUE,
show.summary = TRUE,
statistics = "auto",
drop.empty = TRUE,
show.na = FALSE,
)}
ver_tabla(6)| var1]] | var2]] | Total | |
|---|---|---|---|
| femenino | masculino | ||
| Álvaro Obregón |
21 41.2 % 6.5 % 3 % |
30 58.8 % 7.8 % 4.2 % |
51 100 % 7.2 % 7.2 % |
| Azcapotzalco |
14 45.2 % 4.4 % 2 % |
17 54.8 % 4.4 % 2.4 % |
31 100 % 4.4 % 4.4 % |
| Benito Juárez |
9 27.3 % 2.8 % 1.3 % |
24 72.7 % 6.2 % 3.4 % |
33 100 % 4.7 % 4.7 % |
| Coyoacán |
28 51.9 % 8.7 % 4 % |
26 48.1 % 6.8 % 3.7 % |
54 100 % 7.6 % 7.7 % |
|
Cuajimalpa de Morelos |
5 41.7 % 1.6 % 0.7 % |
7 58.3 % 1.8 % 1 % |
12 100 % 1.7 % 1.7 % |
| Cuauhtémoc |
18 41.9 % 5.6 % 2.5 % |
25 58.1 % 6.5 % 3.5 % |
43 100 % 6.1 % 6 % |
| Gustavo A. Madero |
41 41.8 % 12.8 % 5.8 % |
57 58.2 % 14.8 % 8.1 % |
98 100 % 13.9 % 13.9 % |
| Iztacalco |
15 48.4 % 4.7 % 2.1 % |
16 51.6 % 4.2 % 2.3 % |
31 100 % 4.4 % 4.4 % |
| Iztapalapa |
75 47.2 % 23.4 % 10.6 % |
84 52.8 % 21.8 % 11.9 % |
159 100 % 22.5 % 22.5 % |
|
La Magdalena Contreras |
9 60 % 2.8 % 1.3 % |
6 40 % 1.6 % 0.8 % |
15 100 % 2.1 % 2.1 % |
| Miguel Hidalgo |
16 55.2 % 5 % 2.3 % |
13 44.8 % 3.4 % 1.8 % |
29 100 % 4.1 % 4.1 % |
| Milpa Alta |
3 50 % 0.9 % 0.4 % |
3 50 % 0.8 % 0.4 % |
6 100 % 0.8 % 0.8 % |
| Tláhuac |
10 35.7 % 3.1 % 1.4 % |
18 64.3 % 4.7 % 2.5 % |
28 100 % 4 % 3.9 % |
| Tlalpan |
34 59.6 % 10.6 % 4.8 % |
23 40.4 % 6 % 3.3 % |
57 100 % 8.1 % 8.1 % |
| Venustiano Carranza |
9 33.3 % 2.8 % 1.3 % |
18 66.7 % 4.7 % 2.5 % |
27 100 % 3.8 % 3.8 % |
| Xochimilco |
14 43.8 % 4.4 % 2 % |
18 56.2 % 4.7 % 2.5 % |
32 100 % 4.5 % 4.5 % |
| Total |
321 45.5 % 100 % 45.5 % |
385 54.5 % 100 % 54.5 % |
706 100 % 100 % 100 % |
χ2=16.551 · df=15 · Cramer’s V=0.153 · Fisher’s p=0.359 |
Extra: test de independencia & efecto tamaño
Test de independencia En caso de que se significativa la asociación, se puede explorar de forma individualizada con el siguiente codigo
library(gmodels)
#CrossTable(cdmx$DERECHOHAB, cdmx$GENERO,
# prop.t=TRUE,
# prop.r=TRUE,
# prop.c=TRUE,
# format="SPSS")Test de efecto tamaño
Comprobación de las relaciones cuando el chi-square se encuentre muy cercano al nivel de confianza, y con ello descartar todo posible sesgo en las métricas
# nivel genero-derechohabiente
#gen_med<-ftable(xtabs(~ DERECHOHAB + GENERO, data = cdmx))
# Cohen's w
#n <- sum(gen_med)
#k <- min(dim(gen_med)) - 1
#X2 <- chisq.test(gen_med,
# simulate.p.value = TRUE)$statistic
#w <- sqrt(X2 / (n * k))
#wBLOQUE 2: TS
Arreglos
Originalmente, se consideraban la desagregación a nivel mensual, pero solo es aplicable en terminos de defunciones registradas. Por lo que este apartado seguirá conservandose para su referencia y estimación del pronóstico como ejercició adicional, pero ello no significa que se incluya para la TBM.
1. Valores faltantes
Dado que se trabajará con series de tiempo, es necesario verificar las fechas faltantes a lo largo de la serie con la cual se pretende trabajar. Por ello, es necesario realizar una tabla calendario, y sobre esta vincular las fechas faltantes del conjunto de datos.
library(data.table)
# tabla calendario mensual: inicio y fin
calendario <- data.table(fecha_mes = seq.Date(as.Date("1998-01-01"),
as.Date("2024-12-01"),
by = "month")
)2. Construcción del conjunto de datos.
Construir el conjunto de observaciones
# extraer las fechas del conjunto de datos original (conteo de eventos)
cdmx_mensual<- with(cdmx, table(fecha_mes))
# resumir el conteo de eventos
cdmx_lla<-as.data.frame(cdmx_mensual)
# construir el data frame
cdmx_lla$fecha_mes<-as.Date(cdmx_lla$fecha_mes)
# columna de conteo de eventos a nivel mensual
names(cdmx_lla)[names(cdmx_lla) == "Freq"] <- "Casos"3. Dato faltantes: fechas y registros
- fechas faltantes
library(padr)
pad(cdmx_lla)## fecha_mes Casos
## 1 1998-01-01 2
## 2 1998-02-01 2
## 3 1998-03-01 2
## 4 1998-04-01 2
## 5 1998-05-01 1
## 6 1998-06-01 2
## 7 1998-07-01 1
## 8 1998-08-01 3
## 9 1998-09-01 NA
## 10 1998-10-01 2
## 11 1998-11-01 3
## 12 1998-12-01 4
## 13 1999-01-01 2
## 14 1999-02-01 1
## 15 1999-03-01 1
## 16 1999-04-01 1
## 17 1999-05-01 3
## 18 1999-06-01 2
## 19 1999-07-01 4
## 20 1999-08-01 7
## 21 1999-09-01 1
## 22 1999-10-01 NA
## 23 1999-11-01 5
## 24 1999-12-01 6
## 25 2000-01-01 3
## 26 2000-02-01 2
## 27 2000-03-01 6
## 28 2000-04-01 3
## 29 2000-05-01 4
## 30 2000-06-01 1
## 31 2000-07-01 2
## 32 2000-08-01 2
## 33 2000-09-01 1
## 34 2000-10-01 2
## 35 2000-11-01 3
## 36 2000-12-01 5
## 37 2001-01-01 3
## 38 2001-02-01 NA
## 39 2001-03-01 2
## 40 2001-04-01 2
## 41 2001-05-01 4
## 42 2001-06-01 3
## 43 2001-07-01 4
## 44 2001-08-01 NA
## 45 2001-09-01 1
## 46 2001-10-01 4
## 47 2001-11-01 3
## 48 2001-12-01 1
## 49 2002-01-01 4
## 50 2002-02-01 1
## 51 2002-03-01 1
## 52 2002-04-01 3
## 53 2002-05-01 9
## 54 2002-06-01 2
## 55 2002-07-01 2
## 56 2002-08-01 1
## 57 2002-09-01 3
## 58 2002-10-01 NA
## 59 2002-11-01 3
## 60 2002-12-01 1
## 61 2003-01-01 NA
## 62 2003-02-01 1
## 63 2003-03-01 1
## 64 2003-04-01 3
## 65 2003-05-01 1
## 66 2003-06-01 2
## 67 2003-07-01 NA
## 68 2003-08-01 1
## 69 2003-09-01 2
## 70 2003-10-01 3
## 71 2003-11-01 3
## 72 2003-12-01 1
## 73 2004-01-01 1
## 74 2004-02-01 4
## 75 2004-03-01 3
## 76 2004-04-01 2
## 77 2004-05-01 1
## 78 2004-06-01 3
## 79 2004-07-01 1
## 80 2004-08-01 3
## 81 2004-09-01 4
## 82 2004-10-01 2
## 83 2004-11-01 4
## 84 2004-12-01 NA
## 85 2005-01-01 4
## 86 2005-02-01 1
## 87 2005-03-01 3
## 88 2005-04-01 1
## 89 2005-05-01 3
## 90 2005-06-01 4
## 91 2005-07-01 7
## 92 2005-08-01 4
## 93 2005-09-01 3
## 94 2005-10-01 3
## 95 2005-11-01 4
## 96 2005-12-01 3
## 97 2006-01-01 2
## 98 2006-02-01 1
## 99 2006-03-01 NA
## 100 2006-04-01 6
## 101 2006-05-01 4
## 102 2006-06-01 4
## 103 2006-07-01 NA
## 104 2006-08-01 NA
## 105 2006-09-01 4
## 106 2006-10-01 2
## 107 2006-11-01 1
## 108 2006-12-01 2
## 109 2007-01-01 3
## 110 2007-02-01 3
## 111 2007-03-01 1
## 112 2007-04-01 2
## 113 2007-05-01 2
## 114 2007-06-01 7
## 115 2007-07-01 5
## 116 2007-08-01 3
## 117 2007-09-01 1
## 118 2007-10-01 3
## 119 2007-11-01 NA
## 120 2007-12-01 NA
## 121 2008-01-01 1
## 122 2008-02-01 4
## 123 2008-03-01 2
## 124 2008-04-01 2
## 125 2008-05-01 3
## 126 2008-06-01 1
## 127 2008-07-01 1
## 128 2008-08-01 3
## 129 2008-09-01 2
## 130 2008-10-01 NA
## 131 2008-11-01 2
## 132 2008-12-01 3
## 133 2009-01-01 3
## 134 2009-02-01 2
## 135 2009-03-01 2
## 136 2009-04-01 1
## 137 2009-05-01 4
## 138 2009-06-01 NA
## 139 2009-07-01 1
## 140 2009-08-01 1
## 141 2009-09-01 1
## 142 2009-10-01 1
## 143 2009-11-01 NA
## 144 2009-12-01 3
## 145 2010-01-01 2
## 146 2010-02-01 NA
## 147 2010-03-01 4
## 148 2010-04-01 5
## 149 2010-05-01 3
## 150 2010-06-01 4
## 151 2010-07-01 3
## 152 2010-08-01 5
## 153 2010-09-01 2
## 154 2010-10-01 1
## 155 2010-11-01 3
## 156 2010-12-01 2
## 157 2011-01-01 2
## 158 2011-02-01 4
## 159 2011-03-01 1
## 160 2011-04-01 1
## 161 2011-05-01 2
## 162 2011-06-01 NA
## 163 2011-07-01 3
## 164 2011-08-01 NA
## 165 2011-09-01 4
## 166 2011-10-01 3
## 167 2011-11-01 3
## 168 2011-12-01 1
## 169 2012-01-01 3
## 170 2012-02-01 1
## 171 2012-03-01 2
## 172 2012-04-01 2
## 173 2012-05-01 NA
## 174 2012-06-01 2
## 175 2012-07-01 2
## 176 2012-08-01 1
## 177 2012-09-01 4
## 178 2012-10-01 1
## 179 2012-11-01 4
## 180 2012-12-01 4
## 181 2013-01-01 NA
## 182 2013-02-01 4
## 183 2013-03-01 1
## 184 2013-04-01 4
## 185 2013-05-01 4
## 186 2013-06-01 NA
## 187 2013-07-01 3
## 188 2013-08-01 4
## 189 2013-09-01 3
## 190 2013-10-01 1
## 191 2013-11-01 2
## 192 2013-12-01 3
## 193 2014-01-01 4
## 194 2014-02-01 NA
## 195 2014-03-01 2
## 196 2014-04-01 1
## 197 2014-05-01 NA
## 198 2014-06-01 3
## 199 2014-07-01 4
## 200 2014-08-01 NA
## 201 2014-09-01 4
## 202 2014-10-01 3
## 203 2014-11-01 3
## 204 2014-12-01 1
## 205 2015-01-01 2
## 206 2015-02-01 1
## 207 2015-03-01 1
## 208 2015-04-01 3
## 209 2015-05-01 2
## 210 2015-06-01 3
## 211 2015-07-01 2
## 212 2015-08-01 1
## 213 2015-09-01 1
## 214 2015-10-01 NA
## 215 2015-11-01 6
## 216 2015-12-01 1
## 217 2016-01-01 3
## 218 2016-02-01 2
## 219 2016-03-01 1
## 220 2016-04-01 4
## 221 2016-05-01 3
## 222 2016-06-01 1
## 223 2016-07-01 4
## 224 2016-08-01 2
## 225 2016-09-01 NA
## 226 2016-10-01 3
## 227 2016-11-01 4
## 228 2016-12-01 1
## 229 2017-01-01 2
## 230 2017-02-01 6
## 231 2017-03-01 3
## 232 2017-04-01 6
## 233 2017-05-01 4
## 234 2017-06-01 1
## 235 2017-07-01 1
## 236 2017-08-01 NA
## 237 2017-09-01 1
## 238 2017-10-01 NA
## 239 2017-11-01 4
## 240 2017-12-01 1
## 241 2018-01-01 NA
## 242 2018-02-01 1
## 243 2018-03-01 NA
## 244 2018-04-01 1
## 245 2018-05-01 2
## 246 2018-06-01 NA
## 247 2018-07-01 1
## 248 2018-08-01 1
## 249 2018-09-01 5
## 250 2018-10-01 3
## 251 2018-11-01 2
## 252 2018-12-01 2
## 253 2019-01-01 2
## 254 2019-02-01 5
## 255 2019-03-01 3
## 256 2019-04-01 4
## 257 2019-05-01 1
## 258 2019-06-01 1
## 259 2019-07-01 6
## 260 2019-08-01 2
## 261 2019-09-01 NA
## 262 2019-10-01 NA
## 263 2019-11-01 2
## 264 2019-12-01 4
## 265 2020-01-01 1
## 266 2020-02-01 2
## 267 2020-03-01 1
## 268 2020-04-01 2
## 269 2020-05-01 NA
## 270 2020-06-01 1
## 271 2020-07-01 NA
## 272 2020-08-01 1
## 273 2020-09-01 2
## 274 2020-10-01 NA
## 275 2020-11-01 1
## 276 2020-12-01 2
## 277 2021-01-01 3
## 278 2021-02-01 2
## 279 2021-03-01 4
## 280 2021-04-01 2
## 281 2021-05-01 4
## 282 2021-06-01 2
## 283 2021-07-01 2
## 284 2021-08-01 3
## 285 2021-09-01 2
## 286 2021-10-01 2
## 287 2021-11-01 1
## 288 2021-12-01 2
## 289 2022-01-01 3
## 290 2022-02-01 2
## 291 2022-03-01 1
## 292 2022-04-01 3
## 293 2022-05-01 3
## 294 2022-06-01 2
## 295 2022-07-01 NA
## 296 2022-08-01 NA
## 297 2022-09-01 2
## 298 2022-10-01 1
## 299 2022-11-01 1
## 300 2022-12-01 2
## 301 2023-01-01 3
## 302 2023-02-01 1
## 303 2023-03-01 1
## 304 2023-04-01 1
## 305 2023-05-01 1
## 306 2023-06-01 3
## 307 2023-07-01 4
## 308 2023-08-01 4
## 309 2023-09-01 4
## 310 2023-10-01 3
## 311 2023-11-01 NA
## 312 2023-12-01 1
## 313 2024-01-01 1
## 314 2024-02-01 4
## 315 2024-03-01 1
## 316 2024-04-01 2
## 317 2024-05-01 1
## 318 2024-06-01 2
## 319 2024-07-01 3
## 320 2024-08-01 4
## 321 2024-09-01 1- valores ausentes
Se realiza una unión sobre el conjunto de datos original con respecto a la tabla calendario. Se observa que existen al menos 16 datos faltantes, que presumiblemente pueden corresponder a los registros de fallecimiento con fecha desconocida. (Previamente se habían omitido los registros con fecha desconocida)
# merge de calendario y EDR
merged_df <- merge(cdmx_lla,
calendario,
by = "fecha_mes", all = TRUE)
# conjunto de datos a trabajar
capital<-merged_df
colSums(is.na(capital))## fecha_mes Casos
## 0 433.1 La columna fecha no es un atributo per-se en el conjunto de datos a estudiar, por lo que debe de colocarse como columna indice (nivel de observaciones), y con ello preservar las variables de estudio en su espacio correspondiente
library(tibble)
# indicar la columna que pasa a index
capital <- capital |> column_to_rownames(var = 'fecha_mes')Permitiendo así, observar los casos faltantes sobre la serie de tiempo
library(imputeTS)
# momentos con datos faltantes
ggplot_na_distribution(capital)
# valores consecuentes faltantes en la serie
ggplot_na_gapsize(capital)
4. Impputación de valores faltantes
Se procede a realizar una imputación de datos, bajo el metodo de suavizamiento de Kalman bajo el modelo estructural de maxima verosimilitud.
# imputación de kalman
capital <- na_kalman(capital,
model = "StructTS")
# redondear los valores estimados
capital$Casos<-round(capital$Casos,
digits = 0)
colSums(is.na(capital))## Casos
## 0TBM
Originalmente el conjunto de datos sólo indica la frecuencia/conteo de casos por fallecimiento asociado por LLA para la CDMX. En ese sentido, para construir la tasa cruda de mortalidad asociada por dicha enfermedad, se recurre a la siguiente formula (Week, R. J. 2005. Encyclopedia of Social Measure):
TBM = [d (numero de muertes ocurridas en el año de interés) / p (población total en el año de referencia)] *1 000
La razón por la cual se trabaja con la TBM, es por las siguientes cuestiones:
- no se encuentra disponible la información de carácter público sobre el mecanismo del FPCG a nivel individuo, y por lo tanto no se puede conocer las características de las personas que recibieron tratamiento.
- la TBM agrupa por igual a los eventos ocurridos, sin importar cuestiones socio demográficas. Es una tasa general de mortalidad.
- la EDR desagrega las cuestiones socio demográficas de cada defunción y aunque es posible realizar una vinculación de la tasa de mortalidad a nivel género; esto conllevaría que ante tasas poblacionales diferentes, resultaría que el mecanismo de FPCG fue exitoso en uno u otro grupo. Por eso, la literatura recomienda que ante comparaciones de grupos (que no es el caso de la investigación), se realice métodos experimentales (grupos de control y de intervención).
2. Calculo de la TBM
La población estimada a mitad de año se encuentra expresada en millones de habitantes, por lo que conviene convertirlo en una razon de 10 mil habitantes con el objeto de que la cifra se manejable
# Convertir dataframe mensual a anual directamente
casos_lla_cdmx<-cdmx[, .N, by = .(ANIO_OCUR)]
setnames(casos_lla_cdmx, "ANIO_OCUR", "año")
setnames(casos_lla_cdmx, "N", "def_rgstr_cdmx")
# union de datasets
cdmx_lla <- merge(casos_lla_cdmx,
pob_total_cdmx,
by = "año",
all.x = TRUE)
# poblacion expresada en unidades de cien mil habiantes
cdmx_lla$pob_cien_miel <-cdmx_lla$poblacion_cdmx / 10000
cdmx_lla## Key: <año>
## año def_rgstr_cdmx poblacion_cdmx pob_cien_miel
## <int> <int> <int> <num>
## 1: 1998 24 8674660 867.4660
## 2: 1999 33 8682629 868.2629
## 3: 2000 34 8706087 870.6087
## 4: 2001 27 8758667 875.8667
## 5: 2002 30 8819709 881.9709
## 6: 2003 18 8878329 887.8329
## 7: 2004 28 8931633 893.1633
## 8: 2005 40 8977190 897.7190
## 9: 2006 26 9001741 900.1741
## 10: 2007 30 9011907 901.1907
## 11: 2008 24 9022695 902.2695
## 12: 2009 19 9029228 902.9228
## 13: 2010 34 9044151 904.4151
## 14: 2011 24 9074286 907.4286
## 15: 2012 26 9102171 910.2171
## 16: 2013 29 9121489 912.1489
## 17: 2014 25 9131261 913.1261
## 18: 2015 23 9152227 915.2227
## 19: 2016 28 9196836 919.6836
## 20: 2017 29 9246758 924.6758
## 21: 2018 18 9296268 929.6268
## 22: 2019 30 9343071 934.3071
## 23: 2020 13 9332593 933.2593
## 24: 2021 29 9272649 927.2649
## 25: 2022 20 9237644 923.7644
## 26: 2023 26 9221637 922.1637
## 27: 2024 19 9204018 920.4018
## año def_rgstr_cdmx poblacion_cdmx pob_cien_miel
## <int> <int> <int> <num>Luego, el calculo de la tasa se multiplica por 1mil para que el resultado se pueda expresarse en una relacion de miles de habitantes por cada caso de defuncion registrada mensualmente (como lo sugiere la formula de la TBM)
cdmx_tbm = copy(cdmx_lla)
# tasa calculada por formula
cdmx_tbm$TBM <- (cdmx_tbm$def_rgstr_cdmx / cdmx_tbm$pob_cien_miel) * 1000
# casos completos
cdmx_tbm<-cdmx_tbm[complete.cases(cdmx_tbm),]
cdmx_tbm## Key: <año>
## año def_rgstr_cdmx poblacion_cdmx pob_cien_miel TBM
## <int> <int> <int> <num> <num>
## 1: 1998 24 8674660 867.4660 27.66679
## 2: 1999 33 8682629 868.2629 38.00692
## 3: 2000 34 8706087 870.6087 39.05314
## 4: 2001 27 8758667 875.8667 30.82661
## 5: 2002 30 8819709 881.9709 34.01473
## 6: 2003 18 8878329 887.8329 20.27409
## 7: 2004 28 8931633 893.1633 31.34925
## 8: 2005 40 8977190 897.7190 44.55737
## 9: 2006 26 9001741 900.1741 28.88330
## 10: 2007 30 9011907 901.1907 33.28929
## 11: 2008 24 9022695 902.2695 26.59959
## 12: 2009 19 9029228 902.9228 21.04277
## 13: 2010 34 9044151 904.4151 37.59336
## 14: 2011 24 9074286 907.4286 26.44836
## 15: 2012 26 9102171 910.2171 28.56461
## 16: 2013 29 9121489 912.1489 31.79305
## 17: 2014 25 9131261 913.1261 27.37847
## 18: 2015 23 9152227 915.2227 25.13050
## 19: 2016 28 9196836 919.6836 30.44525
## 20: 2017 29 9246758 924.6758 31.36234
## 21: 2018 18 9296268 929.6268 19.36261
## 22: 2019 30 9343071 934.3071 32.10936
## 23: 2020 13 9332593 933.2593 13.92968
## 24: 2021 29 9272649 927.2649 31.27477
## 25: 2022 20 9237644 923.7644 21.65054
## 26: 2023 26 9221637 922.1637 28.19456
## 27: 2024 19 9204018 920.4018 20.64316
## año def_rgstr_cdmx poblacion_cdmx pob_cien_miel TBM
## <int> <int> <int> <num> <num>Seleccionamos solo la variable a trabajar en las series de tiempo
# selección por posición df
cdmx_tbm<-copy(cdmx_tbm[,c(1,5)])
# construccion de ts del conjunto de datos original
cdmx_tbm_ts <- ts(cdmx_tbm$TBM,
start = c(1998,1),
frequency=1)Validación de la serie
La serie de la TBM debe de reunir requisitos mínimos a fin de que los resultados sean fiables y optimos. Guia de validación
1. Linealidad
1.1 Test de linealidad
Una serie de tiempo lineal significaría que el conjunto de datos muestra una estructura regular en el tiempo; una serie de tiempo no-lineal se referirá a que los patrones de cambio desvían una linea de tendencia. Cada uno por su parte, contempla estrategias diferenciadas para su análisis.
El test de linealidad vía Nomralización Bispectrum reproduce un gráfico de probabilidades donde se asume que la serie se encuentra bajo un proceso linea con innovaciones de iid. Donde:
- H0: el bispectrum es constante sobre toda las frecuencias, ergo, los datos provienen de procesos lineal con innovaciones de i.i.d
- H1: contrario a H0
library(astsa)
# test de linearidad
test.linear(cdmx_tbm_ts,
color = TRUE,
detrend = FALSE,
main = NULL)
> Conclusiones: > - la serie de tiempo es lineal a lo largo del
horizonte de tiempo
1.2. Tests de no-linealidad
Los tests de no-linealidad asumen que la serie de tiempo no es lineal, cuando:
+H0: serie lineal (p-value > 0.05) +H1: serie no lineal (p-value < 0.05)
library(nonlinearTseries)
# función de no-linearidad
nonlinearityTest(cdmx_tbm_ts, verbose = TRUE)## ** Teraesvirta's neural network test **
## Null hypothesis: Linearity in "mean"
## X-squared = 0.9130956 df = 2 p-value = 0.6334667
##
## ** White neural network test **
## Null hypothesis: Linearity in "mean"
## X-squared = 1.350713 df = 2 p-value = 0.508975
##
## ** Keenan's one-degree test for nonlinearity **
## Null hypothesis: The time series follows some AR process
## F-stat = 0.03535562 p-value = 0.8523703
##
## ** McLeod-Li test **
## Null hypothesis: The time series follows some ARIMA process
## Maximum p-value = 0.9077316
##
## ** Tsay's Test for nonlinearity **
## Null hypothesis: The time series follows some AR process
## F-stat = 0.6596094 p-value = 0.4250245
##
## ** Likelihood ratio test for threshold nonlinearity **
## Null hypothesis: The time series follows some AR process
## Alternativce hypothesis: The time series follows some TAR process
## X-squared = 0.4518465 p-value = 0.3219993## $Terasvirta
##
## Teraesvirta Neural Network Test
##
## data: ts(time.series)
## X-squared = 0.9131, df = 2, p-value = 0.6335
##
##
## $White
##
## White Neural Network Test
##
## data: ts(time.series)
## X-squared = 1.3507, df = 2, p-value = 0.509
##
##
## $Keenan
## $Keenan$test.stat
## [1] 0.03535562
##
## $Keenan$df1
## [1] 1
##
## $Keenan$df2
## [1] 25
##
## $Keenan$p.value
## [1] 0.8523703
##
## $Keenan$order
## [1] 0
##
##
## $McLeodLi
## $McLeodLi$p.values
## [1] 0.6799542 0.7536351 0.8694800 0.9077316 0.3947171 0.2515985 0.3151309
## [8] 0.4127917 0.5125134 0.5347459 0.5357992 0.6208451 0.6183194 0.6929277
##
##
## $Tsay
## $Tsay$test.stat
## [1] 0.6596094
##
## $Tsay$p.value
## [1] 0.4250245
##
## $Tsay$order
## [1] 1
##
##
## $TarTest
## $TarTest$percentiles
## [1] 23.07692 76.92308
##
## $TarTest$test.statistic
## [1] 0.4518465
##
## $TarTest$p.value
## [1] 0.3219993Conclusiones: - los tests de no-linearidad confrman que la serie de tiempo es lineal
1.3. Datos subrogados
En esencia, el test genera severos datos sustitutos para luego después ser comparados con valores de una estadística discriminante entre los datos originales y subrogados
- H0: la series es un proceso lineal gausiano
- H1: la serie no es un proceso lineal gausiano
# función de datos subrrogados
surrogateTest(cdmx_tbm_ts,
significance = 0.05,
one.sided = F,
FUN = timeAsymmetry,
do.plot= TRUE,
verbose = TRUE)
## $surrogates.statistics
## [1] 69.364546 -18.160322 -543.418869 160.941848 213.342197 255.114640
## [7] 213.074381 -547.514866 238.135535 -305.033398 -138.630958 48.219895
## [13] 74.061969 -459.526583 7.757824 26.712998 154.441373 -358.650245
## [19] 129.488434 -92.253066 127.015680 76.767753 -129.313834 26.515849
## [25] -313.612310 204.100567 476.325665 187.802784 -168.756814 -28.201281
## [31] -226.598947 -24.288865 -145.167922 -280.645541 173.435421 -50.616699
## [37] 459.716160 396.290447 196.081483
##
## $data.statistic
## [1] -148.1081
##
## attr(,"class")
## [1] "surrogateTest"Conclusion: - Data comes from a linear stochastic process
2. Descomposición
- Descomposición clásica.
Tanto una descomposición de orden aditivo o multiplicativo, se realiza mediante medias móviles para cada unidad de tiempo sobre todo los periodos. Esto conlleva que se omitan unidades de tiempo tanto al inicio y final del periodo por la venta de tiempo simétrica con igualdad de pesos.
1.1 Funciones de componentes Aditivas/Multiplicativas
library(TSstudio)
#descomposición multiplicatica/aditiva
if(FALSE) {
ts_decompose(cdmx_tbm_ts, type = "both", showline = TRUE)
# componentes de la ts en data.frame
additive_compnt_cdmx<-decompose(cdmx_tbm_ts,
type = c("additive"),
filter = NULL)
# componentes de la ts en data.frame
multiplicative_compnt_cdmx<-decompose(cdmx_tbm_ts,
type = c("multiplicative"),
filter = NULL)
}- Descomposición por método STL
Para efectos de la investigación, utilizaremos el metodo conocido como “Seasonal and Trend decomposition using Loess”. Las ventajas de esta descomposición son las siguientes (Hyndman, R.J., & Athanasopoulos, G., 2021):
- maneja cualquier tipo de estacionalidad
- el componente estacional puede definirse por el usuario
- la uniformidad de la tendencia-ciclo puede ser controlada por el usuario
- los outliers usualmente no afectan la tendencia-ciclo y la estacionalidad.
- no se omiten unidades de tiempo al llevar a cabo las estimaciones de descomposción.
# Verificar que tienes al menos 2 periodos
n_periodos <- length(cdmx_tbm_ts) / frequency(cdmx_tbm_ts)
cat("Número de periodos:", n_periodos, "\n")## Número de periodos: 27# funcion stl
if(FALSE) {
stl_decomp <- stl(cdmx_tbm_ts,
s.window = 7, # Ventana móvil
t.window = 13, # Suavizado de tendencia
robust = TRUE)
# informe
summary(stl_decomp)
}if(FALSE) {
plot(stl_decomp)
stl_decomp$time.series
}Conclusiones:
- la serie de tiempo no exhibe una tendencia definida, es aleatoria su comportamiento.
- las variaciones atribuidas al componente estocástico de residual/reminder son considerablemente altas.
3. Distribución de los datos
hist(cdmx_tbm_ts)
Con el fin de evaluar el grado de ajuste del método
(normalidad), se realizará la evaluación sobre la serie irregular (que
elimina la tendencia y estacionalidad) resultante de la descomposición
del método STL con el objeto de saber si esta serie es realmente
aleatoria
Distribución de normalidad:
- H0: distribución normal (p-vlue > 0.05)
- H1: distribución no normal (p-value < 0.05)
- Test’s
library(nortest)
# Kolmogorov-Smirnov test
ks.test(cdmx_tbm_ts, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: cdmx_tbm_ts
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided# Cramer-von Mises normality test
cvm.test(cdmx_tbm_ts)##
## Cramer-von Mises normality test
##
## data: cdmx_tbm_ts
## W = 0.050952, p-value = 0.484# Jarque Bera Test
library(tseries)
jarque.bera.test(cdmx_tbm_ts)##
## Jarque Bera Test
##
## data: cdmx_tbm_ts
## X-squared = 0.0011315, df = 2, p-value = 0.9994# Anderson-Darling normality test
ad.test(cdmx_tbm_ts)##
## Anderson-Darling normality test
##
## data: cdmx_tbm_ts
## A = 0.29918, p-value = 0.5597# Shapiro-Francia normality test
sf.test(cdmx_tbm_ts)##
## Shapiro-Francia normality test
##
## data: cdmx_tbm_ts
## W = 0.97595, p-value = 0.6704- Q-Q Plot
qqnorm(cdmx_tbm_ts)
qqline(cdmx_tbm_ts, col = "red")
- Conocer el tipo de distribución de la serie
library(fitdistrplus)
descdist(cdmx_tbm$TBM, discrete = FALSE, boot = 500)
## summary statistics
## ------
## min: 13.92968 max: 44.55737
## median: 28.8833
## mean: 28.94239
## estimated sd: 6.783751
## estimated skewness: 0.01265694
## estimated kurtosis: 3.285313Conclusiones: + se rechaza H0 en todo los tests de normalidad + qqplot muestra una distribución no aceptable de la normalidad de los datos
Descomposición temporal de la serie
Corrección en la distribución.
Antes de continuar, es necesario realizar una corrección sobre la serie de tiempo con el objeto de obtener mejores estimaciones en las evaluaciones previas a la construcción del pronóstico y en el diseño del modelo ITS.
Para ello, se recurre a una serie de transformaciones mediante estadística de normalidad a fin de evaluar la mejor transformación posible en la serie. Lista de transformaciones Best normalizer package Normalizing Transformation Functions
library(bestNormalize)
# función
BN_obj<-bestNormalize(cdmx_tbm_ts,
loo = TRUE)
BN_obj## Best Normalizing transformation with 27 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 1.2963
## - Box-Cox: 1.1778
## - Center+scale: 1.1778
## - Double Reversed Log_b(x+a): 0.7077
## - Exp(x): 34.7185
## - Log_b(x+a): 1.2963
## - orderNorm (ORQ): 0.1111
## - sqrt(x + a): 1.4148
## - Yeo-Johnson: 1.1778
## Estimation method: Out-of-sample via leave-one-out CV
##
## Based off these, bestNormalize chose:
## orderNorm Transformation with 27 nonmissing obs and no ties
## - Original quantiles:
## 0% 25% 50% 75% 100%
## 13.930 25.789 28.883 31.951 44.557# histograma
hist(BN_obj$x.t,main="Transformación OQR",xlab="",ylab=" ")
De acuerdo con las estadísticas de normalidad, la mejor transformación
de los datos originales corresponde del tipo orderNorm (ORQ).
Transformaciones y back-transformation de la TBM
# Perform transformation
cdmx_tbm_ts_oqr <- predict(BN_obj)
head(cdmx_tbm_ts_oqr)## Time Series:
## Start = 1998
## End = 2003
## Frequency = 1
## [1] -0.2822161 1.3249577 1.5932188 0.1867561 0.9674216 -1.3249577# Perform BACK-TRANSFORMATION
head(predict(BN_obj, newdata = cdmx_tbm_ts_oqr, inverse = TRUE))## Time Series:
## Start = 1998
## End = 2003
## Frequency = 1
## [1] 27.66679 38.00692 39.05314 30.82661 34.01473 20.27409Por procedimiento, se valida la serie en términos de distribución
Descomposición por STL
# funcion STL
if(FALSE){
stl_decomp_oqr<- stl(cdmx_tbm_ts_oqr,
s.window = "periodic",
robust = TRUE)
plot(stl_decomp_oqr)
}Test de distribución
# Kolmogorov-Smirnov test
ks.test(cdmx_tbm_ts_oqr, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: cdmx_tbm_ts_oqr
## D = 0.018519, p-value = 1
## alternative hypothesis: two-sided# Cramer-von Mises normality test
cvm.test(cdmx_tbm_ts_oqr)##
## Cramer-von Mises normality test
##
## data: cdmx_tbm_ts_oqr
## W = 0.003104, p-value = 1# Jarque Bera Test
library(tseries)
jarque.bera.test(cdmx_tbm_ts_oqr)##
## Jarque Bera Test
##
## data: cdmx_tbm_ts_oqr
## X-squared = 0.18758, df = 2, p-value = 0.9105# Anderson-Darling normality test
ad.test(cdmx_tbm_ts_oqr)##
## Anderson-Darling normality test
##
## data: cdmx_tbm_ts_oqr
## A = 0.03465, p-value = 1# Shapiro-Francia normality test
sf.test(cdmx_tbm_ts_oqr)##
## Shapiro-Francia normality test
##
## data: cdmx_tbm_ts_oqr
## W = 0.99982, p-value = 1QQ-plot
qqnorm(cdmx_tbm_ts_oqr)
qqline(cdmx_tbm_ts_oqr, col = "red")
Confirmación sobre el tipo de distribución
library(healthyR.ts)
if(FALSE){
# time serie transformed to data frame
x<- ts_to_tbl(cdmx_tbm_ts_oqr)
# función de distribución
descdist(x$value,
discrete = FALSE,
boot = 500)
}Conclusiones:
- se validan H0 en todo los tests de normalidad
- qqplot muestra una normalidad excelente de los datos transformados Una vez realizado la transformación de la serie de tiempo por el método de OQR, ésta pasará a ser la nueva serie de tiempo a trabajar
4. Outliers
library(tidyverse)
library(anomalize)
library(fpp2)Para la detección de outliers potenciales se estima mediante dos vias: supsum para series no estacionales, y una descomposición periódica de stl para series estacionales.En caso de hallarse ouliers, se deberá de utilizar la función tsclean()
# función de outliers
tsoutliers(cdmx_tbm_ts_oqr)## $index
## integer(0)
##
## $replacements
## numeric(0)Conclusión:la serie de tiempo no presenta outliers.
5. Tendencia
Con el objeto de conocer si la serie de tiempo presenta una tendencia a lo largo de los años, se recurrirá a una serie de tests no-parametricos.
Test for trends in time series
- Student test
- H0: no existe tendencia (p-value > 0.05)
- H1: tendencia lineal (p-value < 0.05)
library(funtimes)
# Student's t-test
notrend_test(cdmx_tbm_ts_oqr,
B = 1000,
ic = "BIC",
test = "t")##
## Sieve-bootstrap Student's t-test for a linear trend
##
## data: cdmx_tbm_ts_oqr
## Student's t value = -2.5375, p-value = 0.017
## alternative hypothesis: linear trend.
## sample estimates:
## $AR_order
## [1] 0
##
## $AR_coefficients
## numeric(0)- Mann-Kendall
- HO: No hay tendencia monótona (p-value > 0.05)
- H1: Hay tendenica monótona (p-value < 0.05)
# Mann–Kendall test
notrend_test(cdmx_tbm_ts_oqr,
B = 1000,
ic = "BIC",
test = "MK")##
## Sieve-bootstrap Mann--Kendall's trend test
##
## data: cdmx_tbm_ts_oqr
## Mann--Kendall's tau = -0.30484, p-value = 0.025
## alternative hypothesis: monotonic trend.
## sample estimates:
## $AR_order
## [1] 0
##
## $AR_coefficients
## numeric(0)- WAVK
- HO: Existe tendencia monótona (p-value > 0.05)
- H1: No tendencia monótona (p-value < 0.05)
# WAVK test
notrend_test(cdmx_tbm_ts_oqr,
B = 1000,
ic = "BIC",
test = "WAVK")##
## Sieve-bootstrap WAVK trend test
##
## data: cdmx_tbm_ts_oqr
## WAVK test statistic = -0.23766, moving window = 3, p-value = 0.838
## alternative hypothesis: (non-)monotonic trend.
## sample estimates:
## $AR_order
## [1] 0
##
## $AR_coefficients
## numeric(0)Conclusiones: + la serie de tiempo no presenta una tendencia lineal (test de Student) + la serie de tiempo presentaría una tendencia monotonica (Mann-Kendall, WAVK)
6. Estacionalidad
library(TSstudio)
library(fpp2)
library (seastests)- Test de Ollech y Webel:
Este test combina los resultados del QS-test y el kw-test, ambos calculados por el residual no-estacional del modelo ARIMA.
El resultado esta en función de p-value en QS-test por debajo de 0.01 o el p-value en kw-test por debajo 0.002, por lo que el WO-test indicará si la serie de tiempo es estacional.
# test Ollech-Welch
if(FALSE){
summary(combined_test(cdmx_tbm_ts_oqr,
freq = NA))
}- Test seasonal
La evaluación de la temporalidad de una serie puede ser mediante la combinación del QS test, Friedman test,Kruskall-Wallis test, F-test, Welch test, dando como resultado una evaluación Boleana.
# función estacional
if(FALSE){
isSeasonal(cdmx_tbm_ts_oqr,
test = "combined",
freq = 12)
}#nsdiffs(cdmx_tbm_ts_oqr)- Peridograma
Una opción de análisis espectral para el análisis de series de tiempo puede ser el periodograma acumulativo. Donde, el eje de las abscisas representa la frecuencia de ciclos de ocurrencia, y donde los intervalos de confianza delimitan tanto la freuencia como los periodos de ocurrencia (eje de las ordenadas)
require(graphics)
# cpgram
cpgram(cdmx_tbm_ts_oqr,
taper = 0.1,
main = "Cumulative periodogram")
# auto correlacion serial
acf(cdmx_tbm_ts_oqr)
- Canova-Hansen test
Para conocer patrones temporales de la serie se recurre al test de Canova-Hansen donde
- H0: existen patrones temporales estables (p-value < 0.05)
- H1: no existen patrones temporales estables (p-value > 0.05)
if(FALSE){
library(uroot)
# funcion Canova-Hansen
ch<-ch.test(cdmx_tbm_ts_oqr, type = "dummy" )
summary(ch$fitted.model)
}Conclusión: no existe evidencia para considerar que la serie de tiempo presente estacionalidades.
6. Crecimientos interanuales
Comportamiento mensuales en terminos absolutos de casos registrados sobre LLA
Las tasas de crecimiento resultantes pueden ayudar a descubrir tendencias que podrían no ser evidentes en los datos originales. Usar una tasa de crecimiento logarítmicamente diferenciada es especialmente útil para capturar el cambio porcentual, lo que facilita la interpretación de los datos
if(FALSE){
# subseries estacionales
ggmonthplot(cdmx_ts_oqr,
main = "Crecimientos temporales por defunciones de LLA en CDMX",
ylab = "TBM") +
theme_classic()
library(healthyR.ts)
# % crecimientos mensuales
ts_growth_rate_vec(cdmx_ts,.log_diff = TRUE)|> head(12)
# gráfico de crecimientos mensuales
plot(ts(ts_growth_rate_vec(cdmx_ts, .log_diff = TRUE)))
}# boxplot con ciclos de serie
boxplot(cdmx_tbm_ts_oqr ~ cycle(cdmx_tbm_ts_oqr))
Conclusiones: el crecimiento de la TBM a lo largo de los años muestra cierta tedencia bajista, por lo que se presenta una tdencia acumulada hacia la baja de forma acelarada
7. Estacionariedad
- Phillip-Perron
- H0: Existe raíz unitaria, la serie no es estacionaria (p-value > 0.05)
- H1: No existe raíz unitaria, la serie es estacionaria (p-value < 0.05)
pp.test(cdmx_tbm_ts_oqr)##
## Phillips-Perron Unit Root Test
##
## data: cdmx_tbm_ts_oqr
## Dickey-Fuller Z(alpha) = -35.095, Truncation lag parameter = 2, p-value
## = 0.01
## alternative hypothesis: stationary- Augmented Dickey-Fuller:
- H0: Existe raíz unitaria, la serie no es estacionaria (p-value > 0.05)
- H1: No existe raíz unitaria, la serie es estacionaria (p-value < 0.05)
library(tseries)
adf.test(cdmx_tbm_ts_oqr)##
## Augmented Dickey-Fuller Test
##
## data: cdmx_tbm_ts_oqr
## Dickey-Fuller = -4.0844, Lag order = 2, p-value = 0.02007
## alternative hypothesis: stationary- Kwiatkowski-Phillips-Schmidt-Shin (KPSS)
El test Kwiatkowski-Phillips-Schmidt-Shin (KPSS) determina si una serie de tiempo es estacionaria alrededor de una tendencia determinista frente a la raíz unitaria.
- H0: la serie es estacionaria y no tiene una raíz unitaria (tendencia estocástica) (p-alfa > 0.05)
- H1: la serie no es estacionaria y presenta una raíz unitaria (tendencia no estocástica) (p-alfa < 0.05)
library(tseries)
kpss.test(cdmx_tbm_ts_oqr, null="Trend" )##
## KPSS Test for Trend Stationarity
##
## data: cdmx_tbm_ts_oqr
## KPSS Trend = 0.043453, Truncation lag parameter = 2, p-value = 0.1Conclusiones: la serie de tiempo es estacionaria conforme a los Test de prueba, por lo que los posibles resultados pueden ser del siguiente orden:
- ambos tests concluyen que la serie de tiempo es no estacionaria,
- ambos tests contluyen que la serie de tiempo es estacionaria
- si KPSS indica estacionariedad y ADF indica no estacionariedad, entonces se debe eliminar la tendencia porque la serie es estacionaria en tendencia
- si KPSS indica no estacionariedad y ADF indica estacionariedad, entonces se debe diferenciar la serie porque esta presenta estacionarieda en diferencia.
https://www.statsmodels.org/dev/examples/notebooks/generated/stationarity_detrending_adf_kpss.html
Correccion de estacionariedad (extra, no es necesario)
En caso de que la serie de tiempo no fuera estacionaria, se puede realizar un procesamiento rápido de auto-estacionarizacion con la paqueteria auto_stationarize()
library(healthyR.ts)
head(auto_stationarize(cdmx_tbm_ts_oqr))## The time series is already stationary via ts_adf_test().## Time Series:
## Start = 1998
## End = 2003
## Frequency = 1
## [1] -0.2822161 1.3249577 1.5932188 0.1867561 0.9674216 -1.3249577Diferencias requeridas para estacionarizar la serie
Numero de diferencias requeridos para convertir la serie en estacionaria, tanto a nivel como a tendencia (en caso de se requerido)
ndiffs(cdmx_tbm_ts_oqr,
alpha = 0.05,
test = c("kpss", "adf", "pp"),
type = c("level", "trend"),
max.d = 2)## [1] 1Si el test ADF y Phillip-Perron sugieren una estacionariedad, en tanto que el test KPSS sugiere una no-estacionariedad, entonces para reconciliar las hipotesis debemos de hacer una serie de ajuste en la serie de tiempo sobre una primera diferencia.
# función de 1st diferencia
diff_cdmx_ts_oqr<-diff(cdmx_tbm_ts_oqr,
differences = 2,
lag = 1)Tests de validacion
# Dickie-Fuller
adf.test(diff_cdmx_ts_oqr)##
## Augmented Dickey-Fuller Test
##
## data: diff_cdmx_ts_oqr
## Dickey-Fuller = -4.6446, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationary# Phillip-Perron
pp.test(diff_cdmx_ts_oqr)##
## Phillips-Perron Unit Root Test
##
## data: diff_cdmx_ts_oqr
## Dickey-Fuller Z(alpha) = -38.595, Truncation lag parameter = 2, p-value
## = 0.01
## alternative hypothesis: stationary# KPSS
library(urca)
summary(ur.kpss(diff_cdmx_ts_oqr))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.0733
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.7398. Autocorrelación
- ACF & PACF
ACF: patrón estacional sobre la serie (MA) PACF: patrón autorregresivo de la serie (AR)
Interpreting ACF and PACF Plots for Time Series Forecasting
library(forecast)
# correlograma
ggtsdisplay(cdmx_tbm_ts_oqr,
main = "ACF & PACF para la serie de tiempo CDMX_TS",
smooth = TRUE,
theme = theme_minimal())
# coeficientes de autocorrelacion
acf(cdmx_tbm_ts_oqr,
pl=F)##
## Autocorrelations of series 'cdmx_tbm_ts_oqr', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 -0.165 0.185 0.073 -0.058 0.238 0.263 -0.084 0.024 0.021 -0.139
## 11 12 13 14
## 0.090 0.041 -0.141 -0.003auto.arima(cdmx_tbm_ts_oqr)## Series: cdmx_tbm_ts_oqr
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -1.6398 0.8273
## s.e. 0.2349 0.2401
##
## sigma^2 = 0.7551: log likelihood = -34.17
## AIC=74.33 AICc=75.42 BIC=78.11Conclusiones:
- la serie presenta una patron autorregresivo de orden 1 (AR) y un patron esatcional de 2 (MA)
- la serie no muestra una autocorrelación serial significativa, dado que presenta un ACF, PACF y frecuencia de 0
9. Análisis de frecuencias (ciclos)
Frecuencia dominante en la serie de tiempo:
- datos estacionales: regresa los periodos estacionales
- datos ciclicos: regresa el promedio del largo del ciclo
findfrequency(cdmx_tbm_ts_oqr)## [1] 1Frecuencia no parametrica
library(seewave)
# 1st option
spec_np<- spectrum(cdmx_tbm_ts_oqr,
spans=c(3,5,3),
taper = 0.1,
pad = 0,
fast = TRUE,
demean = FALSE,
detrend = TRUE,
main = "Smoothed periodogram")
# 2nd option
spec.pgram(cdmx_tbm_ts_oqr)
Interpretación del peridograma
- Axes:
Horizontal Axis (Frequency): This axis represents the frequencies present in the time series data. It typically ranges from 0 to 0.5 cycles per unit of time (often normalized, where 0.5 is the Nyquist frequency, representing the highest frequency that can be resolved given the sampling rate). Lower frequencies correspond to longer cycles or trends, while higher frequencies correspond to shorter, more rapid cycles.
Vertical Axis (Power/Spectral Density): This axis represents the power or intensity of each frequency component. Higher values on this axis indicate a stronger presence of that particular frequency in the original time series. This is often displayed on a logarithmic scale (e.g., in decibels) to better visualize variations across a wide range of power values.
Frecuencia parametrica basada en modelos AR [AR(1)][https://sjmiller8182.github.io/LearningTS/build/autoregressive_models.html]
spec_pm<-spec.ar(cdmx_tbm_ts_oqr)
plot(spec_pm$freq, spec_pm$spec,
type = "l", log = "y",
ylab = "spectrum", xlab = "frequency",
bty = "l")
lines(spec_np$freq, spec_np$spec, lty = 2)
legend("topright",
c("parametric", "nonparametric"),
lty = 1:2, bty = "n")
10. Ruido blanco
Se espera que cada termino de autocorrelación sea cercano a cero. Para un white noise, se espera que el 95% de los spikes en el ACF se encuenren dentro del limite de confianza del ACF de ± 1.96; pero si uno o más spikes se encuentran fuera de los limites, entonces es probable que la serie no sea ruido blanco. Para ello, los supuestos son los siguientes:
- media cero
- desviación estandar es constante en el tiempo
- correlación entre lags es igual a cero
- Test de Box-Pierce o Ljung-Box:
- H0: las observaciones son distribuidos independientemente (p-value > 0.5)
- H1: los observaciones no son distribuidos independientemente y exhiben una autocorrelacion serial significativa (p-value < 0.05)
library(zoo)
# Test portmanteau
Box.test(cdmx_tbm_ts_oqr,
type = "Ljung-Box")##
## Box-Ljung test
##
## data: cdmx_tbm_ts_oqr
## X-squared = 0.8229, df = 1, p-value = 0.3643Box.test(cdmx_tbm_ts_oqr,
type = "Box-Pierce")##
## Box-Pierce test
##
## data: cdmx_tbm_ts_oqr
## X-squared = 0.73777, df = 1, p-value = 0.3904- Bartlett
- HO: The variances of the groups (or the time series segments) are equal (homoscedasticity) where p-value > 0.05
- H1: At least one group’s variance is significantly different from the others (heteroscedasticity) where p-value < 0.05
library(hwwntest)
bartlettB.test(cdmx_tbm_ts_oqr, plot.it = F)##
## Bartlett B Test for white noise
##
## data:
## = 0.64605, p-value = 0.7981Retorno Verdadero si el objeto testeado contiene valores numericos con no variabilidad.
is.constant(cdmx_tbm_ts_oqr)## [1] FALSEConclusiones: la serie de tiempo es ruido blanco, y se confirma complementariamente con los correlogramas que confirman que la serie no presenta autocorrelacion.
11. Cambio estructural
Se ponen a prueba varios métodos analíticos de cambio estructural con el objetivo de evaluar su perfomance.
comparasion of changue point detection methods
De acuerdo con la información pública existente, se conoce que el FPCG entro en operaciones en 2004, pero de acuerdo con la programación presupuestaria de la SHCP para el fideicomico del FPCG ocurrió en 2005. Por lo tanto, un potencial cambio estructural de la TBM se econtraría en ese intervalo de tiempo.
Nota: “If the purpose is to look for turning points in a series, and interpret any changes in direction, then it is better to use the trend-cycle component rather than the seasonally adjusted data” (Hyndman, R.J., & Athanasopoulos, G., 2018)
strucchange: changes in (time series) regression model
Fstats(): test from the F test framework, which assume that there is one (un-known) breakpoint
Test de hipotesis nula donde la TBM permanece constante sobre los años
How to detect and quantify a structural break in time-series
library(strucchange)
# funcion
fs.cdmx<- Fstats(cdmx_tbm_ts_oqr~ 1)
sctest(fs.cdmx)##
## supF test
##
## data: fs.cdmx
## sup.F = 6.4266, p-value = 0.1343# punto de quiebre
breakpoints(fs.cdmx)##
## Optimal 2-segment partition:
##
## Call:
## breakpoints.Fstats(obj = fs.cdmx)
##
## Breakpoints at observation number:
## 10
##
## Corresponding to breakdates:
## 2007# gráfico
plot(fs.cdmx,
pval = TRUE,
alpha = 0.01)
lines(breakpoints(fs.cdmx))
library(strucchange)
# Test supF
fs.cdmx <- Fstats(cdmx_tbm_ts_oqr ~ 1)
test_result <- sctest(fs.cdmx)
cat("\n=== RESULTADOS DEL TEST ===\n\n")##
## === RESULTADOS DEL TEST ===cat("Estadístico sup.F:", test_result$statistic, "\n")## Estadístico sup.F: 6.42658cat("P-valor:", test_result$p.value, "\n\n")## P-valor: 0.134294# Interpretación del p-valor
if (test_result$p.value < 0.01) {
cat("✅ SIGNIFICATIVO al 1%: Fuerte evidencia de breakpoint\n")
} else if (test_result$p.value < 0.05) {
cat("✅ SIGNIFICATIVO al 5%: Evidencia de breakpoint\n")
} else if (test_result$p.value < 0.10) {
cat("⚠️ MARGINALMENTE SIGNIFICATIVO al 10%: Señal débil\n")
} else {
cat("❌ NO SIGNIFICATIVO: No hay evidencia suficiente de breakpoint\n")
cat(" (Aunque p=", round(test_result$p.value, 3), " sugiere una tendencia)\n")
}## ❌ NO SIGNIFICATIVO: No hay evidencia suficiente de breakpoint
## (Aunque p= 0.134 sugiere una tendencia)# Breakpoints
bp <- breakpoints(fs.cdmx)
cat("\n=== BREAKPOINT DETECTADO ===\n\n")##
## === BREAKPOINT DETECTADO ===cat("Posición:", bp$breakpoints, "\n")## Posición: 10# Convertir a fecha
bp_año <- time(cdmx_tbm_ts_oqr)[bp$breakpoints]
cat("Año:", bp_año, "\n\n")## Año: 2007# Estadísticas por segmento
cat("=== ESTADÍSTICAS POR SEGMENTO ===\n\n")## === ESTADÍSTICAS POR SEGMENTO ===seg1 <- window(cdmx_tbm_ts_oqr, end = bp_año - 1/frequency(cdmx_tbm_ts_oqr))
seg2 <- window(cdmx_tbm_ts_oqr, start = bp_año)
cat("Segmento 1 (", start(seg1), "-", end(seg1), "):\n", sep = "")## Segmento 1 (19981-20061):cat(" Observaciones:", length(seg1), "\n")## Observaciones: 9cat(" Media:", mean(seg1), "\n")## Media: 0.5478735cat(" SD:", sd(seg1), "\n\n")## SD: 1.054257cat("Segmento 2 (", start(seg2), "-", end(seg2), "):\n", sep = "")## Segmento 2 (20071-20241):cat(" Observaciones:", length(seg2), "\n")## Observaciones: 18cat(" Media:", mean(seg2), "\n")## Media: -0.2739368cat(" SD:", sd(seg2), "\n\n")## SD: 0.8683267# Cambio
cambio_media <- mean(seg2) - mean(seg1)
cambio_pct <- (cambio_media / mean(seg1)) * 100
cat("Cambio en la media:", round(cambio_media, 4),
"(", round(cambio_pct, 2), "%)\n")## Cambio en la media: -0.8218 ( -150 %)if (cambio_media > 0) {
cat("→ INCREMENTO después de 2007\n")
} else {
cat("→ DISMINUCIÓN después de 2007\n")
}## → DISMINUCIÓN después de 2007gefp(): computes an empirical M-fluctuation process from the scores of a fitted model.
# funcion
gefp<- gefp(cdmx_tbm_ts_oqr~ 1,
fit = glm,
scores = estfun)
sctest(gefp)##
## M-fluctuation test
##
## data: gefp
## f(efp) = 1.1347, p-value = 0.1522# gráfico
plot(gefp,aggregate = FALSE)
breakpoints(): give a number of breaks the function computes the optimal breakpoints.
El criterio de información bayesina (BIC) sugiere elegir los momentos como puntos de cambio estructural, acompañado del residual de la suma de los cuadrados (RSS) que cuantifica la variabilidad de los datos que no es explicada dentro del modelo de regresión
# función de breakpoint
bp.cdmx<-breakpoints(cdmx_tbm_ts_oqr~1,
h = 0.15)
# summary
summary(bp.cdmx)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = cdmx_tbm_ts_oqr ~ 1, h = 0.15)
##
## Breakpoints at observation number:
##
## m = 1 10
## m = 2 10 20
## m = 3 6 10 20
## m = 4 6 10 14 20
## m = 5 6 10 14 18 22
##
## Corresponding to breakdates:
##
## m = 1 2007
## m = 2 2007 2017
## m = 3 2003 2007 2017
## m = 4 2003 2007 2011 2017
## m = 5 2003 2007 2011 2015 2019
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 25.76 20.49 19.02 18.61 18.49 18.60
## BIC 81.95 82.36 86.94 92.94 99.36 106.11# criterios BIC & RSS
plot(bp.cdmx)
Coeficientes de los momentos del cambio estructural
# coeficientes por cambio estructural
coef(bp.cdmx)## (Intercept)
## 1998 - 2024 6.409876e-17Representación gráfica del breakpoint
# intervalo de confianza conforme al BIC
interconf<- confint(bp.cdmx,
breaks = 1, # cambia conforme al BIC
het.err = FALSE,
level = 0.95,
method= "delta")
interconf##
## Confidence intervals for breakpoints
## of optimal 2-segment partition:
##
## Call:
## confint.breakpointsfull(object = bp.cdmx, level = 0.95, breaks = 1,
## het.err = FALSE, method = "delta")
##
## Breakpoints at observation number:
## 2.5 % breakpoints 97.5 %
## 1 -1 10 21
##
## Corresponding to breakdates:## 2.5 % breakpoints 97.5 %
## 1 1996 2007 2018# grafico
plot(cdmx_tbm_ts_oqr)
lines(bp.cdmx)
lines(interconf)
Generalized fluctuation test (GFT)
efp():empirical fluctuation process
# 0. Test recursive CUSUM process
efp_result <- efp(cdmx_tbm_ts_oqr~1,
type = "OLS-CUSUM", # "Score-CUSUM"
dynamic = TRUE)
# Plot the process
plot(efp_result,
alpha = 0.05,
boundary = TRUE,
functional = NULL)
# 1. Test de CUSUM
cusum_test <- efp(cdmx_tbm_ts_oqr ~ 1, type = "Rec-CUSUM")
plot(cusum_test)
sctest(cusum_test)##
## Recursive CUSUM test
##
## data: cusum_test
## S = 0.61228, p-value = 0.3902# 2. Test de MOSUM
mosum_test <- efp(cdmx_tbm_ts_oqr ~ 1, type = "OLS-MOSUM")
plot(mosum_test)
sctest(mosum_test)##
## OLS-based MOSUM test
##
## data: mosum_test
## M0 = 0.78734, p-value = 0.3679sctest(): generic function for performing structural change tests
# Perform a significance test
sctest(efp_result)##
## OLS-based CUSUM test
##
## data: efp_result
## S0 = 1.2963, p-value = 0.0694Conclusiones: no se presenta un cambio estructural en la serie con fecha de 2024
2. Changue-point
Se busca identificar los periodos específicos de tiempo que se relaciona ante un cambio en la probabilidad de distribución de un proceso estocástico. Para ello se realiza un ejercicio sobre el cambio en la media y en la varianza a lo largo de la serie. Change point methods
Convertir la serie de tiempo en formato anual
if(FALSE){
# serie de tiempo de forma anual
chp_tbm<-ts_to_tbl(cdmx_ts)
chp_tbm$año<- year(chp_tbm$date_col)
chp_tbm <- chp_tbm %>%
group_by(año) %>%
summarize(tbm=sum(value))
chp_tbm$fecha<-as.Date(with(chp_tbm,paste(año,1,1,sep="-")),"%Y-%m-%d")
chp_tbm<-chp_tbm[,c(2,3)]
chp_tbm<-ts(chp_tbm$tbm,
start = c(1998,1),
frequency=1)
}- changepoint: changes in mean and variance of time series.
if(FALSE){
library(changepoint)
mean<-cpt.mean(cdmx_ts,
method = "AMOC")
print(mean)
plot(mean, diagnostic = TRUE)
}if(FALSE){
# función por media y por varianza
mean_var <- cpt.meanvar(cdmx_ts,
penalty = "AIC",
test.stat="Normal",
method = "AMOC",
class=TRUE)
mean_var
# gráfico
plot(mean_var,
type = "l",
cpt.col = "red",
xlab = "Periodo",
cpt.width = 2,
diagnostic = TRUE)
}Conclusión:
- Segment
if(FALSE){
library(tidychangepoint)
tidy_cpt <- segment(cdmx_ts, method = "pelt")
changepoints(tidy_cpt)
augment(tidy_cpt)
tidy(tidy_cpt)
plot(tidy_cpt)
}Data frame con los valores originales, asignadole los puntos de cambio a cada observación de la serie
if(FALSE){
library(tidyr)
ts_to_tbl(cdmx_ts) %>%
mutate(observation = 1:n()) %>%
mutate(Breakpoint = case_when(
observation == 2 ~ "1bp",
observation == 7 ~ "2bp",
observation == 13 ~ "3bp",
observation == 16 ~ "4bp",
observation == 18 ~ "5bp",
observation == 54 ~ "6bp",
observation == 56 ~ "7bp",
observation == 59 ~ "10bp",
observation == 127 ~ "11bp",
observation == 132 ~ "12bp",
observation == 301 ~ "13bp",
observation == 305 ~ "14bp",
observation == 310 ~ "15bp")) %>%
fill(Breakpoint, .direction="down")
}- Time-to-event analysis
El análisis solo funciona con series de tiempo estacionales, para el caso objeto de estudio no reúne tal requisito.
if(FALSE){
library(datasets)
# time series to dataframe
df <- ts_to_tbl(cdmx_ts) %>%
select(-index) # %>%
filter (date_col <= '2004-05-01')
# event-analysis
event_tbl <- ts_time_event_analysis_tbl(
.data = df,
.horizon = 6,
.percent_change = 0.05,
.date_col = date_col,
.value_col = value,
.direction = "both",
.filter_non_event_groups = T
)
glimpse(event_tbl)
}if(FALSE){
ts_event_analysis_plot(
.data = event_tbl,
.plot_type = "individual",
.interactive = TRUE
)
}