Índice de Calidad de Vida Municipal de Nuevo León
Yaqing Fan & Osiris Figón Barajas
2026-06-11
Importar base de datos
library(readxl)
DB <- read_excel("Base_Multi.xlsx", sheet = "DB")
Glosario <- read_excel("Base_Multi.xlsx", sheet = "Glosario")La base traía 52 filas pero la última está vacía, por lo que se elimina. Quedan los 51 municipios de Nuevo León.
DB <- DB[!is.na(DB$CVE_MUN), ]
dim(DB)## [1] 51 32Dejar solo las variables para el índice
Se quitan los identificadores (CVE_MUN y NOM_MUN) y se deja las 30 variables de las cuatro dimensiones (educación, salud, vivienda e inclusión económica y digital).
D <- DB[, !(names(DB) %in% c("CVE_MUN","NOM_MUN"))]
dim(D)## [1] 51 30Describir la base de datos
summary(D)## E_ALF E_ESC E_ASISS E_ASISQ
## Min. :88.36 Min. : 6.591 Min. :86.51 Min. :20.97
## 1st Qu.:95.40 1st Qu.: 7.949 1st Qu.:92.70 1st Qu.:31.78
## Median :96.78 Median : 8.845 Median :94.14 Median :36.99
## Mean :96.25 Mean : 8.933 Mean :93.83 Mean :37.49
## 3rd Qu.:98.01 3rd Qu.: 9.647 3rd Qu.:95.54 3rd Qu.:43.28
## Max. :98.91 Max. :13.157 Max. :96.99 Max. :58.00
##
## E_SEC E_PREPA E_UNI E_DOC
## Min. :39.22 Min. :10.26 Min. : 0.9393 Min. : 91.55
## 1st Qu.:50.86 1st Qu.:18.41 1st Qu.: 4.6711 1st Qu.:221.67
## Median :57.30 Median :24.03 Median : 5.9226 Median :286.53
## Mean :58.50 Mean :25.52 Mean : 9.0486 Mean :308.18
## 3rd Qu.:66.85 3rd Qu.:28.25 3rd Qu.:10.5171 3rd Qu.:383.73
## Max. :79.64 Max. :60.40 Max. :43.9834 Max. :728.21
##
## S_AFIL S_IMSS S_ISSSTE S_SPOP
## Min. :66.43 Min. : 2.518 Min. : 0.9552 Min. : 3.75
## 1st Qu.:80.93 1st Qu.:26.054 1st Qu.: 2.9215 1st Qu.:12.50
## Median :85.50 Median :48.648 Median : 4.1453 Median :41.68
## Mean :85.00 Mean :48.245 Mean : 4.7252 Mean :41.83
## 3rd Qu.:89.04 3rd Qu.:77.169 3rd Qu.: 6.1189 3rd Qu.:64.98
## Max. :94.75 Max. :88.745 Max. :11.2405 Max. :95.75
##
## S_SUPERV VSPT VCE VCD
## Min. :91.54 Min. : 76.94 Min. :96.22 Min. :45.61
## 1st Qu.:93.83 1st Qu.: 98.69 1st Qu.:98.64 1st Qu.:94.12
## Median :94.42 Median : 99.25 Median :99.19 Median :97.38
## Mean :94.44 Mean : 98.15 Mean :98.96 Mean :94.42
## 3rd Qu.:95.47 3rd Qu.: 99.54 3rd Qu.:99.66 3rd Qu.:99.10
## Max. :96.25 Max. :100.00 Max. :99.93 Max. :99.91
## NA's :15
## VCEX VCAE VTMR VCR
## Min. : 94.11 Min. :28.87 Min. :23.61 Min. :70.93
## 1st Qu.: 98.25 1st Qu.:89.76 1st Qu.:68.17 1st Qu.:93.39
## Median : 99.13 Median :94.24 Median :87.06 Median :95.72
## Mean : 98.70 Mean :87.71 Mean :77.55 Mean :94.09
## 3rd Qu.: 99.67 3rd Qu.:97.80 3rd Qu.:95.70 3rd Qu.:96.91
## Max. :100.00 Max. :99.70 Max. :99.22 Max. :98.91
##
## V_TIN V_CIS B_OCP B_INT
## Min. : 4.630 Min. : 0.3755 Min. :92.06 Min. : 4.41
## 1st Qu.: 9.741 1st Qu.: 1.2386 1st Qu.:97.85 1st Qu.:28.33
## Median :16.006 Median : 3.0748 Median :98.23 Median :38.96
## Mean :16.967 Mean : 4.0762 Mean :97.92 Mean :42.32
## 3rd Qu.:21.333 3rd Qu.: 4.5042 3rd Qu.:98.69 3rd Qu.:56.06
## Max. :54.239 Max. :27.8734 Max. :99.23 Max. :86.90
##
## B_COMP B_TELFIX B_CEL B_LAV
## Min. : 5.263 Min. : 8.434 Min. :62.39 Min. :67.72
## 1st Qu.:20.140 1st Qu.:18.342 1st Qu.:87.90 1st Qu.:80.05
## Median :25.952 Median :31.033 Median :90.74 Median :86.22
## Mean :29.513 Mean :33.362 Mean :88.46 Mean :84.18
## 3rd Qu.:35.946 3rd Qu.:39.458 3rd Qu.:92.46 3rd Qu.:88.50
## Max. :75.899 Max. :80.380 Max. :96.27 Max. :93.38
##
## B_AUTO B_MOTO
## Min. :30.03 Min. : 3.471
## 1st Qu.:53.51 1st Qu.: 4.960
## Median :61.84 Median : 7.404
## Mean :60.99 Mean : 8.642
## 3rd Qu.:71.39 3rd Qu.:10.676
## Max. :82.65 Max. :21.138
## colSums(is.na(D))## E_ALF E_ESC E_ASISS E_ASISQ E_SEC E_PREPA E_UNI E_DOC
## 0 0 0 0 0 0 0 0
## S_AFIL S_IMSS S_ISSSTE S_SPOP S_SUPERV VSPT VCE VCD
## 0 0 0 0 15 0 0 0
## VCEX VCAE VTMR VCR V_TIN V_CIS B_OCP B_INT
## 0 0 0 0 0 0 0 0
## B_COMP B_TELFIX B_CEL B_LAV B_AUTO B_MOTO
## 0 0 0 0 0 0La variable S_SUPERV (supervivencia de hijas e hijos nacidos vivos) tiene 15 valores faltantes de 51 municipios, el 29% de la base. Imputar casi una tercera parte con tan pocas observaciones sesgaría el índice, así que se elimina la variable. Es una pérdida para la dimensión salud que se reporta como limitación.
D$S_SUPERV <- NULLEn el resto no se observan valores faltantes. Sí hay municipios extremos (Mier y Noriega y Doctor Arroyo en el sur, San Pedro y San Nicolás en la zona metropolitana), pero como queremos crear estrategias para municipios con rezago, nos sirve la identificación de observaciones extremas y no se eliminan. La heterogeneidad es enorme, por ejemplo el agua entubada (VCAE) va de 28.9% a 99.7% y el internet (B_INT) de 4.4% a 86.9%. Notar que E_DOC está en otra escala (docentes por mil estudiantes),pero se estandarizará.
Sobre la dirección de la estrategia; como el índice es de calidad de vida, todas las variables deben ir en el mismo sentido en este caso más alto = mejores condiciones. Más adelante el alfa de Cronbach con check.keys nos dirá si alguna en realidad se mueve en sentido inverso.
Validar linealidad
Las curvas le pegan a la correlación y queremos que haya correlación.
library(car)## Warning: package 'car' was built under R version 4.4.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.4.3scatterplotMatrix(D[, c("E_ALF","E_ESC","E_PREPA","S_IMSS",
"VCD","VCAE","B_INT","B_AUTO")], smooth = FALSE)
No se observan curvas, se observa dispersión entre variables pero no curvas, así que es adecuada.Corroboramos con prueba de Bartlett.
Correlación a través de la prueba de Bartlett
H0: R = I (no se debe de utilizar la técnica de análisis de factores) Ha: R diferente I (sí se puede aplicar)
R = Matriz de correlación I = Es la matriz identidad
library(psych)## Warning: package 'psych' was built under R version 4.4.3##
## Attaching package: 'psych'## The following object is masked from 'package:car':
##
## logitR <- cor(D)
cortest.bartlett(R, n = nrow(D))## $chisq
## [1] 2077.779
##
## $p.value
## [1] 5.640558e-222
##
## $df
## [1] 406El p-value es de 5.6e-222 y como alfa = 0.05 > p-value, rechazamos la H0 a favor de Ha, esto significa que sí existe correlación entre las variables. Ahora hay que evaluar la calidad de dicha correlación.
Evaluación de la calidad con el índice Kaiser-Meyer-Olkin
KMO(R)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = R)
## Overall MSA = 0.71
## MSA for each item =
## E_ALF E_ESC E_ASISS E_ASISQ E_SEC E_PREPA E_UNI E_DOC
## 0.75 0.80 0.46 0.58 0.52 0.56 0.57 0.65
## S_AFIL S_IMSS S_ISSSTE S_SPOP VSPT VCE VCD VCEX
## 0.58 0.81 0.46 0.76 0.65 0.67 0.86 0.76
## VCAE VTMR VCR V_TIN V_CIS B_OCP B_INT B_COMP
## 0.83 0.72 0.71 0.78 0.39 0.16 0.89 0.84
## B_TELFIX B_CEL B_LAV B_AUTO B_MOTO
## 0.76 0.71 0.73 0.62 0.64|Criterio |Evaluación | |MSA ≥ 0.9 |Excelente | |0.8 ≤ MSA < 0.9 |Bueno | |0.7 ≤ MSA < 0.8 |Aceptable | |0.6 ≤ MSA < 0.7 |Regular | |0.5 ≤ MSA < 0.6 |Bajo | |MSA < 0.5 |Inaceptable|
El índice global es de 0.71, aceptable, pero aquí sí hay variables con MSA menor a 0.5, que deberán quitarse:
B_OCP con 0.16: la ocupación municipal casi no varía entre municipios (va de 92% a 99%), no comparte varianza con el resto. V_CIS con 0.39: tener cisterna no acompaña al resto de los indicadores; además refleja escasez o intermitencia del agua, no bienestar. E_ASISS con 0.46: la asistencia de 6 a 14 años es casi universal (87% a 97%), casi no tiene varianza que aportar. S_ISSSTE con 0.46: la derechohabiencia ISSSTE depende del empleo público federal/estatal y no se relaciona con el resto del bienestar municipal.
Quitar variables y volver a correr el índice KMO
D$B_OCP <- NULL
D$V_CIS <- NULL
D$E_ASISS <- NULL
D$S_ISSSTE <- NULL
R <- cor(D)
KMO(R)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = R)
## Overall MSA = 0.78
## MSA for each item =
## E_ALF E_ESC E_ASISQ E_SEC E_PREPA E_UNI E_DOC S_AFIL
## 0.88 0.80 0.58 0.60 0.62 0.61 0.70 0.76
## S_IMSS S_SPOP VSPT VCE VCD VCEX VCAE VTMR
## 0.81 0.80 0.69 0.63 0.83 0.80 0.83 0.81
## VCR V_TIN B_INT B_COMP B_TELFIX B_CEL B_LAV B_AUTO
## 0.73 0.73 0.93 0.89 0.90 0.89 0.74 0.56
## B_MOTO
## 0.67Al quitar las cuatro, el KMO global sube de 0.71 a 0.78 (aceptable, cerca de bueno) y ya ninguna variable tiene MSA menor a 0.5 (la más pequeña es B_AUTO con 0.56). Ya no es posible mejorar el resultado quitando variables por este criterio.
Vamos a reafirmar con alfa de Cronbach
psych::alpha(data.frame(scale(D)), check.keys = TRUE)## Warning in response.frequencies(x, max = max): response.frequency has been
## deprecated and replaced with responseFrequecy. Please fix your call## Number of categories should be increased in order to count frequencies.## Warning in psych::alpha(data.frame(scale(D)), check.keys = TRUE): Some items were negatively correlated with the first principal component and were automatically reversed.
## This is indicated by a negative sign for the variable name.##
## Reliability analysis
## Call: psych::alpha(x = data.frame(scale(D)), check.keys = TRUE)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.94 0.94 0.99 0.37 15 0.013 -0.25 0.63 0.36
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.91 0.94 0.96
## Duhachek 0.91 0.94 0.96
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## E_ALF 0.93 0.93 0.99 0.36 13 0.014 0.073 0.35
## E_ESC 0.93 0.93 0.99 0.36 13 0.014 0.073 0.36
## E_ASISQ 0.94 0.94 0.99 0.38 15 0.013 0.074 0.37
## E_SEC 0.94 0.94 0.99 0.38 15 0.013 0.073 0.38
## E_PREPA 0.94 0.94 0.99 0.38 15 0.013 0.071 0.38
## E_UNI 0.94 0.94 0.99 0.38 15 0.013 0.072 0.38
## E_DOC 0.94 0.94 0.99 0.38 15 0.013 0.073 0.37
## S_AFIL- 0.94 0.94 0.99 0.38 15 0.013 0.076 0.37
## S_IMSS 0.93 0.93 0.99 0.36 14 0.014 0.072 0.35
## S_SPOP- 0.93 0.93 0.99 0.36 14 0.014 0.072 0.36
## VSPT 0.93 0.93 0.99 0.37 14 0.013 0.074 0.36
## VCE 0.94 0.94 0.99 0.38 14 0.013 0.076 0.36
## VCD 0.93 0.93 0.99 0.36 14 0.014 0.076 0.35
## VCEX 0.93 0.93 0.99 0.36 14 0.014 0.074 0.35
## VCAE 0.93 0.93 0.99 0.36 13 0.014 0.072 0.35
## VTMR 0.93 0.93 0.99 0.37 14 0.014 0.075 0.36
## VCR 0.93 0.93 0.99 0.36 14 0.014 0.074 0.35
## V_TIN- 0.93 0.93 0.99 0.37 14 0.013 0.077 0.37
## B_INT 0.93 0.93 0.99 0.35 13 0.014 0.073 0.35
## B_COMP 0.93 0.93 0.99 0.36 13 0.014 0.074 0.35
## B_TELFIX 0.93 0.93 0.99 0.36 13 0.014 0.075 0.35
## B_CEL 0.93 0.93 0.99 0.36 13 0.014 0.074 0.35
## B_LAV 0.93 0.93 0.99 0.37 14 0.014 0.076 0.36
## B_AUTO 0.94 0.94 0.99 0.39 15 0.012 0.069 0.38
## B_MOTO- 0.94 0.94 0.99 0.38 14 0.013 0.077 0.37
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## E_ALF 51 0.85 0.85 0.85 0.83 -6.1e-15 1
## E_ESC 51 0.83 0.83 0.83 0.80 -5.1e-16 1
## E_ASISQ 51 0.42 0.42 0.42 0.37 -6.5e-17 1
## E_SEC 51 0.37 0.37 0.36 0.31 2.8e-16 1
## E_PREPA 51 0.39 0.39 0.39 0.34 -4.4e-16 1
## E_UNI 51 0.42 0.42 0.42 0.36 -5.6e-17 1
## E_DOC 51 0.45 0.45 0.45 0.40 -3.9e-16 1
## S_AFIL- 51 0.41 0.41 0.38 0.36 -1.6e+00 1
## S_IMSS 51 0.72 0.72 0.73 0.69 6.0e-16 1
## S_SPOP- 51 0.78 0.78 0.78 0.75 -1.6e+00 1
## VSPT 51 0.58 0.58 0.58 0.53 7.3e-17 1
## VCE 51 0.50 0.50 0.49 0.45 -2.5e-14 1
## VCD 51 0.72 0.72 0.72 0.69 -9.8e-16 1
## VCEX 51 0.72 0.72 0.71 0.68 9.5e-15 1
## VCAE 51 0.83 0.83 0.83 0.81 -7.7e-16 1
## VTMR 51 0.69 0.69 0.68 0.65 7.2e-16 1
## VCR 51 0.74 0.74 0.75 0.71 -1.8e-15 1
## V_TIN- 51 0.55 0.55 0.54 0.50 -1.6e+00 1
## B_INT 51 0.88 0.88 0.88 0.86 -3.2e-16 1
## B_COMP 51 0.84 0.84 0.84 0.82 2.8e-16 1
## B_TELFIX 51 0.81 0.81 0.82 0.79 -7.4e-17 1
## B_CEL 51 0.81 0.81 0.80 0.78 -2.7e-15 1
## B_LAV 51 0.67 0.67 0.67 0.63 -8.1e-16 1
## B_AUTO 51 0.23 0.23 0.23 0.17 9.7e-17 1
## B_MOTO- 51 0.49 0.49 0.47 0.44 -1.6e+00 1Revisar el raw_alpha: - Mínimo > 0.60 - Ideal > 0.80
El alfa es de 0.94, muy arriba del ideal. PERO al revisar la hoja de las variables, sí hay variables con un “-” al lado del nombre: S_SPOP, S_AFIL, V_TIN y B_MOTO. El check.keys las invirtió automáticamente, lo que significa que esas cuatro se mueven en sentido INVERSO al resto: a más Seguro Popular, más afiliación general, más tinacos y más motos, PEOR es el resto de los indicadores de bienestar.
Entonces, ¿por qué se mueven al revés, si sus nombres suenan positivos?
cor(D$S_IMSS, D$S_SPOP)## [1] -0.9790126cor(D$S_AFIL, D$B_INT)## [1] -0.3950387cor(D$V_TIN, D$B_INT)## [1] -0.3260777cor(D$B_MOTO, D$B_INT)## [1] -0.3590838S_SPOP correlaciona -0.98 con S_IMSS: son espejo casi perfecto. En México el Seguro Popular afiliaba a la población SIN seguridad social laboral, así que marca informalidad, no acceso pleno a salud. Tenerla junto con S_IMSS además mete colinealidad extrema al modelo.
S_AFIL (% afiliada a servicios de salud) sale más ALTA en los municipios más rezagados (Mier y Noriega 92.9%, Doctor Arroyo 89.6%) que en San Pedro (82.9%) o Monterrey (79.2%), por el peso del Seguro Popular en la afiliación rural.
V_TIN: el tinaco es la respuesta de los hogares a un servicio de agua irregular; abunda donde el agua es intermitente.
B_MOTO: la moto sustituye al auto en contextos de menor ingreso.
Se decide quitar estas variables, ya que invertirlas no tiene sustento conceptual: “% no afiliada a salud” o “% sin tinaco” no son indicadores defendibles de calidad de vida; su relación inversa con el bienestar es un artefacto del contexto (informalidad, agua intermitente), no una carencia que se pueda voltear con sentido. La decisión es teórica, no solo estadística.
D$S_SPOP <- NULL
D$S_AFIL <- NULL
D$V_TIN <- NULL
D$B_MOTO <- NULLLa dimensión salud queda representada solo por S_IMSS (derechohabiencia formal); se reconoce como limitación.
Reconfirmar las validaciones para las variables finales
Quedan 21 variables, todas en la misma dirección. Se vuelven a correr las tres verificaciones para el conjunto definitivo:
dim(D)## [1] 51 21R <- cor(D)
cortest.bartlett(R, n = nrow(D))## $chisq
## [1] 1614.505
##
## $p.value
## [1] 6.176902e-215
##
## $df
## [1] 210KMO(R)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = R)
## Overall MSA = 0.77
## MSA for each item =
## E_ALF E_ESC E_ASISQ E_SEC E_PREPA E_UNI E_DOC S_IMSS
## 0.88 0.83 0.55 0.62 0.63 0.61 0.71 0.79
## VSPT VCE VCD VCEX VCAE VTMR VCR B_INT
## 0.64 0.57 0.90 0.77 0.80 0.92 0.72 0.94
## B_COMP B_TELFIX B_CEL B_LAV B_AUTO
## 0.84 0.86 0.87 0.78 0.50psych::alpha(data.frame(scale(D)), check.keys = TRUE)## Warning in response.frequencies(x, max = max): response.frequency has been
## deprecated and replaced with responseFrequecy. Please fix your call## Number of categories should be increased in order to count frequencies.##
## Reliability analysis
## Call: psych::alpha(x = data.frame(scale(D)), check.keys = TRUE)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.93 0.93 0.99 0.39 13 0.015 -1.3e-15 0.65 0.39
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.9 0.93 0.96
## Duhachek 0.9 0.93 0.96
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## E_ALF 0.92 0.92 0.99 0.38 12 0.016 0.080 0.37
## E_ESC 0.92 0.92 0.99 0.38 12 0.016 0.081 0.38
## E_ASISQ 0.93 0.93 0.99 0.40 13 0.015 0.083 0.43
## E_SEC 0.93 0.93 0.99 0.41 14 0.014 0.080 0.43
## E_PREPA 0.93 0.93 0.99 0.40 14 0.014 0.077 0.43
## E_UNI 0.93 0.93 0.99 0.40 13 0.014 0.079 0.43
## E_DOC 0.93 0.93 0.99 0.40 13 0.014 0.080 0.43
## S_IMSS 0.93 0.93 0.99 0.39 13 0.015 0.080 0.38
## VSPT 0.93 0.93 0.99 0.39 13 0.015 0.080 0.40
## VCE 0.93 0.93 0.99 0.40 13 0.015 0.083 0.41
## VCD 0.93 0.93 0.99 0.39 13 0.015 0.084 0.38
## VCEX 0.93 0.93 0.99 0.39 13 0.015 0.080 0.38
## VCAE 0.92 0.92 0.99 0.38 12 0.016 0.079 0.37
## VTMR 0.93 0.93 0.99 0.39 13 0.015 0.083 0.38
## VCR 0.92 0.92 0.99 0.38 12 0.016 0.081 0.37
## B_INT 0.92 0.92 0.99 0.37 12 0.016 0.081 0.37
## B_COMP 0.92 0.92 0.99 0.37 12 0.016 0.082 0.38
## B_TELFIX 0.92 0.92 0.99 0.38 12 0.016 0.084 0.37
## B_CEL 0.92 0.92 0.99 0.38 12 0.016 0.081 0.37
## B_LAV 0.93 0.93 0.99 0.38 12 0.015 0.085 0.38
## B_AUTO 0.93 0.93 0.99 0.41 14 0.014 0.077 0.43
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## E_ALF 51 0.82 0.82 0.82 0.79 -6.1e-15 1
## E_ESC 51 0.82 0.82 0.82 0.80 -5.1e-16 1
## E_ASISQ 51 0.48 0.48 0.47 0.42 -6.5e-17 1
## E_SEC 51 0.40 0.40 0.40 0.33 2.8e-16 1
## E_PREPA 51 0.43 0.43 0.43 0.36 -4.4e-16 1
## E_UNI 51 0.45 0.45 0.45 0.39 -5.6e-17 1
## E_DOC 51 0.48 0.48 0.47 0.42 -3.9e-16 1
## S_IMSS 51 0.69 0.69 0.68 0.64 6.0e-16 1
## VSPT 51 0.58 0.58 0.58 0.52 7.3e-17 1
## VCE 51 0.50 0.50 0.49 0.44 -2.5e-14 1
## VCD 51 0.69 0.69 0.68 0.65 -9.8e-16 1
## VCEX 51 0.69 0.69 0.69 0.65 9.5e-15 1
## VCAE 51 0.80 0.80 0.80 0.77 -7.7e-16 1
## VTMR 51 0.67 0.67 0.65 0.62 7.2e-16 1
## VCR 51 0.76 0.76 0.76 0.73 -1.8e-15 1
## B_INT 51 0.88 0.88 0.88 0.86 -3.2e-16 1
## B_COMP 51 0.85 0.85 0.85 0.83 2.8e-16 1
## B_TELFIX 51 0.83 0.83 0.83 0.80 -7.4e-17 1
## B_CEL 51 0.77 0.77 0.76 0.74 -2.7e-15 1
## B_LAV 51 0.71 0.71 0.72 0.68 -8.1e-16 1
## B_AUTO 51 0.29 0.29 0.29 0.22 9.7e-17 1Prueba de Bartlett con un p-value de 6.8e-210, se sigue rechazando H0.
KMO global de 0.77, aceptable, con la mínima en B_AUTO (0.51), ninguna debajo de 0.5.
Alfa de 0.93, arriba del ideal, y ya ninguna variable tiene “-” al lado, todas van al mismo sentido.
Ahoira ya se puede usar los datos para aplicar análisis de factores con componentes principales.
Estimar análisis de factores
Se reescalan las variables estandarizándolas.
Dz <- data.frame(scale(D))
modelo1 <- principal(Dz, nfactors = 21, rotate = "none")
modelo1## Principal Components Analysis
## Call: principal(r = Dz, nfactors = 21, rotate = "none")
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12
## E_ALF 0.87 -0.16 -0.23 0.14 -0.27 0.00 0.11 -0.05 0.08 0.11 0.03 -0.03
## E_ESC 0.84 0.27 -0.29 -0.33 -0.02 0.03 0.01 -0.03 -0.05 -0.04 0.03 -0.06
## E_ASISQ 0.41 0.55 0.17 -0.50 -0.09 0.31 0.16 0.27 -0.10 0.01 -0.19 0.00
## E_SEC 0.26 0.72 0.12 0.40 0.06 -0.02 0.43 -0.07 -0.14 0.01 0.09 -0.03
## E_PREPA 0.27 0.85 0.06 0.39 0.09 -0.05 0.07 -0.03 0.01 -0.06 0.03 0.08
## E_UNI 0.31 0.84 0.01 0.35 0.08 -0.06 -0.19 0.08 0.06 -0.05 -0.03 0.09
## E_DOC 0.35 0.80 -0.06 0.27 0.01 -0.05 -0.33 0.04 0.04 0.05 -0.10 -0.16
## S_IMSS 0.76 -0.13 -0.54 -0.06 0.06 0.00 -0.05 0.06 0.04 -0.27 0.04 0.12
## VSPT 0.63 -0.47 0.35 0.06 0.39 -0.21 -0.01 0.04 -0.17 -0.10 -0.04 -0.09
## VCE 0.53 -0.34 -0.10 0.45 0.20 0.57 -0.07 0.00 0.06 0.03 0.12 -0.06
## VCD 0.72 -0.17 -0.04 0.28 -0.53 -0.04 -0.10 0.12 -0.17 0.04 0.08 0.02
## VCEX 0.74 -0.43 0.05 0.34 0.10 0.08 0.03 -0.15 -0.01 0.11 -0.26 0.14
## VCAE 0.85 -0.34 0.10 0.19 -0.19 -0.14 -0.06 0.02 -0.11 0.04 -0.01 0.01
## VTMR 0.72 -0.06 -0.38 -0.15 0.34 -0.19 0.06 0.26 0.10 0.28 0.10 0.05
## VCR 0.79 -0.35 0.44 0.08 0.16 0.00 -0.03 0.10 -0.10 -0.08 0.04 -0.02
## B_INT 0.90 0.13 -0.21 -0.25 0.07 -0.01 0.01 -0.19 0.00 0.01 0.02 -0.02
## B_COMP 0.85 0.26 -0.04 -0.42 -0.02 0.02 -0.03 -0.09 0.04 -0.09 0.08 0.01
## B_TELFIX 0.82 0.28 0.04 -0.37 0.03 -0.01 -0.10 -0.24 -0.07 0.13 -0.03 0.00
## B_CEL 0.82 -0.27 -0.09 0.17 -0.13 -0.13 0.22 0.05 0.29 -0.12 -0.13 -0.11
## B_LAV 0.68 0.00 0.69 -0.13 -0.03 0.09 -0.06 0.03 0.06 -0.04 0.07 0.05
## B_AUTO 0.22 0.09 0.92 -0.12 -0.09 -0.05 0.00 -0.03 0.21 0.05 0.07 0.03
## PC13 PC14 PC15 PC16 PC17 PC18 PC19 PC20 PC21 h2 u2 com
## E_ALF -0.13 0.11 0.08 -0.02 -0.03 0.00 -0.01 0.02 0.00 1 8.9e-16 1.7
## E_ESC -0.02 -0.02 0.08 0.04 0.07 -0.03 0.02 -0.05 0.00 1 -4.4e-16 1.9
## E_ASISQ -0.01 0.00 0.00 -0.02 -0.01 0.02 0.01 0.01 0.00 1 2.6e-15 5.0
## E_SEC 0.01 0.03 -0.05 0.05 0.00 0.02 -0.02 -0.01 0.01 1 1.7e-15 2.9
## E_PREPA -0.01 -0.04 0.04 -0.04 0.00 -0.04 0.05 0.02 -0.03 1 -8.9e-16 1.8
## E_UNI -0.01 -0.04 0.06 -0.03 -0.01 0.02 -0.04 -0.02 0.03 1 4.4e-16 1.9
## E_DOC -0.01 0.07 -0.06 0.04 0.01 0.00 0.01 0.01 -0.01 1 5.6e-16 2.2
## S_IMSS -0.01 0.07 -0.04 0.05 -0.02 0.04 0.03 0.00 0.01 1 -1.1e-15 2.4
## VSPT 0.03 0.05 0.06 -0.02 0.00 -0.01 0.01 0.02 0.03 1 1.2e-15 4.1
## VCE 0.02 -0.04 0.01 0.00 0.00 0.01 0.01 0.01 0.01 1 -6.7e-16 4.2
## VCD 0.14 0.02 0.01 -0.01 0.00 -0.01 0.00 0.00 0.00 1 -1.3e-15 2.7
## VCEX 0.03 0.04 0.00 0.03 0.04 -0.01 -0.01 0.00 -0.01 1 -1.8e-15 2.8
## VCAE -0.12 -0.14 -0.05 0.01 0.03 0.02 0.02 0.01 0.01 1 -2.2e-15 1.8
## VTMR 0.02 -0.01 -0.02 0.01 0.00 -0.01 0.00 0.00 0.00 1 0.0e+00 3.3
## VCR -0.03 0.00 0.01 -0.01 -0.03 0.04 -0.02 -0.03 -0.05 1 6.7e-16 2.3
## B_INT 0.02 0.03 -0.07 -0.12 0.04 0.03 0.00 -0.01 0.00 1 -6.7e-16 1.5
## B_COMP 0.03 -0.03 0.00 0.03 0.03 -0.02 -0.05 0.05 -0.01 1 0.0e+00 1.8
## B_TELFIX 0.04 -0.05 0.02 0.03 -0.10 0.01 0.02 -0.01 0.01 1 3.3e-16 2.1
## B_CEL 0.06 -0.07 -0.01 -0.01 -0.03 -0.01 0.00 -0.01 0.00 1 -1.8e-15 2.1
## B_LAV -0.05 0.04 -0.08 -0.01 -0.02 -0.08 -0.01 -0.02 0.02 1 1.1e-15 2.3
## B_AUTO 0.04 0.03 0.04 0.03 0.04 0.06 0.03 0.00 0.00 1 7.8e-16 1.4
##
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11
## SS loadings 9.58 4.05 2.34 1.78 0.79 0.57 0.46 0.32 0.28 0.25 0.19
## Proportion Var 0.46 0.19 0.11 0.08 0.04 0.03 0.02 0.02 0.01 0.01 0.01
## Cumulative Var 0.46 0.65 0.76 0.85 0.88 0.91 0.93 0.95 0.96 0.97 0.98
## Proportion Explained 0.46 0.19 0.11 0.08 0.04 0.03 0.02 0.02 0.01 0.01 0.01
## Cumulative Proportion 0.46 0.65 0.76 0.85 0.88 0.91 0.93 0.95 0.96 0.97 0.98
## PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19 PC20 PC21
## SS loadings 0.11 0.06 0.06 0.05 0.03 0.02 0.02 0.01 0.01 0.01
## Proportion Var 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## Cumulative Var 0.99 0.99 0.99 1.00 1.00 1.00 1.00 1.00 1.00 1.00
## Proportion Explained 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## Cumulative Proportion 0.99 0.99 0.99 1.00 1.00 1.00 1.00 1.00 1.00 1.00
##
## Mean item complexity = 2.5
## Test of the hypothesis that 21 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0
## with the empirical chi square 0 with prob < NA
##
## Fit based upon off diagonal values = 1El método del porcentaje (Cumulative Var, al menos el 60%) sugiere 2 componentes, porque dos explican el 64.9%.
El método de las raíces unitarias (raíces de 1 o más, primer renglón SS loadings) sugiere 4 componentes: 9.58, 4.05, 2.34 y 1.78.
Y el mean item complexity indica 2.5 componentes, es decir, entre 2 y 3.
Sabemos entonces que hay entre 2 y 4 componentes. Vamos a estimar las tres sugerencias rotadas y comparar antes de decidir.
Rotar e interpretar los componentes
Varimax es el método de rotación que usaremos: las correlaciones altas las hace grandes y las chiquitas más chiquitas, para poder ponerle nombre a cada índice. Y 0.5 es el valor de corte para la prioridad 1.
modelo_2f <- principal(Dz, nfactors = 2, rotate = "varimax")
print(modelo_2f$loadings, cut = 0.5, sort = TRUE)##
## Loadings:
## RC1 RC2
## E_ALF 0.869
## E_ESC 0.669 0.572
## S_IMSS 0.753
## VSPT 0.759
## VCE 0.621
## VCD 0.736
## VCEX 0.852
## VCAE 0.917
## VTMR 0.686
## VCR 0.860
## B_INT 0.787
## B_COMP 0.691 0.561
## B_TELFIX 0.655 0.570
## B_CEL 0.859
## B_LAV 0.633
## E_ASISQ 0.662
## E_SEC 0.769
## E_PREPA 0.894
## E_UNI 0.893
## E_DOC 0.874
## B_AUTO
##
## RC1 RC2
## SS loadings 8.780 4.853
## Proportion Var 0.418 0.231
## Cumulative Var 0.418 0.649Con 2 componentes se mezclan en un solo factor los servicios de vivienda, la salud y la conectividad, y B_AUTO se queda fuera. Demasiado agregado.
modelo_4f <- principal(Dz, nfactors = 4, rotate = "varimax")
print(modelo_4f$loadings, cut = 0.5, sort = TRUE)##
## Loadings:
## RC1 RC4 RC2 RC3
## E_ALF 0.779
## VSPT 0.682
## VCE 0.775
## VCD 0.740
## VCEX 0.903
## VCAE 0.856
## VCR 0.740 0.575
## B_CEL 0.799
## E_ESC 0.912
## E_ASISQ 0.683
## S_IMSS 0.593 0.654
## VTMR 0.652
## B_INT 0.836
## B_COMP 0.896
## B_TELFIX 0.829
## E_SEC 0.870
## E_PREPA 0.972
## E_UNI 0.943
## E_DOC 0.875
## B_LAV 0.837
## B_AUTO 0.956
##
## RC1 RC4 RC2 RC3
## SS loadings 6.168 5.185 3.709 2.694
## Proportion Var 0.294 0.247 0.177 0.128
## Cumulative Var 0.294 0.541 0.717 0.846Con 4 componentes aparecen cargas cruzadas fuertes (E_ALF, S_IMSS, VTMR y VCR cargan arriba de 0.4 en dos factores a la vez) y un componente mezcla educación básica con TIC e IMSS, difícil de nombrar.
modelo3 <- principal(Dz, nfactors = 3, rotate = "varimax")
print(modelo3$loadings, cut = 0.5, sort = TRUE)##
## Loadings:
## RC1 RC2 RC3
## E_ALF 0.900
## E_ESC 0.771 0.516
## S_IMSS 0.912
## VCE 0.604
## VCD 0.701
## VCEX 0.757
## VCAE 0.809
## VTMR 0.791
## B_INT 0.840
## B_COMP 0.693 0.515
## B_TELFIX 0.631 0.531
## B_CEL 0.831
## E_ASISQ 0.658
## E_SEC 0.775
## E_PREPA 0.899
## E_UNI 0.893
## E_DOC 0.868
## VSPT 0.554 0.614
## VCR 0.626 0.734
## B_LAV 0.871
## B_AUTO 0.914
##
## RC1 RC2 RC3
## SS loadings 8.100 4.627 3.246
## Proportion Var 0.386 0.220 0.155
## Cumulative Var 0.386 0.606 0.761La estimación con 3 componentes es la mejor. Explica el 76.1% (arriba del 60%), respeta tres de las cuatro raíces unitarias y prácticamente cada variable carga con un poco más de claridad en un solo componente.
Índice 1. Mide servicios básicos, salud y conectividad del hogar (38.6% de la varianza): VCAE % de viviendas con agua entubada dentro de la vivienda VCEX % de viviendas con excusado VCD % de viviendas con drenaje VCE % de viviendas con electricidad VTMR % de viviendas con techos resistentes S_IMSS % derechohabiente en el IMSS E_ALF % de población alfabetizada E_ESC % de personas con escolaridad B_INT % de viviendas con internet B_CEL % de viviendas con celular B_COMP % de viviendas con computadora B_TELFIX % de viviendas con teléfono fijo
Es el piso material del bienestar: servicios de la vivienda, formalidad en salud, alfabetización y conectividad.
Índice 2. Mide capital educativo medio y superior (22%): E_PREPA % de población con preparatoria E_UNI % de población con licenciatura E_DOC Docentes por cada mil estudiantes E_SEC % de población con secundaria completa E_ASISQ % de asistencia de 15 a 24 años
Índice 3. Mide equipamiento y patrimonio del hogar (15.5%): B_AUTO % de viviendas con automóvil o camioneta B_LAV % de viviendas con lavadora VCR % de viviendas con refrigerador VSPT % de viviendas sin piso de tierra
Si una variable aparece en dos índices, para su interpretación se tiene que elegir en dónde la dejaremos. En este caso VCR y VSPT cargan en los índices 1 y 3, pero por interpretación las dejamos en el 3.
Estimación de índices
Ahora usaremos el método de regresión para obtener los puntajes. Por lo general los scores están en escala Z y se tienen que reescalar a 0-100 con la ecuación de la recta.
modelo_final <- principal(Dz, nfactors = 3, rotate = "varimax",
scores = TRUE, method = "regression")Vamos por los puntajes para guardarlos en la base de datos.
DB$PC1 <- modelo_final$scores[, 1]
DB$PC2 <- modelo_final$scores[, 2]
DB$PC3 <- modelo_final$scores[, 3]Vamos a convertir los índices en escala 0 a 100.
DB$Indice1 <- with(DB, 100 * (PC1 - min(PC1)) / (max(PC1) - min(PC1)))
DB$Indice2 <- with(DB, 100 * (PC2 - min(PC2)) / (max(PC2) - min(PC2)))
DB$Indice3 <- with(DB, 100 * (PC3 - min(PC3)) / (max(PC3) - min(PC3)))El signo de cada componente es arbitrario, hay que verificar que en cada índice 100 sea mejor calidad de vida. Verificar por si algún índice sale invertido, hay que invertir primero los índices parciales y luego calcular de nuevo el índice ponderado; si solo se invierte el global al final, la descripción de los clusters no tiene sentido.
cor(DB$Indice1, DB$VCAE) ## [1] 0.8086104#debe ser POSITIVA (servicios)
cor(DB$Indice2, DB$E_PREPA) ## [1] 0.8986256#debe ser POSITIVA (educación)
cor(DB$Indice3, DB$B_AUTO) ## [1] 0.9143253#debe ser POSITIVA (equipamiento)Las tres correlaciones salen positivas, así que con esta base no es necesario invertir.
Índice global
Para sacar el índice global ponderamos por el porcentaje explicado de la varianza:
Proportion Var 0.386 0.220 0.155
w = ponderador w1 = 0.386/(0.386+0.220+0.155) w2 = 0.220/(0.386+0.220+0.155) w3 = 0.155/(0.386+0.220+0.155)
pv <- modelo_final$Vaccounted["Proportion Var", ]
w <- pv / sum(pv)
round(w, 4) ## RC1 RC2 RC3
## 0.5071 0.2897 0.2032DB$Indice_Global <- w[1]*DB$Indice1 + w[2]*DB$Indice2 + w[3]*DB$Indice3
#Verificación final de dirección
cor(DB$Indice_Global, DB$B_INT) ## [1] 0.920721#positiva:más alto = mejorlibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.4.3##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union#Municipios de menor a mayor calidad de vida
DB %>% arrange(Indice_Global) %>%
select(CVE_MUN, NOM_MUN, Indice1, Indice2, Indice3, Indice_Global) %>%
mutate(across(where(is.numeric), ~round(., 1))) %>% print(n = 51)## # A tibble: 51 × 6
## CVE_MUN NOM_MUN Indice1 Indice2 Indice3 Indice_Global
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 19036 Mier y Noriega 0 30.2 44.6 17.8
## 2 19043 Rayones 11.6 11.3 58.2 21
## 3 19014 Doctor Arroyo 15.7 27.9 51.5 26.5
## 4 19024 General Zaragoza 5.4 84.1 0 27.1
## 5 19030 Iturbide 23.8 33 31 27.9
## 6 19007 Aramberri 17.5 43.6 35.6 28.7
## 7 19017 Galeana 23.8 48.4 54.6 37.2
## 8 19050 Vallecillo 46.1 0 81.3 39.9
## 9 19015 Doctor Coss 37 22.6 100 45.6
## 10 19037 Mina 64.5 7.1 53.9 45.7
## 11 19013 China 50.5 15.3 79.6 46.2
## 12 19022 General Terán 50.8 18.6 81.1 47.6
## 13 19016 Doctor González 57.2 18 67.1 47.8
## 14 19003 Los Aldamas 41.4 22.8 99.6 47.9
## 15 19042 Los Ramones 46.2 24.3 87.5 48.2
## 16 19032 Lampazos de Naranjo 55.1 10.3 85.2 48.2
## 17 19028 Higueras 58 25.2 57.3 48.4
## 18 19027 Los Herreras 46.4 32.1 84.5 50
## 19 19020 General Bravo 54.2 28.5 78.7 51.7
## 20 19002 Agualeguas 56.5 22.3 85 52.4
## 21 19023 General Treviño 55.4 29.5 81.5 53.2
## 22 19051 Villaldama 59.3 33 73.4 54.5
## 23 19011 Cerralvo 59.3 28.4 83.1 55.2
## 24 19008 Bustamante 59.3 38.7 71.5 55.8
## 25 19035 Melchor Ocampo 60.6 28 92.8 57.7
## 26 19012 Ciénega de Flores 97.3 8.8 29.9 58
## 27 19033 Linares 76.3 21.4 65.6 58.2
## 28 19005 Anáhuac 46.8 64.9 77.4 58.2
## 29 19034 Marín 78.7 14.3 69.9 58.3
## 30 19010 El Carmen 96.8 11.7 31 58.8
## 31 19045 Salinas Victoria 89.5 20.7 37.3 59
## 32 19041 Pesquería 97.4 14.7 26 59
## 33 19001 Abasolo 71.4 36.3 67.1 60.4
## 34 19029 Hualahuises 63.2 43 78.5 60.5
## 35 19009 Cadereyta Jiménez 83.9 22.5 58.2 60.9
## 36 19038 Montemorelos 74.1 33.5 71.7 61.9
## 37 19040 Parás 49.2 71.6 84.2 62.8
## 38 19004 Allende 77.1 30.1 84.5 65
## 39 19025 General Zuazua 93.8 34.4 39.2 65.5
## 40 19018 García 97.5 23.6 46.4 65.7
## 41 19044 Sabinas Hidalgo 71.8 51.2 81.7 67.8
## 42 19031 Juárez 92.6 40 48.4 68.4
## 43 19049 Santiago 67 67.4 77.4 69.2
## 44 19021 General Escobedo 95.2 31.6 61.3 69.9
## 45 19047 Hidalgo 66.8 94.2 65.5 74.5
## 46 19026 Guadalupe 94.6 43.2 73.7 75.4
## 47 19006 Apodaca 100 37 73.5 76.4
## 48 19039 Monterrey 87.8 62.6 67.6 76.4
## 49 19048 Santa Catarina 85.2 71.1 65.7 77.2
## 50 19019 San Pedro Garza García 71.2 90.6 76.8 78
## 51 19046 San Nicolás de los Garza 87.7 100 75.8 88.8#Municipios de mayor a menor calidad de vida
DB %>% arrange(desc(Indice_Global)) %>%
select(CVE_MUN, NOM_MUN, Indice1, Indice2, Indice3, Indice_Global) %>%
mutate(across(where(is.numeric), ~round(., 1))) %>% print(n = 51)## # A tibble: 51 × 6
## CVE_MUN NOM_MUN Indice1 Indice2 Indice3 Indice_Global
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 19046 San Nicolás de los Garza 87.7 100 75.8 88.8
## 2 19019 San Pedro Garza García 71.2 90.6 76.8 78
## 3 19048 Santa Catarina 85.2 71.1 65.7 77.2
## 4 19039 Monterrey 87.8 62.6 67.6 76.4
## 5 19006 Apodaca 100 37 73.5 76.4
## 6 19026 Guadalupe 94.6 43.2 73.7 75.4
## 7 19047 Hidalgo 66.8 94.2 65.5 74.5
## 8 19021 General Escobedo 95.2 31.6 61.3 69.9
## 9 19049 Santiago 67 67.4 77.4 69.2
## 10 19031 Juárez 92.6 40 48.4 68.4
## 11 19044 Sabinas Hidalgo 71.8 51.2 81.7 67.8
## 12 19018 García 97.5 23.6 46.4 65.7
## 13 19025 General Zuazua 93.8 34.4 39.2 65.5
## 14 19004 Allende 77.1 30.1 84.5 65
## 15 19040 Parás 49.2 71.6 84.2 62.8
## 16 19038 Montemorelos 74.1 33.5 71.7 61.9
## 17 19009 Cadereyta Jiménez 83.9 22.5 58.2 60.9
## 18 19029 Hualahuises 63.2 43 78.5 60.5
## 19 19001 Abasolo 71.4 36.3 67.1 60.4
## 20 19041 Pesquería 97.4 14.7 26 59
## 21 19045 Salinas Victoria 89.5 20.7 37.3 59
## 22 19010 El Carmen 96.8 11.7 31 58.8
## 23 19034 Marín 78.7 14.3 69.9 58.3
## 24 19005 Anáhuac 46.8 64.9 77.4 58.2
## 25 19033 Linares 76.3 21.4 65.6 58.2
## 26 19012 Ciénega de Flores 97.3 8.8 29.9 58
## 27 19035 Melchor Ocampo 60.6 28 92.8 57.7
## 28 19008 Bustamante 59.3 38.7 71.5 55.8
## 29 19011 Cerralvo 59.3 28.4 83.1 55.2
## 30 19051 Villaldama 59.3 33 73.4 54.5
## 31 19023 General Treviño 55.4 29.5 81.5 53.2
## 32 19002 Agualeguas 56.5 22.3 85 52.4
## 33 19020 General Bravo 54.2 28.5 78.7 51.7
## 34 19027 Los Herreras 46.4 32.1 84.5 50
## 35 19028 Higueras 58 25.2 57.3 48.4
## 36 19032 Lampazos de Naranjo 55.1 10.3 85.2 48.2
## 37 19042 Los Ramones 46.2 24.3 87.5 48.2
## 38 19003 Los Aldamas 41.4 22.8 99.6 47.9
## 39 19016 Doctor González 57.2 18 67.1 47.8
## 40 19022 General Terán 50.8 18.6 81.1 47.6
## 41 19013 China 50.5 15.3 79.6 46.2
## 42 19037 Mina 64.5 7.1 53.9 45.7
## 43 19015 Doctor Coss 37 22.6 100 45.6
## 44 19050 Vallecillo 46.1 0 81.3 39.9
## 45 19017 Galeana 23.8 48.4 54.6 37.2
## 46 19007 Aramberri 17.5 43.6 35.6 28.7
## 47 19030 Iturbide 23.8 33 31 27.9
## 48 19024 General Zaragoza 5.4 84.1 0 27.1
## 49 19014 Doctor Arroyo 15.7 27.9 51.5 26.5
## 50 19043 Rayones 11.6 11.3 58.2 21
## 51 19036 Mier y Noriega 0 30.2 44.6 17.8Arriba quedan San Nicolás de los Garza (88.8), San Pedro Garza García (78), Santa Catarina, Monterrey y Apodaca. Abajo quedan Mier y Noriega (17.8), Rayones, Doctor Arroyo, General Zaragoza e Iturbide, es decir, el sur rural del estado, exactamente los municipios que CONAPO clasifica con mayor marginación. Esa coincidencia es una validación del índice.
Hay municipios con desempeño mixto: General Zaragoza tiene educación de 84.1 pero equipamiento de 0; Vallecillo tiene equipamiento de 81.3 pero educación de 0; la periferia metropolitana (García, Juárez, Pesquería) tiene servicios casi perfectos pero con capital educativo bajo.
Procedemos al análisis cluster # Análisis Cluster
Se agrupa con los índices parciales y el índice global, no con las 21 variables originales, los índices ya resumen la información en dimensiones interpretables y en la misma escala 0 a 100, por eso la distancia euclidiana tiene sentido. Ward une en cada paso los grupos que minimizan el aumento de la varianza interna, lo que da grupos compactos, ideales para tipologías de política pública.
Estimar la matriz de distancia con la distancia euclidiana
d <- dist(DB[, c("Indice1","Indice2","Indice3","Indice_Global")],
method = "euclidean")Obtener cada uno de los pasos de agrupaciones usando el método de aglomeración de Ward, para posteriormente calcular el dendograma:
library(factoextra)## Warning: package 'factoextra' was built under R version 4.4.3## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.4.3##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alpha## Welcome to factoextra!## Want to learn more? See two factoextra-related books at https://www.datanovia.com/en/product/practical-guide-to-principal-component-methods-in-r/HC <- hcut(d, hc_method = "ward.D")
HC##
## Call:
## stats::hclust(d = x, method = hc_method)
##
## Cluster method : ward.D
## Distance : euclidean
## Number of objects: 51Dendograma
fviz_dend(HC)
Gráfico del codo
hcut_ward <- function(x, k) {
if (k == 1) {
list(cluster = rep(1, nrow(x)))
} else {
hcut(x, k = k, hc_method = "ward.D")
}
}
fviz_nbclust(as.data.frame(DB[, c("Indice1","Indice2","Indice3","Indice_Global")]),
FUNcluster = hcut_ward, method = "wss")
A partir de cinco clusters se observa que la disminución de la varianza ya es marginal El codo sugiere entre 5 y 6 clusters. Elegimos 5.
Estimación de clusters
Vamos a estimar los 5 clusters:
HC <- hcut(d, k = 5, hc_method = "ward.D")
DB$ClusterHC <- HC$cluster
fviz_dend(HC, k = 5, rect = TRUE)
Vamos a hacer una tabla para describir los clusters con los índices. # Tabla descripción de clusters
library(dplyr)
#Tabla con el total de observaciones
Tabla1 <- DB %>%
summarize("n" = n(),
"Índice1" = round(mean(Indice1), 1),
"Índice2" = round(mean(Indice2), 1),
"Índice3" = round(mean(Indice3), 1),
"ÍndiceG" = round(mean(Indice_Global), 1))
#Tabla por cada cluster
Tabla2 <- DB %>%
group_by(ClusterHC) %>%
summarize("n" = n(),
"Índice1" = round(mean(Indice1), 1),
"Índice2" = round(mean(Indice2), 1),
"Índice3" = round(mean(Indice3), 1),
"ÍndiceG" = round(mean(Indice_Global), 1))
#Tabla agregada
Tabla <- bind_rows(Tabla2, Tabla1)
Tabla## # A tibble: 6 × 6
## ClusterHC n Índice1 Índice2 Índice3 ÍndiceG
## <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 15 74.2 29 68.5 59.9
## 2 2 13 50.7 21.7 86.1 49.5
## 3 3 9 70.4 74.8 74.7 72.5
## 4 4 7 14 39.8 39.3 26.6
## 5 5 7 95 22 36.9 62
## 6 NA 51 62.1 35.8 65.7 55.2#Municipios por cluster
split(DB$NOM_MUN, DB$ClusterHC)## $`1`
## [1] "Abasolo" "Allende" "Apodaca"
## [4] "Bustamante" "Cadereyta Jiménez" "Doctor González"
## [7] "General Escobedo" "Guadalupe" "Higueras"
## [10] "Hualahuises" "Linares" "Marín"
## [13] "Mina" "Montemorelos" "Villaldama"
##
## $`2`
## [1] "Agualeguas" "Los Aldamas" "Cerralvo"
## [4] "China" "Doctor Coss" "General Bravo"
## [7] "General Terán" "General Treviño" "Los Herreras"
## [10] "Lampazos de Naranjo" "Melchor Ocampo" "Los Ramones"
## [13] "Vallecillo"
##
## $`3`
## [1] "Anáhuac" "San Pedro Garza García"
## [3] "Monterrey" "Parás"
## [5] "Sabinas Hidalgo" "San Nicolás de los Garza"
## [7] "Hidalgo" "Santa Catarina"
## [9] "Santiago"
##
## $`4`
## [1] "Aramberri" "Doctor Arroyo" "Galeana" "General Zaragoza"
## [5] "Iturbide" "Mier y Noriega" "Rayones"
##
## $`5`
## [1] "El Carmen" "Ciénega de Flores" "García"
## [4] "General Zuazua" "Juárez" "Pesquería"
## [7] "Salinas Victoria"n = municipios de cada cluster, y lo que sacó cada cluster en cada índice. Los cinco perfiles:
Alta calidad de vida integral: San Pedro, San Nicolás, Monterrey, Santa Catarina, Santiago con buen desempeño en las tres dimensiones.
Periferia metropolitana en expansión: García, Juárez, Pesquería, Zuazua, El Carmen, Ciénega de Flores y Salinas Victoria. Servicios casi universales pero capital educativo y equipamiento bajos = vivienda nueva con servicios, población joven sin acumulación educativa ni patrimonial.
Calidad de vida media: Apodaca, Guadalupe, Escobedo, Allende, Linares, etc. Sin carencias graves, débiles en educación post-básica; candidatos a becas y absorción en media superior.
Rural norte equipado pero con rezago en servicios y educación: Agualeguas, Los Aldamas, China, Doctor Coss, Vallecillo, con alto equipamiento (86.9) pero rezago en agua y drenaje y educación.
Rezago multidimensional del sur del estado: Mier y Noriega, Doctor Arroyo, Aramberri, Galeana, General Zaragoza, Iturbide y Rayones.Con carencias simultáneas en servicios y equipamiento; coinciden con ser los municipios más marginados en el índice de la CONAPO.
Cluster No Jerárquico, K-Medias
El cluster de k-medias hace menos operaciones que el jerárquico: elige aleatoriamente observaciones como centros y agrupa cada observación donde la distancia sea menor. Por eso se usa en big data; con muchos datos el jerárquico se puede ciclar haciendo cálculos. Aun haciendo menos operaciones, ambos métodos arrojan resultados muy similares.
Gráfico del codo
fviz_nbclust(DB[, c("Indice1","Indice2","Indice3","Indice_Global")],
kmeans, method = "wss")
#en este gráfico cambiamos hcut por kmeansAl igual que con el jerárquico, la varianza comienza a decrecer de manera marginal a partir de 4 ó 5 clusters.
Estimación de los clusters con k-medias
Se establece set.seed(61) para que a todos nos dé la misma solución, dado que se eligen aleatoriamente las observaciones que sirven como centros.
set.seed(61)
KM <- kmeans(DB[, c("Indice1","Indice2","Indice3","Indice_Global")],
centers = 5)
DB$ClusterKM <- KM$clusterHacer la tabla para interpretar los clusters de k-medias con los índices:
#Tabla con el total de observaciones
Tabla1 <- DB %>%
summarize("n" = n(),
"Índice1" = round(mean(Indice1), 1),
"Índice2" = round(mean(Indice2), 1),
"Índice3" = round(mean(Indice3), 1),
"ÍndiceG" = round(mean(Indice_Global), 1))
#Tabla por cada cluster
Tabla2 <- DB %>%
group_by(ClusterKM) %>% #Se cambia ClusterHC por ClusterKM
summarize("n" = n(),
"Índice1" = round(mean(Indice1), 1),
"Índice2" = round(mean(Indice2), 1),
"Índice3" = round(mean(Indice3), 1),
"ÍndiceG" = round(mean(Indice_Global), 1))
#Tabla agregada
Tabla <- bind_rows(Tabla2, Tabla1)
Tabla## # A tibble: 6 × 6
## ClusterKM n Índice1 Índice2 Índice3 ÍndiceG
## <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 11 94.4 26.2 47.7 65.2
## 2 2 7 14 39.8 39.3 26.6
## 3 3 11 67.2 27.3 69.1 56
## 4 4 13 50.7 21.7 86.1 49.5
## 5 5 9 70.4 74.8 74.7 72.5
## 6 NA 51 62.1 35.8 65.7 55.2#Comparación de asignaciones entre métodos
table(Jerarquico = DB$ClusterHC, KMedias = DB$ClusterKM)## KMedias
## Jerarquico 1 2 3 4 5
## 1 4 0 11 0 0
## 2 0 0 0 13 0
## 3 0 0 0 0 9
## 4 0 7 0 0 0
## 5 7 0 0 0 0Ambos métodos coinciden en aproximadamente 46 de 51 municipios más o menos el 90%; las diferencias se concentran en municipios de nivel medio que están en la frontera entre dos grupos, y los cinco perfiles aparecen en los dos.
Exportar los datos a la base original
library(writexl)
Base_Final <- DB %>%
mutate(Ranking = rank(-Indice_Global)) %>%
arrange(Ranking)
write_xlsx(Base_Final, "Base_Final_Indice_Calidad_Vida.xlsx")