Objetivo de la práctica: Aprender a analizar series temporales biomédicas de principio a fin: desde la gestión de fechas y la exploración visual, hasta el modelado ARIMA y la predicción. Trabajaremos con datos reales de ingresos hospitalarios por COVID-19 y mortalidad por cáncer de mama.

1 Instalación y carga de paquetes

Antes de comenzar cualquier análisis en R, debemos cargar los paquetes necesarios. Cada uno tiene un rol específico:

Paquete	Función principal
`lubridate`	Manejo y manipulación de fechas
`tseries`	Tests estadísticos de series temporales (ADF)
`forecast`	Modelado ARIMA, predicción y diagnóstico
`ggplot2`	Visualización estática de alta calidad
`zoo`	Series temporales con fechas reales como índice
`xts`	Extensión de zoo, especialmente para datos financieros y temporales
`dygraphs`	Gráficos interactivos (solo en HTML)

# Descomenta la siguiente línea si es la primera vez que ejecutas esto:
# install.packages(c("lubridate", "tseries", "forecast", "ggplot2", "zoo", "xts", "dygraphs"))

library(lubridate)
library(tseries)
library(forecast)
library(ggplot2)
library(zoo)
library(xts)
library(dygraphs)

2 Bloque 1: Manejo de fechas con `lubridate`

2.1 ¿Por qué es tan importante el manejo de fechas?

En una serie temporal, el tiempo es el eje fundamental. Sin fechas correctamente formateadas, R no puede ordenar los datos cronológicamente ni detectar patrones temporales. El paquete lubridate facilita enormemente este proceso.

2.2 Funciones principales de `lubridate`

Las tres funciones más usadas leen fechas en distintos formatos según el orden de los componentes:

ymd() → formato año-mes-día (ISO internacional, ej: "2020-03-14")
dmy() → formato día-mes-año (europeo, ej: "14-03-2020")
mdy() → formato mes-día-año (americano, ej: "03-14-2020")

# Creamos la fecha del inicio del estado de alarma por COVID-19 en España
inicio_pandemia <- ymd("2020-03-14")
inicio_pandemia

#> [1] "2020-03-14"

class(inicio_pandemia)  # R la reconoce como objeto de tipo "Date"

#> [1] "Date"

¿Qué significa class(inicio_pandemia) = "Date"? R guarda internamente las fechas como el número de días transcurridos desde el 1 de enero de 1970. Al convertirlas a clase Date, podemos hacer aritmética temporal directamente (sumar días, calcular diferencias, etc.).

# Podemos extraer cualquier componente de una fecha
year(inicio_pandemia)                              # Año

#> [1] 2020

month(inicio_pandemia)                             # Mes (número)

#> [1] 3

day(inicio_pandemia)                               # Día del mes

#> [1] 14

wday(inicio_pandemia, label = TRUE, abbr = FALSE)  # Día de la semana

#> [1] Saturday
#> 7 Levels: Sunday < Monday < Tuesday < Wednesday < Thursday < ... < Saturday

Interpretación: El 14 de marzo de 2020 fue sábado, el día que el Gobierno de España declaró el estado de alarma por COVID-19.

# Calculamos cuántos días duró el estado de alarma
fin_estado_alarma <- ymd("2020-06-21")
dias_alarma       <- as.numeric(fin_estado_alarma - inicio_pandemia)

cat("El estado de alarma duró", dias_alarma, "días\n")

#> El estado de alarma duró 99 días

cat("Inicio:", format(inicio_pandemia, "%d/%m/%Y"), "\n")

#> Inicio: 14/03/2020

cat("Fin:   ", format(fin_estado_alarma,  "%d/%m/%Y"), "\n")

#> Fin:    21/06/2020

Nota: La resta entre dos objetos Date devuelve un objeto de clase difftime. Lo convertimos a número con as.numeric() para obtener directamente los días.

También podemos incluir hora y zona horaria cuando necesitamos precisión horaria:

# Ejemplo: primer caso reportado a la OMS (31 diciembre 2019, Wuhan)
primer_caso_oms <- ymd_hm("2019-12-31 12:00", tz = "Asia/Shanghai")
primer_caso_oms

#> [1] "2019-12-31 12:00:00 CST"

3 Bloque 2: Importar y preparar la serie temporal

3.1 Conceptos clave: objetos `zoo` y `xts`

El objeto ts base de R es útil, pero tiene una limitación importante: solo trabaja con frecuencias regulares y no guarda la fecha real (solo posición). Los objetos zoo y xts resuelven esto al asociar cada valor a su fecha exacta, respetando el orden cronológico y permitiendo huecos o fechas irregulares.

3.2 Carga y exploración del dataset de ingresos COVID

# Cargamos el dataset de ingresos hospitalarios
admissions <- read.csv("admissions.csv")

# Vista de las primeras filas
head(admissions)

# Estructura del dataset
str(admissions)

#> 'data.frame':    9263 obs. of  5 variables:
#>  $ HospitalAdmitDate  : chr  "2020-03-07" "2020-03-10" "2020-03-11" "2020-03-12" ...
#>  $ HospitalID         : chr  "C" "B" "B" "B" ...
#>  $ ZipCD_short        : chr  "02537" "02111" "02114" "02356" ...
#>  $ COVID_Encounter_Ref: chr  "covid_second" "covid_prim" "covid_prim" "covid_second" ...
#>  $ NumberOfAdmissions : int  1 1 1 1 1 1 1 1 1 1 ...

# Dimensiones
cat("Filas:", nrow(admissions), " | Columnas:", ncol(admissions), "\n")

#> Filas: 9263  | Columnas: 5

# Hospitales disponibles en el dataset
cat("Hospitales:", paste(unique(admissions$HospitalID), collapse = ", "), "\n")

#> Hospitales: C, B, E, A, H, I, D, F, J, G

Descripción del dataset: - HospitalAdmitDate: fecha de ingreso - HospitalID: identificador del hospital (A, B o C) - ZipCD_short: código postal del paciente - COVID_Encounter_Ref: tipo de encuentro COVID (primario o secundario) - NumberOfAdmissions: número de ingresos ese día en ese código postal

3.3 Filtrado y agregación del Hospital A

Un hospital puede registrar varios códigos postales el mismo día. Para construir la serie temporal necesitamos sumar todos los ingresos por fecha:

# Filtramos solo el Hospital A
admissions_A <- admissions[admissions$HospitalID == "A", ]

# Agregamos los ingresos por fecha (sumamos todos los códigos postales)
data_A <- aggregate(admissions_A$NumberOfAdmissions,
                    by  = list(admissions_A$HospitalAdmitDate),
                    FUN = sum)
names(data_A) <- c("Dates", "Admissions")

# Convertimos la columna de fechas al tipo Date de R
data_A$Dates <- as.Date(data_A$Dates, format = "%Y-%m-%d")

head(data_A, 10)

3.4 Creación del objeto `zoo`

# Creamos la serie temporal con fechas reales como índice
ts_A <- zoo(x = data_A$Admissions, order.by = data_A$Dates)

# Información básica de la serie
cat("Inicio de la serie: ", format(start(ts_A), "%d/%m/%Y"), "\n")

#> Inicio de la serie:  15/03/2020

cat("Fin de la serie:    ", format(end(ts_A),   "%d/%m/%Y"), "\n")

#> Fin de la serie:     11/04/2021

cat("Total observaciones:", length(ts_A), "días\n")

#> Total observaciones: 336 días

Interpretación: La serie del Hospital A abarca todo el período principal de la pandemia COVID-19 en España, con una observación por cada día.

4 Bloque 3: Exploración visual de la serie

4.1 Por qué siempre visualizamos primero

Antes de aplicar cualquier modelo estadístico, la visualización nos permite detectar:

Tendencia: ¿Los valores suben, bajan o se mantienen estables a lo largo del tiempo?
Estacionalidad: ¿Existen patrones que se repiten con cierta periodicidad (semanal, mensual)?
Valores atípicos (outliers): ¿Hay picos o caídas abruptas que no siguen el patrón?

4.2 Gráfico de la serie temporal

# Gráfico estático con plot base
plot(ts_A,
     col  = "steelblue",
     lwd  = 2,
     main = "Ingresos diarios por COVID-19 — Hospital A",
     ylab = "Número de ingresos",
     xlab = "Fecha")
grid(col = "gray90")

Interpretación visual:

Se distinguen claramente varias olas pandémicas, con picos pronunciados de ingresos seguidos de periodos de descenso.

La serie no es estacionaria: la media y la varianza cambian a lo largo del tiempo, lo que es esperable en una epidemia con distintas fases.

Se observan oscilaciones de alta frecuencia (probablemente variación semanal) superpuestas a la tendencia de cada ola.

4.3 Gráfico interactivo con `dygraphs`

(Disponible solo en la versión HTML del documento)

dygraph(ts_A, main = "Ingresos diarios COVID-19 — Hospital A") %>%
  dyAxis("y", label = "Número de ingresos") %>%
  dyOptions(strokeWidth = 1.5, colors = "steelblue") %>%
  dyRangeSelector()

4.4 Descomposición de la serie temporal

La función decompose() separa la serie en tres componentes aditivos:

\[Y_t = T_t + S_t + R_t\]

\(T_t\) Tendencia (Trend): dirección general a largo plazo
\(S_t\) Componente estacional (Seasonal): patrón cíclico que se repite con la misma frecuencia
\(R_t\) Residual: lo que no explican los anteriores (idealmente debe ser ruido aleatorio)

Con datos diarios, usamos frequency = 7 para capturar el patrón semanal:

# Convertimos a objeto ts con frecuencia semanal (7 días)
ts_A_weekly <- ts(data_A$Admissions, frequency = 7)

# Descomposición
decomp_A <- decompose(ts_A_weekly)

# Visualización con ggplot2
autoplot(decomp_A) +
  ggtitle("Descomposición de la serie — Hospital A (frecuencia semanal)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

Interpretación de la descomposición:

Tendencia (trend): Muestra claramente las distintas olas de la pandemia. Se observan picos (momentos de máxima presión hospitalaria) y valles (descensos entre olas). La tendencia no es lineal, lo que confirma la naturaleza epidémica de los datos.

Estacionalidad (seasonal): Existe un patrón semanal repetitivo. Los lunes y martes tienden a registrar más ingresos, mientras que los fines de semana los ingresos son menores. Esto tiene sentido clínico: las admisiones del fin de semana a menudo se registran el lunes.

Residual: Los residuos deberían ser aleatorios (ruido blanco). Si tienen estructura (picos o patrones), significa que el modelo no ha capturado toda la información de la serie.

4.5 Estacionalidad semanal: ¿qué día hay más ingresos?

# Añadimos el día de la semana al dataset
data_A$DiaSemana <- wday(data_A$Dates, label = TRUE, abbr = FALSE,
                          week_start = 1)

# Media de ingresos por día de la semana
media_por_dia <- aggregate(data_A$Admissions,
                            by  = list(data_A$DiaSemana),
                            FUN = mean)
names(media_por_dia) <- c("Dia", "MediaIngresos")

# Ordenamos para ver el ranking
media_por_dia[order(media_por_dia$MediaIngresos, decreasing = TRUE), ]

ggplot(media_por_dia,
       aes(x = Dia, y = MediaIngresos, fill = MediaIngresos)) +
  geom_col() +
  scale_fill_gradient(low = "steelblue", high = "tomato") +
  ggtitle("Media de ingresos COVID por día de la semana — Hospital A") +
  ylab("Media de ingresos") + xlab("") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title      = element_text(hjust = 0.5, face = "bold"))

Interpretación clínica: El día de la semana con mayor media de ingresos es clave para la planificación de recursos hospitalarios (personal, camas UCI). Este patrón sistemático refleja comportamientos tanto de la enfermedad como del sistema sanitario (retrasos en el registro del fin de semana).

4.6 Comparación Hospital A vs Hospital C

# Preparamos los datos del Hospital C de la misma manera
admissions_C <- admissions[admissions$HospitalID == "C", ]
data_C <- aggregate(admissions_C$NumberOfAdmissions,
                    by  = list(admissions_C$HospitalAdmitDate),
                    FUN = sum)
names(data_C) <- c("Dates", "Admissions")
data_C$Dates  <- as.Date(data_C$Dates)
ts_C <- zoo(x = data_C$Admissions, order.by = data_C$Dates)

# Totales acumulados
total_A <- sum(coredata(ts_A), na.rm = TRUE)
total_C <- sum(coredata(ts_C), na.rm = TRUE)

cat("Total ingresos Hospital A:", total_A, "\n")

#> Total ingresos Hospital A: 1341

cat("Total ingresos Hospital C:", total_C, "\n")

#> Total ingresos Hospital C: 1952

cat("Diferencia (A - C):       ", total_A - total_C, "\n")

#> Diferencia (A - C):        -611

# Gráfico comparativo interactivo
dygraph(cbind(ts_A, ts_C),
        main = "Ingresos COVID-19: Hospital A vs Hospital C") %>%
  dySeries("ts_A", label = "Hospital A", color = "steelblue") %>%
  dySeries("ts_C", label = "Hospital C", color = "tomato")  %>%
  dyOptions(strokeWidth = 1.5) %>%
  dyRangeSelector()

Interpretación: Ambos hospitales experimentaron las mismas olas pandémicas de forma sincronizada, lo que confirma que el patrón es poblacional y no específico de un centro. Las diferencias en magnitud reflejan distintas capacidades hospitalarias o distintas zonas de influencia.

5 Bloque 4: Datos faltantes y limpieza

5.1 ¿Por qué son problemáticos los datos faltantes en series temporales?

Las series temporales requieren periodos iguales y continuos. Un NA rompe esa continuidad porque los modelos ARIMA y de descomposición asumen que los datos están igualmente espaciados en el tiempo. Hay que imputar (rellenar) los valores faltantes antes de modelar.

5.2 Carga y detección de valores faltantes

mydata <- read.csv("Rmissing.csv")
myts   <- ts(mydata$mydata)

# ¿Cuántos NA tiene la serie?
n_na <- sum(is.na(myts))
cat("Número de valores faltantes (NA):", n_na, "\n")

#> Número de valores faltantes (NA): 5

cat("Porcentaje de datos faltantes:   ", round(n_na/length(myts)*100, 1), "%\n")

#> Porcentaje de datos faltantes:    2 %

plot(myts, main = "Serie con valores faltantes y outliers",
     col = "steelblue", lwd = 1.5, ylab = "Valor")
grid(col = "gray90")

5.3 Los tres métodos de imputación

5.3.1 Método 1: LOCF — Last Observation Carried Forward

Copia el último valor conocido hacia adelante hasta el siguiente valor disponible.

Ventaja: Muy simple.
Desventaja: Puede distorsionar la varianza y crear “escalones” artificiales en la serie.

myts_locf <- na.locf(myts)

5.3.2 Método 2: Interpolación lineal ✅ (Recomendado)

Estima los NA trazando una línea recta entre el valor anterior y el posterior conocido.

Ventaja: Respeta la continuidad temporal, produce transiciones suaves.
Desventaja: Asume linealidad entre los puntos (no siempre realista para series muy volátiles).

myts_interp <- na.interp(myts)

5.3.3 Método 3: Relleno con la media global

Sustituye todos los NA por la media de toda la serie.

Ventaja: Neutro respecto a la tendencia global.
Desventaja: Ignora completamente la estructura temporal local; introduce discontinuidades.

media_serie <- mean(myts, na.rm = TRUE)
myts_fill   <- na.fill(myts, media_serie)
cat("Media global de la serie:", round(media_serie, 2), "\n")

#> Media global de la serie: 50.71

5.4 Comparación visual de los tres métodos

par(mfrow = c(3, 1), mar = c(3, 4, 2.5, 1))
plot(myts_locf,   main = "Método 1: LOCF (último valor hacia adelante)",
     col = "steelblue",  lwd = 1.5, ylab = "Valor")
plot(myts_interp, main = "Método 2: Interpolación lineal (RECOMENDADO)",
     col = "darkorange", lwd = 1.5, ylab = "Valor")
plot(myts_fill,   main = "Método 3: Media global",
     col = "darkgreen",  lwd = 1.5, ylab = "Valor")

par(mfrow = c(1, 1))

¿Cuándo usar cada método?

Interpolación lineal: La opción más adecuada para series temporales biomédicas donde los valores faltantes son pocos y la serie no tiene cambios bruscos. Preserva la dinámica temporal.

LOCF: Útil cuando los datos son de tipo “estado” (ej: dosis de medicamento que se mantiene hasta el siguiente cambio).

Media global: Solo si los NA son muy pocos, aislados, y la serie es estacionaria sin tendencia visible.

5.5 Limpieza completa con `tsclean()`

La función tsclean() hace dos cosas a la vez: elimina outliers (valores extremos) e imputa NA, usando interpolación robusta basada en descomposición STL.

myts_clean <- tsclean(myts)

plot(myts_clean,
     main = "Serie limpia con tsclean() — outliers e NA tratados",
     col  = "purple", lwd = 2, ylab = "Valor")
grid(col = "gray90")

¿Cuándo usar tsclean()?
Cuando la serie tiene tanto valores faltantes como outliers. Es especialmente útil en datos biomédicos donde pueden ocurrir errores de registro o lecturas extremas de equipos. La desventaja es que modifica automáticamente más datos, por lo que hay que usarla con criterio clínico.

6 Bloque 5: Estacionariedad y test de Dickey-Fuller

6.1 ¿Qué es la estacionariedad y por qué importa?

Una serie es estacionaria cuando sus propiedades estadísticas (media, varianza, autocorrelación) son constantes en el tiempo. Es decir, no tiene tendencia ni cambios sistemáticos en la variabilidad.

La mayoría de los modelos ARIMA requieren que la serie sea estacionaria. Si no lo es, hay que transformarla (diferenciando).

\[\text{Una serie } Y_t \text{ es estacionaria si:} \quad E[Y_t] = \mu \quad \forall t \quad \text{y} \quad \text{Var}(Y_t) = \sigma^2 \quad \forall t\]

6.2 Test de Dickey-Fuller Aumentado (ADF)

El test ADF contrasta formalmente si una serie tiene raíz unitaria (no estacionaria).

\[H_0: \text{La serie NO es estacionaria (tiene raíz unitaria)}\] \[H_1: \text{La serie ES estacionaria}\]

Regla de decisión:

p-valor	Decisión	Conclusión
p < 0.05	Rechazamos H₀	La serie es estacionaria
p ≥ 0.05	No rechazamos H₀	La serie NO es estacionaria

6.3 Aplicación al Hospital A

adf_resultado <- adf.test(data_A$Admissions)
adf_resultado

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  data_A$Admissions
#> Dickey-Fuller = -2.409, Lag order = 6, p-value = 0.4041
#> alternative hypothesis: stationary

Interpretación del test ADF — Serie original:

El p-valor obtenido (0.4041) nos permite tomar una decisión estadística:

Si p < 0.05: la serie es estacionaria (no necesita diferenciación, d = 0)

Si p ≥ 0.05: la serie NO es estacionaria (necesita diferenciación)

Desde el punto de vista clínico, que la serie no sea estacionaria tiene sentido: los ingresos por COVID-19 no fluctúan alrededor de un nivel fijo, sino que siguen la dinámica de las olas pandémicas, con medias muy distintas en cada periodo.

6.4 Diferenciación para conseguir estacionariedad

La diferenciación elimina la tendencia calculando la diferencia entre observaciones consecutivas:

\[Y'_t = Y_t - Y_{t-1}\]

El parámetro d en ARIMA(p, d, q) indica cuántas diferenciaciones se aplicaron.

ts_A_diff <- diff(data_A$Admissions)

# Visualización comparativa
par(mfrow = c(2, 1), mar = c(3, 4, 2.5, 1))
plot(data_A$Admissions, type = "l",
     main = "Serie original — Hospital A (con tendencia)", 
     col = "steelblue", lwd = 1.5, ylab = "Ingresos")
plot(ts_A_diff, type = "l",
     main = "Serie diferenciada (d = 1) — sin tendencia",
     col = "tomato", lwd = 1.5, ylab = "Cambio en ingresos")
abline(h = 0, lty = 2, col = "gray40")

par(mfrow = c(1, 1))

adf_diff <- adf.test(ts_A_diff)
adf_diff

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  ts_A_diff
#> Dickey-Fuller = -9.5667, Lag order = 6, p-value = 0.01
#> alternative hypothesis: stationary

Interpretación tras diferenciación:

Después de aplicar una diferenciación (d = 1), el p-valor cae a 0.01.

Si p < 0.05: ✅ La serie diferenciada ya es estacionaria → usaremos d = 1 en ARIMA

Esto significa que los cambios diarios en ingresos (cuánto varía de un día al siguiente) sí son estables, aunque el nivel absoluto no lo sea.

7 Bloque 6: ACF y PACF — Identificación del modelo

7.1 ¿Para qué sirven ACF y PACF?

Estas dos funciones nos ayudan a identificar los parámetros p y q del modelo ARIMA:

Función	Qué mide	Parámetro que informa
ACF (Función de Autocorrelación)	Correlación de la serie con sus propios retardos pasados	q (parte MA)
PACF (Función de Autocorrelación Parcial)	Correlación directa con cada retardo, eliminando el efecto de los intermedios	p (parte AR)

Las líneas azules discontinuas marcan el umbral de significación estadística (95%). Un retardo es significativo si su barra supera ese umbral.

Reglas empíricas: - Si la ACF decae gradualmente y la PACF corta abruptamente → modelo AR(p) - Si la PACF decae gradualmente y la ACF corta abruptamente → modelo MA(q) - Si ambas decaen gradualmente → modelo ARMA(p,q)

par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))
acf(ts_A_diff,  main = "ACF — Serie diferenciada Hospital A",  lag.max = 30)
pacf(ts_A_diff, main = "PACF — Serie diferenciada Hospital A", lag.max = 30)

par(mfrow = c(1, 1))

Interpretación de ACF y PACF:

ACF: Observamos los retardos que superan las líneas azules. Esto nos da una estimación del parámetro q (número de términos de media móvil).

PACF: Observamos los retardos significativos. Esto nos da una estimación del parámetro p (número de términos autorregresivos).

Aunque podemos identificar p y q manualmente, en la práctica usamos auto.arima() que prueba muchas combinaciones automáticamente y selecciona la mejor por criterio estadístico.

8 Bloque 7: Modelado ARIMA y predicción

8.1 ¿Qué es un modelo ARIMA?

ARIMA(p, d, q) es el modelo estándar para series temporales:

\[\underbrace{AR(p)}_{\text{Autorregresivo}} + \underbrace{I(d)}_{\text{Integrado}} + \underbrace{MA(q)}_{\text{Media Móvil}}\]

p: Número de retardos autorregresivos (la serie depende de sus p valores anteriores)
d: Orden de diferenciación (cuántas veces se diferencia para lograr estacionariedad)
q: Número de términos de media móvil (la serie depende de q errores pasados)

8.2 Auto ARIMA: selección automática del mejor modelo

auto.arima() prueba múltiples combinaciones de (p, d, q) y selecciona la que minimiza el Criterio de Información de Akaike (AIC):

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

donde \(\mathcal{L}\) es la verosimilitud del modelo y \(k\) el número de parámetros. A menor AIC, mejor modelo (equilibrio entre ajuste y parsimonia).

# Creamos objeto xts para usar con forecast
df_A.ts <- xts(data_A$Admissions, order.by = data_A$Dates)
colnames(df_A.ts) <- "Admissions"

# Auto ARIMA con búsqueda exhaustiva (más lento pero más preciso)
modelo_arima <- auto.arima(df_A.ts,
                            stepwise      = FALSE,
                            approximation = FALSE,
                            trace         = FALSE)   # trace=TRUE para ver el proceso

# Resumen del modelo seleccionado
summary(modelo_arima)

#> Series: df_A.ts 
#> ARIMA(1,1,2) 
#> 
#> Coefficients:
#>          ar1      ma1     ma2
#>       0.8614  -1.6913  0.7427
#> s.e.  0.0866   0.0871  0.0713
#> 
#> sigma^2 = 5.831:  log likelihood = -769.71
#> AIC=1547.42   AICc=1547.55   BIC=1562.68
#> 
#> Training set error measures:
#>                      ME     RMSE     MAE      MPE     MAPE      MASE       ACF1
#> Training set 0.01226008 2.400334 1.69382 -35.2989 63.43883 0.7958342 0.03133127

cat("╔══════════════════════════════════════╗\n")

#> ╔══════════════════════════════════════╗

cat("  Modelo seleccionado: ARIMA",
    "(",modelo_arima$arma[1],",",
    modelo_arima$arma[6],",",
    modelo_arima$arma[2],")\n")

#>   Modelo seleccionado: ARIMA ( 1 , 1 , 2 )

cat("  AIC:", round(modelo_arima$aic, 2), "\n")

#>   AIC: 1547.42

cat("  AICc:", round(modelo_arima$aicc, 2), "\n")

#>   AICc: 1547.55

cat("  BIC:", round(modelo_arima$bic, 2), "\n")

#>   BIC: 1562.68

cat("╚══════════════════════════════════════╝\n")

#> ╚══════════════════════════════════════╝

Interpretación del modelo ARIMA seleccionado:

Los valores de p, d, q que aparecen en ARIMA(p,d,q) significan:

p = número de retardos autorregresivos: el modelo usa los últimos p días para predecir el siguiente

d = orden de diferenciación: confirma cuántas veces hubo que diferenciar para lograr estacionariedad

q = términos de media móvil: el modelo también incorpora los errores de predicción recientes

El AIC más bajo indica el mejor compromiso entre ajuste a los datos y simplicidad del modelo.

8.3 Diagnóstico de residuos

Un buen modelo ARIMA debe dejar residuos que se comporten como ruido blanco: aleatorios, sin estructura, con media cero y sin autocorrelación. Si los residuos tienen estructura, significa que el modelo no ha capturado toda la información de la serie.

El test de Ljung-Box contrasta formalmente si hay autocorrelación en los residuos:

\[H_0: \text{Los residuos no tienen autocorrelación (son ruido blanco)}\] \[H_1: \text{Los residuos sí tienen autocorrelación}\]

checkresiduals(modelo_arima)

#> 
#>  Ljung-Box test
#> 
#> data:  Residuals from ARIMA(1,1,2)
#> Q* = 8.9529, df = 7, p-value = 0.2561
#> 
#> Model df: 3.   Total lags used: 10

lb_test <- Box.test(residuals(modelo_arima), type = "Ljung-Box", lag = 20)
cat("Test de Ljung-Box:\n")

#> Test de Ljung-Box:

cat("  Estadístico Q:", round(lb_test$statistic, 4), "\n")

#>   Estadístico Q: 12.6552

cat("  p-valor:      ", round(lb_test$p.value, 4), "\n")

#>   p-valor:       0.8917

cat("  Conclusión:   ",
    ifelse(lb_test$p.value > 0.05,
           "✅ Residuos son ruido blanco (buen modelo)",
           "⚠️  Residuos tienen autocorrelación (modelo mejorable)"), "\n")

#>   Conclusión:    ✅ Residuos son ruido blanco (buen modelo)

Interpretación del diagnóstico de residuos:

Analizamos los cuatro paneles de checkresiduals():

Residuos en el tiempo: Deben fluctuar aleatoriamente alrededor de cero, sin tendencia ni heteroscedasticidad (varianza cambiante). Si hay patrones, el modelo no es suficiente.

ACF de los residuos: La mayoría de las barras deben estar dentro de las líneas azules. Si hay barras significativas, hay autocorrelación residual y el modelo necesita parámetros adicionales.

Histograma: Los residuos deben seguir aproximadamente una distribución normal centrada en cero.

Test de Ljung-Box: Si p > 0.05, no rechazamos H₀ → los residuos son ruido blanco → el modelo es adecuado.

8.4 Predicción a 30 días

# Predicción para los próximos 30 días
prediccion_arima <- forecast(modelo_arima, h = 30)

autoplot(prediccion_arima) +
  ggtitle("Predicción a 30 días — Hospital A (Auto ARIMA)") +
  ylab("Ingresos diarios") +
  xlab("Tiempo") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Valores numéricos de la predicción para los primeros 7 días
pred_df <- data.frame(
  Dia     = 1:7,
  Pred    = round(prediccion_arima$mean[1:7], 1),
  IC_80_L = round(prediccion_arima$lower[1:7, 1], 1),
  IC_80_U = round(prediccion_arima$upper[1:7, 1], 1),
  IC_95_L = round(prediccion_arima$lower[1:7, 2], 1),
  IC_95_U = round(prediccion_arima$upper[1:7, 2], 1)
)
print(pred_df)

#>   Dia Pred IC_80_L IC_80_U IC_95_L IC_95_U
#> 1   1  2.9    -0.2     6.0    -1.9     7.6
#> 2   2  2.9    -0.3     6.0    -1.9     7.7
#> 3   3  2.9    -0.3     6.1    -2.0     7.8
#> 4   4  2.9    -0.3     6.2    -2.1     7.9
#> 5   5  3.0    -0.4     6.3    -2.2     8.1
#> 6   6  3.0    -0.5     6.4    -2.3     8.3
#> 7   7  3.0    -0.6     6.5    -2.5     8.4

Interpretación de la predicción:

El gráfico muestra:

Línea azul oscuro: Valor predicho (punto estimado)

Banda azul media: Intervalo de confianza al 80% — hay un 80% de probabilidad de que el valor real caiga dentro de esta banda

Banda azul clara: Intervalo de confianza al 95% — hay un 95% de probabilidad de que el valor real caiga aquí

¿Intervalos amplios o estrechos? En series pandémicas, los intervalos tienden a ensancharse rápidamente con el horizonte de predicción. Esto refleja la alta incertidumbre inherente al comportamiento epidémico, que depende de factores externos (variantes virales, medidas de control, comportamiento social) que el modelo no puede anticipar.

Implicación clínica: Para planificación hospitalaria a corto plazo (3-7 días), el modelo puede ser útil. Para predicciones a más de 2-3 semanas, la incertidumbre es tan alta que los intervalos de confianza pierden valor práctico.

8.5 Ajuste del modelo: valores reales vs estimados

dygraph(cbind(Real     = prediccion_arima$x,
              Estimado = fitted(prediccion_arima)),
        main = "Ingresos reales vs estimados por el modelo — Hospital A") %>%
  dySeries("Real",     label = "Real",     color = "steelblue") %>%
  dySeries("Estimado", label = "Estimado", color = "tomato") %>%
  dyOptions(strokeWidth = 1.5) %>%
  dyRangeSelector()

Interpretación del ajuste:

Comparamos la serie real (azul) con los valores que el modelo habría predicho en ese mismo periodo (rojo/naranja). Un buen modelo debe:

Seguir los picos y valles de la serie real

No tener errores sistemáticos (siempre sobreestimar o siempre subestimar)

Los periodos donde el modelo sobreestima (estimado > real) podrían corresponder a momentos donde medidas de control redujeron bruscamente los ingresos. Los periodos donde subestima suelen coincidir con el inicio de nuevas olas (el modelo no anticipa el cambio de tendencia hasta que ya ha ocurrido).

9 Bloque 8: Ejercicio integrador — Cáncer de mama

9.1 Preparación de los datos

El dataset de mortalidad por cáncer de mama tiene particularidades importantes:

España: Valores numéricos directos (miles de muertes)
Francia: Valores con sufijo "k" (ej: "12.5k") que hay que transformar

df_cancer <- read.csv("breast_cancer_number_of_female_deaths.csv",
                      check.names = FALSE,
                      row.names    = 1)

# España: valores numéricos normales
ts_spain <- ts(as.numeric(df_cancer["Spain", ]),
               start = 1990, frequency = 1)

# Francia: eliminar "k" y multiplicar × 1000
france_raw <- as.character(unlist(df_cancer["France", ]))
france_num <- as.numeric(gsub("k", "", france_raw)) * 1000
ts_france  <- ts(france_num, start = 1990, frequency = 1)

# Verificación de calidad
cat("NAs en España: ", sum(is.na(ts_spain)),  "\n")

#> NAs en España:  0

cat("NAs en Francia:", sum(is.na(ts_france)), "\n")

#> NAs en Francia: 0

cat("\nPrimeras observaciones España:\n"); print(head(ts_spain))

#> 
#> Primeras observaciones España:

#> Time Series:
#> Start = 1990 
#> End = 1995 
#> Frequency = 1 
#> [1] 6480 6620 6650 6760 6820 6840

cat("\nPrimeras observaciones Francia (en muertes reales):\n"); print(head(ts_france))

#> 
#> Primeras observaciones Francia (en muertes reales):

#> Time Series:
#> Start = 1990 
#> End = 1995 
#> Frequency = 1 
#> [1] 12500 12700 12800 13100 13100 13300

9.2 Análisis para España

9.2.1 Visualización

autoplot(ts_spain) +
  ggtitle("Muertes anuales por cáncer de mama — España (1990–2019)") +
  ylab("Número de muertes") +
  xlab("Año") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

Interpretación visual — España:
La serie muestra una tendencia que incluye dos fases diferenciadas: un crecimiento inicial desde 1990 hasta aproximadamente el año 2000 (posiblemente por aumento de la incidencia y mejora en el diagnóstico), seguido de una relativa estabilización o ligero ascenso, aunque el número absoluto de muertes refleja también el crecimiento y envejecimiento de la población femenina española.

9.2.2 Test de estacionariedad

adf_spain <- adf.test(ts_spain)
adf_spain

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  ts_spain
#> Dickey-Fuller = -1.1619, Lag order = 3, p-value = 0.8947
#> alternative hypothesis: stationary

Interpretación ADF — España:
p-valor = 0.8947

Si p ≥ 0.05: No rechazamos H₀ → la serie NO es estacionaria → auto.arima() necesitará aplicar diferenciación (d ≥ 1).

Esto es esperable: las muertes por cáncer de mama en España siguen una tendencia temporal, no fluctúan aleatoriamente alrededor de un valor fijo.

9.2.3 Auto ARIMA para España

modelo_spain <- auto.arima(ts_spain,
                            stepwise      = FALSE,
                            approximation = FALSE,
                            trace         = FALSE)

summary(modelo_spain)

#> Series: ts_spain 
#> ARIMA(0,1,0) with drift 
#> 
#> Coefficients:
#>         drift
#>       51.7241
#> s.e.  13.0859
#> 
#> sigma^2 = 5145:  log likelihood = -164.55
#> AIC=333.1   AICc=333.56   BIC=335.83
#> 
#> Training set error measures:
#>                     ME     RMSE      MAE         MPE      MAPE      MASE
#> Training set 0.2142758 69.29529 53.61657 -0.01706144 0.7592384 0.7547964
#>                   ACF1
#> Training set 0.1536229

cat("\nModelo seleccionado para España: ARIMA(",
    modelo_spain$arma[1],",",
    modelo_spain$arma[6],",",
    modelo_spain$arma[2],")\n")

#> 
#> Modelo seleccionado para España: ARIMA( 0 , 1 , 0 )

cat("AIC:", round(modelo_spain$aic, 2), "\n")

#> AIC: 333.1

9.2.4 Diagnóstico de residuos

checkresiduals(modelo_spain)

#> 
#>  Ljung-Box test
#> 
#> data:  Residuals from ARIMA(0,1,0) with drift
#> Q* = 7.6472, df = 6, p-value = 0.2651
#> 
#> Model df: 0.   Total lags used: 6

lb_spain <- Box.test(residuals(modelo_spain), type = "Ljung-Box")
cat("Test de Ljung-Box (España):\n")

#> Test de Ljung-Box (España):

cat("  p-valor:", round(lb_spain$p.value, 4), "\n")

#>   p-valor: 0.3768

cat("  Conclusión:",
    ifelse(lb_spain$p.value > 0.05,
           "✅ Residuos son ruido blanco",
           "⚠️  Residuos con autocorrelación"), "\n")

#>   Conclusión: ✅ Residuos son ruido blanco

9.2.5 Predicción para España: próximos 5 años

pred_spain <- forecast(modelo_spain, h = 5)

autoplot(pred_spain) +
  ggtitle("Predicción 5 años — Muertes por cáncer de mama en España") +
  ylab("Número de muertes") +
  xlab("Año") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

cat("Predicción para España (próximos 5 años):\n")

#> Predicción para España (próximos 5 años):

pred_spain_df <- data.frame(
  Año    = time(pred_spain$mean),
  Pred   = round(pred_spain$mean, 0),
  IC_95_L = round(pred_spain$lower[, 2], 0),
  IC_95_U = round(pred_spain$upper[, 2], 0)
)
print(pred_spain_df)

#>    Año Pred IC_95_L IC_95_U
#> 1 2020 8032    7891    8172
#> 2 2021 8083    7885    8282
#> 3 2022 8135    7892    8379
#> 4 2023 8187    7906    8468
#> 5 2024 8239    7924    8553

Interpretación de la predicción para España:

El modelo proyecta la tendencia histórica hacia los próximos años. Los intervalos de confianza al 95% nos indican el rango de incertidumbre de esa predicción.

Implicación sanitaria: Una tendencia al alza en números absolutos no necesariamente implica fracaso de las políticas de salud pública. El envejecimiento de la población aumenta la incidencia de cáncer de mama. Para evaluar correctamente el impacto de programas de cribado (mamografía) y tratamientos, debería analizarse la tasa de mortalidad por 100.000 mujeres ajustada por edad, no el número absoluto.

9.3 Comparación España vs Francia

# Combinamos las dos series
ts_ambos <- cbind(ts_spain, ts_france)

# Totales acumulados
total_spain  <- sum(coredata(ts_spain),  na.rm = TRUE)
total_france <- sum(coredata(ts_france), na.rm = TRUE)

cat("Total muertes acumuladas 1990–2019:\n")

#> Total muertes acumuladas 1990–2019:

cat("  España: ", format(total_spain,  big.mark = "."), "\n")

#>   España:  212.150

cat("  Francia:", format(total_france, big.mark = "."), "\n")

#>   Francia: 424.200

cat("  Ratio Francia/España:", round(total_france/total_spain, 2), "\n")

#>   Ratio Francia/España: 2

dygraph(ts_ambos,
        main = "Muertes por cáncer de mama: España vs Francia (1990–2019)") %>%
  dySeries("ts_spain",  label = "España",  color = "steelblue") %>%
  dySeries("ts_france", label = "Francia", color = "tomato")   %>%
  dyOptions(strokeWidth = 2) %>%
  dyRangeSelector()

Interpretación comparativa España vs Francia:

Francia registra consistentemente más muertes absolutas que España a lo largo de toda la serie. Sin embargo, esta comparación directa es estadísticamente incompleta por dos razones fundamentales:

Diferencia poblacional: Francia tiene aproximadamente el doble de población que España (~67 millones vs ~47 millones de habitantes). Por tanto, es esperable que tenga más muertes absolutas.

Estructura demográfica: La distribución por edad de la población femenina es diferente entre ambos países, y el cáncer de mama es mucho más frecuente en mujeres mayores.

Conclusión metodológica: Para una comparación válida entre países se deben usar tasas estandarizadas por edad por 100.000 mujeres, que permiten comparar la mortalidad eliminando el efecto del tamaño y la estructura de la población. Los números absolutos solo son informativos para la planificación de recursos sanitarios (oncólogos, radioterapia, quimioterapia) dentro de cada país.

10 Resumen y conclusiones

10.1 Conceptos clave aprendidos

Resumen de conceptos de la práctica
Concepto	Descripción
Objeto zoo/xts	Serie temporal con fechas reales como índice
Descomposición	Separa tendencia, estacionalidad y residuos
Datos faltantes	Imputar antes de modelar (preferible interpolación lineal)
Test ADF	Contrasta estacionariedad: p<0.05 → estacionaria
Diferenciación	Elimina tendencia; el número de diferenciaciones = d en ARIMA
ACF y PACF	ACF → q (MA); PACF → p (AR)
Auto ARIMA	Selecciona ARIMA(p,d,q) minimizando AIC
Ljung-Box	Contrasta si residuos son ruido blanco: p>0.05 → OK
Predicción	La incertidumbre crece con el horizonte temporal

10.2 Flujo de trabajo en análisis de series temporales

Datos brutos
     ↓
1. Formatear fechas (lubridate)
     ↓
2. Crear serie temporal (zoo/xts)
     ↓
3. Visualizar y descomponer (plot, decompose)
     ↓
4. Tratar datos faltantes (na.interp / tsclean)
     ↓
5. Test de estacionariedad (adf.test)
     ↓ (si no es estacionaria)
6. Diferenciar (diff) hasta p < 0.05
     ↓
7. Identificar p y q (ACF, PACF)
     ↓
8. Ajustar modelo (auto.arima)
     ↓
9. Diagnosticar residuos (checkresiduals, Ljung-Box)
     ↓ (si los residuos son ruido blanco)
10. Predecir (forecast)

# Información de la sesión para reproducibilidad
sessionInfo()

#> R version 4.5.1 (2025-06-13 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: Europe/Madrid
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] dygraphs_1.1.1.6 xts_0.14.2       zoo_1.8-15       ggplot2_4.0.2   
#> [5] forecast_9.0.2   tseries_0.10-61  lubridate_1.9.5 
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10        generics_0.1.4     lattice_0.22-7     digest_0.6.37     
#>  [5] magrittr_2.0.3     evaluate_1.0.5     grid_4.5.1         timechange_0.4.0  
#>  [9] RColorBrewer_1.1-3 fastmap_1.2.0      jsonlite_2.0.0     scales_1.4.0      
#> [13] jquerylib_0.1.4    cli_3.6.5          rlang_1.1.6        withr_3.0.2       
#> [17] cachem_1.1.0       yaml_2.3.10        tools_4.5.1        parallel_4.5.1    
#> [21] dplyr_1.1.4        colorspace_2.1-2   curl_7.0.0         vctrs_0.6.5       
#> [25] R6_2.6.1           lifecycle_1.0.5    htmlwidgets_1.6.4  pkgconfig_2.0.3   
#> [29] urca_1.3-4         bslib_0.9.0        pillar_1.11.0      gtable_0.3.6      
#> [33] glue_1.8.0         quantmod_0.4.28    Rcpp_1.1.0         xfun_0.56         
#> [37] tibble_3.3.0       tidyselect_1.2.1   rstudioapi_0.18.0  knitr_1.50        
#> [41] farver_2.1.2       nlme_3.1-168       htmltools_0.5.8.1  labeling_0.4.3    
#> [45] rmarkdown_2.30     timeDate_4051.111  fracdiff_1.5-3     compiler_4.5.1    
#> [49] quadprog_1.5-8     S7_0.2.0           TTR_0.24.4

Documento generado con R Markdown — Máster en Bioinformática, Asignatura de Series Temporales

Práctica: series temporales

Análisis completo con R - COVID-19, Cáncer de mama y datos faltantes

Series temporales

16 de April de 2026