library(tidyverse)
library(broom) # para tidy()
library(NobBS) # trae los datos de dengue
library(scoringutils) # evaluación de pronósticos
library(forecast) # ARIMA
library(bsts) # modelos estructurales bayesianos
library(WaveletComp) # escalograma14 Series de tiempo
Dengue en Puerto Rico: de la regresión lineal al ARIMA y al bsts
NotaResumen
Usamos veinte años de casos semanales de dengue para entender por qué una regresión lineal no sirve para pronosticar series de tiempo. Descomponemos la estacionalidad con un escalograma, y vemos cómo el ARIMA y el bsts resuelven los dos problemas que la regresión lineal no puede: la dependencia entre observaciones y el crecimiento de la incertidumbre. Comparamos los tres con el WIS.
AdvertenciaDos detalles de instalación que te van a ahorrar un mal rato
El orden de las librerías importa aquí. scoringutils y forecast definen los dos una función para imprimir pronósticos, y la que se carga al último gana. Si pones scoringutils después de forecast, los pronósticos del ARIMA dejan de imprimirse bien. Cárgalo antes, como arriba.
Y hay que apagar el paralelismo de data.table:
data.table::setDTthreads(1)scoringutils usa data.table por dentro, y en algunas computadoras el paralelismo de data.table choca con el de los paquetes bayesianos (bsts trae su propio motor en C++) y tumba la sesión de R sin dar error. Ponerlo en un solo hilo lo resuelve y no lo vas a notar: aquí no hay nada tan pesado como para necesitar varios núcleos.
Tip
Este capítulo es la continuación natural de Modelos y de Regresiones (parte 1). Aquí vas a ver el flujo de trabajo de diez pasos aplicado de principio a fin a un problema real.
14.1 Los datos
Usaremos los datos de dengue de Puerto Rico que vienen dentro del paquete NobBS (McGough et al. 2026). La ventaja de que vengan en un paquete es que no dependen de que un sitio web siga vivo: se instalan con install.packages("NobBS") y ya.
McGough, Sarah, Nicolas Menzies, Marc Lipsitch, Michael Johansson, Rami Yaari, y Rodrigo Zepeda Tello. 2026. NobBS: Nowcasting by Bayesian Smoothing. https://doi.org/10.32614/CRAN.package.NobBS.
data(denguedat, package = "NobBS")
head(denguedat) onset_week report_week gender
1 1990-01-01 1990-01-01 Male
2 1990-01-01 1990-01-01 Female
3 1990-01-01 1990-01-01 Female
4 1990-01-01 1990-01-08 Female
5 1990-01-01 1990-01-08 Male
6 1990-01-01 1990-01-15 FemaleNota que no es una serie de tiempo todavía: es una lista de casos, un renglón por persona enferma, con la semana en que empezaron los síntomas (onset_week) y la semana en que el caso se reportó (report_week).
nrow(denguedat)[1] 52987Para tener una serie de tiempo hay que contar cuántos casos hubo por semana:
dengue <- denguedat %>%
count(onset_week, name = "casos") %>%
rename(semana = onset_week) %>%
arrange(semana)
head(dengue) semana casos
1 1990-01-01 61
2 1990-01-08 50
3 1990-01-15 44
4 1990-01-22 46
5 1990-01-29 39
6 1990-02-05 3414.1.1 El retraso de reporte y por qué recortamos el final
Los casos no se reportan el mismo día que la persona se enferma. Veamos cuánto tardan:
retraso <- as.numeric(denguedat$report_week - denguedat$onset_week) / 7
c(mediana = median(retraso), p95 = quantile(retraso, 0.95, names = FALSE))mediana p95
1 4
AdvertenciaLas últimas semanas siempre parecen tener menos casos
La mitad de los casos se reportan en la misma semana o la siguiente, pero el 5% tarda 4 semanas o más. Eso significa que las últimas semanas de cualquier base de vigilancia están incompletas: no es que haya menos enfermos, es que sus casos todavía no llegan.
Si no lo tomas en cuenta, vas a concluir que un brote está cediendo justo cuando más está creciendo. Éste es el error más peligroso de la vigilancia epidemiológica en tiempo real, y es tan importante que hay paquetes enteros dedicados a corregirlo (a eso se dedica NobBS: se llama nowcasting).
Como aquí queremos modelar y no hacer nowcasting, simplemente quitamos las últimas 6 semanas, que son las que están claramente truncadas.
dengue <- head(dengue, nrow(dengue) - 6)
nrow(dengue)[1] 1085range(dengue$semana)[1] "1990-01-01" "2010-10-18"14.1.2 Guardar los datos en csv
Ya que la tenemos armada, guardemos la serie para no depender del paquete:
write_csv(dengue, "datasets/dengue_puerto_rico_semanal.csv")Y así se vería para volver a leerla después:
dengue <- read_csv("datasets/dengue_puerto_rico_semanal.csv")14.2 Mirar la serie antes de modelar
Nunca modeles sin graficar primero (es el paso 2 del flujo de trabajo).
ggplot(dengue) +
geom_line(aes(x = semana, y = casos), color = "#003f5c", linewidth = 0.4) +
labs(
x = "", y = "Casos semanales",
title = "Dengue en Puerto Rico",
caption = "Fuente: datos denguedat del paquete NobBS"
) +
scale_y_continuous(labels = scales::label_comma()) +
theme_bw()
A ojo se ven dos cosas: hay picos que se repiten todos los años y hay epidemias grandes cada varios años. Vamos a separar esas dos señales.
14.3 Deconstruir la estacionalidad
14.3.1 El efecto de la semana del año
La forma más directa: ¿en qué semanas del año hay más dengue?
dengue %>%
mutate(semana_anio = as.numeric(format(semana, "%V"))) %>%
ggplot() +
geom_boxplot(aes(x = semana_anio, y = casos, group = semana_anio),
fill = "#bc5090", alpha = 0.6, outlier.size = 0.5
) +
labs(
x = "Semana del año", y = "Casos",
title = "El dengue tiene una estación bien marcada"
) +
theme_bw()
El patrón es clarísimo: los casos suben en la segunda mitad del año. Esto es estacionalidad: una parte de la serie que se repite con periodo fijo.
14.3.2 Separar tendencia, estacionalidad y resto
STL descompone la serie en tres pedazos que suman la serie original:
descomposicion <- stl(ts(dengue$casos, frequency = 52), s.window = "periodic")
# Las tres partes regresan al data.frame como columnas nuevas
dengue <- dengue |>
mutate(
tendencia = as.numeric(descomposicion$time.series[, "trend"]),
estacional = as.numeric(descomposicion$time.series[, "seasonal"]),
resto = as.numeric(descomposicion$time.series[, "remainder"])
)
# Y como son columnas, se grafican como cualquier otra cosa
glimpse(dengue)Rows: 1,085
Columns: 5
$ semana <date> 1990-01-01, 1990-01-08, 1990-01-15, 1990-01-22, 1990-01-29…
$ casos <int> 61, 50, 44, 46, 39, 34, 24, 17, 17, 16, 17, 22, 15, 19, 4, …
$ tendencia <dbl> 36.18225, 36.19914, 36.21603, 36.23292, 36.24981, 36.26670,…
$ estacional <dbl> 0.5226546, -0.9363916, -0.1097235, -0.7592459, -4.6944824, …
$ resto <dbl> 24.29509602, 14.73725289, 7.89369550, 10.52632848, 7.444675…
NotaLa única función de este capítulo que pide un
ts
stl() es de las pocas que exige un objeto ts (el formato viejo de series de tiempo en R). En vez de arrastrar ese objeto por todo el capítulo, lo creamos dentro de la llamada y de inmediato guardamos los resultados como columnas del data.frame.
De aquí en adelante todo sale del data.frame, incluidos el ARIMA y el bsts.
dengue |>
select(semana, casos, tendencia, estacional, resto) |>
pivot_longer(-semana, names_to = "parte", values_to = "valor") |>
mutate(parte = factor(parte,
levels = c("casos", "tendencia", "estacional", "resto")
)) |>
ggplot(aes(x = semana, y = valor)) +
geom_line(color = "#003f5c", linewidth = 0.3) +
facet_wrap(~parte, ncol = 1, scales = "free_y") +
labs(x = "", y = "", title = "Descomposición de la serie de dengue") +
theme_bw()
Y comprobamos que de verdad suman la serie original:
dengue |>
summarise(
diferencia_maxima = max(abs(
(tendencia + estacional + resto) - casos
))
) diferencia_maxima
1 2.842171e-14
Nota
Fíjate en las escalas de cada panel: son distintas. La barra gris del lado derecho te dice el tamaño relativo de cada componente. Si la barra de un panel es grande, ese componente aporta poco a la serie total.
14.3.3 El escalograma
El boxplot de arriba supone que el ciclo dura exactamente un año. ¿Y si hay otros ciclos? ¿Y si la estacionalidad cambia con el tiempo?
Para eso sirve un escalograma (o wavelet power spectrum): en lugar de preguntar “¿qué tan fuerte es el ciclo anual?”, pregunta “¿qué ciclos hay, y en qué épocas?”. El eje horizontal es el tiempo, el vertical es la duración del ciclo, y el color es qué tan fuerte está ese ciclo en ese momento.
onda <- analyze.wavelet(
data.frame(x = dengue$casos), "x",
loess.span = 0, # no quitar tendencia
dt = 1, # cada observación es 1 semana
dj = 1 / 12,
lowerPeriod = 8, # de 8 semanas...
upperPeriod = 512, # ...a ~10 años
make.pval = TRUE,
n.sim = 10,
verbose = FALSE
)wt.image(onda,
main = "¿Qué ciclos tiene el dengue y cuándo?",
periodlab = "Periodo del ciclo (semanas)",
timelab = "Semanas desde 1990",
legend.params = list(lab = "Potencia")
)
¿Cuál es el ciclo dominante?
tibble(periodo_semanas = onda$Period, potencia = onda$Power.avg) %>%
slice_max(potencia, n = 3) %>%
mutate(periodo_anios = periodo_semanas / 52)# A tibble: 3 × 3
periodo_semanas potencia periodo_anios
<dbl> <dbl> <dbl>
1 50.8 0.391 0.977
2 53.8 0.369 1.03
3 47.9 0.348 0.922
ImportanteQué nos dice el escalograma
El periodo dominante está alrededor de las 51 semanas, o sea un año. Eso confirma numéricamente lo que vimos en el boxplot.
Pero el escalograma dice algo que el boxplot no puede: la banda anual no es igual de intensa todo el tiempo. Hay épocas donde el ciclo anual es muy fuerte y otras donde casi se apaga, y además se alcanzan a ver franjas en periodos más largos (de varios años), que son las epidemias grandes.
Esto importa para modelar: un modelo que suponga que la estacionalidad es siempre la misma va a quedar corto. Guárdate esta idea, porque es justamente la diferencia entre el ARIMA que vamos a ajustar y el bsts.
14.3.4 Términos de Fourier: la estacionalidad como columnas
Antes de partir la base agregamos las columnas que van a representar la estacionalidad. Suena mucho más complicado de lo que es: son senos y cosenos con periodo de un año, y se agregan como columnas normales.
K <- 3 # cuántos pares de seno y coseno
periodo <- 52 # una vuelta completa = 52 semanas
dengue <- dengue |>
mutate(t = row_number()) |>
mutate(
sin1 = sin(2 * pi * 1 * t / periodo),
cos1 = cos(2 * pi * 1 * t / periodo),
sin2 = sin(2 * pi * 2 * t / periodo),
cos2 = cos(2 * pi * 2 * t / periodo),
sin3 = sin(2 * pi * 3 * t / periodo),
cos3 = cos(2 * pi * 3 * t / periodo)
)
terminos_fourier <- c("sin1", "cos1", "sin2", "cos2", "sin3", "cos3")
dengue |>
select(semana, casos, all_of(terminos_fourier)) |>
head() semana casos sin1 cos1 sin2 cos2 sin3 cos3
1 1990-01-01 61 0.1205367 0.9927089 0.2393157 0.9709418 0.3546049 0.9350162
2 1990-01-08 50 0.2393157 0.9709418 0.4647232 0.8854560 0.6631227 0.7485107
3 1990-01-15 44 0.3546049 0.9350162 0.6631227 0.7485107 0.8854560 0.4647232
4 1990-01-22 46 0.4647232 0.8854560 0.8229839 0.5680647 0.9927089 0.1205367
5 1990-01-29 39 0.5680647 0.8229839 0.9350162 0.3546049 0.9709418 -0.2393157
6 1990-02-05 34 0.6631227 0.7485107 0.9927089 0.1205367 0.8229839 -0.5680647
TipQué son esas columnas
sin1 y cos1 completan una vuelta al año: capturan la temporada alta y la baja. sin2 y cos2 dan dos vueltas al año y sin3 y cos3 tres, así que van agregando detalle a la forma de la curva estacional.
Entre más pares (más K), más flexible la estacionalidad; con K demasiado grande empiezas a ajustar ruido. Es exactamente el mismo compromiso que ya viste con el ancho de banda de un histograma.
Lo importante: son columnas comunes y corrientes. No hace falta convertir nada a un formato especial de series de tiempo.
14.4 Partir los datos
Para poder evaluar honestamente (paso 7 del flujo de trabajo), apartamos el último año y no lo dejamos ver a ningún modelo:
H <- 52 # vamos a pronosticar un año
entrena <- head(dengue, nrow(dengue) - H)
prueba <- tail(dengue, H)
# el tiempo como número, para la regresión lineal
entrena$t <- seq_len(nrow(entrena))
prueba$t <- nrow(entrena) + seq_len(H)
c(entrena = nrow(entrena), prueba = nrow(prueba))entrena prueba
1033 52 14.5 Primer intento: una regresión lineal
Empecemos con lo más simple que conocemos (paso 2: el modelo más simple que podría servir). Una recta contra el tiempo:
\[ \text{casos}_t = \beta_0 + \beta_1 \cdot t + \varepsilon_t \]
modelo_lm <- lm(casos ~ t, data = entrena)
tidy(modelo_lm)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 61.6 2.75 22.4 1.29e-90
2 t -0.0343 0.00462 -7.44 2.12e-13entrena %>%
mutate(ajustado = fitted(modelo_lm)) %>%
ggplot() +
geom_line(aes(x = semana, y = casos), color = "gray60", linewidth = 0.3) +
geom_line(aes(x = semana, y = ajustado), color = "#ff6361", linewidth = 1) +
labs(x = "", y = "Casos", title = "La recta no capta absolutamente nada") +
theme_bw()
La recta es una línea casi plana en medio de la nube. Pero el problema de fondo no es que se vea fea: son dos problemas más serios y menos obvios.
14.5.1 Problema 1: la dependencia
La regresión lineal supone que los errores son independientes entre sí: que saber cuánto se equivocó el modelo esta semana no te dice nada sobre cuánto se va a equivocar la siguiente.
En una serie de tiempo eso es falso y es fácil de comprobar. Veamos el residual de cada semana contra el de la semana anterior:
res_lm <- residuals(modelo_lm)
tibble(anterior = head(res_lm, -1), actual = tail(res_lm, -1)) %>%
ggplot(aes(x = anterior, y = actual)) +
geom_point(color = "#003f5c", alpha = 0.4, size = 0.8) +
geom_smooth(method = "lm", formula = "y ~ x", color = "#ff6361") +
labs(
x = "Residual de la semana anterior", y = "Residual de esta semana",
title = "Los errores NO son independientes"
) +
theme_bw()
Es casi una recta perfecta. Podemos medirlo con la autocorrelación:
acf(res_lm, main = "ACF de los residuales de la regresión lineal")
acf_lm <- acf(res_lm, plot = FALSE, lag.max = 3)$acf[2:4]
round(acf_lm, 3)[1] 0.946 0.897 0.845
ImportanteÉsta es la idea central de las series de tiempo
La autocorrelación en el primer rezago es de 0.95. Es decir: si esta semana el modelo se quedó corto por 100 casos, la próxima semana se va a quedar corto por unos 95 casos.
Eso es información valiosísima que la regresión lineal está tirando a la basura. Y no sólo la desperdicia: los errores estándar, los valores p y los intervalos que reporta están todos calculados suponiendo independencia, así que están mal.
Cuando los datos vienen ordenados en el tiempo, lo que pasó ayer te dice algo de lo que va a pasar hoy. Un modelo que no usa eso, está incompleto.
14.5.2 Problema 2: la varianza no crece
El segundo problema es más sutil y, para pronosticar, más grave. Pidámosle a la regresión que pronostique el año que apartamos:
pred_lm <- predict(modelo_lm, newdata = prueba, se.fit = TRUE)
# desviación de la predicción (incluye el error del modelo y el ruido)
sd_lm <- sqrt(pred_lm$se.fit^2 + summary(modelo_lm)$sigma^2)
# ancho del intervalo al 95%, en la primera y la última semana pronosticada
ancho_lm <- 2 * qt(0.975, df = modelo_lm$df.residual) * sd_lm
round(c(semana_1 = ancho_lm[1], semana_52 = ancho_lm[H]), 1) semana_1.1034 semana_52.1085
174 174
AdvertenciaEl intervalo mide lo mismo dentro de una semana que dentro de un año
Ése es el problema. La regresión lineal dice estar igual de segura de lo que va a pasar la semana que entra que de lo que va a pasar dentro de doce meses.
Piénsalo con sentido común: ¿tú confiarías igual en un pronóstico de dengue a una semana que a un año? Claro que no. Pero el modelo sí, porque en una regresión lineal la incertidumbre depende de qué tan lejos estás del promedio de las x, no de qué tan lejos estás en el futuro.
Para una regresión normal esto está bien: predecir el peso de una persona de 1.90 m no es “más al futuro” que predecir el de una de 1.70 m. Pero en una serie de tiempo, el futuro sí se vuelve más incierto conforme te alejas, y el modelo tiene que saberlo.
14.6 El ARIMA: modelar la dependencia
Un ARIMA ataca justo el problema 1: en vez de suponer que los errores son independientes, modela explícitamente cómo cada observación depende de las anteriores. Su nombre son sus tres piezas:
- AR (autorregresivo): el valor de hoy depende de los valores pasados.
- I (integrado): en lugar de modelar la serie, se modelan sus cambios (la diferencia entre una semana y la anterior). Es lo que le quita la tendencia.
- MA (medias móviles): el valor de hoy depende de los errores pasados.
Para la estacionalidad usaremos los términos de Fourier que ya agregamos a la base. Es la forma recomendada cuando hay muchas estaciones (52 semanas son demasiadas para un ARIMA estacional clásico).
Y ajustamos. auto.arima acepta el vector de casos directamente, y los términos de Fourier van en xreg como una matriz:
modelo_arima <- auto.arima(
entrena$casos,
xreg = as.matrix(entrena[, terminos_fourier]),
seasonal = FALSE
)
modelo_arimaSeries: entrena$casos
Regression with ARIMA(0,1,1) errors
Coefficients:
ma1 sin1 cos1 sin2 cos2 sin3 cos3
-0.0958 -29.4416 14.6667 -1.2218 -5.5083 -0.5122 -2.1219
s.e. 0.0316 4.6552 4.6313 2.3323 2.3312 1.5655 1.5660
sigma^2 = 199.4: log likelihood = -4193.13
AIC=8402.27 AICc=8402.41 BIC=8441.78¿Se arregló la dependencia? Veamos los residuales:
acf(residuals(modelo_arima), main = "ACF de los residuales del ARIMA")
acf(residuals(modelo_arima), plot = FALSE, lag.max = 3)$acf[2:4][1] 0.001773791 -0.019940648 0.013940173
TipEsto es lo que significa que un modelo “esté bien”
La autocorrelación pasó de 0.95 a prácticamente cero. Los residuales ya no tienen estructura: son ruido.
Eso es exactamente lo que uno quiere. Si sobra estructura en los residuales, significa que el modelo dejó información sin usar. Cuando los residuales parecen ruido, ya le exprimiste a los datos lo que se podía.
Y ahora el pronóstico:
pron_arima <- forecast(
modelo_arima,
xreg = as.matrix(prueba[, terminos_fourier]),
level = c(50, 80, 90, 95)
)
# ancho del intervalo al 95% en la primera y la última semana
ancho_arima <- c(
pron_arima$upper[1, 4] - pron_arima$lower[1, 4],
pron_arima$upper[H, 4] - pron_arima$lower[H, 4]
)
round(ancho_arima, 1) 95% 95%
55.3 361.7 # El pronóstico también se guarda como data.frame, con sus fechas reales
pronostico_df <- prueba |>
select(semana, observado = casos) |>
mutate(
ajuste = as.numeric(pron_arima$mean),
bajo = as.numeric(pron_arima$lower[, 4]),
alto = as.numeric(pron_arima$upper[, 4])
)
ggplot(pronostico_df, aes(x = semana)) +
geom_ribbon(aes(ymin = bajo, ymax = alto), fill = "#003f5c", alpha = 0.2) +
geom_line(aes(y = ajuste), color = "#003f5c", linewidth = 1) +
geom_line(aes(y = observado), color = "#ff6361") +
labs(
x = "", y = "Casos",
title = "Pronóstico ARIMA a un año",
subtitle = "En rojo lo observado; en azul el pronóstico con su intervalo al 95%"
) +
theme_bw()
ImportanteEl ARIMA sí resuelve los dos problemas
- Dependencia: la autocorrelación de los residuales se fue a cero.
- Varianza: el intervalo pasa de medir ~55 casos en la primera semana a ~362 en la semana 52. El modelo admite que sabe menos conforme se aleja, que es justo lo que la regresión lineal no hacía.
Ese abanico que se abre es la firma visual de un modelo de series de tiempo bien planteado.
14.7 El bsts: estructura y honestidad bayesiana
El bsts (Bayesian Structural Time Series) parte de otra idea: en vez de modelar la dependencia con rezagos, arma la serie por pedazos que tienen significado, y deja que cada pedazo evolucione con el tiempo:
\[ \text{casos}_t = \underbrace{\text{nivel}_t}_{\text{dónde está}} + \underbrace{\text{estacionalidad}_t}_{\text{en qué época del año}} + \text{ruido}_t \]
La diferencia clave con el ARIMA es que aquí el nivel y la estacionalidad cambian de una semana a otra. Justo lo que vimos en el escalograma: que el ciclo anual no es igual de fuerte en todas las épocas.
# La semilla va ANTES de construir la estructura: esos comandos también usan
# números aleatorios, así que si la pones después el resultado cambia entre
# corridas. El argumento seed = de bsts fija además el muestreo interno.
set.seed(97531)
estructura <- AddLocalLevel(list(), entrena$casos)
estructura <- AddSeasonal(estructura, entrena$casos, nseasons = 52)
modelo_bsts <- bsts(entrena$casos,
state.specification = estructura,
niter = 2000,
seed = 97531,
ping = 0
)
Nota
El bsts es bayesiano: en vez de un solo resultado, produce mil versiones posibles del modelo (eso es el niter = 1000). Las primeras iteraciones se descartan porque el algoritmo todavía se está acomodando; a eso se le llama burn-in.
Por eso tarda más que el ARIMA: alrededor de un minuto y medio contra menos de un segundo.
Podemos ver los componentes que estimó:
plot(modelo_bsts, "components", burn = 200)
Y el pronóstico:
pron_bsts <- predict(modelo_bsts, horizon = H, burn = 200)
# ancho del intervalo al 95%
q_bsts <- t(apply(pron_bsts$distribution, 2, quantile,
probs = c(0.025, 0.975)
))
ancho_bsts <- c(q_bsts[1, 2] - q_bsts[1, 1], q_bsts[H, 2] - q_bsts[H, 1])
round(ancho_bsts, 1)97.5% 97.5%
55.7 362.2
TipDos caminos distintos, conclusión parecida
Su intervalo pasa de ~56 casos a ~362 en un año: prácticamente lo mismo que el ARIMA (~55 y ~362).
Eso es tranquilizador, no aburrido. Son dos modelos con filosofías muy distintas —uno modela rezagos, el otro componentes que evolucionan— y llegan por separado a la misma conclusión sobre cuánta incertidumbre hay. Cuando dos métodos independientes coinciden, uno confía más en los dos.
Lo que sigue es ver si esa incertidumbre estaba bien calibrada.
AdvertenciaUna advertencia sobre el
bsts que te va a ahorrar horas
Aquí usamos AddLocalLevel, que supone que el nivel de la serie se mueve sin rumbo fijo. Existe también AddLocalLinearTrend, que además le estima una pendiente que también se mueve.
Suena mejor, pero para pronosticar a un año es una trampa: el modelo extrapola esa pendiente 52 semanas hacia adelante y los intervalos explotan. Probándolo con dos semillas distintas nos dio coberturas del 100% y del 81% para el mismo intervalo al 95%: o sea, el resultado dependía del azar y no de los datos.
Si tu modelo bayesiano cambia de conclusión al cambiar la semilla, el problema no es la semilla: es que el modelo no está bien identificado con los datos que tienes. Cámbialo por uno más simple.
14.8 ¿Cuál predice mejor? El WIS
Ya tenemos tres pronósticos del mismo año. Necesitamos una forma justa de compararlos.
AdvertenciaNota sobre el nombre
La medida se llama WIS, de Weighted Interval Score (puntaje de intervalo ponderado). A veces se le dice “WIC” por confusión con los criterios de información tipo AIC, pero no tiene ninguna relación con ellos: el AIC compara qué tan bien un modelo ajusta los datos que ya vio, mientras que el WIS evalúa qué tan buenos son los pronósticos sobre datos que no vio.
14.8.1 Qué es el WIS
Comparar pronósticos con el error promedio (por ejemplo, la diferencia entre lo predicho y lo observado) tiene un problema grave: ignora la incertidumbre. Un modelo que dice “van a ser 100 casos, seguro” y un modelo que dice “van a ser entre 20 y 300” pueden tener el mismo error en el centro, pero no son igual de buenos.
El WIS evalúa el intervalo completo. Para un intervalo dado, cobra tres cosas:
Qué tan ancho es el intervalo. Los intervalos anchos son poco útiles, así que pagan una multa por su ancho. Esto evita el truco de decir “entre 0 y un millón” para nunca fallar.
Qué tanto se quedó corto: si lo observado cayó arriba del intervalo, se cobra la distancia.
Qué tanto se pasó: si lo observado cayó abajo del intervalo, también se cobra la distancia.
Esto se calcula para varios intervalos a la vez (50%, 80%, 90%, 95%…) y se promedia. Entre más chico el WIS, mejor, igual que un error.
TipLo bonito del
WIS: se descompone
El WIS total se separa exactamente en esas tres partes, y esa descomposición te dice por qué un modelo es malo:
- Mucha dispersión (
dispersion) → tus intervalos son demasiado anchos: el modelo no se compromete. - Mucha subestimación (
underprediction) → los casos reales salieron por arriba de tus intervalos: te quedaste corto. - Mucha sobreestimación (
overprediction) → predijiste más de lo que pasó.
Un modelo puede tener buen WIS por ser preciso, o mal WIS por dos razones opuestas. La descomposición te dice cuál te tocó.
14.8.2 Preparar los pronósticos
scoringutils necesita los pronósticos en formato largo: un renglón por modelo, por semana y por cuantil. Un cuantil de 0.025 es el piso del intervalo al 95%, 0.5 es la mediana, y así.
niveles <- c(0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975)De la regresión lineal los sacamos con la distribución t:
q_lm <- sapply(niveles, function(q) {
pred_lm$fit + qt(q, df = modelo_lm$df.residual) * sd_lm
})Del ARIMA, los intervalos que ya calculamos:
q_arima <- cbind(
pron_arima$lower[, 4], pron_arima$lower[, 3],
pron_arima$lower[, 2], pron_arima$lower[, 1],
pron_arima$mean,
pron_arima$upper[, 1], pron_arima$upper[, 2],
pron_arima$upper[, 3], pron_arima$upper[, 4]
)Y del bsts, directo de sus mil simulaciones:
q_bsts_todos <- t(apply(pron_bsts$distribution, 2, quantile, probs = niveles))Los juntamos:
armar <- function(cuantiles, nombre) {
tibble(
modelo = nombre,
semana = rep(prueba$semana, times = length(niveles)),
observado = rep(prueba$casos, times = length(niveles)),
nivel = rep(niveles, each = H),
prediccion = as.vector(cuantiles)
)
}
pronosticos <- bind_rows(
armar(q_lm, "Regresion lineal"),
armar(q_arima, "ARIMA"),
armar(q_bsts_todos, "bsts")
)
head(pronosticos)# A tibble: 6 × 5
modelo semana observado nivel prediccion
<chr> <date> <int> <dbl> <dbl>
1 Regresion lineal 2009-10-26 68 0.025 -60.9
2 Regresion lineal 2009-11-02 80 0.025 -60.9
3 Regresion lineal 2009-11-09 89 0.025 -61.0
4 Regresion lineal 2009-11-16 85 0.025 -61.0
5 Regresion lineal 2009-11-23 110 0.025 -61.0
6 Regresion lineal 2009-11-30 98 0.025 -61.114.8.3 Calcular el WIS
evaluacion <- pronosticos %>%
as_forecast_quantile(
forecast_unit = c("modelo", "semana"),
observed = "observado",
predicted = "prediccion",
quantile_level = "nivel"
) %>%
score() %>%
summarise_scores(by = "modelo")
evaluacion %>%
select(modelo, wis, dispersion, underprediction, overprediction) %>%
arrange(wis) modelo wis dispersion underprediction overprediction
<char> <num> <num> <num> <num>
1: bsts 51.73490 12.152973 39.57070 0.01122503
2: ARIMA 54.06773 11.901665 42.15374 0.01232778
3: Regresion lineal 79.64174 8.438959 71.20278 0.00000000
ImportanteCómo se lee esta tabla
El bsts gana (el WIS más chico), seguido de cerca por el ARIMA, y la regresión lineal queda muy atrás. Fíjate en las distancias: entre bsts y ARIMA hay unos pocos puntos, mientras que de ahí a la regresión lineal hay un abismo.
La lección no es “usa bsts”, es “usa un modelo de series de tiempo”. Los dos que sí modelan la dependencia terminan parecidos; el que no la modela pierde por mucho, sin importar cuál de los otros dos elijas.
Y la descomposición dice por qué pierde:
La regresión lineal tiene la dispersión más baja (los intervalos más angostos de los tres) y a la vez la subestimación más alta. Traducido: estuvo muy segura y muy equivocada. Es el peor escenario posible en salud pública, porque un pronóstico angosto se lee como confiable.
Los otros dos “pagan” un poco más de dispersión —sus intervalos son más anchos— y a cambio se quedan mucho menos cortos. Ése es el intercambio que el
WISmide: vale la pena admitir incertidumbre si eso te evita fallar feo.
Nota que la sobreestimación es casi cero en los tres: ninguno se pasó. Todos se quedaron cortos, unos más que otros.
14.8.4 La prueba definitiva: la cobertura
Hay una forma todavía más directa de ver quién decía la verdad. Un intervalo al 50% debería contener el valor real el 50% de las veces. Uno al 90%, el 90%. Veamos:
evaluacion %>%
select(modelo, interval_coverage_50, interval_coverage_90) %>%
arrange(desc(interval_coverage_90)) modelo interval_coverage_50 interval_coverage_90
<char> <num> <num>
1: bsts 0.4038462 0.7115385
2: ARIMA 0.2692308 0.6730769
3: Regresion lineal 0.1730769 0.5769231
AdvertenciaEl orden se mantiene… pero los tres se quedan cortos
La regresión lineal es otra vez la peor: sus intervalos al 50% atinaron sólo el 17% de las veces y los del 90% atinaron el 58%. Cuando decía “estoy 90% seguro”, se equivocaba cuatro de cada diez veces.
Pero mira los otros dos con cuidado, porque aquí viene lo incómodo: el bsts, que es el mejor, apenas llega a 40% y 71%. También se queda corto de lo que promete.
Ninguno de los tres modelos está bien calibrado. Y eso no es un error de programación: es que el año que apartamos tuvo una epidemia más grande que casi todo lo que había en los veinte años previos. Ningún modelo que sólo mire el pasado de la serie podía anticiparla, porque la información no estaba en la serie: estaba en el clima, en el serotipo circulante, en la inmunidad de la población.
ImportanteLa conclusión honesta de este capítulo
Es tentador terminar con “y entonces usa bsts”. La conclusión de verdad es más útil y menos cómoda:
- Modelar la dependencia sí sirve: bajó el
WISde 80 a ~52 y subió la cobertura del 58% al 71%. - Aun así, el mejor modelo sigue estando mal calibrado. Si reportaras estos intervalos como si fueran del 90%, estarías prometiendo de más.
- La forma de saberlo no fue mirar las gráficas, que se ven bonitas en los tres casos. Fue medir la cobertura contra datos que los modelos no vieron.
Un intervalo que no cubre lo que dice cubrir no es conservador ni cauteloso: es incorrecto. Y si alguien planea camas de hospital con él, el error se paga en la vida real.
Ése es el paso 10 del flujo de trabajo: reportar el resultado junto con sus límites. Aquí el límite es que estos modelos sirven para semanas, no para anticipar una epidemia a un año.
14.8.5 Verlo todo junto
bandas <- pronosticos %>%
filter(nivel %in% c(0.025, 0.5, 0.975)) %>%
mutate(cual = case_when(
nivel == 0.025 ~ "bajo",
nivel == 0.5 ~ "medio",
TRUE ~ "alto"
)) %>%
select(-nivel) %>%
pivot_wider(names_from = cual, values_from = prediccion)
ggplot(bandas) +
geom_ribbon(aes(x = semana, ymin = bajo, ymax = alto),
fill = "#003f5c", alpha = 0.25
) +
geom_line(aes(x = semana, y = medio), color = "#003f5c") +
geom_line(aes(x = semana, y = observado), color = "#ff6361", linewidth = 0.8) +
facet_wrap(~modelo, ncol = 1, scales = "free_y") +
labs(
x = "", y = "Casos",
title = "Pronóstico a un año contra lo observado",
subtitle = "En rojo lo que pasó; en azul el pronóstico con su intervalo al 95%"
) +
theme_bw()
14.9 Ejercicios
14.9.1 A mano
Antes de correr nada: si la autocorrelación de los residuales en el primer rezago es
0.95, y esta semana el modelo se quedó corto por 200 casos, ¿por cuánto esperas que se quede corto la semana que entra? ¿Y dentro de dos semanas, si el efecto se multiplica?Dibuja en tu cuaderno cómo se vería el intervalo de pronóstico (el abanico) de los tres modelos, uno junto a otro, a 52 semanas. Después compáralo con la última gráfica del capítulo.
Un modelo dice “entre 0 y 1,000,000 de casos” todas las semanas. Nunca falla: su cobertura al 95% sería del 100%. ¿Por qué el
WISsí lo castiga? ¿Cuál de los tres componentes se le dispararía?
14.9.2 Encuentra el error
# a) Queremos el ACF de los residuales del modelo lineal
acf(modelo_lm)# b) Queremos pronosticar 52 semanas con el ARIMA
forecast(modelo_arima, h = 52)# c) Queremos evaluar qué tan bueno fue el ajuste del ARIMA
pron <- forecast(modelo_arima, h = 52)
mean(abs(pron$mean - entrena$casos))
NotaRespuestas
acfnecesita un vector de números, no el objeto del modelo. Correcto:acf(residuals(modelo_lm)).El modelo se ajustó con variables externas (
xreg, los términos de Fourier), así que para pronosticar hay que darle los valores futuros de esas variables. Correcto:forecast(modelo_arima, xreg = as.matrix(prueba[, terminos_fourier])). Nota que entonces ya no hace falta lah: el horizonte lo define cuántos renglones tenga elxreg.Dos errores. Uno técnico: compara el pronóstico contra
entrena, que son los datos con los que se construyó el modelo, en vez de contraprueba. Y uno conceptual: evaluar un pronóstico contra los datos de entrenamiento no es evaluarlo. Correcto:mean(abs(pron$mean - prueba$casos)).
14.9.3 Con apoyo de una IA
Pídele a una IA que te explique qué es un
ARIMA(0,1,1)“como si nunca hubiera visto una serie de tiempo”. Contrasta su explicación con las tres piezas (AR, I, MA) de este capítulo. ¿Explicó bien qué significa el1de en medio?Pásale la tabla de descomposición del
WISy pregúntale cuál modelo es mejor y por qué. Verifica si menciona la cobertura por su cuenta. Si no lo hace, es una omisión importante: pregúntale luego por qué la cobertura importa.Pídele que te proponga un cuarto modelo para esta serie. Antes de correr nada, revisa si su propuesta puede producir casos negativos (como nos pasó con el
ARIMAen el capítulo de regresiones). Ése es el paso 3 del flujo de trabajo y casi ninguna IA lo hace sola.
14.9.4 Para practicar
Ajusta el
ARIMAconK = 1y conK = 6términos de Fourier en lugar deK = 3. ¿Cambia elWIS? ¿Qué representa cadaKen términos de la forma de la curva estacional?Agrega un cuarto modelo a la comparación: un
ARIMAsin los términos de Fourier (es decir, sin estacionalidad). Calcula suWISy su cobertura. ¿Cuánto aportaba la estacionalidad?Cambia el horizonte de
H = 52aH = 4semanas y vuelve a evaluar los tres modelos. ¿Se mantiene el mismo orden? Piensa por qué la regresión lineal es menos mala en horizontes cortos.
14.10 Sistema
sessioninfo::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.5.3 (2026-03-11)
os macOS Sequoia 15.7.3
system x86_64, darwin20
ui X11
language (EN)
collate es_ES.UTF-8
ctype es_ES.UTF-8
tz America/Mexico_City
date 2026-08-07
pandoc 3.10.1 @ /usr/local/bin/ (via rmarkdown)
quarto 1.10.18 @ /usr/local/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
backports 1.5.1 2026-04-03 [1] CRAN (R 4.5.2)
bit 4.6.0 2025-03-06 [1] CRAN (R 4.5.0)
bit64 4.8.2 2026-05-19 [1] CRAN (R 4.5.2)
Boom * 0.9.16 2025-09-02 [1] CRAN (R 4.5.1)
BoomSpikeSlab * 1.2.7 2025-09-04 [1] CRAN (R 4.5.1)
broom * 1.0.12 2026-01-27 [1] CRAN (R 4.5.1)
bsts * 0.9.11 2025-09-05 [1] CRAN (R 4.5.1)
checkmate 2.3.4 2026-02-03 [1] CRAN (R 4.5.1)
cli 3.6.6 2026-04-09 [1] CRAN (R 4.5.2)
coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.5.0)
colorspace 2.1-2 2025-09-22 [1] CRAN (R 4.5.1)
crayon 1.5.3 2024-06-20 [1] CRAN (R 4.5.0)
data.table 1.18.4 2026-05-06 [1] CRAN (R 4.5.2)
digest 0.6.39 2025-11-19 [1] CRAN (R 4.5.1)
dplyr * 1.2.1 2026-04-03 [1] CRAN (R 4.5.2)
evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.5.1)
farver 2.1.2 2024-05-13 [1] CRAN (R 4.5.0)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.5.0)
forcats * 1.0.1 2025-09-25 [1] CRAN (R 4.5.1)
forecast * 9.0.2 2026-03-18 [1] CRAN (R 4.5.2)
fracdiff 1.5-3 2024-02-01 [1] CRAN (R 4.5.0)
generics 0.1.4 2025-05-09 [1] CRAN (R 4.5.0)
ggplot2 * 4.0.3 2026-04-22 [1] CRAN (R 4.5.2)
glue 1.8.1 2026-04-17 [1] CRAN (R 4.5.2)
gtable 0.3.6 2024-10-25 [1] CRAN (R 4.5.0)
hms 1.1.4 2025-10-17 [1] CRAN (R 4.5.1)
htmltools 0.5.9 2025-12-04 [1] CRAN (R 4.5.1)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.5.0)
jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.5.0)
knitr 1.51 2025-12-20 [1] CRAN (R 4.5.1)
labeling 0.4.3 2023-08-29 [1] CRAN (R 4.5.0)
lattice 0.22-9 2026-02-09 [1] CRAN (R 4.5.3)
lifecycle 1.0.5 2026-01-08 [1] CRAN (R 4.5.1)
lubridate * 1.9.5 2026-02-04 [1] CRAN (R 4.5.1)
magrittr 2.0.5 2026-04-04 [1] CRAN (R 4.5.2)
MASS 7.3-65 2025-02-28 [1] CRAN (R 4.5.3)
Matrix 1.7-5 2026-03-21 [1] CRAN (R 4.5.2)
mgcv 1.9-4 2025-11-07 [1] CRAN (R 4.5.3)
nlme 3.1-168 2025-03-31 [1] CRAN (R 4.5.3)
NobBS * 1.1.1 2026-06-01 [1] CRAN (R 4.5.2)
otel 0.2.0 2025-08-29 [1] CRAN (R 4.5.1)
pillar 1.11.1 2025-09-17 [1] CRAN (R 4.5.1)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.5.0)
purrr * 1.2.2 2026-04-10 [1] CRAN (R 4.5.2)
R6 2.6.1 2025-02-15 [1] CRAN (R 4.5.0)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.5.0)
Rcpp 1.1.2 2026-07-05 [1] CRAN (R 4.5.2)
readr * 2.2.0 2026-02-19 [1] CRAN (R 4.5.2)
rjags 4-17 2025-03-24 [1] CRAN (R 4.5.0)
rlang 1.3.0 2026-07-05 [1] CRAN (R 4.5.2)
rmarkdown 2.31 2026-03-26 [1] CRAN (R 4.5.2)
S7 0.2.2 2026-04-22 [1] CRAN (R 4.5.2)
scales 1.4.0 2025-04-24 [1] CRAN (R 4.5.0)
scoringRules 1.1.3 2024-09-18 [1] CRAN (R 4.5.0)
scoringutils * 2.2.0 2026-04-05 [1] CRAN (R 4.5.2)
sessioninfo 1.2.3 2025-02-05 [1] CRAN (R 4.5.0)
stringi 1.8.7 2025-03-27 [1] CRAN (R 4.5.0)
stringr * 1.6.0 2025-11-04 [1] CRAN (R 4.5.1)
tibble * 3.3.1 2026-01-11 [1] CRAN (R 4.5.1)
tidyr * 1.3.2 2025-12-19 [1] CRAN (R 4.5.1)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.5.0)
tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.5.0)
timechange 0.4.0 2026-01-29 [1] CRAN (R 4.5.1)
timeDate 4052.112 2026-01-28 [1] CRAN (R 4.5.1)
tzdb 0.5.0 2025-03-15 [1] CRAN (R 4.5.0)
urca 1.3-4 2024-05-27 [1] CRAN (R 4.5.0)
utf8 1.2.6 2025-06-08 [1] CRAN (R 4.5.0)
vctrs 0.7.3 2026-04-11 [1] CRAN (R 4.5.2)
vroom 1.7.1 2026-03-31 [1] CRAN (R 4.5.2)
WaveletComp * 1.2 2025-04-29 [1] CRAN (R 4.5.0)
withr 3.0.3 2026-06-19 [1] CRAN (R 4.5.2)
xfun 0.60 2026-07-09 [1] CRAN (R 4.5.2)
xts * 0.14.2 2026-02-28 [1] CRAN (R 4.5.2)
yaml 2.3.12 2025-12-10 [1] CRAN (R 4.5.1)
zoo * 1.8-15 2025-12-15 [1] CRAN (R 4.5.1)
[1] /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────