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 11, 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. Los datos provienen del archivo Consumo_Final.csv, que consolida el registro mensual de energía eléctrica de la Universidad Surcolombiana. Para su construcción, se totalizaron todas las cuentas de la universidad por período, obteniendo la energía activa y reactiva agregada de cada mes. A partir de estos totales se calculó la relación Reactiva/Activa y se determinó si cada período cumple o no cumple con el límite normativo establecido (Reactiva/Activa ≤ 0.48).

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_consumo <- dt_consumo[dt_consumo$Relacion > 0, ]   # registros con consumo real
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 = ",")
}

cat("Registros válidos para el análisis:", nrow(dt), "\n")
Registros válidos para el análisis: 292
Ver código 
cat("Período:", min(dt$Ano), "–", max(dt$Ano), "\n")
Período: 2001 – 2025
Ver código 
cat("Cuentas con incumplimiento normativo:", sum(dt$Incumple), "\n\n")
Cuentas con incumplimiento normativo: 1
Ver código 
head(dt[, c("Ano","Mes","Activa","Reactiva","Relacion","Cumple","Pago_Total","Tiempo")], 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 <- if(n <= 5000) shapiro.test(x_s) else { set.seed(42); shapiro.test(sample(x_s, 5000)) }
sw_y <- if(n <= 5000) shapiro.test(y_s) else { set.seed(42); shapiro.test(sample(y_s, 5000)) }

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",
         nota      = "Criterio: p > 0.05 → No se rechaza H₀ de normalidad.")
Prueba de Normalidad — Shapiro-Wilk
Variable W p-valor Resultado Coef. Recomendado
Tiempo (X) 0.9546 0 ✗ No Normal Spearman / Kendall
Pago Total (Y) 0.8308 0 ✗ No Normal Spearman / Kendall
Nota. Criterio: p > 0.05 → No se rechaza H₀ de normalidad.
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 (monótona, por rangos) como principal y Kendall como confirmación robusta.

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)),
  `R²`           = 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 R² Interpretación
Pearson (r) 0.7372 54.35% FUERTE positiva
Spearman (ρ) 0.8006 — FUERTE positiva
Kendall (τ) 0.6287 — MODERADA positiva
TipPrueba de Significancia — Pearson
  • t = 18.5815 | p-valor = 2.630615e-51
  • IC 95%: [0.68, 0.7855]
  • Conclusión: Se rechaza H₀ → Correlación estadísticamente SIGNIFICATIVA entre Tiempo y Pago Total.

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",
         nota      = paste0("Modelo estimado con N = ", nrow(dt), " observaciones."))
Coeficientes del Modelo de Regresión Lineal Simple
Parámetro Estimador (COP) Error Estándar t p-valor
β₀ (Intercepto) 9090202.1 5554378.40 1.6366 1.028023e-01
β₁ (Pendiente: Tiempo) 566537.9 30489.35 18.5815 2.630615e-51
Nota. Modelo estimado con N = 292 observaciones.

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"),
          nota      = "Suma de cuadrados tipo I (secuencial).")
Prueba Omnibus ANOVA — Modelo Simple
Fuente df Suma de Cuadrados Media Cuadrática F p
Tiempo 1 6.659172e+17 6.659172e+17 345.3 <.001
Residuals 290 5.593152e+17 1.928673e+15
Nota. Suma de cuadrados tipo I (secuencial).

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"),
         nota      = paste0("Models estimated using sample size of N = ", nrow(dt), "."))
Medidas de Ajuste del Modelo
Indicador Valor
R 0.7372
R² 0.5435
R² Ajustado 0.5419
F-estadístico 345.27
p-valor (F) 2.630615e-51
RMSE 43766001.91
Nota. Models estimated using sample size of N = 292.
ImportantEcuación del Modelo Simple

\[\hat{Y} = 9.090.202 + 566.537,9 \cdot Tiempo\]

Interpretación de los parámetros:

  • \(\hat{\beta}_0\) = 9.090.202 COP → Pago estimado en el período inicial (t = 0). Valor teórico base del modelo.
  • \(\hat{\beta}_1\) = 566.537,9 COP → Por cada mes adicional transcurrido, el pago total de la USCO aumenta en promedio 566.537,9 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: Permite verificar linealidad y homocedasticidad. Se espera una nube aleatoria centrada en cero sin patrones sistemáticos y con dispersión constante. Si la línea LOWESS se aleja de e = 0 o la varianza aumenta con Ŷ, hay indicios de no-linealidad o heterocedasticidad.

Q-Q Plot de Residuos: Contrasta los cuantiles de los residuos con los de una distribución normal teórica. Si los puntos siguen la línea diagonal, los errores son aproximadamente normales. Desviaciones en las colas indican asimetría o colas pesadas.

Histograma de Residuos: Muestra la distribución empírica de los errores. La curva roja superpuesta es la densidad normal con la misma media y desviación estándar. Una forma aproximadamente simétrica y acampanada es consistente con el supuesto de normalidad.

Residuos en el Tiempo: Evalúa el supuesto de independencia. Los residuos deben dispersarse aleatoriamente alrededor de cero sin tendencias, ciclos ni aumento progresivo de la varianza. Patrones sistemáticos sugieren autocorrelación (confirmable con Durbin-Watson).

Ver código 
# Pruebas
t_test_res <- t.test(res, mu = 0)
bp_test    <- bptest(modelo_simple)
dw_test    <- dwtest(modelo_simple)
sw_res     <- if(length(res) <= 5000) shapiro.test(res) else {
                set.seed(42); shapiro.test(sample(res, 5000))
              }

tabla_sup <- data.frame(
  Supuesto    = c("Linealidad", "Media cero del error",
                  "Homocedasticidad", "Independencia", "Normalidad de residuos"),
  Prueba      = c("Inspección visual (Res. vs Ŷ)", "Prueba t (μ = 0)",
                  "Breusch-Pagan", "Durbin-Watson", "Shapiro-Wilk"),
  Estadístico = c("Ver gráfico",
                  paste0("t = ", round(t_test_res$statistic, 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(t_test_res$p.value, 6),
                  round(bp_test$p.value,    6),
                  round(dw_test$p.value,    6),
                  round(sw_res$p.value,     6)),
  Resultado   = c("Inspección visual",
    ifelse(t_test_res$p.value > 0.05, "✓ Cumple", "✗ No cumple"),
    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"),
         nota      = "Criterio de decisión: p > 0.05 → No se rechaza H₀.") |>
  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 Prueba t (μ = 0) t = 0 1 ✓ Cumple
Homocedasticidad Breusch-Pagan BP = 56.1137 0 ✗ No cumple
Independencia Durbin-Watson DW = 0.669 0 ✗ No cumple
Normalidad de residuos Shapiro-Wilk W = 0.9087 0 ✗ No cumple
Nota. Criterio de decisión: p > 0.05 → No se rechaza H₀.
WarningTratamiento de Supuestos Incumplidos
  • Heterocedasticidad → Usar errores estándar robustos de White (HC1) con el paquete sandwich.
  • Autocorrelación → Dado que los datos son series de tiempo, es esperable. Solución: estimador de Newey-West o mínimos cuadrados generalizados (MCG).
  • No normalidad → Considerar transformación logarítmica: log(Pago_Total), ya que el pago tiene distribución asimétrica positiva típica de datos económicos.

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",
         nota      = paste0("Modelo estimado con N = ", nrow(dt), " observaciones."))
Coeficientes del Modelo de Regresión Lineal Múltiple
Parámetro Estimador (COP) Error Estándar t p-valor
β₀ (Intercepto) -37410095.47 8643732.14 -4.3280 2.075891e-05
β₁ (Activa) 664.62 39.48 16.8357 8.509317e-45
β₂ (Incumple) -73656560.20 46288625.95 -1.5912 1.126476e-01
Nota. Modelo estimado con N = 292 observaciones.

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"),
          nota      = "Suma de cuadrados tipo I (secuencial).")
Prueba Omnibus ANOVA — Modelo Múltiple
Fuente df Suma de Cuadrados Media Cuadrática F p
Activa 1 6.029682e+17 6.029682e+17 282.5 <.001
Incumple 1 5.404590e+15 5.404590e+15 2.5 0.1126
Residuals 289 6.168595e+17 2.134462e+15
Nota. Suma de cuadrados tipo I (secuencial).
ImportantEcuación del Modelo Múltiple

\[\hat{Y} = -37.410.095 + 664,62 \cdot Activa + -73.656.560 \cdot Incumple\]

Interpretación de los parámetros:

  • \(\hat{\beta}_0\) = -37.410.095 COP → Pago base cuando la energía activa es 0 y se cumple la norma. Valor teórico.
  • \(\hat{\beta}_1\) = 664,62 COP/kWh → Por cada kWh adicional de energía activa consumida, el pago total aumenta en promedio 664,62 pesos, manteniendo constante el cumplimiento normativo.
  • \(\hat{\beta}_2\) = -73.656.560 COP → El incumplimiento de la norma reactiva genera un sobrecosto adicional de -73.656.560 pesos por período, en comparación con los meses donde la USCO cumple la norma (manteniendo el consumo activo constante).
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)
R² 5.435000e-01 4.965000e-01
R² Ajustado 5.419000e-01 4.931000e-01
AIC 1.110977e+04 1.114037e+04
BIC 1.112080e+04 1.115507e+04
RMSE 4.376600e+07 4.596230e+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"),
         nota      = "Criterio: VIF > 10 indica multicolinealidad severa.")
Factor de Inflación de la Varianza (VIF)
Variable VIF Diagnóstico
Activa 1.0004 ✓ Aceptable
Incumple 1.0004 ✓ Aceptable
Nota. Criterio: VIF > 10 indica multicolinealidad severa.

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 la verosimilitud maximizada. 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),
  "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"),
         nota      = "Modelos ordenados de menor a mayor AIC. Menor AIC = mejor modelo.")
Comparación de Modelos por AIC
Modelo Parámetros R² R² corregida AIC BIC RMSE Ranking
M5: Pago ~ Activa + Tiempo + Incumple 4 0.7661 0.7637 10918.50 10936.88 31327141 1
M1: Pago ~ Tiempo 2 0.5435 0.5419 11109.77 11120.80 43766002 2
M4: Pago ~ Tiempo + Incumple 3 0.5440 0.5408 11111.48 11126.19 43744524 3
M3: Pago ~ Activa + Incumple 3 0.4965 0.4931 11140.37 11155.07 45962296 4
M2: Pago ~ Activa 2 0.4921 0.4904 11140.91 11151.94 46163206 5
Nota. Modelos ordenados de menor a mayor AIC. Menor AIC = mejor modelo.
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 : 10916.55
Ver código 
cat(sprintf("R²        : %.4f\n", summary(modelo_step)$r.squared))
R²        : 0.7661
Ver código 
cat(sprintf("R² Ajust. : %.4f\n", summary(modelo_step)$adj.r.squared))
R² Ajust. : 0.7644
TipModelo Recomendado

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

8 Conclusiones

1. Correlación (Tiempo → Pago Total): La correlación de Pearson es \(r\) = 0.7372 (\(R^2\) = 54.35%), 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}\) = 9.090.202 + 566.537,9 · Tiempo, con \(R^2\) = 0.5435. Por cada mes adicional, el pago mensual de la USCO aumenta en promedio 566.537,9 pesos.

3. Modelo Múltiple (Activa + Incumple → Pago): El modelo múltiple captura dos efectos simultáneos: (i) el costo estructural asociado al nivel de consumo activo (\(\beta_1\)), y (ii) el sobrecosto por incumplimiento normativo (\(\beta_2\) = -73.656.560 COP). El \(R^2\) = 0.4965 representa una variación respecto al modelo simple.

4. Selección AIC: El modelo con menor AIC (1.091655^{4}) y mejor balance entre ajuste y parsimonia es el recomendado por stepwise: Pago_Total ~ Activa + Tiempo.

5. Cumplimiento Normativo: El 99.66% de los períodos registrados cumplen con la relación Reactiva/Activa ≤ 0.48. El único período de incumplimiento registra un pago notablemente superior al esperado por su nivel de consumo activo, lo cual es consistente con la lógica del recargo tarifario.

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. El residuo de un mes predice el del siguiente, violando el supuesto de independencia.

Quiebre estructural — 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.

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)