Taller: Modelo de Regresión Logística

Consumo Energético – Universidad Surcolombiana (USCO) – Especialización en Estadística

Author

Sergio Andres Beltran, Juan Pablo Donato

Published

April 9, 2026

1 Descripción del Problema e Hipótesis

La Sede Central de la Universidad Surcolombiana (USCO) cuenta con cinco cuentas eléctricas activas (Central 1 a Central 5). La normativa eléctrica colombiana (CREG 015 de 2018) establece que la relación entre energía reactiva y energía activa no debe superar 0.48. El incumplimiento de esta norma genera recargos económicos directos en la factura eléctrica.

El análisis se aborda en dos niveles:

  • Nivel Sede: ¿El consumo total de energía activa de la sede predice si alguna cuenta incumplirá la norma ese mes?
  • Nivel Cuenta: ¿El consumo de energía activa de una cuenta individual predice si esa cuenta cumplirá o no la norma?
NoteHipótesis del Estudio

Modelo Simple — Nivel Sede:

  • \(H_0: \beta_1 = 0\) → La energía activa total NO predice el incumplimiento normativo de la sede.
  • \(H_1: \beta_1 \neq 0\) → La energía activa total SÍ predice el incumplimiento normativo.

Modelo Múltiple — Nivel Sede:

  • \(H_0: \beta_1 = \beta_2 = 0\) → Ni la energía activa ni el mes del año explican el incumplimiento.
  • \(H_1\): Al menos un \(\beta_k \neq 0\) → Alguna variable predice el incumplimiento.

Modelos por Cuenta Individual:

  • \(H_0: \beta_1 = 0\) → El consumo activo de la cuenta NO predice su propio cumplimiento.
  • \(H_1: \beta_1 \neq 0\) → El consumo activo SÍ predice el cumplimiento de la cuenta.

2 Descripción de Variables

Ver código 
# ── Ruta base — ajusta según tu equipo ───────────────────────────────────────
RUTA_BASE <- "C:/Users/ESTUDIANTES/Documents/Working Holydays/"

# ── Cargar las dos tablas exportadas por el script de unificación ─────────────
dt_cuenta <- read.csv(paste0(RUTA_BASE, "Consumo_Central.csv"),
                      header = TRUE, sep = ",")
dt_sede   <- read.csv(paste0(RUTA_BASE, "Resumen_Sede.csv"),
                      header = TRUE, sep = ",")

# ── Conversiones de tipo necesarias para los modelos ─────────────────────────
dt_cuenta$Cuenta       <- factor(dt_cuenta$Cuenta)
dt_cuenta$Nombre_Cuenta<- factor(dt_cuenta$Nombre_Cuenta)
dt_cuenta$Cumple       <- factor(dt_cuenta$Cumple, levels = c("Sí", "No"))
dt_cuenta$Cumple_bin   <- as.integer(dt_cuenta$Cumple_bin)
dt_cuenta$Semestre     <- factor(dt_cuenta$Semestre,  labels = c("Ene–Jun", "Jul–Dic"))
dt_cuenta$Trimestre    <- factor(dt_cuenta$Trimestre,
                                  labels = c("T1 (Ene-Mar)", "T2 (Abr-Jun)",
                                             "T3 (Jul-Sep)", "T4 (Oct-Dic)"))

dt_sede$Sede_Incumple  <- as.integer(dt_sede$Sede_Incumple)

# ── Tabla resumen de datos cargados ─────────────────────────────────────
resumen_datos <- data.frame(
  Concepto = c("Registros cuenta-período (Tabla 1)",
               "Registros sede-período (Tabla 2)",
               "Cuentas activas",
               "Período cubierto",
               "Períodos con incumplimiento (sede)"),
  Valor = c(as.character(nrow(dt_cuenta)),
            as.character(nrow(dt_sede)),
            as.character(length(unique(dt_cuenta$Cuenta))),
            paste(min(dt_cuenta$Ano), "–", max(dt_cuenta$Ano)),
            as.character(sum(dt_sede$Sede_Incumple)))
)

apa_table(resumen_datos,
          col.names = c("Concepto", "Valor"),
          caption   = "Resumen de los Datos Cargados",
          align     = c("l", "r"))
Resumen de los Datos Cargados
Concepto Valor
Registros cuenta-período (Tabla 1) 1112
Registros sede-período (Tabla 2) 291
Cuentas activas 5
Período cubierto 2001 – 2025
Períodos con incumplimiento (sede) 142
Ver código 
apa_table(
  head(dt_cuenta[, c("Ano", "Mes", "Cuenta", "Nombre_Cuenta",
                     "Activa", "Reactiva", "Relacion", "Cumple", "Pago")], 8),
  caption = "Primeros 8 registros — Consumo por Cuenta (Tabla 1)",
  align   = c("c", "c", "c", "l", "r", "r", "r", "c", "r")
)
Primeros 8 registros — Consumo por Cuenta (Tabla 1)
Ano Mes Cuenta Nombre_Cuenta Activa Reactiva Relacion Cumple Pago
2001 5 167382131 Central 1 26000 960 0.0369231 Sí 5689170
2001 6 167382131 Central 1 13840 160 0.0115607 Sí 3054810
2001 7 167382131 Central 1 23520 400 0.0170068 Sí 5378590
2001 8 167382131 Central 1 22960 640 0.0278746 Sí 5263020
2001 9 167382131 Central 1 27600 1200 0.0434783 Sí 6507590
2001 10 167382131 Central 1 23840 720 0.0302013 Sí 5639490
2001 11 167382131 Central 1 22000 1040 0.0472727 Sí 4978020
2001 12 167382131 Central 1 16480 560 0.0339806 Sí 3663650
Ver código 
apa_table(
  head(dt_sede[, c("Ano", "Mes", "Tiempo", "Activa_Total",
                   "Cuentas_Pagadas", "Cuentas_Incumplen", "Sede_Incumple")], 8),
  caption = "Primeros 8 registros — Resumen Sede (Tabla 2)",
  align   = c("c", "c", "c", "r", "c", "c", "c")
)
Primeros 8 registros — Resumen Sede (Tabla 2)
Ano Mes Tiempo Activa_Total Cuentas_Pagadas Cuentas_Incumplen Sede_Incumple
2001 5 1 110240 3 2 1
2001 6 2 59080 3 1 1
2001 7 3 98400 3 1 1
2001 8 4 103720 3 0 0
2001 9 5 125360 3 0 0
2001 10 6 100560 3 1 1
2001 11 7 104320 3 0 0
2001 12 8 62000 3 0 0
Ver código 
tabla_vars <- data.frame(
  Tabla     = c("Tabla 2", "Tabla 2", "Tabla 2", "Tabla 2",
                "Tabla 1", "Tabla 1", "Tabla 1", "Tabla 1", "Tabla 1"),
  Variable  = c("Sede_Incumple", "Activa_Total", "Cuentas_Incumplen",
                "Mes / Semestre",
                "Cumple_bin", "Activa", "Cuenta",
                "Semestre", "Pago"),
  Tipo      = c("Binaria (0/1)", "Continua", "Discreta (0–5)",
                "Ordinal / Factor",
                "Binaria (0/1)", "Continua", "Categórica (5 niveles)",
                "Factor (2 niveles)", "Continua"),
  Rol       = c("**Y — Modelo Sede**", "X principal — Modelo Sede",
                "Descriptiva", "Covariable — Modelo Sede",
                "**Y — Modelo por Cuenta**", "X principal — Modelo Cuenta",
                "Covariable — Modelo Unificado",
                "Covariable — Modelo Cuenta", "Descriptiva"),
  Descripcion = c(
    "1 = al menos una cuenta incumple ese mes; 0 = todas cumplen",
    "Suma de energía activa de las 5 cuentas (kWh)",
    "Número de cuentas con Reactiva/Activa > 0.48 ese mes",
    "Mes del año (1–12) o agrupado en semestres",
    "1 = cuenta cumple norma ese mes; 0 = incumple",
    "Energía activa individual de la cuenta (kWh)",
    "Identificador de cuenta (Central 1 a Central 5)",
    "1er semestre (Ene–Jun) vs 2do semestre (Jul–Dic)",
    "Valor pagado por la cuenta ese período (COP)"
  )
)

apa_table(tabla_vars,
         col.names = c("Tabla", "Variable", "Tipo", "Rol", "Descripción"),
         caption   = "Diccionario de Variables del Estudio",
         align     = c("l", "l", "l", "l", "l"))
Diccionario de Variables del Estudio
Tabla Variable Tipo Rol Descripción
Tabla 2 Sede_Incumple Binaria (0/1) **Y — Modelo Sede** 1 = al menos una cuenta incumple ese mes; 0 = todas cumplen
Tabla 2 Activa_Total Continua X principal — Modelo Sede Suma de energía activa de las 5 cuentas (kWh)
Tabla 2 Cuentas_Incumplen Discreta (0–5) Descriptiva Número de cuentas con Reactiva/Activa > 0.48 ese mes
Tabla 2 Mes / Semestre Ordinal / Factor Covariable — Modelo Sede Mes del año (1–12) o agrupado en semestres
Tabla 1 Cumple_bin Binaria (0/1) **Y — Modelo por Cuenta** 1 = cuenta cumple norma ese mes; 0 = incumple
Tabla 1 Activa Continua X principal — Modelo Cuenta Energía activa individual de la cuenta (kWh)
Tabla 1 Cuenta Categórica (5 niveles) Covariable — Modelo Unificado Identificador de cuenta (Central 1 a Central 5)
Tabla 1 Semestre Factor (2 niveles) Covariable — Modelo Cuenta 1er semestre (Ene–Jun) vs 2do semestre (Jul–Dic)
Tabla 1 Pago Continua Descriptiva Valor pagado por la cuenta ese período (COP)
WarningVariables excluidas del modelo — Razón metodológica

Relacion y Reactiva no se incluyen como predictoras porque Cumple se define directamente como Relacion ≤ 0.48, donde Relacion = Reactiva / Activa. Incluirlas sería tautológico: el modelo estaría explicando Y con la misma información usada para construir Y, generando coeficientes artificialmente perfectos sin valor predictivo real.

3 Exploración de Datos

3.1 Distribución del Cumplimiento Normativo

Ver código 
# ── Período máximo de referencia ──────────────────────────────────────────────
periodo_max <- max(dt_cuenta$Tiempo, na.rm = TRUE)  # 291

# ── Tabla EPV por cuenta con cobertura temporal ──────────────────────────────
epv_tabla <- dt_cuenta %>%
  group_by(Cuenta, Nombre_Cuenta) %>%
  summarise(
    Inicio      = min(Tiempo, na.rm = TRUE),
    Períodos    = n(),
    Esperados   = periodo_max - Inicio + 1,
    Faltantes   = Esperados - Períodos,
    Cumple      = sum(Cumple == "Sí", na.rm = TRUE),
    No_Cumple   = sum(Cumple == "No", na.rm = TRUE),
    Pct_Cumple  = round(Cumple  / Períodos * 100, 1),
    EPV_max_vars = floor(No_Cumple / 10),
    .groups = "drop"
  )

# ── Fila resumen sede ─────────────────────────────────────────────────────────
epv_sede <- dt_sede %>%
  summarise(
    Inicio       = min(Tiempo, na.rm = TRUE),
    Períodos     = n(),
    Esperados    = periodo_max - Inicio + 1,
    Faltantes    = Esperados - Períodos,
    Cumple       = sum(Sede_Incumple == 0),
    No_Cumple    = sum(Sede_Incumple == 1),
    Pct_Cumple   = round(Cumple / Períodos * 100, 1),
    EPV_max_vars = floor(No_Cumple / 10)
  ) %>%
  mutate(Cuenta = "SEDE", Nombre_Cuenta = "Sede Central (todas)") %>%
  select(Cuenta, Nombre_Cuenta, everything())

epv_completa <- bind_rows(epv_tabla, epv_sede)

apa_table(epv_completa,
         col.names = c("Cuenta", "Nombre", "Inicio", "Períodos",
                       "Esperados", "Faltantes", "Cumple",
                       "Incumple", "% Cumple", "Máx. vars (EPV)"),
         caption   = "Distribución del Cumplimiento Normativo, Cobertura Temporal y Capacidad del Modelo (EPV)",
         align     = c("l", "l", "c", "c", "c", "c", "r", "r", "r", "c"),
         nota      = paste0("Inicio: período secuencial en que la cuenta registra su primera medición. ",
                            "Esperados: períodos desde el inicio hasta el período ", periodo_max,
                            ". Faltantes: períodos sin registro desde el inicio (huecos en la serie)."))
Distribución del Cumplimiento Normativo, Cobertura Temporal y Capacidad del Modelo (EPV)
Cuenta Nombre Inicio Períodos Esperados Faltantes Cumple Incumple % Cumple Máx. vars (EPV)
167382131 Central 1 1 286 291 5 282 0 98.6 0
167383918 Central 4 1 277 291 14 274 0 98.9 0
167385482 Central 2 1 286 291 5 147 0 51.4 0
357154485 Central 5 32 239 260 21 215 0 90.0 0
847747377 Central 3 265 24 27 3 17 0 70.8 0
SEDE Sede Central (todas) 1 291 291 0 149 142 51.2 14
Nota. Inicio: período secuencial en que la cuenta registra su primera medición. Esperados: períodos desde el inicio hasta el período 291. Faltantes: períodos sin registro desde el inicio (huecos en la serie).
NoteRegla de los 10 Eventos por Variable (EPV)

La regla EPV del módulo establece: \(\text{Máx. variables} = \frac{\text{Total incumplimientos}}{10}\). Para el modelo a nivel sede (142 incumplimientos) el modelo soporta hasta 14 variables. Para los modelos individuales se usará Semestre en lugar de Mes (12 dummies) para respetar el EPV de cada cuenta.

3.2 Gráfico Principal: Pago, Cuentas Pagadas e Incumplimiento en el Tiempo

Ver código 
# Escala de transformación para el eje secundario (0–5 → escala del pago)
pago_max  <- max(dt_sede$Pago_Total_Sede, na.rm = TRUE)
escala_k  <- pago_max / 5          # factor para llevar 0–5 a la escala del pago

ggplot(dt_sede, aes(x = Tiempo)) +

  # ── Área: Cuentas pagadas (fondo) ──────────────────────────────────────────
  geom_area(aes(y = Cuentas_Pagadas * escala_k),
            fill = "#AED6F1", alpha = 0.4) +

  # ── Barras: Cuentas que incumplen ──────────────────────────────────────────
  geom_col(aes(y = Cuentas_Incumplen * escala_k),
           fill  = "#E74C3C", alpha = 0.6, width = 0.8) +

  # ── Línea: Pago Total Sede ─────────────────────────────────────────────────
  geom_line(aes(y = Pago_Total_Sede),
            color = "#2C3E50", linewidth = 1.1) +
  geom_point(aes(y = Pago_Total_Sede,
                 color = factor(Sede_Incumple)),
             size = 1.8, alpha = 0.8) +
  scale_color_manual(
    values = c("0" = "#27AE60", "1" = "#E74C3C"),
    labels = c("0" = "Sede cumple", "1" = "Sede incumple"),
    name   = "Estado normativo"
  ) +

  # ── Eje secundario (derecho): escala 0–5 ──────────────────────────────────
  scale_y_continuous(
    name   = "Pago Total Sede (COP)",
    labels = label_comma(big.mark = ".", decimal.mark = ","),
    sec.axis = sec_axis(
      transform = ~ . / escala_k,
      name      = "Número de cuentas (máx. 5)",
      breaks    = 0:5
    )
  ) +

  labs(
    title    = "Pago Total, Cuentas Pagadas y Cuentas en Incumplimiento — Sede Central USCO",
    subtitle = "Línea: Pago total mensual | Área azul: cuentas pagadas | Barras rojas: cuentas que incumplen",
    x        = "Período (meses)",
    caption  = "Fuente: USCO — Consumo Sede Central 2000–2025"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title       = element_text(face = "bold", color = "#2C3E50", size = 13),
    plot.subtitle    = element_text(color = "#555555", size = 10),
    axis.title.y     = element_text(color = "#2C3E50"),
    axis.title.y.right = element_text(color = "#E74C3C"),
    legend.position  = "bottom",
    panel.grid.minor = element_blank()
  )

Interpretación: Los períodos donde las barras rojas son más altas (más cuentas incumpliendo) tienden a coincidir visualmente con picos en el valor del pago, lo que sugiere que el incumplimiento normativo tiene un efecto económico observable. Esta relación motiva la construcción del modelo logístico.

3.3 Consumo Activo por Cuenta en el Tiempo

Ver código 
ggplot(dt_cuenta, aes(x = Tiempo, y = Activa / 1000,
                       color = Nombre_Cuenta, group = Nombre_Cuenta)) +
  geom_line(alpha = 0.8, linewidth = 0.9) +
  geom_point(aes(shape = Cumple), size = 1.5, alpha = 0.6) +
  scale_shape_manual(values = c("Sí" = 19, "No" = 4),
                     name   = "Cumple norma") +
  scale_color_brewer(palette = "Set1", name = "Cuenta") +
  labs(
    title    = "Energía Activa por Cuenta — Sede Central USCO",
    subtitle = "Punto ✕ = período de incumplimiento normativo",
    x        = "Período (meses)",
    y        = "Energía Activa (miles de kWh)",
    caption  = "Fuente: USCO — Consumo Sede Central 2000–2025"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", color = "#2C3E50"),
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

4 Modelo Logístico Simple — Nivel Sede

Variable Dependiente (\(Y\)): Sede_Incumple (0 = todas cumplen, 1 = al menos una incumple) Variable Predictora (\(X\)): Activa_Total (kWh totales de la sede ese mes)

4.1 La Función Logística

El modelo no estima directamente si la sede incumplirá, sino la probabilidad de que ocurra el incumplimiento, acotada estrictamente entre 0 y 1 mediante la función sigmoide:

\[P(Sede\_Incumple = 1 \mid X) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 \cdot Activa\_Total)}}\]

4.2 Estimación del Modelo

Ver código 
modelo_sede_simple <- glm(Sede_Incumple ~ Activa_Total,
                          data   = dt_sede,
                          family = binomial(link = "logit"))

resumen_s <- summary(modelo_sede_simple)
b0_s <- coef(modelo_sede_simple)[1]
b1_s <- coef(modelo_sede_simple)[2]

# ── Tabla de coeficientes ──────────────────────────────────────────────────────
coef_tabla_s <- data.frame(
  Parámetro   = c("β₀ (Intercepto)", "β₁ (Activa_Total)"),
  Estimador   = round(resumen_s$coefficients[, 1], 6),
  Error_Std   = round(resumen_s$coefficients[, 2], 6),
  z           = round(resumen_s$coefficients[, 3], 4),
  p_valor     = format(resumen_s$coefficients[, 4], scientific = TRUE)
)

apa_table(coef_tabla_s,
         col.names = c("Parámetro", "Estimador", "Error Estándar", "z", "p-valor"),
         caption   = "Coeficientes del Modelo Logístico Simple — Nivel Sede",
         align     = c("l", "r", "r", "r", "r"))
Coeficientes del Modelo Logístico Simple — Nivel Sede
Parámetro Estimador Error Estándar z p-valor
β₀ (Intercepto) -2.535903 0.388575 -6.5262 6.747827e-11
β₁ (Activa_Total) 0.000019 0.000003 6.6734 2.498793e-11
Ver código 
# ── Indicadores de ajuste ─────────────────────────────────────────────────────
ajuste_s <- data.frame(
  Indicador = c("Devianza nula", "Devianza residual", "AIC", "N observaciones"),
  Valor     = c(round(resumen_s$null.deviance, 2),
                round(resumen_s$deviance, 2),
                round(resumen_s$aic, 2),
                nrow(dt_sede))
)

apa_table(ajuste_s,
         col.names = c("Indicador", "Valor"),
         caption   = "Indicadores de Ajuste — Modelo Logístico Simple (Sede)",
         align     = c("l", "r"))
Indicadores de Ajuste — Modelo Logístico Simple (Sede)
Indicador Valor
Devianza nula 403.24
Devianza residual 345.72
AIC 349.72
N observaciones 291.00

4.3 Ecuación del Modelo

ImportantEcuación del Modelo Logístico Simple

Escala Log-Odds (Logit):

\[\ln\left(\frac{P}{1-P}\right) = -2.5359 + 0.00001930 \cdot Activa\_Total\]

Escala de Probabilidad (Sigmoide):

\[P(Incumple) = \frac{1}{1 + e^{-(-2.5359 + 0.00001930 \cdot Activa\_Total)}}\]

Interpretación del coeficiente \(\beta_1\):

  • Por cada kWh adicional de Activa_Total, el log-odds de incumplimiento cambia en 0.00001930 unidades.
  • En términos del Odds Ratio: \(OR = e^{\beta_1} = e^{0.00001930} = 1.000019\)

4.4 Interpretación: Odds Ratios e Intervalos de Confianza

Ver código 
or_sede_s <- exp(cbind(OR = coef(modelo_sede_simple),
                        confint(modelo_sede_simple)))

or_tabla <- as.data.frame(round(or_sede_s, 6))
or_tabla$Variable    <- rownames(or_tabla)
or_tabla$Significativo <- c("—",
  ifelse(resumen_s$coefficients[2, 4] < 0.05, "✓ Sí (p < 0.05)", "✗ No"))
or_tabla$Interpretacion <- c(
  "Punto de anclaje matemático (sin interpretación directa)",
  ifelse(exp(b1_s) > 1,
         "Factor de RIESGO: mayor activa → mayor prob. de incumplir",
         "Factor PROTECTOR: mayor activa → menor prob. de incumplir")
)

apa_table(or_tabla[, c("Variable", "OR", "2.5 %", "97.5 %",
                       "Significativo", "Interpretacion")],
         col.names = c("Variable", "OR", "IC 2.5%", "IC 97.5%",
                       "Significativo", "Interpretación"),
         caption   = "Odds Ratios — Modelo Logístico Simple (Sede)",
         align     = c("l", "r", "r", "r", "c", "l"))
Odds Ratios — Modelo Logístico Simple (Sede)
Variable OR IC 2.5% IC 97.5% Significativo Interpretación
(Intercept) 0.079190 0.035904 0.165371 — Punto de anclaje matemático (sin interpretación directa)
Activa_Total 1.000019 1.000014 1.000025 ✓ Sí (p < 0.05) Factor de RIESGO: mayor activa → mayor prob. de incumplir

4.5 Curva Logística Ajustada

Ver código 
rango_activa <- data.frame(
  Activa_Total = seq(min(dt_sede$Activa_Total, na.rm = TRUE),
                     max(dt_sede$Activa_Total, na.rm = TRUE),
                     length.out = 300)
)
rango_activa$Prob <- predict(modelo_sede_simple,
                              newdata = rango_activa, type = "response")

ggplot() +
  geom_jitter(data = dt_sede,
              aes(x = Activa_Total / 1000, y = Sede_Incumple,
                  color = factor(Sede_Incumple)),
              height = 0.02, alpha = 0.4, size = 1.5) +
  geom_line(data = rango_activa,
            aes(x = Activa_Total / 1000, y = Prob),
            color = "#2C3E50", linewidth = 1.4) +
  geom_hline(yintercept = 0.5, color = "#E74C3C",
             linetype = "dashed", linewidth = 0.8) +
  scale_color_manual(values = c("0" = "#27AE60", "1" = "#E74C3C"),
                     labels = c("0" = "Cumple", "1" = "Incumple"),
                     name   = "Estado real") +
  labs(
    title    = "Curva Logística: Activa_Total → Probabilidad de Incumplimiento",
    subtitle = "Línea discontinua roja: umbral de decisión p = 0.5",
    x        = "Energía Activa Total Sede (miles de kWh)",
    y        = "P(Sede Incumple)"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#2C3E50"),
        legend.position = "bottom")

4.6 Matriz de Confusión y Métricas

Ver código 
prob_s   <- predict(modelo_sede_simple, type = "response")
pred_s   <- ifelse(prob_s > 0.5, 1, 0)
mc_s     <- table(Prediccion = pred_s, Real = dt_sede$Sede_Incumple)

VN_s <- mc_s["0", "0"]; FN_s <- mc_s["0", "1"]
FP_s <- mc_s["1", "0"]; VP_s <- mc_s["1", "1"]

acc_s  <- (VP_s + VN_s) / sum(mc_s)
sens_s <- VP_s / (VP_s + FN_s)
esp_s  <- VN_s / (VN_s + FP_s)

# ── Gráfico: Matriz de Confusión estilo heatmap ──────────────────────────────
datos_mc_s <- data.frame(
  Predicho   = factor(c("Cumple (0)", "Incumple (1)", "Cumple (0)", "Incumple (1)"),
                      levels = c("Cumple (0)", "Incumple (1)")),
  Real       = factor(c("Cumple (0)", "Cumple (0)", "Incumple (1)", "Incumple (1)"),
                      levels = c("Incumple (1)", "Cumple (0)")),
  Frecuencia = c(VN_s, FP_s, FN_s, VP_s)
)

ggplot(datos_mc_s, aes(x = Predicho, y = Real, fill = Frecuencia)) +
  geom_tile(color = "black", linewidth = 1) +
  geom_text(aes(label = Frecuencia), size = 8, fontface = "bold") +
  scale_fill_gradient(low = "white", high = "#E74C3C") +
  labs(title    = "Matriz de Confusión: Modelo Logístico Simple — Sede",
       subtitle = paste0("Evaluación empírica sobre n = ", sum(mc_s), " períodos"),
       x        = "Predicción del Modelo",
       y        = "Realidad (Observada)") +
  theme_minimal() +
  theme(axis.text  = element_text(size = 12, face = "bold"),
        title      = element_text(size = 14, face = "bold"),
        legend.position = "none")

Ver código 
metricas_s <- data.frame(
  Métrica = c("Exactitud global (Accuracy)",
              "Sensibilidad (Recall — detectar incumplimientos)",
              "Especificidad (descartar cumplimientos)"),
  Valor   = paste0(round(c(acc_s, sens_s, esp_s) * 100, 2), "%"),
  Interpretacion = c(
    "Proporción total de clasificaciones correctas",
    "De todos los meses que SÍ hubo incumplimiento, ¿cuántos detectó el modelo?",
    "De todos los meses donde se cumplió la norma, ¿cuántos identificó correctamente?"
  )
)

apa_table(metricas_s,
         col.names = c("Métrica", "Valor", "Interpretación"),
         caption   = "Métricas de Desempeño — Modelo Simple Sede",
         align     = c("l", "r", "l"))
Métricas de Desempeño — Modelo Simple Sede
Métrica Valor Interpretación
Exactitud global (Accuracy) 68.04% Proporción total de clasificaciones correctas
Sensibilidad (Recall — detectar incumplimientos) 56.34% De todos los meses que SÍ hubo incumplimiento, ¿cuántos detectó el modelo?
Especificidad (descartar cumplimientos) 79.19% De todos los meses donde se cumplió la norma, ¿cuántos identificó correctamente?

5 Modelo Logístico Múltiple — Nivel Sede

Se incorpora Mes como covariable categórica (factor con 12 niveles) para evaluar si la estacionalidad mensual — calendarío académico, vacaciones, variaciones climáticas — influye en la probabilidad de incumplimiento, controlando por el nivel de consumo activo. A nivel sede (142 incumplimientos), la regla EPV permite hasta 14 variables, por lo que incluir 11 dummies de mes es estadísticamente viable.

\[\ln\left(\frac{P}{1-P}\right) = \beta_0 + \beta_1 \cdot Activa\_Total + \sum_{k=2}^{12} \beta_k \cdot \mathbb{1}_{(Mes=k)}\]

5.1 Estimación del Modelo

Ver código 
modelo_sede_mult <- glm(Sede_Incumple ~ Activa_Total + factor(Mes),
                         data   = dt_sede,
                         family = binomial(link = "logit"))

resumen_m <- summary(modelo_sede_mult)
coef_m    <- coef(modelo_sede_mult)

# ── Tabla de coeficientes (dinámica para todos los parámetros) ───────────────
nombres_param <- rownames(resumen_m$coefficients)
nombres_param[1] <- "β₀ (Intercepto)"
nombres_param[2] <- "β₁ (Activa_Total)"
nombres_param <- gsub("factor\\(Mes\\)", "Mes ", nombres_param)

coef_tabla_m <- data.frame(
  Parámetro   = nombres_param,
  Estimador   = round(resumen_m$coefficients[, 1], 6),
  Error_Std   = round(resumen_m$coefficients[, 2], 6),
  z           = round(resumen_m$coefficients[, 3], 4),
  p_valor     = format(resumen_m$coefficients[, 4], scientific = TRUE)
)

apa_table(coef_tabla_m,
         col.names = c("Parámetro", "Estimador", "Error Estándar", "z", "p-valor"),
         caption   = "Coeficientes del Modelo Logístico Múltiple — Nivel Sede (Activa + Mes)",
         align     = c("l", "r", "r", "r", "r"),
         nota      = "Categoría de referencia: Mes 1 (Enero).")
Coeficientes del Modelo Logístico Múltiple — Nivel Sede (Activa + Mes)
Parámetro Estimador Error Estándar z p-valor
β₀ (Intercepto) -2.471869 0.542779 -4.5541 5.260993e-06
β₁ (Activa_Total) 0.000030 0.000004 7.4816 7.344313e-14
Mes 2 -1.493035 0.689637 -2.1650 3.039085e-02
Mes 3 -1.964731 0.692669 -2.8365 4.561588e-03
Mes 4 -1.773867 0.689402 -2.5731 1.008065e-02
Mes 5 -1.883626 0.712352 -2.6442 8.187610e-03
Mes 6 -0.865573 0.617690 -1.4013 1.611229e-01
Mes 7 -0.644114 0.615318 -1.0468 2.951923e-01
Mes 8 -1.839221 0.688139 -2.6727 7.523277e-03
Mes 9 -2.847058 0.754432 -3.7738 1.607944e-04
Mes 10 -2.128035 0.711433 -2.9912 2.778890e-03
Mes 11 -1.458182 0.672816 -2.1673 3.021332e-02
Mes 12 -0.116694 0.611622 -0.1908 8.486870e-01
Nota. Categoría de referencia: Mes 1 (Enero).
Ver código 
# ── Indicadores de ajuste ─────────────────────────────────────────────────────
ajuste_m <- data.frame(
  Indicador = c("Devianza nula", "Devianza residual", "AIC",
                "N observaciones", "Parámetros estimados"),
  Valor     = c(round(resumen_m$null.deviance, 2),
                round(resumen_m$deviance, 2),
                round(resumen_m$aic, 2),
                nrow(dt_sede),
                length(coef_m))
)

apa_table(ajuste_m,
         col.names = c("Indicador", "Valor"),
         caption   = "Indicadores de Ajuste — Modelo Logístico Múltiple (Sede)",
         align     = c("l", "r"))
Indicadores de Ajuste — Modelo Logístico Múltiple (Sede)
Indicador Valor
Devianza nula 403.24
Devianza residual 319.77
AIC 345.77
N observaciones 291.00
Parámetros estimados 13.00

5.2 Ecuación e Interpretación

ImportantEcuación del Modelo Logístico Múltiple

Escala Log-Odds:

\[\ln\left(\frac{P}{1-P}\right) = -2.4719 + 0.00002957 \cdot Activa\_Total + \sum_{k=2}^{12} \hat{\beta}_k \cdot \mathbb{1}_{(Mes=k)}\]

Interpretación de los coeficientes:

  • \(\hat{\beta}_1\) = 0.00002957 → Por cada kWh adicional de energía activa total, el log-odds de incumplimiento cambia en 0.00002957 unidades, manteniendo constante el mes.
  • Cada \(\hat{\beta}_k\) (Mes \(k\)) mide la diferencia en log-odds respecto a Enero (categoría de referencia), controlando por consumo activo.
  • Meses estadísticamente significativos (p < 0.05): Mes 2, Mes 3, Mes 4, Mes 5, Mes 8, Mes 9, Mes 10, Mes 11.
Ver código 
or_mult <- exp(cbind(OR = coef(modelo_sede_mult),
                      confint(modelo_sede_mult)))

or_mult_df <- as.data.frame(round(or_mult, 4))
nombres_or <- rownames(or_mult_df)
nombres_or[1] <- "Intercepto"
nombres_or[2] <- "Activa_Total"
nombres_or <- gsub("factor\\(Mes\\)", "Mes ", nombres_or)

or_mult_df$Variable <- nombres_or
or_mult_df$p_valor  <- round(resumen_m$coefficients[, 4], 6)
or_mult_df$Signif   <- ifelse(or_mult_df$p_valor < 0.05, "✓ Sí", "✗ No")

apa_table(or_mult_df[, c("Variable", "OR", "2.5 %", "97.5 %",
                         "p_valor", "Signif")],
         col.names = c("Variable", "OR", "IC 2.5%", "IC 97.5%",
                       "p-valor", "Significativo"),
         caption   = "Odds Ratios — Modelo Logístico Múltiple (Sede: Activa + Mes)",
         align     = c("l", "r", "r", "r", "r", "c"),
         nota      = "Categoría de referencia: Mes 1 (Enero). OR > 1 = mayor riesgo que Enero; OR < 1 = menor riesgo.")
Odds Ratios — Modelo Logístico Múltiple (Sede: Activa + Mes)
Variable OR IC 2.5% IC 97.5% p-valor Significativo
Intercepto 0.0844 0.0283 0.2409 0.000005 ✓ Sí
Activa_Total 1.0000 1.0000 1.0000 0.000000 ✓ Sí
Mes 2 0.2247 0.0567 0.8594 0.030391 ✓ Sí
Mes 3 0.1402 0.0350 0.5361 0.004562 ✓ Sí
Mes 4 0.1697 0.0424 0.6422 0.010081 ✓ Sí
Mes 5 0.1520 0.0363 0.6012 0.008188 ✓ Sí
Mes 6 0.4208 0.1234 1.4079 0.161123 ✗ No
Mes 7 0.5251 0.1548 1.7525 0.295192 ✗ No
Mes 8 0.1589 0.0398 0.5993 0.007523 ✓ Sí
Mes 9 0.0580 0.0125 0.2453 0.000161 ✓ Sí
Mes 10 0.1191 0.0283 0.4683 0.002779 ✓ Sí
Mes 11 0.2327 0.0605 0.8587 0.030213 ✓ Sí
Mes 12 0.8899 0.2663 2.9690 0.848687 ✗ No
Nota. Categoría de referencia: Mes 1 (Enero). OR > 1 = mayor riesgo que Enero; OR < 1 = menor riesgo.

5.3 Comparación de Modelos por AIC — Nivel Sede

Ver código 
# Crear variable Semestre para comparación
dt_sede$Semestre <- factor(ifelse(dt_sede$Mes <= 6, "Ene-Jun", "Jul-Dic"),
                            levels = c("Ene-Jun", "Jul-Dic"))

modelos_sede <- list(
  "M1: ~ Activa_Total"              = glm(Sede_Incumple ~ Activa_Total,
                                           data = dt_sede, family = "binomial"),
  "M2: ~ Semestre"                  = glm(Sede_Incumple ~ Semestre,
                                           data = dt_sede, family = "binomial"),
  "M3: ~ Activa_Total + Semestre"   = glm(Sede_Incumple ~ Activa_Total + Semestre,
                                           data = dt_sede, family = "binomial"),
  "M4: ~ Activa_Total + Mes"        = modelo_sede_mult,
  "M5: ~ Mes"                       = glm(Sede_Incumple ~ factor(Mes),
                                           data = dt_sede, family = "binomial")
)

tabla_aic_sede <- data.frame(
  Modelo     = names(modelos_sede),
  Parámetros = sapply(modelos_sede, function(m) length(coef(m))),
  AIC        = round(sapply(modelos_sede, AIC), 2),
  Devianza   = round(sapply(modelos_sede, function(m) m$deviance), 2)
)
tabla_aic_sede <- tabla_aic_sede[order(tabla_aic_sede$AIC), ]
tabla_aic_sede$Ranking <- 1:nrow(tabla_aic_sede)

apa_table(tabla_aic_sede,
         col.names = c("Modelo", "Parámetros", "AIC", "Devianza", "Ranking"),
         caption   = "Comparación de Modelos por AIC — Nivel Sede (menor AIC = mejor)",
         align     = c("l", "c", "r", "r", "c"))
Comparación de Modelos por AIC — Nivel Sede (menor AIC = mejor)
Modelo Parámetros AIC Devianza Ranking
M4: ~ Activa_Total + Mes 13 345.77 319.77 1
M1: ~ Activa_Total 2 349.72 345.72 2
M3: ~ Activa_Total + Semestre 3 351.46 345.46 3
M2: ~ Semestre 2 407.21 403.21 4
M5: ~ Mes 12 426.30 402.30 5

5.4 Matriz de Confusión — Modelo Múltiple Sede

Ver código 
prob_m  <- predict(modelo_sede_mult, type = "response")
pred_m  <- ifelse(prob_m > 0.5, 1, 0)
mc_m    <- table(Prediccion = pred_m, Real = dt_sede$Sede_Incumple)

VN_m <- mc_m["0","0"]; FN_m <- mc_m["0","1"]
FP_m <- mc_m["1","0"]; VP_m <- mc_m["1","1"]

acc_m  <- (VP_m + VN_m) / sum(mc_m)
sens_m <- VP_m / (VP_m + FN_m)
esp_m  <- VN_m / (VN_m + FP_m)

# ── Gráfico: Matriz de Confusión estilo heatmap ──────────────────────────────
datos_mc_m <- data.frame(
  Predicho   = factor(c("Cumple (0)", "Incumple (1)", "Cumple (0)", "Incumple (1)"),
                      levels = c("Cumple (0)", "Incumple (1)")),
  Real       = factor(c("Cumple (0)", "Cumple (0)", "Incumple (1)", "Incumple (1)"),
                      levels = c("Incumple (1)", "Cumple (0)")),
  Frecuencia = c(VN_m, FP_m, FN_m, VP_m)
)

ggplot(datos_mc_m, aes(x = Predicho, y = Real, fill = Frecuencia)) +
  geom_tile(color = "black", linewidth = 1) +
  geom_text(aes(label = Frecuencia), size = 8, fontface = "bold") +
  scale_fill_gradient(low = "white", high = "#E74C3C") +
  labs(title    = "Matriz de Confusión: Modelo Logístico Múltiple — Sede",
       subtitle = paste0("Evaluación empírica sobre n = ", sum(mc_m), " períodos"),
       x        = "Predicción del Modelo",
       y        = "Realidad (Observada)") +
  theme_minimal() +
  theme(axis.text  = element_text(size = 12, face = "bold"),
        title      = element_text(size = 14, face = "bold"),
        legend.position = "none")

Ver código 
comp_metricas <- data.frame(
  Métrica        = c("Accuracy", "Sensibilidad", "Especificidad", "AIC"),
  Modelo_Simple  = c(paste0(round(acc_s  * 100, 2), "%"),
                     paste0(round(sens_s * 100, 2), "%"),
                     paste0(round(esp_s  * 100, 2), "%"),
                     round(AIC(modelo_sede_simple), 2)),
  Modelo_Multiple = c(paste0(round(acc_m  * 100, 2), "%"),
                      paste0(round(sens_m * 100, 2), "%"),
                      paste0(round(esp_m  * 100, 2), "%"),
                      round(AIC(modelo_sede_mult), 2))
)

apa_table(comp_metricas,
         col.names = c("Métrica", "Modelo Simple", "Modelo Múltiple"),
         caption   = "Comparación de Desempeño — Modelos Nivel Sede",
         align     = c("l", "r", "r"))
Comparación de Desempeño — Modelos Nivel Sede
Métrica Modelo Simple Modelo Múltiple
Accuracy 68.04% 81.1%
Sensibilidad 56.34% 73.94%
Especificidad 79.19% 87.92%
AIC 349.72 345.77

6 Modelos por Cuenta Individual

6.1 Análisis EPV y Viabilidad por Cuenta

Ver código 
apa_table(epv_tabla[, c("Nombre_Cuenta", "Períodos", "No_Cumple",
                        "Pct_Cumple", "EPV_max_vars")],
         col.names = c("Cuenta", "Períodos", "Incumplimientos",
                       "% Cumple", "Máx. vars EPV"),
         caption   = "Capacidad del Modelo por Cuenta (EPV)",
         align     = c("l", "c", "r", "r", "c"))
Capacidad del Modelo por Cuenta (EPV)
Cuenta Períodos Incumplimientos % Cumple Máx. vars EPV
Central 1 286 0 98.6 0
Central 4 277 0 98.9 0
Central 2 286 0 51.4 0
Central 5 239 0 90.0 0
Central 3 24 0 70.8 0
ImportantDecisión metodológica por cuenta
  • Cuentas con EPV ≥ 2: Se ajusta modelo simple (Activa) y modelo múltiple (Activa + Semestre).
  • Cuentas con EPV = 1: Solo modelo simple — agregar Semestre superaría el EPV y los estimadores serían inestables.
  • Cuentas con EPV = 0: No se modela — no hay eventos suficientes para estimar ningún coeficiente.

6.2 Modelos Simples Individuales — Activa → Cumple

Ver código 
cuentas <- levels(dt_cuenta$Nombre_Cuenta)
resultados_ind <- list()

for (cta in cuentas) {
  dt_sub <- dt_cuenta[dt_cuenta$Nombre_Cuenta == cta, ]
  n_inc  <- sum(dt_sub$Cumple_bin == 0, na.rm = TRUE)

  if (n_inc < 1) {
    cat(sprintf("  %-12s → Sin incumplimientos. Modelo no ajustado.\n", cta))
    next
  }

  modelo_i <- tryCatch(
    glm(Cumple_bin ~ Activa, data = dt_sub, family = "binomial"),
    error = function(e) NULL
  )

  if (is.null(modelo_i)) {
    cat(sprintf("  %-12s → Error al ajustar el modelo.\n", cta))
    next
  }

  coef_i <- coef(modelo_i)
  pval_i <- summary(modelo_i)$coefficients[2, 4]
  or_i   <- exp(coef_i[2])
  aic_i  <- AIC(modelo_i)

  resultados_ind[[cta]] <- data.frame(
    Cuenta   = cta,
    b0       = round(coef_i[1], 4),
    b1       = round(coef_i[2], 6),
    OR       = round(or_i, 4),
    p_valor  = round(pval_i, 6),
    AIC      = round(aic_i, 2),
    Signif   = ifelse(pval_i < 0.05, "✓ Sí", "✗ No"),
    Efecto   = ifelse(or_i > 1,
                      "Riesgo (+Activa → menos cumplimiento)",
                      "Protector (+Activa → más cumplimiento)")
  )
}

tabla_ind <- do.call(rbind, resultados_ind)

apa_table(tabla_ind,
         col.names = c("Cuenta", "β₀", "β₁ (Activa)", "OR",
                       "p-valor", "AIC", "Significativo", "Dirección del efecto"),
         caption   = "Resumen de Modelos Logísticos Simples por Cuenta",
         align     = c("l", "r", "r", "r", "r", "r", "c", "l"))
Resumen de Modelos Logísticos Simples por Cuenta
Cuenta β₀ β₁ (Activa) OR p-valor AIC Significativo Dirección del efecto
Central 1 -4.1310 0.000803 1.0008 0.034879 21.43 ✓ Sí Riesgo (+Activa → menos cumplimiento)
Central 2 3.5379 -0.000046 1.0000 0.000000 274.56 ✓ Sí Protector (+Activa → más cumplimiento)
Central 3 -313.4424 0.055895 1.0575 0.996625 4.00 ✗ No Riesgo (+Activa → menos cumplimiento)
Central 4 3.6513 0.000036 1.0000 0.483543 36.64 ✗ No Riesgo (+Activa → menos cumplimiento)
Central 5 1.6161 0.000564 1.0006 0.026553 153.45 ✓ Sí Riesgo (+Activa → menos cumplimiento)

6.3 Modelo Logístico Unificado — Todas las Cuentas

El modelo unificado evalúa si el efecto de la energía activa sobre el cumplimiento es homogéneo entre cuentas, o si hay cuentas estructuralmente más propensas a incumplir.

6.3.1 Modelo Simple Unificado

Ver código 
modelo_unif_s <- glm(Cumple_bin ~ Activa + Cuenta,
                     data   = dt_cuenta,
                     family = binomial(link = "logit"))

resumen_us <- summary(modelo_unif_s)

coef_tabla_us <- data.frame(
  Parámetro   = rownames(resumen_us$coefficients),
  Estimador   = round(resumen_us$coefficients[, 1], 6),
  Error_Std   = round(resumen_us$coefficients[, 2], 6),
  z           = round(resumen_us$coefficients[, 3], 4),
  p_valor     = format(resumen_us$coefficients[, 4], scientific = TRUE)
)

apa_table(coef_tabla_us,
         col.names = c("Parámetro", "Estimador", "Error Estándar", "z", "p-valor"),
         caption   = "Coeficientes — Modelo Unificado Simple (Activa + Cuenta)",
         align     = c("l", "r", "r", "r", "r"),
         nota      = paste0("Devianza nula: ", round(resumen_us$null.deviance, 2),
                            " | Devianza residual: ", round(resumen_us$deviance, 2),
                            " | AIC: ", round(resumen_us$aic, 2)))
Coeficientes — Modelo Unificado Simple (Activa + Cuenta)
Parámetro Estimador Error Estándar z p-valor
(Intercept) 5.327696 0.522354 10.1994 1.995232e-24
Activa -0.000039 0.000005 -8.2993 1.047705e-16
Cuenta167383918 0.309544 0.769190 0.4024 6.873693e-01
Cuenta167385482 -2.245916 0.564166 -3.9810 6.864038e-05
Cuenta357154485 -3.083124 0.563449 -5.4719 4.452884e-08
Cuenta847747377 -3.783688 0.683062 -5.5393 3.036692e-08
Nota. Devianza nula: 974.77 | Devianza residual: 549.29 | AIC: 561.29

6.3.2 Modelo Múltiple Unificado

Ver código 
modelo_unif_m <- glm(Cumple_bin ~ Activa + Cuenta + Semestre,
                     data   = dt_cuenta,
                     family = binomial(link = "logit"))

resumen_um <- summary(modelo_unif_m)

coef_tabla_um <- data.frame(
  Parámetro   = rownames(resumen_um$coefficients),
  Estimador   = round(resumen_um$coefficients[, 1], 6),
  Error_Std   = round(resumen_um$coefficients[, 2], 6),
  z           = round(resumen_um$coefficients[, 3], 4),
  p_valor     = format(resumen_um$coefficients[, 4], scientific = TRUE)
)

apa_table(coef_tabla_um,
         col.names = c("Parámetro", "Estimador", "Error Estándar", "z", "p-valor"),
         caption   = "Coeficientes — Modelo Unificado Múltiple (Activa + Cuenta + Semestre)",
         align     = c("l", "r", "r", "r", "r"),
         nota      = paste0("Devianza nula: ", round(resumen_um$null.deviance, 2),
                            " | Devianza residual: ", round(resumen_um$deviance, 2),
                            " | AIC: ", round(resumen_um$aic, 2)))
Coeficientes — Modelo Unificado Múltiple (Activa + Cuenta + Semestre)
Parámetro Estimador Error Estándar z p-valor
(Intercept) 5.297766 0.534193 9.9173 3.499978e-23
Activa -0.000039 0.000005 -8.3001 1.040532e-16
Cuenta167383918 0.309507 0.769177 0.4024 6.873993e-01
Cuenta167385482 -2.245514 0.564175 -3.9802 6.886468e-05
Cuenta357154485 -3.081095 0.563439 -5.4684 4.541924e-08
Cuenta847747377 -3.779377 0.683266 -5.5313 3.177919e-08
SemestreJul–Dic 0.058084 0.222018 0.2616 7.936175e-01
Nota. Devianza nula: 974.77 | Devianza residual: 549.22 | AIC: 563.22

6.3.3 Odds Ratios — Modelo Unificado Múltiple

Ver código 
or_unif <- exp(cbind(OR = coef(modelo_unif_m), confint(modelo_unif_m)))
or_unif_df <- as.data.frame(round(or_unif, 4))
or_unif_df$p_valor <- round(summary(modelo_unif_m)$coefficients[, 4], 6)
or_unif_df$Signif  <- ifelse(or_unif_df$p_valor < 0.05, "✓ Sí", "✗ No")

apa_table(or_unif_df,
         col.names = c("OR", "IC 2.5%", "IC 97.5%", "p-valor", "Significativo"),
         caption   = "Odds Ratios — Modelo Logístico Unificado Múltiple",
         align     = c("r", "r", "r", "r", "c"))
Odds Ratios — Modelo Logístico Unificado Múltiple
OR IC 2.5% IC 97.5% p-valor Significativo
199.8897 79.2163 676.9462 0.000000 ✓ Sí
1.0000 1.0000 1.0000 0.000000 ✓ Sí
1.3628 0.2974 6.9794 0.687399 ✗ No
0.1059 0.0300 0.2899 0.000069 ✓ Sí
0.0459 0.0130 0.1253 0.000000 ✓ Sí
0.0228 0.0054 0.0843 0.000000 ✓ Sí
1.0598 0.6856 1.6401 0.793618 ✗ No

6.3.4 Comparación por AIC — Modelos Unificados

Ver código 
modelos_unif <- list(
  "M1: ~ Activa"                  = glm(Cumple_bin ~ Activa,
                                         data = dt_cuenta, family = "binomial"),
  "M2: ~ Activa + Cuenta"         = modelo_unif_s,
  "M3: ~ Activa + Semestre"       = glm(Cumple_bin ~ Activa + Semestre,
                                         data = dt_cuenta, family = "binomial"),
  "M4: ~ Activa + Cuenta + Semestre" = modelo_unif_m
)

tabla_aic_unif <- data.frame(
  Modelo     = names(modelos_unif),
  Parámetros = sapply(modelos_unif, function(m) length(coef(m))),
  AIC        = round(sapply(modelos_unif, AIC), 2),
  Devianza   = round(sapply(modelos_unif, function(m) m$deviance), 2)
)
tabla_aic_unif <- tabla_aic_unif[order(tabla_aic_unif$AIC), ]
tabla_aic_unif$Ranking <- 1:nrow(tabla_aic_unif)

apa_table(tabla_aic_unif,
         col.names = c("Modelo", "Parámetros", "AIC", "Devianza", "Ranking"),
         caption   = "Comparación por AIC — Modelos Unificados (menor AIC = mejor)",
         align     = c("l", "c", "r", "r", "c"))
Comparación por AIC — Modelos Unificados (menor AIC = mejor)
Modelo Parámetros AIC Devianza Ranking
M2: ~ Activa + Cuenta 6 561.29 549.29 1
M4: ~ Activa + Cuenta + Semestre 7 563.22 549.22 2
M1: ~ Activa 2 654.10 650.10 3
M3: ~ Activa + Semestre 3 655.84 649.84 4

6.4 Matrices de Confusión — Modelos Unificados

Ver código 
evaluar_modelo <- function(modelo, datos, y_col, nombre) {
  prob  <- predict(modelo, type = "response")
  pred  <- ifelse(prob > 0.5, 1, 0)
  real  <- datos[[y_col]]
  mc    <- table(Prediccion = pred, Real = real)

  VN <- mc["0","0"]; FN <- mc["0","1"]
  FP <- mc["1","0"]; VP <- mc["1","1"]

  data.frame(
    Modelo        = nombre,
    Accuracy      = paste0(round((VP+VN)/sum(mc)*100, 2), "%"),
    Sensibilidad  = paste0(round(VP/(VP+FN)*100, 2), "%"),
    Especificidad = paste0(round(VN/(VN+FP)*100, 2), "%"),
    AIC           = round(AIC(modelo), 2)
  )
}

res_unif <- rbind(
  evaluar_modelo(modelo_unif_s, dt_cuenta, "Cumple_bin", "Simple: Activa + Cuenta"),
  evaluar_modelo(modelo_unif_m, dt_cuenta, "Cumple_bin", "Múltiple: Activa + Cuenta + Semestre")
)

apa_table(res_unif,
         col.names = c("Modelo", "Accuracy", "Sensibilidad", "Especificidad", "AIC"),
         caption   = "Comparación de Desempeño — Modelos Unificados",
         align     = c("l", "r", "r", "r", "r"))
Comparación de Desempeño — Modelos Unificados
Modelo Accuracy Sensibilidad Especificidad AIC
Simple: Activa + Cuenta 91.28% 98.18% 54.8% 561.29
Múltiple: Activa + Cuenta + Semestre 91.37% 98.29% 54.8% 563.22

7 Conclusiones

1. Balance de los datos: El 48.8% de los períodos registran al menos una cuenta en incumplimiento. Este balance casi perfecto (cercano a 50/50) a nivel sede evita el problema de desbalance de clases que advierte el módulo, haciendo al modelo estadísticamente confiable.

2. Modelo Simple — Nivel Sede: La energía activa total es un predictor estadísticamente significativo (p < 0.05) del incumplimiento normativo de la sede. El OR indica que a mayor consumo activo, mayor es la probabilidad de que alguna cuenta incumpla la norma..

3. Selección por AIC — Nivel Sede: El modelo recomendado es M4: ~ Activa_Total + Mes por tener el menor AIC, balanceando ajuste y parsimonia.

4. Modelos Individuales por Cuenta: Los modelos individuales revelan si el efecto de la energía activa es homogéneo entre cuentas o si algunas tienen un comportamiento diferenciado. La tabla de OR comparativa permite identificar cuáles cuentas son estructuralmente más propensas al incumplimiento.

5. Modelo Unificado: El modelo seleccionado por AIC es M2: ~ Activa + Cuenta. Si Cuenta resulta significativa, confirma que existen diferencias estructurales entre cuentas que van más allá del nivel de consumo activo — posiblemente relacionadas con el tipo de carga eléctrica instalada (motores, aires acondicionados, iluminación) en cada edificio.

Análisis desarrollado para la Universidad Surcolombiana (USCO) — Sede Central Herramienta: R R version 4.5.3 (2026-03-11 ucrt)