Taller: Modelo de Regresión Lineal

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 Universidad Surcolombiana (USCO) registra mensualmente su consumo eléctrico y el valor total de pago de sus cuentas desde el año 2000. El análisis integra dos fuentes:

  • Consumo_Final.csv: Sumatoria mensual de todas las cuentas (energía activa, reactiva y su relación normativa).
  • Pago_Cuentas.xls: Valor del pago mensual por cuenta. El pago incluye el costo de la energía activa y, cuando la universidad incumple la relación Reactiva/Activa > 0.48, también el recargo por energía reactiva.

La normativa eléctrica colombiana establece que la relación Reactiva / Activa no debe superar 0.48. El incumplimiento genera penalizaciones que se reflejan directamente en el valor del pago.

NoteHipótesis Principal

“El valor total de pago mensual de la USCO puede ser explicado por el tiempo transcurrido (tendencia de crecimiento) y, en el modelo múltiple, por el nivel de energía activa consumida y por el cumplimiento o incumplimiento de la norma reactiva.”

Hipótesis Estadísticas — Regresión Simple (Tiempo → Pago):

  • \(H_0: \beta_1 = 0\) → El tiempo NO tiene efecto lineal significativo sobre el valor de pago.
  • \(H_1: \beta_1 \neq 0\) → El tiempo SÍ tiene efecto lineal significativo sobre el valor de pago (tendencia creciente).

Hipótesis Estadísticas — Regresión Múltiple (Activa + Cumple → Pago):

  • \(H_0: \beta_1 = \beta_2 = 0\) → Ni la energía activa ni el cumplimiento normativo explican el pago.
  • \(H_1\): Al menos un \(\beta_k \neq 0\) → Alguna variable predictora sí explica el pago.

Hipótesis de Normalidad (Shapiro-Wilk):

  • \(H_0\): Los datos siguen distribución normal.
  • \(H_1\): Los datos no siguen distribución normal.

2 Carga de Datos y Descripción de Variables

Ver código
# ── 1. Cargar consumo mensual ──────────────────────────────────────────────────
dt_consumo <- read.csv("Consumo_Final.csv", header = TRUE, sep = ",")
dt <- dt_consumo

# ── 4. Variable dummy: Cumple (Sí = 0, No = 1) ────────────────────────────────
# "No cumple" significa que la relación > 0.48 y genera recargo en el pago
dt$Incumple <- ifelse(dt$Cumple == "No", 1, 0)

# ── 5. Formato auxiliar ────────────────────────────────────────────────────────
fmt <- function(x, dec = 0) {
  formatC(round(x, dec), format = "f", digits = dec,
          big.mark = ".", decimal.mark = ",")
}

resumen_datos <- data.frame(
  Concepto = c("Registros válidos", "Período", 
               "Períodos con incumplimiento"),
  Valor = c(as.character(nrow(dt)),
            paste(min(dt$Ano), "–", max(dt$Ano)),
            as.character(sum(dt$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 válidos 308
Período 2000 – 2025
Períodos con incumplimiento 1
Ver código
apa_table(
  head(dt[, c("Ano","Mes","Activa","Reactiva",
              "Relacion","Cumple","Pago_Total","Tiempo")], 8),
  caption = "Primeros 8 registros del conjunto de datos",
  align   = c("c","c","r","r","r","c","r","c")
)
Primeros 8 registros del conjunto de datos
Ano Mes Activa Reactiva Relacion Cumple Pago_Total Tiempo
2000 2 97995 0 0 795340 1
2000 3 96381 0 0 806660 2
2000 4 409 0 0 710750 3
2000 5 388 0 0 836190 4
2000 6 376 0 0 807600 5
2000 7 380 0 0 824930 6
2000 8 331 0 0 937200 7
2000 9 361 0 0 830800 8
Ver código
tabla_vars <- data.frame(
  Variable    = c("Tiempo", "Pago_Total", "Activa", "Reactiva",
                  "Relacion", "Cumple / Incumple", "Ano / Mes"),
  Tipo        = c("Numérica discreta", "Numérica continua", "Numérica continua",
                  "Numérica continua", "Numérica continua",
                  "Categórica / Dummy (0-1)", "Numérica"),
  Rol         = c("Independiente (X) — Modelo Simple",
                  "**Dependiente (Y) — ambos modelos**",
                  "Independiente (X₁) — Modelo Múltiple",
                  "Variable de contexto",
                  "Variable de control",
                  "Independiente (X₂) — Modelo Múltiple",
                  "Temporal"),
  Descripción = c(
    "Período secuencial mensual (1 al N)",
    "Suma total del valor pagado mensualmente por todas las cuentas (COP)",
    "Energía activa mensual total (kWh) — incluida en el pago base",
    "Energía reactiva mensual total (kVArh) — genera recargo si incumple",
    "Razón Reactiva/Activa — límite normativo: 0.48",
    "Cumple: relación ≤ 0.48; Incumple (dummy=1): relación > 0.48 → recargo",
    "Año y mes del registro"
  )
)

apa_table(tabla_vars,
         col.names = c("Variable", "Tipo", "Rol en el Modelo", "Descripción"),
         caption   = "Descripción de las Variables del Estudio",
         align     = c("l", "l", "l", "l"))
Descripción de las Variables del Estudio
Variable Tipo Rol en el Modelo Descripción
Tiempo Numérica discreta Independiente (X) — Modelo Simple Período secuencial mensual (1 al N)
Pago_Total Numérica continua **Dependiente (Y) — ambos modelos** Suma total del valor pagado mensualmente por todas las cuentas (COP)
Activa Numérica continua Independiente (X₁) — Modelo Múltiple Energía activa mensual total (kWh) — incluida en el pago base
Reactiva Numérica continua Variable de contexto Energía reactiva mensual total (kVArh) — genera recargo si incumple
Relacion Numérica continua Variable de control Razón Reactiva/Activa — límite normativo: 0.48
Cumple / Incumple Categórica / Dummy (0-1) Independiente (X₂) — Modelo Múltiple Cumple: relación ≤ 0.48; Incumple (dummy=1): relación > 0.48 → recargo
Ano / Mes Numérica Temporal Año y mes del registro
ImportantLógica económica del Pago_Total

El Pago_Total mensual refleja el costo de la energía activa consumida. Cuando la universidad incumple la norma (Reactiva/Activa > 0.48), el operador de red agrega un recargo por reactiva, elevando el valor de la factura. Así, el incumplimiento normativo es una variable que encarece directamente el pago.


3 Análisis Exploratorio y Correlación

3.1 Gráficos Exploratorios

Ver código
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

# Boxplot Tiempo
boxplot(dt$Tiempo,
        main   = "Distribución: Tiempo (X)",
        col    = "#AED6F1", border = "#1A5276",
        ylab   = "Período (meses)")

# Boxplot Pago
boxplot(dt$Pago_Total / 1e6,
        main   = "Distribución: Pago Total (Y)",
        col    = "#A9DFBF", border = "#1D8348",
        ylab   = "Millones de COP")

# Dispersión Tiempo vs Pago
plot(dt$Tiempo, dt$Pago_Total / 1e6,
     main = "Dispersión: Tiempo vs. Pago Total",
     xlab = "Período (meses)", ylab = "Pago Total (millones COP)",
     col  = "#2980B9", pch = 19, cex = 0.7)
abline(lm(Pago_Total / 1e6 ~ Tiempo, data = dt),
       col = "#E74C3C", lwd = 2)
legend("topleft", legend = "Regresión OLS",
       col = "#E74C3C", lty = 1, lwd = 2, bty = "n", cex = 0.8)

# Evolución temporal del pago con cumplimiento coloreado
colores_c <- ifelse(dt$Incumple == 1, "#E74C3C", "#27AE60")
plot(dt$Tiempo, dt$Pago_Total / 1e6,
     type = "p", pch = 19, cex = 0.7, col = colores_c,
     main = "Pago Total en el Tiempo\n(verde: cumple | rojo: incumple)",
     xlab = "Período (meses)", ylab = "Pago Total (millones COP)")
lines(dt$Tiempo, dt$Pago_Total / 1e6, col = "gray70", lwd = 0.8)
legend("topleft", legend = c("Cumple norma", "Incumple norma"),
       col = c("#27AE60", "#E74C3C"), pch = 19, bty = "n", cex = 0.8)

Ver código
par(mfrow = c(1, 1))

Observación sobre valores atípicos — Pago Total: El boxplot del Pago Total revela la presencia de observaciones alejadas de la tendencia central, incluyendo valores negativos y montos excepcionalmente altos que no corresponden al comportamiento típico de la serie. Estos outliers tienen explicaciones operativas concretas:

  • Correcciones retroactivas de factura: El operador de red puede emitir notas crédito o ajustes de períodos anteriores que aparecen como valores negativos o inusualmente bajos en el registro mensual.
  • Errores de lectura del medidor: Fallas en los equipos de medición o en la transcripción de datos generan registros que no reflejan el consumo real del período.
  • Períodos con cuentas consolidadas o divididas: En algunos meses se agrupan o separan cuentas de distintos medidores, alterando artificialmente el total registrado.

Estos puntos atípicos tienen implicaciones directas sobre el análisis de regresión: al minimizar errores cuadráticos, el método OLS es especialmente sensible a valores extremos, lo que puede distorsionar los coeficientes estimados e inflar la varianza de los residuos.

Interpretación: Los gráficos permiten observar la tendencia creciente del pago a lo largo del tiempo (efecto inflacionario y mayor consumo) y cómo los períodos de incumplimiento (puntos rojos) tienden a registrar pagos más elevados por el recargo en energía reactiva.

3.2 Prueba de Normalidad — Shapiro-Wilk

La prueba de Shapiro-Wilk es el criterio formal para determinar qué coeficiente de correlación aplicar.

Ver código
x_s <- dt$Tiempo
y_s <- dt$Pago_Total
n   <- nrow(dt)

sw_x <- shapiro.test(x_s)
sw_y <- shapiro.test(y_s)

tabla_normalidad <- data.frame(
  Variable  = c("Tiempo (X)", "Pago Total (Y)"),
  W         = c(round(sw_x$statistic, 4), round(sw_y$statistic, 4)),
  `p-valor` = c(round(sw_x$p.value,   6), round(sw_y$p.value,   6)),
  Resultado = c(
    ifelse(sw_x$p.value > 0.05, "✓ Normal",    "✗ No Normal"),
    ifelse(sw_y$p.value > 0.05, "✓ Normal",    "✗ No Normal")
  ),
  `Coeficiente recomendado` = c(
    ifelse(sw_x$p.value > 0.05, "Pearson", "Spearman / Kendall"),
    ifelse(sw_y$p.value > 0.05, "Pearson", "Spearman / Kendall")
  )
)

apa_table(tabla_normalidad,
         col.names = c("Variable", "W", "p-valor", "Resultado", "Coef. Recomendado"),
         caption   = "Prueba de Normalidad — Shapiro-Wilk")
Prueba de Normalidad — Shapiro-Wilk
Variable W p-valor Resultado Coef. Recomendado
Tiempo (X) 0.9547 0 ✗ No Normal Spearman / Kendall
Pago Total (Y) 0.8669 0 ✗ No Normal Spearman / Kendall
Ver código
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))

qqnorm(x_s, main = "Q-Q Plot: Tiempo (X)",
       col = "#2980B9", pch = 19, cex = 0.6)
qqline(x_s, col = "#E74C3C", lwd = 2)

qqnorm(y_s, main = "Q-Q Plot: Pago Total (Y)",
       col = "#27AE60", pch = 19, cex = 0.6)
qqline(y_s, col = "#E74C3C", lwd = 2)

Ver código
par(mfrow = c(1, 1))

Q-Q Plot: Tiempo (X): La variable Tiempo es una secuencia numérica uniforme (1, 2, 3, …, N), por lo que su distribución es inherentemente no normal. En el gráfico se observa una curva en S característica: los puntos siguen la línea diagonal en el centro pero se desvían sistemáticamente en ambos extremos, lo cual es consistente con una distribución uniforme. El Shapiro-Wilk confirma formalmente el rechazo de normalidad (p < 0.05).

Q-Q Plot: Pago Total (Y): Los puntos siguen la línea de referencia razonablemente bien en la zona central, pero se alejan hacia arriba en la cola derecha. Esto indica una distribución con sesgo positivo (cola derecha pesada), lo cual es típico en variables de costos y facturas donde la mayoría de los valores son moderados pero existen meses con pagos excepcionalmente altos. El Shapiro-Wilk rechaza la normalidad (p < 0.05), lo que justifica el uso de Spearman como coeficiente de correlación principal.

NoteCriterio de selección del coeficiente
  • Si ambas variables son normales → Se usa Pearson (mide relación lineal).
  • Si alguna no es normal → Se usa Spearman (basada en rangos).

3.3 Coeficientes de Correlación — Tiempo vs. Pago Total

Ver código
r_pearson  <- cor(x_s, y_s, method = "pearson")
r_spearman <- cor(x_s, y_s, method = "spearman")
r_kendall  <- cor(x_s, y_s, method = "kendall")

interpretar_r <- function(r) {
  r_abs <- abs(r)
  dir   <- ifelse(r > 0, "positiva", "negativa")
  if(r_abs > 0.7)       paste0("FUERTE ", dir)
  else if(r_abs > 0.4)  paste0("MODERADA ", dir)
  else if(r_abs > 0.0)  paste0("DÉBIL ", dir)
  else                   "Sin correlación lineal"
}

test_cor <- cor.test(x_s, y_s, method = "pearson")

tabla_cor <- data.frame(
  Método         = c("Pearson (r)", "Spearman (ρ)", "Kendall (τ)"),
  Coeficiente    = c(round(r_pearson,  4),
                     round(r_spearman, 4),
                     round(r_kendall,  4)),
  ``           = c(paste0(round(r_pearson^2  * 100, 2), "%"), "—", "—"),
  Interpretación = c(interpretar_r(r_pearson),
                     interpretar_r(r_spearman),
                     interpretar_r(r_kendall))
)

apa_table(tabla_cor,
         col.names = c("Método", "Coeficiente", "R²", "Interpretación"),
         caption   = "Coeficientes de Correlación — Tiempo vs. Pago Total")
Coeficientes de Correlación — Tiempo vs. Pago Total
Método Coeficiente Interpretación
Pearson (r) 0.7534 56.76% FUERTE positiva
Spearman (ρ) 0.8122 FUERTE positiva
Kendall (τ) 0.6522 MODERADA positiva
TipPrueba de Significancia — Pearson
  • t = 20.0408 | p-valor = 1.189627e-57
  • IC 95%: [0.7006, 0.7979]
  • Conclusión: Como el p-valor (1.189627e-57) es menor que 0.05, se rechaza H₀. Esto significa que la correlación entre Tiempo y Pago Total no es producto del azar — existe una relación lineal positiva real y estadísticamente significativa entre el tiempo transcurrido y el valor del pago mensual de la USCO.

Interpretación económica: Una correlación positiva entre Tiempo y Pago_Total confirma una tendencia de crecimiento estructural en la factura eléctrica de la USCO, explicada por la inflación tarifaria, el aumento de equipos instalados y la expansión de la planta física a lo largo de los años.


4 Modelo de Regresión Lineal Simple

Variable Independiente (\(X\)): Tiempo (período secuencial mensual) Variable Dependiente (\(Y\)): Pago Total mensual (COP)

4.1 Ecuación e Indicadores del Modelo

Ver código
modelo_simple <- lm(Pago_Total ~ Tiempo, data = dt)
resumen       <- summary(modelo_simple)
b0            <- coef(modelo_simple)[1]
b1            <- coef(modelo_simple)[2]

r2     <- resumen$r.squared
r2_adj <- resumen$adj.r.squared
f_stat <- resumen$fstatistic
p_f    <- pf(f_stat[1], f_stat[2], f_stat[3], lower.tail = FALSE)
rmse   <- sqrt(mean(residuals(modelo_simple)^2))

# Tabla de coeficientes
coef_tabla <- data.frame(
  Parámetro    = c("β₀ (Intercepto)", "β₁ (Pendiente: Tiempo)"),
  Estimador    = c(round(b0, 2),       round(b1, 2)),
  `Error Std`  = round(resumen$coefficients[, 2], 2),
  `t`          = round(resumen$coefficients[, 3], 4),
  `p-valor`    = format(resumen$coefficients[, 4], scientific = TRUE)
)

apa_table(coef_tabla,
         col.names = c("Parámetro", "Estimador (COP)", "Error Estándar", "t", "p-valor"),
         caption   = "Coeficientes del Modelo de Regresión Lineal Simple")
Coeficientes del Modelo de Regresión Lineal Simple
Parámetro Estimador (COP) Error Estándar t p-valor
β₀ (Intercepto) 7961642.2 5035465.53 1.5811 1.148851e-01
β₁ (Pendiente: Tiempo) 566119.3 28248.37 20.0408 1.189627e-57

4.2 Análisis de Varianza (ANOVA)

Ver código
# Tabla ANOVA como formato APA
anova_s    <- anova(modelo_simple)
anova_s_df <- as.data.frame(anova_s)
anova_s_df <- cbind(Fuente = rownames(anova_s_df), anova_s_df)
rownames(anova_s_df) <- NULL

# Formatear p-valor
anova_s_df$`Pr(>F)` <- ifelse(is.na(anova_s_df$`Pr(>F)`), "",
                               ifelse(anova_s_df$`Pr(>F)` < 0.001, "<.001",
                                      round(anova_s_df$`Pr(>F)`, 4)))
anova_s_df$`F value` <- ifelse(is.na(anova_s_df$`F value`), "",
                                round(anova_s_df$`F value`, 1))
anova_s_df$`Sum Sq`  <- format(round(anova_s_df$`Sum Sq`, 2),  big.mark = ",")
anova_s_df$`Mean Sq` <- format(round(anova_s_df$`Mean Sq`, 2), big.mark = ",")

apa_table(anova_s_df,
          col.names = c("Fuente", "df", "Suma de Cuadrados", "Media Cuadrática",
                        "F", "p"),
          caption   = "Prueba Omnibus ANOVA — Modelo Simple",
          align     = c("l", "c", "r", "r", "r", "r"))
Prueba Omnibus ANOVA — Modelo Simple
Fuente df Suma de Cuadrados Media Cuadrática F p
Tiempo 1 7.803372e+17 7.803372e+17 401.6 <.001
Residuals 306 5.945311e+17 1.942912e+15

4.3 Medidas de Ajuste del Modelo

Ver código
# Tabla de bondad de ajuste
tabla_ajuste <- data.frame(
  Indicador = c("R", "R²", "R² Ajustado", "F-estadístico",
                "p-valor (F)", "RMSE"),
  Valor     = c(round(sqrt(r2), 4), round(r2, 4), round(r2_adj, 4),
                round(f_stat[1], 2), format(p_f, scientific = TRUE),
                round(rmse, 2))
)

apa_table(tabla_ajuste,
         col.names = c("Indicador", "Valor"),
         caption   = "Medidas de Ajuste del Modelo",
         align     = c("l", "r"))
Medidas de Ajuste del Modelo
Indicador Valor
R 0.7534
0.5676
R² Ajustado 0.5662
F-estadístico 401.63
p-valor (F) 1.189627e-57
RMSE 43935133.28
ImportantEcuación del Modelo Simple

\[\hat{Y} = 7.961.642 + 566.119,3 \cdot Tiempo\]

Interpretación de los parámetros:

  • \(\hat{\beta}_0\) = 7.961.642 COP → Pago estimado en el período inicial (t = 0). Valor teórico base del modelo.
  • \(\hat{\beta}_1\) = 566.119,3 COP → Por cada mes adicional transcurrido, el pago total de la USCO aumenta en promedio 566.119,3 pesos colombianos. Esto refleja la tendencia inflacionaria y de crecimiento institucional.

4.4 Diagrama de Regresión Simple

Ver código
ggplot(dt, aes(x = Tiempo, y = Pago_Total / 1e6)) +
  geom_point(aes(color = Cumple), alpha = 0.7, size = 2.5) +
  scale_color_manual(values = c("Sí" = "#27AE60", "No" = "#E74C3C"),
                     name   = "Cumple norma\n(Reactiva/Activa ≤ 0.48)") +
  geom_smooth(method = "lm", color = "#2C3E50", se = TRUE,
              fill = "#D5D8DC", linewidth = 1.3) +
  labs(
    title    = "Regresión Lineal Simple: Tiempo → Pago Total Mensual",
    subtitle = paste0("Ŷ = ", format(round(b0/1e6, 3), big.mark=".", decimal.mark=","),
                      "M + ", format(round(b1/1e6, 4), big.mark=".", decimal.mark=","),
                      "M·t   |   R² = ", round(r2, 4)),
    x        = "Período (meses)",
    y        = "Pago Total (millones COP)",
    caption  = "Universidad Surcolombiana | Consumo mensual 2000–2025"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title    = element_text(face = "bold", color = "#2C3E50"),
    plot.subtitle = element_text(color = "#555555"),
    panel.grid.minor = element_blank(),
    legend.position  = "right"
  )


5 Verificación de Supuestos de Gauss-Markov

Para que los estimadores OLS sean BLUE (Best Linear Unbiased Estimator), los residuos deben cumplir cinco supuestos fundamentales:

  1. Linealidad: \(E(Y|X) = \beta_0 + \beta_1 X\)
  2. Media cero del error: \(E(\varepsilon_i) = 0\)
  3. Homocedasticidad: \(Var(\varepsilon_i) = \sigma^2\) (constante)
  4. Independencia: \(Cov(\varepsilon_i, \varepsilon_j) = 0\)
  5. Normalidad: \(\varepsilon_i \sim N(0, \sigma^2)\)
Ver código
res <- residuals(modelo_simple)
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

# 1. Residuos vs ajustados (linealidad + homocedasticidad)
plot(fitted(modelo_simple), res,
     main = "Residuos vs. Valores Ajustados\n(Linealidad y Homocedasticidad)",
     xlab = "Valores Ajustados (Ŷ)", ylab = "Residuos (e)",
     col  = "#2980B9", pch = 19, cex = 0.7)
abline(h = 0, col = "#E74C3C", lwd = 2, lty = 2)
lines(lowess(fitted(modelo_simple), res), col = "#27AE60", lwd = 2)
legend("topright", legend = c("e = 0", "LOWESS"),
       col = c("#E74C3C", "#27AE60"), lty = c(2,1), lwd = 2, bty = "n", cex = 0.8)

# 2. Q-Q Plot (normalidad de residuos)
qqnorm(res, main = "Q-Q Plot de Residuos\n(Normalidad)",
       col = "#2980B9", pch = 19, cex = 0.6)
qqline(res, col = "#E74C3C", lwd = 2)

# 3. Histograma de residuos
hist(res, main = "Histograma de Residuos",
     xlab = "Residuos", col = "#AED6F1", border = "#1A5276",
     probability = TRUE)
curve(dnorm(x, mean(res), sd(res)), add = TRUE, col = "#E74C3C", lwd = 2)

# 4. Residuos en el tiempo (independencia)
plot(dt$Tiempo, res,
     main = "Residuos en el Tiempo\n(Independencia)",
     xlab = "Período (meses)", ylab = "Residuos",
     col  = "#8E44AD", pch = 19, cex = 0.7)
abline(h = 0, col = "#E74C3C", lwd = 2, lty = 2)
lines(lowess(dt$Tiempo, res), col = "#27AE60", lwd = 2)

Ver código
par(mfrow = c(1, 1))

Residuos vs. Valores Ajustados: Si los puntos forman una nube aleatoria centrada en cero, el modelo es lineal y la varianza es constante. Patrones o forma de embudo indican problemas.

Q-Q Plot de Residuos: Si los puntos siguen la línea diagonal, los residuos son aproximadamente normales. Desviaciones en las colas indican asimetría.

Histograma de Residuos: Una forma simétrica y acampanada (similar a la curva roja) es consistente con normalidad.

Residuos en el Tiempo: Si se dispersan aleatoriamente alrededor de cero, hay independencia. Patrones o tendencias sugieren autocorrelación.

Ver código
res <- residuals(modelo_simple)
mean(res)
[1] -2.18033e-09

La media de los residuos es prácticamente cero, lo que confirma el cumplimiento de este supuesto. Cuando el modelo incluye intercepto (β₀), el método de mínimos cuadrados garantiza que los residuos sumen exactamente cero.

Conclusión: Supuesto cumplido. El modelo no tiene sesgo sistemático en sus predicciones.

Ver código
# Pruebas formales
bp_test <- bptest(modelo_simple)
dw_test <- dwtest(modelo_simple)
sw_res  <- shapiro.test(res)

tabla_sup <- data.frame(
  Supuesto    = c("Linealidad", "Media cero del error",
                  "Homocedasticidad", "Independencia", "Normalidad de residuos"),
  Prueba      = c("Inspección visual (Res. vs Ŷ)", "Media de residuos",
                  "Breusch-Pagan", "Durbin-Watson", "Shapiro-Wilk"),
  Estadístico = c("Ver gráfico",
                  paste0("ē = ", formatC(mean(res), format="e", digits=4)),
                  paste0("BP = ", round(bp_test$statistic, 4)),
                  paste0("DW = ", round(dw_test$statistic, 4)),
                  paste0("W = ",  round(sw_res$statistic,  4))),
  p_valor     = c("—",
                  "—",
                  round(bp_test$p.value,    6),
                  round(dw_test$p.value,    6),
                  round(sw_res$p.value,     6)),
  Resultado   = c("Inspección visual",
    "✓ Cumple (algebraico)",
    ifelse(bp_test$p.value    > 0.05, "✓ Cumple", "✗ No cumple"),
    ifelse(dw_test$p.value    > 0.05, "✓ Cumple", "✗ No cumple"),
    ifelse(sw_res$p.value     > 0.05, "✓ Cumple", "✗ No cumple"))
)

apa_table(tabla_sup,
         col.names = c("Supuesto", "Prueba", "Estadístico", "p-valor", "Resultado"),
         caption   = "Verificación de Supuestos de Gauss-Markov — Modelo Simple",
         align     = c("l", "l", "r", "r", "c")) |>
  column_spec(5, bold = TRUE)
Verificación de Supuestos de Gauss-Markov — Modelo Simple
Supuesto Prueba Estadístico p-valor Resultado
Linealidad Inspección visual (Res. vs Ŷ) Ver gráfico Inspección visual
Media cero del error Media de residuos ē = -2.1803e-09 ✓ Cumple (algebraico)
Homocedasticidad Breusch-Pagan BP = 58.9213 0 ✗ No cumple
Independencia Durbin-Watson DW = 0.7987 0 ✗ No cumple
Normalidad de residuos Shapiro-Wilk W = 0.904 0 ✗ No cumple

6 Modelo de Regresión Lineal Múltiple

Se incorpora el efecto de la energía activa (costo base de la factura) y el cumplimiento normativo (variable dummy que captura el recargo por reactiva) como predictores del Pago Total.

\[Y_i = \beta_0 + \beta_1 \cdot Activa_i + \beta_2 \cdot Incumple_i + \varepsilon_i\]

Donde \(Incumple_i = 1\) si la relación Reactiva/Activa > 0.48 ese mes (genera recargo), e \(Incumple_i = 0\) si cumple la norma.

6.1 Ajuste y Comparación de Modelos

Ver código
# Ajuste del modelo múltiple
modelo_mult  <- lm(Pago_Total ~ Activa + Incumple, data = dt)
resumen_mult <- summary(modelo_mult)
coef_m       <- coef(modelo_mult)

# Tabla de coeficientes
coef_mult_tabla <- data.frame(
  Parámetro    = c("β₀ (Intercepto)", "β₁ (Activa)", "β₂ (Incumple)"),
  Estimador    = round(resumen_mult$coefficients[, 1], 2),
  `Error Std`  = round(resumen_mult$coefficients[, 2], 2),
  `t`          = round(resumen_mult$coefficients[, 3], 4),
  `p-valor`    = format(resumen_mult$coefficients[, 4], scientific = TRUE)
)

apa_table(coef_mult_tabla,
         col.names = c("Parámetro", "Estimador (COP)", "Error Estándar", "t", "p-valor"),
         caption   = "Coeficientes del Modelo de Regresión Lineal Múltiple")
Coeficientes del Modelo de Regresión Lineal Múltiple
Parámetro Estimador (COP) Error Estándar t p-valor
β₀ (Intercepto) -25779592.61 6936125.39 -3.7167 2.400658e-04
β₁ (Activa) 613.74 32.52 18.8755 3.459999e-53
β₂ (Incumple) -73526310.78 45637535.15 -1.6111 1.081938e-01

6.2 Análisis de Varianza (ANOVA) — Modelo Múltiple

Ver código
# Tabla ANOVA como formato APA
anova_m    <- anova(modelo_mult)
anova_m_df <- as.data.frame(anova_m)
anova_m_df <- cbind(Fuente = rownames(anova_m_df), anova_m_df)
rownames(anova_m_df) <- NULL

anova_m_df$`Pr(>F)` <- ifelse(is.na(anova_m_df$`Pr(>F)`), "",
                               ifelse(anova_m_df$`Pr(>F)` < 0.001, "<.001",
                                      round(anova_m_df$`Pr(>F)`, 4)))
anova_m_df$`F value` <- ifelse(is.na(anova_m_df$`F value`), "",
                                round(anova_m_df$`F value`, 1))
anova_m_df$`Sum Sq`  <- format(round(anova_m_df$`Sum Sq`, 2),  big.mark = ",")
anova_m_df$`Mean Sq` <- format(round(anova_m_df$`Mean Sq`, 2), big.mark = ",")

apa_table(anova_m_df,
          col.names = c("Fuente", "df", "Suma de Cuadrados", "Media Cuadrática",
                        "F", "p"),
          caption   = "Prueba Omnibus ANOVA — Modelo Múltiple",
          align     = c("l", "c", "r", "r", "r", "r"))
Prueba Omnibus ANOVA — Modelo Múltiple
Fuente df Suma de Cuadrados Media Cuadrática F p
Activa 1 7.366541e+17 7.366541e+17 355 <.001
Incumple 1 5.385519e+15 5.385519e+15 2.6 0.1082
Residuals 305 6.328287e+17 2.074848e+15
ImportantEcuación del Modelo Múltiple

\[\hat{Y} = -25.779.593 + 613,74 \cdot Activa - 73.526.311 \cdot Incumple\]

Interpretación de los parámetros:

  • \(\hat{\beta}_0\) = -25.779.593 COP → Pago base cuando la energía activa es 0 y se cumple la norma. Valor teórico.
  • \(\hat{\beta}_1\) = 613,74 COP/kWh → Por cada kWh adicional de energía activa consumida, el pago total aumenta en promedio 613,74 pesos, manteniendo constante el cumplimiento normativo.
  • \(\hat{\beta}_2\) = -73.526.311 COP → El coeficiente negativo indica que los períodos donde la USCO incumple la norma reactiva están asociados a un pago menor en comparación con los períodos de cumplimiento. Sin embargo, este resultado debe interpretarse con cautela: solo existe un registro de incumplimiento en toda la serie, lo que hace que este coeficiente sea estadísticamente poco confiable y su signo posiblemente no refleje el efecto real del recargo tarifario.
Ver código
r2_s  <- summary(modelo_simple)$r.squared
r2_m  <- resumen_mult$r.squared
r2a_s <- summary(modelo_simple)$adj.r.squared
r2a_m <- resumen_mult$adj.r.squared

tabla_comp <- data.frame(
  Indicador = c("R²", "R² Ajustado", "AIC", "BIC", "RMSE"),
  Simple    = c(round(r2_s,  4), round(r2a_s, 4),
                round(AIC(modelo_simple), 2), round(BIC(modelo_simple), 2),
                round(sqrt(mean(residuals(modelo_simple)^2)), 2)),
  Multiple  = c(round(r2_m,  4), round(r2a_m, 4),
                round(AIC(modelo_mult), 2), round(BIC(modelo_mult), 2),
                round(sqrt(mean(residuals(modelo_mult)^2)), 2))
)

apa_table(tabla_comp,
         col.names = c("Indicador", "Modelo Simple (Tiempo → Pago)",
                       "Modelo Múltiple (Activa + Incumple → Pago)"),
         caption   = "Comparación de Modelos: Simple vs. Múltiple",
         align     = c("l", "r", "r"))
Comparación de Modelos: Simple vs. Múltiple
Indicador Modelo Simple (Tiempo → Pago) Modelo Múltiple (Activa + Incumple → Pago)
5.676000e-01 5.397000e-01
R² Ajustado 5.662000e-01 5.367000e-01
AIC 1.172057e+04 1.174180e+04
BIC 1.173176e+04 1.175672e+04
RMSE 4.393513e+07 4.532812e+07

6.3 Diagnóstico de Multicolinealidad (VIF)

Ver código
vif_vals  <- vif(modelo_mult)
tabla_vif <- data.frame(
  Variable  = names(vif_vals),
  VIF       = round(vif_vals, 4),
  Resultado = ifelse(vif_vals > 10, "⚠ Multicolinealidad severa", "✓ Aceptable")
)

apa_table(tabla_vif,
         col.names = c("Variable", "VIF", "Diagnóstico"),
         caption   = "Factor de Inflación de la Varianza (VIF)",
         align     = c("l", "r", "c"))
Factor de Inflación de la Varianza (VIF)
Variable VIF Diagnóstico
Activa 1.0006 ✓ Aceptable
Incumple 1.0006 ✓ Aceptable

Nota: La variable Incumple es una dummy binaria (0/1) con muy pocos casos de valor 1 (solo un período registra incumplimiento). Un VIF cercano a 1 confirma que no hay multicolinealidad entre Activa e Incumple.

6.4 Gráficos del Modelo Múltiple

Ver código
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))

# Predichos vs Reales — Modelo Simple
plot(fitted(modelo_simple) / 1e6, dt$Pago_Total / 1e6,
     main = "Modelo Simple\nPredichos vs. Reales",
     xlab = "Predichos (millones COP)", ylab = "Reales (millones COP)",
     col  = "#2980B9", pch = 19, cex = 0.7)
abline(a = 0, b = 1, col = "#E74C3C", lwd = 2, lty = 2)

# Predichos vs Reales — Modelo Múltiple
colores_c2 <- ifelse(dt$Incumple == 1, "#E74C3C", "#27AE60")
plot(fitted(modelo_mult) / 1e6, dt$Pago_Total / 1e6,
     main = "Modelo Múltiple\nPredichos vs. Reales",
     xlab = "Predichos (millones COP)", ylab = "Reales (millones COP)",
     col  = colores_c2, pch = 19, cex = 0.7)
abline(a = 0, b = 1, col = "#E74C3C", lwd = 2, lty = 2)
legend("topleft", legend = c("Cumple", "Incumple"),
       col = c("#27AE60", "#E74C3C"), pch = 19, bty = "n", cex = 0.8)

Ver código
par(mfrow = c(1, 1))
Ver código
# Visualizar el efecto de Incumple sobre Pago dado nivel de Activa
ggplot(dt, aes(x = Activa / 1000, y = Pago_Total / 1e6, color = factor(Incumple))) +
  geom_point(alpha = 0.7, size = 2.5) +
  geom_smooth(method = "lm", se = TRUE, linewidth = 1.2) +
  scale_color_manual(
    values = c("0" = "#27AE60", "1" = "#E74C3C"),
    labels = c("0" = "Cumple norma", "1" = "Incumple norma"),
    name   = "Estado normativo"
  ) +
  labs(
    title    = "Efecto de Energía Activa y Cumplimiento sobre el Pago Total",
    subtitle = "Líneas paralelas muestran el sobrecosto por incumplimiento (β₂)",
    x        = "Energía Activa (miles de kWh)",
    y        = "Pago Total (millones COP)",
    caption  = "Universidad Surcolombiana | Consumo mensual 2000–2025"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title    = element_text(face = "bold", color = "#2C3E50"),
    plot.subtitle = element_text(color = "#555555"),
    panel.grid.minor = element_blank()
  )


7 Selección del Modelo — Criterio AIC

El Criterio de Información de Akaike (AIC) permite comparar modelos penalizando la complejidad:

\[AIC = 2k - 2\ln(\hat{L})\]

Donde \(k\) es el número de parámetros y \(\hat{L}\) es el ajuste del modelo. El modelo con menor AIC es preferible.

Ver código
# Modelos candidatos
modelos <- list(
  "M1: Pago ~ Tiempo"                   = lm(Pago_Total ~ Tiempo,           data = dt),
  "M2: Pago ~ Activa"                   = lm(Pago_Total ~ Activa,           data = dt),
  "M6: Pago ~ Activa + Tiempo"          = lm(Pago_Total ~ Activa + Tiempo,  data = dt),
  "M3: Pago ~ Activa + Incumple"        = lm(Pago_Total ~ Activa + Incumple, data = dt),
  "M4: Pago ~ Tiempo + Incumple"        = lm(Pago_Total ~ Tiempo + Incumple, data = dt),
  "M5: Pago ~ Activa + Tiempo + Incumple" = lm(Pago_Total ~ Activa + Tiempo + Incumple, data = dt)
)

tabla_aic <- data.frame(
  Modelo = names(modelos),
  Vars   = sapply(modelos, function(m) length(coef(m))),
  R2     = round(sapply(modelos, function(m) summary(m)$r.squared),     4),
  R2_Adj = round(sapply(modelos, function(m) summary(m)$adj.r.squared), 4),
  AIC    = round(sapply(modelos, AIC), 2),
  BIC    = round(sapply(modelos, BIC), 2),
  RMSE   = round(sapply(modelos, function(m) sqrt(mean(residuals(m)^2))), 2)
)
tabla_aic <- tabla_aic[order(tabla_aic$AIC), ]
tabla_aic$Ranking <- 1:nrow(tabla_aic)

apa_table(tabla_aic,
         col.names = c("Modelo", "Parámetros", "R²", "R² corregida",
                       "AIC", "BIC", "RMSE", "Ranking"),
         caption   = "Comparación de Modelos por AIC",
         align     = c("l", "c", "r", "r", "r", "r", "r", "c"))
Comparación de Modelos por AIC
Modelo Parámetros R² corregida AIC BIC RMSE Ranking
M6: Pago ~ Activa + Tiempo 3 0.7601 0.7585 11541.12 11556.04 32725037 1
M5: Pago ~ Activa + Tiempo + Incumple 4 0.7602 0.7578 11542.98 11561.63 32717654 2
M1: Pago ~ Tiempo 2 0.5676 0.5662 11720.57 11731.76 43935133 3
M4: Pago ~ Tiempo + Incumple 3 0.5680 0.5652 11722.26 11737.18 43912876 4
M3: Pago ~ Activa + Incumple 3 0.5397 0.5367 11741.80 11756.72 45328124 5
M2: Pago ~ Activa 2 0.5358 0.5343 11742.41 11753.60 45520591 6
Ver código
# Selección automática por stepwise
modelo_full <- lm(Pago_Total ~ Activa + Tiempo + Incumple, data = dt)
modelo_step <- stepAIC(modelo_full, direction = "both", trace = FALSE)

cat("Modelo seleccionado por stepAIC:\n")
Modelo seleccionado por stepAIC:
Ver código
print(formula(modelo_step))
Pago_Total ~ Activa + Tiempo
Ver código
cat(sprintf("AIC final : %.2f\n", AIC(modelo_step)))
AIC final : 11541.12
Ver código
cat(sprintf("R²        : %.4f\n", summary(modelo_step)$r.squared))
R²        : 0.7601
Ver código
cat(sprintf("R² Ajust. : %.4f\n", summary(modelo_step)$adj.r.squared))
R² Ajust. : 0.7585
TipModelo Recomendado

El modelo seleccionado por stepwise AIC es: Pago_Total ~ Activa + Tiempo con AIC = 1.154112^{4} y R² = 0.7601.


8 Conclusiones

1. Correlación (Tiempo → Pago Total): La correlación de Pearson es \(r\) = 0.7534 (\(R^2\) = 56.76%), lo que indica una correlación FUERTE positiva entre el tiempo transcurrido y el valor del pago. Esto confirma una tendencia estructural de crecimiento en la factura eléctrica de la USCO, atribuible a la inflación tarifaria y al crecimiento de la planta física.

2. Modelo Simple (Tiempo → Pago): La ecuación estimada es \(\hat{Y}\) = 7.961.642 + 566.119,3 · Tiempo, con \(R^2\) = 0.5676. Por cada mes adicional, el pago mensual de la USCO aumenta en promedio 566.119,3 pesos.

3. Modelo Múltiple (Activa + Incumple → Pago): El modelo múltiple incorpora el consumo de energía activa y el incumplimiento normativo como predictores. El coeficiente de Incumple resultó negativo y no significativo (p = 0.1082), por lo que su interpretación debe tomarse con cautela dado que solo existe un registro de incumplimiento en la serie.

4. Selección AIC: El modelo con menor AIC (11541.12) y mejor balance entre ajuste y número de variables es el recomendado por stepwise: Pago_Total ~ Activa + Tiempo.

5. Cumplimiento Normativo: El 99.68% de los períodos registrados cumplen con la relación Reactiva/Activa ≤ 0.48. Solo existe un período de incumplimiento en toda la serie, lo que limita el análisis del efecto del recargo tarifario a nivel agregado.

6. Limitación del nivel de agregación: El análisis se realizó sobre el total mensual de todas las cuentas de la universidad. Al agregar por período, los incumplimientos a nivel de cuenta individual pueden quedar diluidos en el total, subestimando la frecuencia real de penalizaciones por energía reactiva. Un análisis más preciso requeriría trabajar a nivel de sede o cuenta individual, donde el recargo se aplica directamente.


9 Limitaciones del Modelo y Perspectivas

9.1 Por qué los supuestos de Gauss-Markov no se cumplen plenamente

El análisis de regresión lineal es válido como marco metodológico para identificar la tendencia estructural del gasto energético de la USCO. Sin embargo, la verificación de supuestos revela incumplimientos que tienen explicaciones operativas concretas, no simplemente estadísticas:

Heterocedasticidad (Breusch-Pagan, p ≈ 0): La varianza del pago no es constante en el tiempo. A medida que la universidad crece en capacidad instalada y las tarifas eléctricas aumentan, los pagos crecen también en escala — los meses recientes tienen fluctuaciones absolutas mucho mayores que los primeros años de la serie. OLS trata todos los períodos con el mismo peso, lo que subestima la incertidumbre en los valores más recientes.

Autocorrelación (Durbin-Watson = 0.669): El consumo eléctrico de la USCO no es aleatorio entre meses: depende fuertemente del calendario académico. La universidad opera aproximadamente 32 semanas al año (8 meses de clases). Los meses de vacaciones (diciembre–enero y julio) presentan consumos sistemáticamente menores, generando una estacionalidad anual de período s = 12 que el modelo lineal simple no captura.

Pandemia 2020: A partir de marzo de 2020, el cierre físico de la universidad produjo una caída abrupta del consumo durante aproximadamente 6 a 8 períodos, con recuperación gradual hacia 2022. Este quiebre no es ruido aleatorio sino un evento externo identificable. Al no modelarlo, OLS distorsiona la pendiente β₁ y genera los residuos más grandes de toda la serie en ese período.

No-normalidad de residuos (Shapiro-Wilk, p ≈ 0): Consecuencia directa de los dos puntos anteriores: outliers por correcciones de factura y el impacto de la pandemia producen colas pesadas en la distribución de los errores.La prueba de Shapiro-Wilk rechazó la normalidad tanto en las variables originales como en los residuos del modelo. Cuando los datos no cumplen el supuesto de normalidad, una alternativa es el uso de métodos no paramétricos, que no requieren asumir ninguna distribución. Para el análisis de correlación ya se reportaron Spearman y Kendall como alternativas no paramétricas a Pearson.

En conclusión, el modelo lineal simple cumple su propósito pedagógico de cuantificar la tendencia de largo plazo y aplicar la metodología de Gauss-Markov. Sus limitaciones no invalidan el ejercicio — al contrario, los supuestos incumplidos son informativos: señalan exactamente qué estructura del fenómeno real (estacionalidad, pandemia, crecimiento por etapas) el modelo no alcanza a representar y qué extensiones metodológicas serían necesarias para un análisis completo.


Análisis desarrollado para la Universidad Surcolombiana (USCO) – Consumo Energético 2000–2025 Herramienta: R R version 4.5.3 (2026-03-11 ucrt)