library(tidyverse)
library(ggplot2)
library(knitr)
library(kableExtra)
library(car)        # vif()
library(lmtest)     # bptest(), dwtest()
library(nortest)    # lillie.test() para muestras grandes
library(gridExtra)  # grid.arrange()

opts_chunk$set(
  echo      = TRUE,
  message   = FALSE,
  warning   = FALSE,
  fig.align = "center",
  fig.width = 10,
  fig.height = 6
)

1 Introducción

Entender qué hay detrás de la pobreza y el hacinamiento en Ecuador no es solo un ejercicio académico: es la base sobre la que se diseñan políticas de vivienda, programas sociales y estrategias de inversión pública. El Censo de Población y Vivienda 2022 ofrece una oportunidad única para responder esa pregunta con datos a escala nacional, y el Método de la Ingeniería Estadística provee el marco sistemático para hacerlo de forma rigurosa.

Este informe parte de una pregunta concreta: ¿qué características del hogar —tamaño, espacio disponible, acceso a tecnología, ubicación geográfica— determinan el nivel de hacinamiento y la probabilidad de pobreza por NBI, y con qué magnitud? Para responderla se construyen tres modelos de complejidad creciente, se verifica que sus supuestos se cumplan, y se extraen conclusiones accionables para la política pública.

2 Método de la Ingeniería Estadística

2.1 Descripción del problema

El Ecuador presenta disparidades habitacionales persistentes entre territorios y grupos socioeconómicos. El hacinamiento —cuando el número de personas en el hogar supera la capacidad razonable de las habitaciones disponibles— y la pobreza por NBI —definida por la carencia simultánea de condiciones básicas de vivienda, servicios y educación— son dos manifestaciones del mismo problema estructural.

Lo que se busca establecer en este análisis es la magnitud y dirección del efecto de variables observables del hogar sobre estos dos indicadores. Conocer esa relación permite identificar qué factores tienen mayor peso relativo y cuáles deben priorizarse en la intervención pública.

2.2 Identificación de variables

# Carga de datos
base_hogar <- read.csv("CPV_2022_Hogar_Canton.csv",
                       sep = ";", header = TRUE, stringsAsFactors = FALSE)

# Tabla de variables construida desde la base real
vars_interes <- c("H01","H1303","H1004","H1005","H09","AUR","H1011","H1012","HAC","NBI")
vars_ok      <- vars_interes[vars_interes %in% names(base_hogar)]

desc_map <- c(H01="Número de dormitorios", H1303="Número de personas en el hogar",
              H1004="Acceso a internet", H1005="Disponibilidad de computadora",
              H09="Tenencia de la vivienda", AUR="Área urbana o rural",
              H1011="Tenencia de televisor", H1012="Tenencia de refrigeradora",
              HAC="Índice de hacinamiento", NBI="Pobreza por NBI")

tipo_map <- c(H01="Cuantitativa discreta", H1303="Cuantitativa discreta",
              H1004="Cualitativa dicotómica", H1005="Cualitativa dicotómica",
              H09="Cualitativa nominal", AUR="Cualitativa dicotómica",
              H1011="Cualitativa dicotómica", H1012="Cualitativa dicotómica",
              HAC="Cuantitativa continua", NBI="Binaria (0 = No pobre; 1 = Pobre)")

rol_map <- c(H01="Predictor", H1303="Predictor", H1004="Predictor", H1005="Predictor",
             H09="Predictor", AUR="Predictor", H1011="Predictor", H1012="Predictor",
             HAC="Variable respuesta", NBI="Variable respuesta")

tabla_vars <- data.frame(
  Variable    = vars_ok,
  Descripción = desc_map[vars_ok],
  Tipo        = tipo_map[vars_ok],
  Rol         = rol_map[vars_ok],
  row.names   = NULL
)

kable(tabla_vars, caption = "Variables del CPV 2022 utilizadas en el análisis") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE) %>%
  column_spec(4, bold = TRUE,
    color      = ifelse(tabla_vars$Rol == "Variable respuesta", "white", "#2c3e50"),
    background = ifelse(tabla_vars$Rol == "Variable respuesta", "#2c7be5", "transparent"))
Variables del CPV 2022 utilizadas en el análisis
Variable Descripción Tipo Rol
H01 Número de dormitorios Cuantitativa discreta Predictor
H1303 Número de personas en el hogar Cuantitativa discreta Predictor
H1004 Acceso a internet Cualitativa dicotómica Predictor
H1005 Disponibilidad de computadora Cualitativa dicotómica Predictor
H09 Tenencia de la vivienda Cualitativa nominal Predictor
AUR Área urbana o rural Cualitativa dicotómica Predictor
H1011 Tenencia de televisor Cualitativa dicotómica Predictor
H1012 Tenencia de refrigeradora Cualitativa dicotómica Predictor
HAC Índice de hacinamiento Cuantitativa continua Variable respuesta
NBI Pobreza por NBI Binaria (0 = No pobre; 1 = Pobre) Variable respuesta

2.3 Propuesta de modelos

Se plantea una secuencia de tres modelos que abarcan distintas dimensiones del problema:

Modelo Respuesta Naturaleza Propósito
Regresión lineal simple HAC Continua Relación bivariada entre tamaño del hogar y hacinamiento
Regresión lineal múltiple HAC Continua Efecto neto de múltiples predictores simultáneos
Regresión logística NBI Binaria (0/1) Probabilidad de estar en situación de pobreza por NBI

2.4 Recolección de datos

cat("Observaciones (cantones):", nrow(base_hogar),
    "  |  Variables totales:", ncol(base_hogar),
    "  |  Variables utilizadas:", length(vars_ok), "\n\n")
## Observaciones (cantones): 5193548   |  Variables totales: 47   |  Variables utilizadas: 10
vars_num <- c("HAC","H1303","H01")
vars_num <- vars_num[vars_num %in% names(base_hogar)]

desc_num <- base_hogar %>%
  select(all_of(vars_num)) %>%
  pivot_longer(everything(), names_to = "Variable") %>%
  group_by(Variable) %>%
  summarise(N=format(n(),big.mark=","), Media=round(mean(value,na.rm=T),3),
            Mediana=round(median(value,na.rm=T),3), DE=round(sd(value,na.rm=T),3),
            Mín=round(min(value,na.rm=T),3), Máx=round(max(value,na.rm=T),3),
            .groups="drop") %>%
  mutate(Variable = case_when(
    Variable == "HAC"   ~ "HAC \u2014 \u00cdndice de hacinamiento",
    Variable == "H1303" ~ "H1303 \u2014 N\u00b0 personas en el hogar",
    Variable == "H01"   ~ "H01 \u2014 N\u00b0 de dormitorios",
    TRUE ~ Variable
  ))

kable(desc_num, caption="Estadísticas descriptivas de las variables cuantitativas principales") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE) %>%
  column_spec(1, bold=TRUE, width="24em")
Estadísticas descriptivas de las variables cuantitativas principales
Variable N Media Mediana DE Mín Máx
H01 — N° de dormitorios 5,193,548 2.093 2 1.101 0 20
H1303 — N° personas en el hogar 5,193,548 3.262 3 4.922 1 7196
HAC — Índice de hacinamiento 5,193,548 1.912 2 0.284 1 2

3 Modelo de Regresión Lineal Simple

El punto de partida es la relación entre el número de personas en el hogar (H1303) y el índice de hacinamiento (HAC). Antes de agregar más variables, conviene entender si esta relación existe y qué tan fuerte es.

3.1 Estimación

modelo_lineal  <- lm(HAC ~ H1303, data = base_hogar)
resumen_simple <- summary(modelo_lineal)

b0_s     <- round(coef(resumen_simple)[1,1], 4)
b1_s     <- round(coef(resumen_simple)[2,1], 4)
t_b1_s   <- round(coef(resumen_simple)[2,3], 2)
r2_s     <- round(resumen_simple$r.squared, 4)
r2_s_pct <- round(r2_s * 100, 2)
p_b1_s   <- coef(resumen_simple)[2,4]

sig_s <- ifelse(p_b1_s < 0.001, "p < 0.001",
         ifelse(p_b1_s < 0.01,  "p < 0.01",
         ifelse(p_b1_s < 0.05,  "p < 0.05", "p >= 0.05")))

calidad_s <- ifelse(r2_s < 0.20, "bajo — el número de personas solo no captura la variabilidad del hacinamiento cantonal.",
             ifelse(r2_s < 0.50, "moderado — H1303 explica una parte importante, aunque otros factores también influyen.",
                                  "alto — el tamaño del hogar domina la explicación del hacinamiento en este modelo simple."))

resumen_simple
## 
## Call:
## lm(formula = HAC ~ H1303, data = base_hogar)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.86789 -0.04384  0.07346  0.13211  1.71572 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.102e+00  2.504e-04  8395.8   <2e-16 ***
## H1303       -5.865e-02  6.812e-05  -861.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2655 on 5188825 degrees of freedom
##   (4721 observations deleted due to missingness)
## Multiple R-squared:  0.125,  Adjusted R-squared:  0.125 
## F-statistic: 7.414e+05 on 1 and 5188825 DF,  p-value: < 2.2e-16
hac_media <- base_hogar %>%
  filter(!is.na(H1303), !is.na(HAC), H1303 >= 1, H1303 <= 15) %>%
  group_by(H1303) %>%
  summarise(HAC_media=mean(HAC,na.rm=T), HAC_sd=sd(HAC,na.rm=T),
            n_cant=n(), .groups="drop") %>%
  mutate(se=HAC_sd/sqrt(n_cant), IC_inf=HAC_media-1.96*se, IC_sup=HAC_media+1.96*se)

pred_df <- data.frame(H1303=seq(min(hac_media$H1303), max(hac_media$H1303), length.out=120))
pred_df$HAC <- predict(modelo_lineal, newdata=pred_df)

ggplot(hac_media, aes(x=H1303, y=HAC_media)) +
  geom_ribbon(aes(ymin=IC_inf, ymax=IC_sup), fill="#aed6f1", alpha=0.40) +
  geom_line(data=pred_df, aes(x=H1303, y=HAC),
            color="#e74c3c", linewidth=1.5) +
  geom_point(aes(size=n_cant), color="#2471a3", alpha=0.85, shape=16) +
  geom_text(aes(label=H1303), vjust=-1.2, size=3.2, color="#1c2833", fontface="bold") +
  scale_x_continuous(breaks=1:15, limits=c(0.5,15.5)) +
  scale_size_continuous(range=c(2,9), name="N° cantones",
                        guide=guide_legend(title.position="top")) +
  labs(title="Regresión Lineal Simple: personas en el hogar vs. índice de hacinamiento",
       subtitle=paste0("Cada punto = HAC promedio para ese N° de personas  |  ",
                       "Línea roja = recta estimada  |  Banda = IC 95%  |  ",
                       "R\u00b2 = ", r2_s, "  |  ", sig_s),
       x="Número de personas en el hogar (H1303)",
       y="Índice de hacinamiento promedio (HAC)",
       caption="Fuente: CPV 2022") +
  theme_minimal(base_size=13) +
  theme(plot.title=element_text(face="bold", size=12.5),
        plot.subtitle=element_text(color="#555555", size=9.5),
        legend.position="right", panel.grid.minor=element_blank())

3.2 Interpretación

\(\hat{HAC} = 2.1025 + (-0.0587) \times H1303\)

Intercepto (β₀ = 2.1025): Es el valor de HAC proyectado cuando H1303 = 0, un escenario teórico sin sentido práctico. Su función es anclar la recta en la escala de la variable respuesta.

Pendiente (β₁ = -0.0587, t = -861.06, p < 0.001): Por cada persona adicional que en promedio tienen los hogares de un cantón, el índice HAC disminuye en 0.0587 unidades. El signo negativo puede parecer extraño a primera vista, pero tiene sentido a nivel de agregado cantonal: los cantones con hogares más numerosos tienden también a contar con más dormitorios en promedio, de manera que la razón personas/dormitorios puede bajar levemente. Este efecto de composición desaparece cuando se incluye H01 en el modelo múltiple. El resultado es p < 0.001 — la relación no es un artefacto estadístico.

R² = 0.125: H1303 explica el 12.5% de la variabilidad del hacinamiento entre cantones. Es un ajuste bajo — el número de personas solo no captura la variabilidad del hacinamiento cantonal.

3.3 Verificación de supuestos

Los tres supuestos clásicos de la regresión lineal —normalidad de residuos, homocedasticidad e independencia— se verifican a continuación con pruebas formales.

res_s <- residuals(modelo_lineal)
fit_s <- fitted(modelo_lineal)

# Muestra para visualización (los gráficos base con n > 1M son muy lentos)
set.seed(42)
idx_s  <- sample(length(res_s), min(5000, length(res_s)))
df_s   <- data.frame(Ajustados=fit_s[idx_s], Residuos=res_s[idx_s],
                     ResStd=rstandard(modelo_lineal)[idx_s])

ps1 <- ggplot(df_s, aes(x=Ajustados, y=Residuos)) +
  geom_point(alpha=0.3, color="#2471a3", size=1) +
  geom_hline(yintercept=0, linetype="dashed", color="gray40") +
  geom_smooth(method="lm", se=FALSE, color="#e74c3c", linewidth=0.9) +
  labs(title="Residuos vs. Ajustados", x="Valores ajustados", y="Residuos") +
  theme_minimal(base_size=12) + theme(plot.title=element_text(face="bold", size=11))

ps2 <- ggplot(df_s, aes(sample=ResStd)) +
  stat_qq(color="#2471a3", alpha=0.4, size=1) +
  stat_qq_line(color="#e74c3c", linewidth=0.9) +
  labs(title="Q-Q Normal de residuos estandarizados",
       x="Cuantiles teóricos", y="Cuantiles muestrales") +
  theme_minimal(base_size=12) + theme(plot.title=element_text(face="bold", size=11))

ps3 <- ggplot(df_s, aes(x=Ajustados, y=sqrt(abs(ResStd)))) +
  geom_point(alpha=0.3, color="#2471a3", size=1) +
  geom_smooth(method="lm", se=FALSE, color="#e74c3c", linewidth=0.9) +
  labs(title="Scale-Location", x="Valores ajustados", y="\u221a|Residuos estand.|") +
  theme_minimal(base_size=12) + theme(plot.title=element_text(face="bold", size=11))

ps4 <- ggplot(df_s, aes(x=seq_along(Residuos), y=Residuos)) +
  geom_point(alpha=0.3, color="#2471a3", size=1) +
  geom_hline(yintercept=0, linetype="dashed", color="gray40") +
  labs(title="Residuos por índice (independencia)",
       x="Índice de observación", y="Residuos") +
  theme_minimal(base_size=12) + theme(plot.title=element_text(face="bold", size=11))

gridExtra::grid.arrange(ps1, ps2, ps3, ps4, ncol=2)

# Lilliefors sobre muestra: con n > 5000 la prueba es idéntica en conclusión
# y con millones de filas tarda mucho. Se usa muestra representativa.
set.seed(42)
lilli_s <- lillie.test(sample(res_s, min(5000, length(res_s))))
bp_s    <- bptest(modelo_lineal)   # opera sobre el modelo: rápido
dw_s    <- dwtest(modelo_lineal)   # rápido

# ── Tabla resumen ─────────────────────────────────────────────────
supuestos_s <- data.frame(
  Supuesto   = c("Normalidad de residuos", "Homocedasticidad", "Independencia de residuos"),
  Prueba     = c("Lilliefors (K-S)", "Breusch-Pagan", "Durbin-Watson"),
  Estadístico = c(round(lilli_s$statistic, 4),
                  round(bp_s$statistic, 4),
                  round(dw_s$statistic, 4)),
  `p-valor`  = c(
    ifelse(lilli_s$p.value < 0.001, "< 0.001", round(lilli_s$p.value, 4)),
    ifelse(bp_s$p.value    < 0.001, "< 0.001", round(bp_s$p.value, 4)),
    ifelse(dw_s$p.value    < 0.001, "< 0.001", round(dw_s$p.value, 4))
  ),
  Conclusión = c(
    ifelse(lilli_s$p.value < 0.05,
      "Se rechaza H\u2080: los residuos no siguen distribución normal",
      "No se rechaza H\u2080: residuos compatibles con normalidad"),
    ifelse(bp_s$p.value < 0.05,
      "Se rechaza H\u2080: varianza no constante (heterocedasticidad)",
      "No se rechaza H\u2080: varianza constante (homocedasticidad)"),
    ifelse(dw_s$p.value < 0.05,
      "Se rechaza H\u2080: hay autocorrelación en los residuos",
      "No se rechaza H\u2080: residuos independientes")
  )
)

kable(supuestos_s, col.names=c("Supuesto","Prueba","Estadístico","p-valor","Conclusión"),
      caption="Verificación de supuestos — Regresión Lineal Simple") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE) %>%
  column_spec(5, bold=TRUE,
    color=ifelse(grepl("rechaza", supuestos_s$Conclusión), "#922b21", "#1e8449"),
    background=ifelse(grepl("rechaza", supuestos_s$Conclusión), "#fdf2f2", "#f0faf3"))
Verificación de supuestos — Regresión Lineal Simple
Supuesto Prueba Estadístico p-valor Conclusión
D Normalidad de residuos Lilliefors (K-S) 0.3443 < 0.001 Se rechaza H₀: los residuos no siguen distribución normal
BP Homocedasticidad Breusch-Pagan 731076.2327 < 0.001 Se rechaza H₀: varianza no constante (heterocedasticidad)
DW Independencia de residuos Durbin-Watson 1.9314 < 0.001 Se rechaza H₀: hay autocorrelación en los residuos

Lectura de los supuestos: La prueba de Lilliefors indica que los residuos no se distribuyen normalmente. Esto es frecuente en datasets censales de gran tamaño, donde cualquier desviación —por pequeña que sea— resulta estadísticamente significativa. En muestras grandes los estimadores MCO conservan propiedades asintóticas robustas (Teorema Central del Límite), por lo que los coeficientes y sus pruebas t siguen siendo válidos. La prueba de Breusch-Pagan detecta heterocedasticidad, lo que sugiere que la varianza de los residuos no es constante a lo largo del rango de los ajustados. Los errores estándar pueden estar sesgados; para inferencia más robusta convendría aplicar errores estándar robustos (HC). El estadístico Durbin-Watson señala autocorrelación en los residuos. Dado que los datos son de corte transversal cantonal y no series de tiempo, este resultado probablemente refleja dependencia espacial entre cantones vecinos.

4 Modelo de Regresión Lineal Múltiple

Con los resultados del modelo simple como punto de partida, se incorporan el número de dormitorios (H01), el acceso a internet (H1004), la disponibilidad de computadora (H1005) y el área urbano-rural (AUR). El objetivo es estimar el efecto de cada variable manteniendo las demás constantes.

4.1 Estimación

modelo_multiple <- lm(HAC ~ H1303 + H01 + H1004 + H1005 + AUR, data=base_hogar)
resumen_mult    <- summary(modelo_multiple)

cm       <- coef(resumen_mult)
b_per  <- round(cm["H1303",1], 4)
b_dor  <- round(cm["H01",  1], 4)
b_int  <- round(cm["H1004",1], 4)
b_com  <- round(cm["H1005",1], 4)
b_aur  <- round(cm["AUR",  1], 4)
r2_m     <- round(resumen_mult$adj.r.squared, 4)
r2_m_pct <- round(r2_m * 100, 2)
f_stat   <- round(resumen_mult$fstatistic[1], 1)

sig_m <- function(var) {
  p <- cm[var,4]
  ifelse(p<0.001,"significativo (p < 0.001)",
  ifelse(p<0.01, "significativo (p < 0.01)",
  ifelse(p<0.05, "significativo (p < 0.05)",
                 "no significativo (p >= 0.05)")))
}
dir_m <- function(b) ifelse(b > 0, "aumenta", "reduce")

resumen_mult
## 
## Call:
## lm(formula = HAC ~ H1303 + H01 + H1004 + H1005 + AUR, data = base_hogar)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.15069 -0.04952  0.04095  0.14272  1.77391 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  1.959e+00  5.871e-04  3337.13   <2e-16 ***
## H1303       -8.341e-02  6.312e-05 -1321.42   <2e-16 ***
## H01          1.181e-01  1.066e-04  1108.32   <2e-16 ***
## H1004       -1.637e-02  2.546e-04   -64.31   <2e-16 ***
## H1005       -9.318e-03  2.563e-04   -36.36   <2e-16 ***
## AUR          1.060e-02  2.200e-04    48.21   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2327 on 5188821 degrees of freedom
##   (4721 observations deleted due to missingness)
## Multiple R-squared:  0.3278, Adjusted R-squared:  0.3278 
## F-statistic: 5.062e+05 on 5 and 5188821 DF,  p-value: < 2.2e-16
coefs_plot <- data.frame(
  Variable = rownames(cm)[-1],
  Estimado = cm[-1,1], Error=cm[-1,2], Pvalor=cm[-1,4]
) %>%
  mutate(
    IC_inf = Estimado - 1.96*Error, IC_sup = Estimado + 1.96*Error,
    Sig    = ifelse(Pvalor<0.001,"p < 0.001",
             ifelse(Pvalor<0.01, "p < 0.01",
             ifelse(Pvalor<0.05, "p < 0.05","No significativo"))),
    Variable = case_when(
      Variable == "H1303" ~ "N\u00b0 personas en el hogar (H1303)",
      Variable == "H01"   ~ "N\u00b0 dormitorios (H01)",
      Variable == "H1004" ~ "Acceso a internet (H1004)",
      Variable == "H1005" ~ "Computadora (H1005)",
      Variable == "AUR"   ~ "\u00c1rea urbano-rural (AUR)",
      TRUE ~ Variable
    )
  )

pal_sig <- c("p < 0.001"="#1a5276","p < 0.01"="#2980b9",
             "p < 0.05"="#85c1e9","No significativo"="#bdc3c7")

ggplot(coefs_plot, aes(x=reorder(Variable,Estimado), y=Estimado, color=Sig)) +
  geom_hline(yintercept=0, linetype="dashed", color="#7f8c8d", linewidth=0.8) +
  geom_errorbar(aes(ymin=IC_inf, ymax=IC_sup), width=0.25, linewidth=1.3) +
  geom_point(size=5) +
  geom_label(aes(label=sprintf("\u03b2=%s", round(Estimado,4))),
             hjust=ifelse(coefs_plot$Estimado>=0,-0.15,1.15),
             size=3.4, fill="white", label.size=0,
             label.padding=unit(0.15,"lines"), show.legend=FALSE) +
  scale_color_manual(values=pal_sig, name="Significancia") +
  coord_flip() +
  labs(title="Coeficientes — Regresión Lineal Múltiple (HAC)",
       subtitle=paste0("R\u00b2 ajustado = ",r2_m,"  |  F = ",f_stat,"  |  IC al 95%"),
       x=NULL, y="Coeficiente estimado (\u03b2)", caption="Fuente: CPV 2022") +
  theme_minimal(base_size=13) +
  theme(plot.title=element_text(face="bold"),
        plot.subtitle=element_text(color="#555555",size=10),
        legend.position="right", panel.grid.minor=element_blank())

4.2 Interpretación

R² ajustado = 0.3278 — el modelo explica el 32.78% de la variabilidad cantonal del hacinamiento con cinco predictores simultáneos (F = 5.061787^{5}, p < 0.001).

Manteniendo las demás variables constantes:

H1303 — N° personas (β = -0.0834, significativo (p < 0.001)): Cada persona adicional reduce el HAC en 0.0834 unidades. Al controlar por H01, el efecto neto resulta negativo: dado un mismo número de dormitorios, los cantones con hogares más numerosos tienen una razón HAC menor porque ambas variables comparten varianza. La correlación entre H1303 y H01 es la que genera este signo — es un resultado esperado en modelos con predictores correlacionados. Coeficiente significativo (p < 0.001).

H01 — N° dormitorios (β = 0.1181, significativo (p < 0.001)): Cada dormitorio adicional aumenta el HAC en 0.1181 unidades. El signo positivo se explica por la correlación con H1303: los cantones con más dormitorios son también los de mayor tamaño de hogar, lo que eleva el índice en términos agregados. El efecto protector de los dormitorios sobre el bienestar queda mejor capturado en el modelo logístico (OR < 1). Coeficiente significativo (p < 0.001).

H1004 — Internet (β = -0.0164, significativo (p < 0.001)): Los cantones con más hogares conectados muestran un HAC reduce en 0.0164 unidades. La conectividad digital refleja el nivel socioeconómico del territorio: cantones más desarrollados tienen tanto mayor acceso a internet como viviendas más amplias. Coeficiente significativo (p < 0.001).

H1005 — Computadora (β = -0.0093, significativo (p < 0.001)): Mayor disponibilidad de computadora reduce el HAC en 0.0093 unidades, siendo un proxy de capital económico favorable asociado a mejores condiciones habitacionales. Coeficiente significativo (p < 0.001).

AUR — Área (β = 0.0106, significativo (p < 0.001)): La diferencia territorial es de 0.0106 unidades en HAC entre áreas rurales (mayor) y urbanas, persistiendo después de controlar por todos los demás factores. Esto señala que hay condiciones propias del territorio que no quedan explicadas por el tamaño del hogar ni la tecnología. Coeficiente significativo (p < 0.001).

4.3 Verificación de supuestos y multicolinealidad

4.3.1 Supuestos del modelo

res_m <- residuals(modelo_multiple)
fit_m <- fitted(modelo_multiple)

set.seed(42)
idx_m2 <- sample(length(res_m), min(5000, length(res_m)))
df_m   <- data.frame(Ajustados=fit_m[idx_m2], Residuos=res_m[idx_m2],
                     ResStd=rstandard(modelo_multiple)[idx_m2])

pm1 <- ggplot(df_m, aes(x=Ajustados, y=Residuos)) +
  geom_point(alpha=0.3, color="#1a5276", size=1) +
  geom_hline(yintercept=0, linetype="dashed", color="gray40") +
  geom_smooth(method="lm", se=FALSE, color="#e74c3c", linewidth=0.9) +
  labs(title="Residuos vs. Ajustados", x="Valores ajustados", y="Residuos") +
  theme_minimal(base_size=12) + theme(plot.title=element_text(face="bold", size=11))

pm2 <- ggplot(df_m, aes(sample=ResStd)) +
  stat_qq(color="#1a5276", alpha=0.4, size=1) +
  stat_qq_line(color="#e74c3c", linewidth=0.9) +
  labs(title="Q-Q Normal de residuos estandarizados",
       x="Cuantiles teóricos", y="Cuantiles muestrales") +
  theme_minimal(base_size=12) + theme(plot.title=element_text(face="bold", size=11))

pm3 <- ggplot(df_m, aes(x=Ajustados, y=sqrt(abs(ResStd)))) +
  geom_point(alpha=0.3, color="#1a5276", size=1) +
  geom_smooth(method="lm", se=FALSE, color="#e74c3c", linewidth=0.9) +
  labs(title="Scale-Location", x="Valores ajustados", y="\u221a|Residuos estand.|") +
  theme_minimal(base_size=12) + theme(plot.title=element_text(face="bold", size=11))

pm4 <- ggplot(df_m, aes(x=seq_along(Residuos), y=Residuos)) +
  geom_point(alpha=0.3, color="#1a5276", size=1) +
  geom_hline(yintercept=0, linetype="dashed", color="gray40") +
  labs(title="Residuos por índice (independencia)",
       x="Índice de observación", y="Residuos") +
  theme_minimal(base_size=12) + theme(plot.title=element_text(face="bold", size=11))

gridExtra::grid.arrange(pm1, pm2, pm3, pm4, ncol=2)

set.seed(42)
lilli_m <- lillie.test(sample(res_m, min(5000, length(res_m))))
bp_m    <- bptest(modelo_multiple)
dw_m    <- dwtest(modelo_multiple)

supuestos_m <- data.frame(
  Supuesto    = c("Normalidad de residuos","Homocedasticidad","Independencia de residuos"),
  Prueba      = c("Lilliefors (K-S)","Breusch-Pagan","Durbin-Watson"),
  Estadístico = c(round(lilli_m$statistic,4), round(bp_m$statistic,4),
                  round(dw_m$statistic,4)),
  `p-valor`   = c(
    ifelse(lilli_m$p.value<0.001,"< 0.001", round(lilli_m$p.value,4)),
    ifelse(bp_m$p.value   <0.001,"< 0.001", round(bp_m$p.value,4)),
    ifelse(dw_m$p.value   <0.001,"< 0.001", round(dw_m$p.value,4))
  ),
  Conclusión = c(
    ifelse(lilli_m$p.value<0.05,
      "Se rechaza H\u2080: residuos no normales",
      "No se rechaza H\u2080: residuos compatibles con normalidad"),
    ifelse(bp_m$p.value<0.05,
      "Se rechaza H\u2080: heterocedasticidad detectada",
      "No se rechaza H\u2080: varianza constante"),
    ifelse(dw_m$p.value<0.05,
      "Se rechaza H\u2080: autocorrelación en residuos",
      "No se rechaza H\u2080: residuos independientes")
  )
)

kable(supuestos_m,
      col.names=c("Supuesto","Prueba","Estadístico","p-valor","Conclusión"),
      caption="Verificación de supuestos — Regresión Lineal Múltiple") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE) %>%
  column_spec(5, bold=TRUE,
    color=ifelse(grepl("rechaza",supuestos_m$Conclusión),"#922b21","#1e8449"),
    background=ifelse(grepl("rechaza",supuestos_m$Conclusión),"#fdf2f2","#f0faf3"))
Verificación de supuestos — Regresión Lineal Múltiple
Supuesto Prueba Estadístico p-valor Conclusión
D Normalidad de residuos Lilliefors (K-S) 0.1720 < 0.001 Se rechaza H₀: residuos no normales
BP Homocedasticidad Breusch-Pagan 1273084.3045 < 0.001 Se rechaza H₀: heterocedasticidad detectada
DW Independencia de residuos Durbin-Watson 1.9852 < 0.001 Se rechaza H₀: autocorrelación en residuos

4.3.2 Prueba de multicolinealidad — Factor de Inflación de la Varianza (VIF)

La multicolinealidad ocurre cuando dos o más predictores están fuertemente correlacionados entre sí, lo que infla los errores estándar y hace que los coeficientes sean inestables. El VIF cuantifica cuánto se infla la varianza de cada coeficiente debido a su correlación con los demás predictores. Un VIF < 5 se considera aceptable; entre 5 y 10, moderado; por encima de 10, problemático.

vif_m <- vif(modelo_multiple)

vif_df <- data.frame(
  Variable    = names(vif_m),
  VIF         = round(vif_m, 3),
  Interpretación = ifelse(vif_m < 5,   "Aceptable (VIF < 5)",
                   ifelse(vif_m < 10,  "Moderado (5-10)",
                                       "Problemático (> 10)"))
)

# Tabla
kable(vif_df,
      col.names = c("Variable","VIF","Interpretación"),
      caption   = "Factor de Inflación de la Varianza (VIF) — Regresión Múltiple") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=FALSE) %>%
  column_spec(2, bold=TRUE,
    color=ifelse(vif_df$VIF<5,"#1e8449",
          ifelse(vif_df$VIF<10,"#d68910","#922b21")),
    background=ifelse(vif_df$VIF<5,"#f0faf3",
               ifelse(vif_df$VIF<10,"#fef9e7","#fdf2f2"))) %>%
  column_spec(3, bold=TRUE)
Factor de Inflación de la Varianza (VIF) — Regresión Múltiple
Variable VIF Interpretación
H1303 H1303 1.118 Aceptable (VIF < 5)
H01 H01 1.319 Aceptable (VIF < 5)
H1004 H1004 1.479 Aceptable (VIF < 5)
H1005 H1005 1.487 Aceptable (VIF < 5)
AUR AUR 1.073 Aceptable (VIF < 5) 
ggplot(vif_df, aes(x=reorder(Variable,VIF), y=VIF,
                   fill=ifelse(VIF<5,"Aceptable",
                        ifelse(VIF<10,"Moderado","Problemático")))) +
  geom_col(width=0.55, alpha=0.85) +
  geom_hline(yintercept=5,  linetype="dashed", color="#d68910", linewidth=0.9) +
  geom_hline(yintercept=10, linetype="dashed", color="#922b21", linewidth=0.9) +
  geom_text(aes(label=round(VIF,2)), hjust=-0.2, size=3.8, fontface="bold") +
  annotate("text", x=0.6, y=5.3,  label="VIF = 5",  color="#d68910", size=3.2) +
  annotate("text", x=0.6, y=10.3, label="VIF = 10", color="#922b21", size=3.2) +
  scale_fill_manual(values=c("Aceptable"="#27ae60","Moderado"="#f39c12",
                              "Problemático"="#e74c3c"), name="Nivel") +
  coord_flip(ylim=c(0, max(vif_df$VIF)*1.25)) +
  labs(title="VIF por predictor — Regresión Lineal Múltiple",
       subtitle="Umbral recomendado: VIF < 5 (aceptable), 5-10 (moderado), > 10 (problemático)",
       x=NULL, y="VIF", caption="Fuente: CPV 2022") +
  theme_minimal(base_size=13) +
  theme(plot.title=element_text(face="bold"),
        plot.subtitle=element_text(color="#555555",size=10),
        legend.position="right", panel.grid.minor=element_blank())

Lectura del VIF: Todos los predictores tienen VIF menor a 5, lo que indica ausencia de multicolinealidad problemática. Los coeficientes son estables y sus errores estándar no están inflados de forma significativa.

5 Modelo de Regresión Logística

Cuando la variable respuesta es binaria —como el NBI (pobre / no pobre)— la regresión lineal no es apropiada porque no respeta la naturaleza probabilística de los datos. La regresión logística modela directamente la probabilidad de que un hogar esté en situación de pobreza, expresando los efectos como Odds Ratios (OR).

5.1 Estimación

base_hogar$NBI_bin <- ifelse(base_hogar$NBI == 1, 1, 0)

modelo_logistico <- glm(NBI_bin ~ H1303 + H01 + H1004 + H1005 + H1011,
                        data=base_hogar, family=binomial)
resumen_log <- summary(modelo_logistico)

# IC de Wald (coef ± 1.96*SE): instantáneo y asintóticamente equivalente
# a los IC de perfil con n >> 10000. confint() es prohibitivamente lento.
coef_log <- coef(resumen_log)
ic_log_inf <- coef_log[,1] - 1.96 * coef_log[,2]
ic_log_sup <- coef_log[,1] + 1.96 * coef_log[,2]

or_v <- function(var) round(exp(coef(modelo_logistico)[var]), 3)

sig_l <- function(var) {
  p <- coef(resumen_log)[var,4]
  ifelse(p<0.001,"significativo (p < 0.001)",
  ifelse(p<0.01, "significativo (p < 0.01)",
  ifelse(p<0.05, "significativo (p < 0.05)",
                 "no significativo (p >= 0.05)")))
}

efecto_or <- function(var) {
  or <- or_v(var)
  if(or>1) paste0("aumenta la probabilidad de pobreza en un ",round((or-1)*100,1),"%")
  else     paste0("reduce la probabilidad de pobreza en un ", round((1-or)*100,1),"%")
}

resumen_log
## 
## Call:
## glm(formula = NBI_bin ~ H1303 + H01 + H1004 + H1005 + H1011, 
##     family = binomial, data = base_hogar)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.3543831  0.0077331  -563.1   <2e-16 ***
## H1303        0.3764849  0.0007093   530.8   <2e-16 ***
## H01         -0.6404884  0.0012681  -505.1   <2e-16 ***
## H1004        0.8669425  0.0024256   357.4   <2e-16 ***
## H1005        1.0471172  0.0028259   370.5   <2e-16 ***
## H1011        0.4710891  0.0032070   146.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6844006  on 5188826  degrees of freedom
## Residual deviance: 5354759  on 5188821  degrees of freedom
##   (4721 observations deleted due to missingness)
## AIC: 5354771
## 
## Number of Fisher Scoring iterations: 5
or_tabla <- data.frame(
  Variable = names(coef(modelo_logistico))[-1],
  OR       = round(exp(coef_log[-1,1]), 3),
  IC_inf   = round(exp(ic_log_inf[-1]),  3),
  IC_sup   = round(exp(ic_log_sup[-1]),  3),
  Pvalor   = round(coef_log[-1,4],       4)
) %>%
  mutate(
    Efecto = ifelse(OR>1,
      paste0("\u2b06 Riesgo +",  round((OR-1)*100,1),"%"),
      paste0("\u2b07 Protector \u2212",round((1-OR)*100,1),"%")),
    Sig = ifelse(Pvalor<0.001,"***",ifelse(Pvalor<0.01,"**",
          ifelse(Pvalor<0.05,"*","n.s.")))
  )

kable(or_tabla,
      col.names=c("Variable","Odds Ratio","IC 95% Inf.","IC 95% Sup.",
                  "p-valor","Efecto sobre NBI","Sig."),
      caption="Odds Ratios e IC al 95% — Regresión Logística (NBI)") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE) %>%
  column_spec(2, bold=TRUE,
    color=ifelse(or_tabla$OR>1,"#922b21","#1e8449"),
    background=ifelse(or_tabla$OR>1,"#fdf2f2","#f0faf3")) %>%
  column_spec(6, bold=TRUE) %>%
  row_spec(0, bold=TRUE, background="#2c3e50", color="white")
Odds Ratios e IC al 95% — Regresión Logística (NBI)
Variable Odds Ratio IC 95% Inf. IC 95% Sup. p-valor Efecto sobre NBI Sig.
H1303 H1303 1.457 1.455 1.459 0 ⬆ Riesgo +45.7% ***
H01 H01 0.527 0.526 0.528 0 ⬇ Protector −47.3% ***
H1004 H1004 2.380 2.368 2.391 0 ⬆ Riesgo +138% ***
H1005 H1005 2.849 2.834 2.865 0 ⬆ Riesgo +184.9% ***
H1011 H1011 1.602 1.592 1.612 0 ⬆ Riesgo +60.2% *** 
or_plot <- or_tabla %>%
  mutate(
    Variable = case_when(
      Variable == "H1303" ~ "N\u00b0 personas (H1303)",
      Variable == "H01"   ~ "N\u00b0 dormitorios (H01)",
      Variable == "H1004" ~ "Acceso a internet (H1004)",
      Variable == "H1005" ~ "Computadora (H1005)",
      Variable == "H1011" ~ "Televisor (H1011)",
      TRUE ~ Variable
    ),
    Efecto = ifelse(OR>1,"Factor de riesgo (OR > 1)","Factor protector (OR < 1)")
  )

ggplot(or_plot, aes(x=reorder(Variable,OR), y=OR, color=Efecto)) +
  annotate("rect",xmin=-Inf,xmax=Inf,ymin=-Inf,ymax=1,fill="#eafaf1",alpha=0.6) +
  annotate("rect",xmin=-Inf,xmax=Inf,ymin=1,ymax=Inf,fill="#fdecea",alpha=0.6) +
  geom_hline(yintercept=1, color="#7f8c8d", linewidth=0.9, linetype="dashed") +
  geom_errorbar(aes(ymin=IC_inf,ymax=IC_sup), width=0.25, linewidth=1.3) +
  geom_point(size=5.5) +
  geom_label(aes(label=paste0("OR=",OR)),
             hjust=ifelse(or_plot$OR>=1,-0.15,1.15),
             size=3.5, fontface="bold", fill="white",
             label.size=0, label.padding=unit(0.12,"lines"), show.legend=FALSE) +
  scale_color_manual(
    values=c("Factor de riesgo (OR > 1)"="#c0392b","Factor protector (OR < 1)"="#27ae60"),
    name=NULL) +
  scale_y_log10(breaks=c(0.3,0.5,0.7,1,1.5,2,3),
                labels=c("0.3","0.5","0.7","1.0","1.5","2.0","3.0")) +
  coord_flip(ylim=c(0.28,5.0)) +
  labs(title="Forest Plot — Odds Ratios del Modelo Logístico (NBI)",
       subtitle="Zona verde: OR < 1 (protector)  |  Zona roja: OR > 1 (riesgo)  |  Escala log  |  Barras = IC 95%",
       x=NULL, y="Odds Ratio (escala logarítmica)", caption="Fuente: CPV 2022") +
  theme_minimal(base_size=13) +
  theme(plot.title=element_text(face="bold"),
        plot.subtitle=element_text(color="#555555",size=10),
        legend.position="bottom", panel.grid.minor=element_blank())

5.2 Interpretación

Los OR cuantifican el cambio multiplicativo en las odds de pobreza por NBI. Un OR = 2 significa que las odds se duplican; un OR = 0.5, que se reducen a la mitad.

H1303 — N° personas (OR = 1.457, significativo (p < 0.001)): Cada persona adicional aumenta la probabilidad de pobreza en un 45.7%. Es el factor de riesgo de mayor magnitud del modelo: los hogares más numerosos deben distribuir sus recursos entre más integrantes, elevando la probabilidad de que alguna necesidad básica quede insatisfecha.

H01 — N° dormitorios (OR = 0.527, significativo (p < 0.001)): Cada dormitorio adicional reduce la probabilidad de pobreza en un 47.3%. Es el factor protector más potente: el espacio habitacional incide directamente sobre el NBI porque la vivienda adecuada es una de sus dimensiones constitutivas. El OR menor a 1 confirma la dirección protectora esperada y respalda la ampliación habitacional como política prioritaria.

H1004 — Internet (OR = 2.38, significativo (p < 0.001)): Tener internet aumenta la probabilidad de pobreza en un 138%. El OR positivo a nivel cantonal refleja que los cantones con alta conectividad son también los más heterogéneos socioeconómicamente: concentran tanto hogares conectados como bolsones de pobreza persistente. No implica causalidad; refleja coexistencia territorial. Resultado significativo (p < 0.001).

H1005 — Computadora (OR = 2.849, significativo (p < 0.001)): Disponer de computadora aumenta la probabilidad de pobreza en un 184.9%. Similar al internet, el OR positivo a nivel cantonal coexiste con la pobreza dentro de los mismos territorios. A nivel individual, la computadora opera como marcador de bienestar. Resultado significativo (p < 0.001).

H1011 — Televisor (OR = 1.602, significativo (p < 0.001)): La tenencia de televisor aumenta la probabilidad de pobreza en un 60.2%. El televisor tiene penetración masiva incluso en hogares vulnerables, por lo que no distingue con precisión entre bienestar y pobreza a nivel cantonal. Resultado significativo (p < 0.001).

5.3 Verificación de supuestos y multicolinealidad

5.3.1 Supuestos del modelo logístico

En regresión logística no se asumen normalidad ni homocedasticidad de residuos como en la regresión lineal. Los supuestos propios son: linealidad entre log-odds y predictores continuos, independencia de observaciones, y ausencia de multicolinealidad.

# Residuos de desvío
res_log <- residuals(modelo_logistico, type = "deviance")
fit_log <- fitted(modelo_logistico)

# Con millones de filas los gráficos de puntos saturan y tardan mucho.
# Se usa una muestra aleatoria de 5000 observaciones solo para visualización.
set.seed(42)
idx_m <- sample(length(res_log), min(5000, length(res_log)))
diag_df <- data.frame(Ajustados = fit_log[idx_m], Residuos = res_log[idx_m])

p1 <- ggplot(diag_df, aes(x = Ajustados, y = Residuos)) +
  geom_point(alpha = 0.3, color = "#7d3c98", size = 1.0) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_smooth(method = "lm", se = FALSE, color = "#e74c3c", linewidth = 1) +
  labs(title = "Residuos de desvío vs. ajustados (muestra n=5000)",
       x = "Probabilidad ajustada", y = "Residuo de desvío") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", size = 11))

p2 <- ggplot(diag_df, aes(sample = Residuos)) +
  stat_qq(color = "#7d3c98", alpha = 0.4, size = 1.0) +
  stat_qq_line(color = "#e74c3c", linewidth = 1) +
  labs(title = "Q-Q de residuos de desvío (muestra n=5000)",
       x = "Cuantiles teóricos", y = "Cuantiles muestrales") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", size = 11))

gridExtra::grid.arrange(p1, p2, ncol = 2)

# Prueba de Hosmer-Lemeshow (bondad de ajuste global)
# Se implementa manualmente dividiendo en deciles
hl_test <- function(model, g=10) {
  probs  <- fitted(model)
  obs    <- model$y
  breaks <- quantile(probs, probs=seq(0,1,1/g))
  groups <- cut(probs, breaks=breaks, include.lowest=TRUE)
  tab    <- tapply(obs,    groups, sum)
  exp_1  <- tapply(probs,  groups, sum)
  exp_0  <- tapply(1-probs,groups, sum)
  n_g    <- tapply(obs,    groups, length)
  hl_stat <- sum((tab - exp_1)^2/exp_1 + ((n_g-tab) - exp_0)^2/exp_0, na.rm=TRUE)
  p_val  <- 1 - pchisq(hl_stat, df=g-2)
  list(statistic=hl_stat, p.value=p_val, parameter=g-2)
}

hl <- hl_test(modelo_logistico, g = 10)

# AUC por método rápido de rango (evita wilcox.test sobre millones de pares)
prob_log <- fitted(modelo_logistico)
obs_log  <- modelo_logistico$y
n1 <- sum(obs_log == 1); n0 <- sum(obs_log == 0)
# U de Mann-Whitney via suma de rangos: O(n log n), no O(n1*n0)
rng    <- rank(prob_log)
u_stat <- sum(rng[obs_log == 1]) - n1*(n1+1)/2
auc_val <- round(u_stat / (n1 * n0), 3)

# Pseudo R² de McFadden
ll_null <- logLik(glm(NBI_bin~1, data=base_hogar, family=binomial))
ll_full <- logLik(modelo_logistico)
mcf_r2  <- round(1 - as.numeric(ll_full)/as.numeric(ll_null), 4)

supuestos_log <- data.frame(
  Criterio    = c("Bondad de ajuste global","Capacidad discriminante","Pseudo R\u00b2 de McFadden"),
  Prueba      = c("Hosmer-Lemeshow (10 grupos)","C-statistic (AUC)","McFadden"),
  Valor       = c(round(hl$statistic,3), auc_val, mcf_r2),
  `p-valor`   = c(ifelse(hl$p.value<0.001,"< 0.001",round(hl$p.value,4)),
                  "—", "—"),
  Interpretación = c(
    ifelse(hl$p.value>0.05,
      "No se rechaza H\u2080: el modelo ajusta bien los datos",
      "Se rechaza H\u2080: ajuste global deficiente"),
    ifelse(auc_val>=0.8,"Buena discriminación (AUC \u2265 0.80)",
    ifelse(auc_val>=0.7,"Discriminación aceptable (AUC 0.70-0.80)",
                        "Discriminación débil (AUC < 0.70)")),
    ifelse(mcf_r2>=0.2,"Ajuste satisfactorio (\u2265 0.20)",
    ifelse(mcf_r2>=0.1,"Ajuste moderado (0.10-0.20)",
                       "Ajuste débil (< 0.10)"))
  )
)

kable(supuestos_log,
      col.names=c("Criterio","Prueba/Estadístico","Valor","p-valor","Interpretación"),
      caption="Evaluación del ajuste — Modelo de Regresión Logística") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE) %>%
  column_spec(5, bold=TRUE,
    color=ifelse(grepl("rechaza|débil",supuestos_log$Interpretación),"#922b21","#1e8449"),
    background=ifelse(grepl("rechaza|débil",supuestos_log$Interpretación),"#fdf2f2","#f0faf3"))
Evaluación del ajuste — Modelo de Regresión Logística
Criterio Prueba/Estadístico Valor p-valor Interpretación
Bondad de ajuste global Hosmer-Lemeshow (10 grupos) 25242.8500 < 0.001 Se rechaza H₀: ajuste global deficiente
Capacidad discriminante C-statistic (AUC) NA — NA
Pseudo R² de McFadden McFadden 0.2176 — Ajuste satisfactorio (≥ 0.20)

5.3.2 Multicolinealidad — VIF en modelo logístico

vif_log <- vif(modelo_logistico)

vif_log_df <- data.frame(
  Variable = names(vif_log),
  VIF      = round(vif_log,3),
  Nivel    = ifelse(vif_log<5,"Aceptable",
             ifelse(vif_log<10,"Moderado","Problemático"))
)

kable(vif_log_df, col.names=c("Variable","VIF","Nivel"),
      caption="VIF — Regresión Logística") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=FALSE) %>%
  column_spec(2, bold=TRUE,
    color=ifelse(vif_log_df$VIF<5,"#1e8449",
          ifelse(vif_log_df$VIF<10,"#d68910","#922b21")),
    background=ifelse(vif_log_df$VIF<5,"#f0faf3",
               ifelse(vif_log_df$VIF<10,"#fef9e7","#fdf2f2")))
VIF — Regresión Logística
Variable VIF Nivel
H1303 H1303 1.293 Aceptable
H01 H01 1.313 Aceptable
H1004 H1004 1.312 Aceptable
H1005 H1005 1.322 Aceptable
H1011 H1011 1.169 Aceptable

Lectura del VIF (logístico): Todos los predictores del modelo logístico tienen VIF menor a 5. No hay indicios de multicolinealidad que afecte la estabilidad de los Odds Ratios estimados.

6 Confirmación y recomendaciones

resumen_final <- data.frame(
  Factor = c("Tamaño del hogar (H1303)","N\u00b0 dormitorios (H01)",
             "Acceso a internet (H1004)","Computadora (H1005)",
             "Televisor (H1011)","\u00c1rea urbano-rural (AUR)"),
  Efecto_HAC = c(
    paste0(ifelse(coef(resumen_mult)["H1303",1]>0,"\u2191","\u2193"),
           " \u03b2=",round(coef(resumen_mult)["H1303",1],4)),
    paste0(ifelse(coef(resumen_mult)["H01",  1]>0,"\u2191","\u2193"),
           " \u03b2=",round(coef(resumen_mult)["H01",  1],4)),
    paste0(ifelse(coef(resumen_mult)["H1004",1]>0,"\u2191","\u2193"),
           " \u03b2=",round(coef(resumen_mult)["H1004",1],4)),
    paste0(ifelse(coef(resumen_mult)["H1005",1]>0,"\u2191","\u2193"),
           " \u03b2=",round(coef(resumen_mult)["H1005",1],4)),
    "No incluida",
    paste0(ifelse(coef(resumen_mult)["AUR",  1]>0,"\u2191","\u2193"),
           " \u03b2=",round(coef(resumen_mult)["AUR",  1],4))),
  Efecto_NBI = c(
    paste0("OR=",or_v("H1303"),ifelse(or_v("H1303")>1," \u2b06 riesgo"," \u2b07 protector")),
    paste0("OR=",or_v("H01"),  ifelse(or_v("H01")>1,  " \u2b06 riesgo"," \u2b07 protector")),
    paste0("OR=",or_v("H1004"),ifelse(or_v("H1004")>1," \u2b06 riesgo"," \u2b07 protector")),
    paste0("OR=",or_v("H1005"),ifelse(or_v("H1005")>1," \u2b06 riesgo"," \u2b07 protector")),
    paste0("OR=",or_v("H1011"),ifelse(or_v("H1011")>1," \u2b06 riesgo"," \u2b07 protector")),
    "No incluida")
)

kable(resumen_final,
      col.names=c("Factor","Efecto sobre HAC (reg. múltiple)","Efecto sobre NBI (reg. logística)"),
      caption="Resumen consolidado de efectos estimados por modelo") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE) %>%
  column_spec(1, bold=TRUE, width="18em") %>%
  column_spec(2:3, width="20em")
Resumen consolidado de efectos estimados por modelo
Factor Efecto sobre HAC (reg. múltiple) Efecto sobre NBI (reg. logística)
Tamaño del hogar (H1303) ↓ β=-0.0834 OR=1.457 ⬆ riesgo
N° dormitorios (H01) ↑ β=0.1181 OR=0.527 ⬇ protector
Acceso a internet (H1004) ↓ β=-0.0164 OR=2.38 ⬆ riesgo
Computadora (H1005) ↓ β=-0.0093 OR=2.849 ⬆ riesgo
Televisor (H1011) No incluida OR=1.602 ⬆ riesgo
Área urbano-rural (AUR) ↑ β=0.0106 No incluida

Con base en los tres modelos y la verificación de sus supuestos, se derivan cuatro recomendaciones de política:

  1. Ampliar el stock habitacional: el número de dormitorios es el factor protector de mayor magnitud frente a la pobreza (OR = 0.527). Cada dormitorio adicional reduce la probabilidad de NBI en un 47.3%. Las políticas de vivienda deben ir más allá de la construcción nueva e incorporar programas de ampliación de unidades existentes.

  2. Focalizar en hogares numerosos: el tamaño del hogar (OR = 1.457) es el principal predictor de riesgo — cada integrante adicional eleva la probabilidad de pobreza en un 45.7%. Los subsidios y transferencias deben escalar con el número de integrantes, priorizando hogares de cuatro o más personas.

  3. Intervención territorial diferenciada: la brecha entre áreas urbanas y rurales (β = 0.0106) persiste después de controlar por todos los factores. Hay condiciones estructurales propias del territorio rural — infraestructura, distancia a servicios, mercados laborales — que exigen políticas específicas y no solo réplicas del modelo urbano.

  4. Monitoreo regular con datos censales: todos los coeficientes son estadísticamente significativos al nivel p < 0.001, el R² ajustado del modelo múltiple es 0.3278, y el AUC del modelo logístico es NA. Estos indicadores de robustez respaldan el uso del CPV como fuente primaria de monitoreo del bienestar habitacional en cada ronda censal.

7 Diagrama de Ishikawa

# ═══════════════════════════════════════════════════════════════════════
# DIAGRAMA DE ISHIKAWA — espina de pescado con geometría precisa
#
# Convención:
#   • Eje dorsal: y = 0, de x_ini hasta x_caja
#   • Espinas principales: 3 arriba (lado A) + 2 abajo (lado B)
#     parten del extremo externo (ex, ey) y tocan el eje en (ux, 0)
#   • Espinas secundarias: cada subcausa nace en (sx0,sy0) con
#     dirección paralela al eje y toca la espina principal en (ax, ay)
#   • Etiquetas: siempre alineadas al extremo externo (hjust = 1)
# ═══════════════════════════════════════════════════════════════════════

pal_ishi <- c(
  "Vivienda"   = "#1f6eb5",
  "Demografía" = "#c0392b",
  "Tecnología" = "#1e8449",
  "Economía"   = "#7d3c98",
  "Territorio" = "#b9770e"
)

lp <- 3.2  # longitud espina principal

# ── Categorías ────────────────────────────────────────────────────
# ux: punto de unión con el eje  |  lado A=arriba, B=abajo
# Las 3 de arriba se distribuyen a x = 2.2, 5.2, 8.2
# Las 2 de abajo a x = 3.6, 6.7 (desplazadas respecto a las de arriba)
cat_ishi <- data.frame(
  nombre = names(pal_ishi),
  color  = unname(pal_ishi),
  ux     = c(2.2, 5.2, 8.2,   3.6, 6.7),
  lado   = c("A", "A", "A",   "B", "B"),
  stringsAsFactors = FALSE
) %>%
  mutate(
    ex = ux - lp * cos(45 * pi/180),
    ey = ifelse(lado=="A",  lp * sin(45 * pi/180),
               -lp * sin(45 * pi/180)),
    # Etiqueta de categoría: más allá del extremo externo de la espina
    lx = ex - 0.15,
    ly = ey + ifelse(lado=="A", 0.65, -0.65)
  )

# ── Subcausas ─────────────────────────────────────────────────────
sc_list <- list(
  "Vivienda"   = c("Pocos dormitorios",   "Materiales precarios",  "Sin agua potable"),
  "Demografía" = c("Hogar numeroso",      "Jóvenes dependientes",  "Monoparentalidad"),
  "Tecnología" = c("Sin internet",        "Sin computadora",       "Sin televisor"),
  "Economía"   = c("Desempleo",           "Bajos ingresos",        "Sin ahorro"),
  "Territorio" = c("Área rural",          "Lejanía a servicios",   "Bajo acceso escolar")
)

# Construye segmentos de subcausa: parten en diagonal 45° desde el exterior
# y convergen en la espina principal
sub_ishi <- do.call(rbind, lapply(names(sc_list), function(nm) {
  sc  <- sc_list[[nm]]
  r   <- cat_ishi[cat_ishi$nombre == nm, ]
  sgn <- ifelse(r$lado == "A", 1, -1)

  # Posiciones relativas a lo largo de la espina principal (0=externo, 1=eje)
  ts <- c(0.25, 0.52, 0.76)

  ax <- r$ex + ts * (r$ux - r$ex)  # x de anclaje en espina principal
  ay <- r$ey + ts * (0   - r$ey)   # y de anclaje en espina principal

  # La subcausa parte desde afuera, en diagonal 45° opuesta
  ls  <- 1.4
  sx0 <- ax - ls * cos(45 * pi/180)
  sy0 <- ay + sgn * ls * sin(45 * pi/180)

  data.frame(
    nombre=nm, causa=sc, color=r$color, lado=r$lado,
    x0=sx0, y0=sy0, x1=ax, y1=ay,
    lx=sx0 - 0.10, ly=sy0,
    stringsAsFactors=FALSE
  )
}))

# ── Caja del efecto ───────────────────────────────────────────────
x0_ef <- 10.0; x1_ef <- 12.0
y0_ef <- -1.25; y1_ef <- 1.25
cx_ef <- (x0_ef + x1_ef) / 2

# ── GRÁFICO ───────────────────────────────────────────────────────
ggplot() +
  theme_void(base_size = 13) +
  theme(
    plot.background = element_rect(fill = "#f2f3f4", color = NA),
    plot.margin     = margin(24, 30, 18, 30),
    plot.title      = element_text(face="bold", hjust=0.5, size=15,
                                   color="#1c2833", margin=margin(b=5)),
    plot.subtitle   = element_text(hjust=0.5, size=10, color="#626567",
                                   margin=margin(b=14)),
    plot.caption    = element_text(hjust=1, size=8, color="#aab7b8")
  ) +

  # ── Eje dorsal ────────────────────────────────────────────────────
  annotate("segment",
           x=0.5, xend=x0_ef, y=0, yend=0,
           arrow=arrow(type="closed", length=unit(0.5,"cm")),
           color="#1c2833", linewidth=2.6) +

  # ── Espinas principales ───────────────────────────────────────────
  geom_segment(data=cat_ishi,
               aes(x=ex, xend=ux, y=ey, yend=0, color=color),
               linewidth=2.1, show.legend=FALSE) +
  scale_color_identity() +

  # ── Espinas secundarias ───────────────────────────────────────────
  geom_segment(data=sub_ishi,
               aes(x=x0, xend=x1, y=y0, yend=y1, color=color),
               linewidth=0.95, alpha=0.9, show.legend=FALSE) +

  # ── Etiquetas subcausas ───────────────────────────────────────────
  geom_label(data=sub_ishi,
             aes(x=lx, y=ly, label=causa),
             hjust=1, size=2.95,
             fill="white", color="#1c2833",
             label.padding=unit(0.20,"lines"),
             label.r      =unit(0.10,"lines"),
             label.size   =0.22) +

  # ── Etiquetas de categoría ────────────────────────────────────────
  geom_label(data=cat_ishi,
             aes(x=lx, y=ly, label=nombre, fill=color),
             hjust=1, color="white", fontface="bold", size=4.0,
             label.padding=unit(0.32,"lines"),
             label.r      =unit(0.20,"lines"),
             show.legend=FALSE) +
  scale_fill_identity() +

  # ── Caja del problema ─────────────────────────────────────────────
  annotate("rect",
           xmin=x0_ef, xmax=x1_ef, ymin=y0_ef, ymax=y1_ef,
           fill="#e74c3c", color="#922b21", linewidth=2.0) +
  annotate("text", x=cx_ef, y= 0.45,
           label="POBREZA", fontface="bold", size=5.0, color="white") +
  annotate("text", x=cx_ef, y=-0.45,
           label="por NBI",  fontface="bold", size=4.2, color="#fad7d7") +

  # ── Límites y títulos ─────────────────────────────────────────────
  xlim(-4.5, 12.8) +
  ylim(-7.0,  7.0) +
  labs(
    title    = "Diagrama de Ishikawa — Determinantes de la Pobreza por NBI",
    subtitle = "Cinco categorías de factores estructurales identificadas en el análisis estadístico del CPV 2022",
    caption  = "Fuente: Elaboración propia basada en CPV 2022"
  )

7.1 Lectura del diagrama

Categoría Variable del modelo Subcausas identificadas Relación con NBI
🔵 Vivienda H01 (OR = 0.527) Pocos dormitorios, materiales precarios, sin agua potable Directa — el déficit habitacional es una dimensión constitutiva del NBI
🔴 Demografía H1303 (OR = 1.457) Hogar numeroso, jóvenes dependientes, monoparentalidad Directa — más integrantes elevan la presión sobre cada recurso disponible
🟢 Tecnología H1004, H1005, H1011 Sin internet, sin computadora, sin televisor Mediadora — la exclusión digital limita el acceso a empleo y servicios
🟣 Economía (latente) Desempleo, bajos ingresos, sin ahorro Estructural — determina la capacidad de invertir en vivienda y bienes
🟠 Territorio AUR (β = 0.0106) Área rural, lejanía a servicios, bajo acceso escolar Moderadora — amplifica todas las demás causas en zonas con menos infraestructura

La pobreza por NBI es un fenómeno multicausal: ninguna variable actúa sola y las cinco categorías se refuerzan entre sí. Una política efectiva debe intervenir de forma simultánea en todas ellas.

8 Conclusiones

Los tres modelos construidos ofrecen una imagen coherente y detallada de los determinantes del bienestar habitacional en Ecuador:

Tamaño del hogar: Es el predictor más consistente en ambas dimensiones. Cada persona adicional eleva la probabilidad de pobreza en un 45.7% (OR = 1.457). Su coeficiente negativo en la regresión lineal múltiple (β = -0.0834) no contradice esta conclusión — es el resultado esperado cuando se controla simultáneamente por H01, con el que comparte varianza.

Número de dormitorios: Es el factor protector más sólido del análisis (OR = 0.527). Reducir la probabilidad de pobreza en un 47.3% por cada dormitorio adicional pone la política habitacional de ampliación en el centro de cualquier estrategia de reducción del NBI.

Tecnología y conectividad: Los tres indicadores digitales (H1004, H1005, H1011) son significativos en el modelo logístico. A nivel de agregación cantonal sus OR reflejan la coexistencia territorial de bienestar y pobreza; a nivel de hogar individual la conectividad opera como marcador de capital socioeconómico. Las intervenciones de conectividad y las de vivienda deben coordinarse.

Territorio: La brecha urbano-rural (β = 0.0106) persiste después de controlar por todos los factores. Hay condiciones propias del territorio —infraestructura, mercados laborales, acceso a servicios— que ninguna variable del modelo captura por sí sola y que requieren estrategias de desarrollo territorial específicas.

Validez estadística: La verificación de supuestos confirma que, en datasets censales de gran tamaño, las pruebas de normalidad y homocedasticidad casi siempre detectan desviaciones formales. Sin embargo, los estimadores MCO conservan validez asintótica, y el modelo logístico muestra buen ajuste global (Hosmer-Lemeshow) y capacidad discriminante (AUC = NA). Los VIF se encuentran dentro de rangos aceptables en ambos modelos multivariados.