Conforme se ha referido previamente, se emplea el conjunto de datos global informado diariamente, disponible para su descarga en https://covid19.who.int/WHO-COVID-19-global-data.csv. Resulta relevante destacar que la base de datos consultada corresponde al 16 de Enero del 2023, restringiéndose a los datos concernientes exclusivamente a los casos confirmados en Irán entre el 20 de febrero y el 15 de agosto de 2020.
El análisis posterior se ha llevado a cabo empleando el software R versión 4.3.2. Con el propósito de realizar el análisis de los datos y la generación de gráficos, se procedió a convertir los datos al formato ts, lo que permitió su representación como una serie temporal.
# Se crea un objeto 'Date' diario
inds <- seq(as.Date("2020-02-20"), as.Date("2020-08-15"), by = "day")
# Se crea un objeto 'serie de tiempo' de frecuencia diaria
Confirmed_ts <- ts(Confirmed_df[2],
start = c(2020, as.numeric(format(inds[1], "%j"))),
frequency = 365)La gráfica de la Figura 8.1 exhibe la serie temporal derivada de la base de datos, en la cual se evidencia la ausencia de información para los días 27 y 29 de Febrero, así como para el 02 de Marzo y el 05 de Abril de 2020. Para subsanar esta carencia de datos, se llevó a cabo una interpolación promedio a fin de sustituir los valores faltantes. La Tabla 8.1 muestra la base de datos con las modificaciones efectuadas, así como la serie de tiempo (Figura 8.2) resultante de estas correcciones.
Con el propósito de llevar a cabo una auditoría de los datos y al mismo tiempo una descripción preliminar, se ejecuta un estudio de estadística descriptiva que arroja los resultados correspondientes, incluyendo un gráfico Boxplot (Figura 8.3) para representar la información o
Min. 1st Qu. Median Mean 3rd Qu. Max.
3 1304 2181 1923 2529 3574 Los componentes identificados en la serie de tiempo de casos confirmados de COVID-19 en Irán, revelan distintos patrones y características.
En primer lugar, se observa una tendencia discernible en el gráfico de la serie temporal (Figura 8.2). Por ejemplo, entre el 30 de marzo y el 03 de mayo de 2020, se evidencia una tendencia negativa o decreciente, seguida por una tendencia creciente a partir del 03 de mayo en adelante. Estos cambios en la tendencia podrían indicar fluctuaciones significativas en la evolución de los casos confirmados durante esos periodos específicos.
En cuanto a la estacionalidad, aunque no se identifica claramente a simple vista en el periodo observado, la extensión del análisis a un periodo más amplio podría revelar patrones recurrentes o ciclos temporales característicos. Es posible que ciertos patrones estacionales se manifiesten en intervalos más extensos de la serie temporal, lo que implicaría variaciones sistemáticas y repetitivas en los datos en períodos específicos.
Por último, se destacan pequeñas subidas y bajadas en el gráfico que sugieren la presencia de ruido en la serie temporal. Estas fluctuaciones irregulares podrían atribuirse a diversas causas, como posibles errores en la recolección de datos o fluctuaciones aleatorias inherentes al comportamiento de la enfermedad. Es importante considerar estas variaciones no sistemáticas al analizar la serie temporal, ya que podrían influir en la interpretación de los patrones y tendencias observadas.
A continuación, se emplea el test de Dickey-Fuller para examinar la presencia de estacionariedad en la serie temporal. Este test fue utilizado con la finalidad de identificar la existencia de raíces unitarias en la serie, lo cual permite inferir la presencia o ausencia de estacionariedad en los datos analizados.
adf.test(Confirmed_ts, alternative = "stationary")
Augmented Dickey-Fuller Test
data: Confirmed_ts
Dickey-Fuller = -2.9529, Lag order = 5, p-value = 0.1781
alternative hypothesis: stationaryLa hipótesis nula \((H_0)\) asume la presencia de raíces unitarias, lo que indica no estacionariedad en la serie. Al obtener un \(p-\)valor superior al nivel de significancia establecido el cuál es del \(95\%\), no se rechaza la hipótesis nula, sugiriendo la ausencia de estacionariedad en la serie de tiempo de casos confirmados.
Además, se complementa la evaluación de la estacionalidad mediante la inspección de los gráficos de la función de autocorrelación (ACF) y la función de autocorrelación parcial (PACF). Estos gráficos se utilizan para identificar patrones de autocorrelación en la serie temporal, lo que permite visualizar la presencia de estacionalidad, tendencias o ciclos.
La serie de tiempo representada en la Figura 8.2 exhibe un comportamiento característico de deambulación aleatoria. Dado que el valor de la variable \(X_{t+1}\) generalmente se encuentra en proximidad al valor \(X_t\), se evidencia una autocorrelación positiva notablemente marcada entre las variables \(X_t\) y \(X_{t+1}\).
autoplot(acf(Confirmed_ts, plot = FALSE),
main="Autocorrelograma de casos confirmados.")
En la Figura 8.4 se observa que la autocorrelación (vea la ecuación (5.2)) entre \(X_t\) y \(X_{t+k}\) decrece con el incremento del retraso \(k\). Este declive conduce a la constatación de que, a un desfase de \(20\), existe una correlación bastante débil entre \(X_t\) y \(X_{t+20}\). Al analizar el gráfico de la función de autocorrelación (ACF), se aprecia que \(\rho_{20}\approx 0.19\).
La gráfica de la Función de Autocorrelación Parcial (PACF) proporciona información valiosa sobre la estructura de autocorrelación de una serie temporal una vez han sido eliminadas las correlaciones debidas a los intervalos de tiempo intermedios.
ggPacf((Confirmed_ts), main = 'Autocorrelograma parcial de casos confirmados.')
Considerando que los datos se ajustan a un modelo de series de tiempo, la Figura 8.5 indica que el valor de correlación \(\phi_{15}\) es ligeramente superior a \(0.25\), aproximadamente \(\phi_{52}\approx 0.16\), y \(\phi_{76}\approx 0.15\), mientras que para los restantes valores, la correlación parcial no es nula.
Observación. De acuerdo con la gráfica de la Función de Autocorrelación Parcial Figura 8.5, se observa un corte abrupto después del rezago 4, lo cual sugiere que las autocorrelaciones parciales más allá de ese punto no poseen significancia estadística. Por consiguiente, se infiere la posibilidad de ajustar un modelo autoregresivo AR(4) a la base de datos.
library(tswge)
coeff <- est.ar.wge(Confirmed_ts, p=4)
Coefficients of AR polynomial:
0.8656 0.1724 0.0592 -0.1318
AR Factor Table
Factor Roots Abs Recip System Freq
1-0.9605B 1.0411 0.9605 0.0000
1-0.5355B 1.8675 0.5355 0.0000
1+0.6303B+0.2563B^2 -1.2298+-1.5459i 0.5062 0.3570
coeff$phi #coeficientes[1] 0.86563811 0.17237111 0.05917558 -0.13180534coeff$xbar #media[1] 1922.868coeff$avar #varianza finita[1] 47818.04El modelo autoregresivo AR(4) se expresa mediante la siguiente ecuación:
\[ (1-0.865B-0.172B^2-0.059B^3+0.131B^4)(X_t-1922.868)+a_t, \tag{8.1}\]
donde \(\hat{\sigma}_a^2 = 47818.04\).
El análisis del ACF y PACF proporcionó información sobre la relación de los puntos de datos con sus rezagos, permitiendo observar posibles patrones estacionales. La presencia de picos significativos en estos gráficos podría indicar la existencia de estacionalidad en la serie de tiempo.
Se procede a la evaluación del rendimiento de métodos destinados al ajuste y consecuente pronóstico. Específicamente, se contempla el método de suavizamiento exponencial de Holt-Winters y el ajuste mediante un modelo de red neuronal del tipo perceptrón multicapa. Ambos procedimientos requieren la subdivisión de los datos en conjuntos destinados a entrenamiento y prueba. El set inicial, compuesto por el \(70\%\) de los datos, se emplea para el entrenamiento de los modelos, mientras que el \(30\%\) restante se reservará para llevar a cabo las pruebas pertinentes.
Confirmed_ts <- ts(Confirmed_ts,frequency=1)
tsize <- round(0.7 * nrow(Confirmed_df))
train_confirmed <- window(Confirmed_ts,end=tsize)
test_confirmed <- window(Confirmed_ts,start=tsize+1)Con el fin de determinar la descomposición más adecuada para los datos en cuestión, se empleó un criterio elaborado basado en el coeficiente de variación, el cual proporciona una recomendación entre las dos versiones disponibles.
DescRec <- function(x){
n = length(x)
di = rep(0, n-1)
ci = rep(0, n-1)
for (i in 1:n-1) {
di[i] = x[i+1] - x[i]
ci[i] = x[i+1] / x[i]
}
d <- cv(di)
c <- cv(ci) / mean(di)
if(d < c)
print("Se recomienda la descomposición aditiva")
else
print("Se recomienda la descomposición multiplicativa")
}
DescRec(train_confirmed)[1] "Se recomienda la descomposición multiplicativa"De acuerdo con la recomendación observada, se sugiere la utilización de la versión multiplicativa (vea la ecuación (5.19)). En consecuencia, se procede a mostrar la representación gráfica de la descomposición multiplicativa de la serie temporal.
ts_train <- ts(train_confirmed, frequency = 2)
components_ts <- decompose(ts_train, type = 'mult')
plot(components_ts)
Se procede ahora a la aplicación del modelo multiplicativo de Holt-Winters a la serie temporal de los datos de entrenamiento utilizando una frecuencia de dos, con el fin de permitir la aplicabilidad del modelo.
HWc <- HoltWinters(ts_train, seasonal = 'mult')
HWcHolt-Winters exponential smoothing with trend and multiplicative seasonal component.
Call:
HoltWinters(x = ts_train, seasonal = "mult")
Smoothing parameters:
alpha: 0.729859
beta : 0
gamma: 0.6494094
Coefficients:
[,1]
a 2289.283568
b 2.250000
s1 1.119185
s2 1.113547Finalmente, utilizando el modelo de entrenamiento desarrollado en la fase previa, se lleva a cabo la proyección con un horizonte de predicción igual en extensión a los datos de prueba, acompañado de un intervalo de confianza que oscila entre el \(80\%\) y el \(95\%\).
HWc_for <- forecast(HWc, h=length(test_confirmed))Las funciones aplicadas en esta sección son parte de la librería stats (2023) de R.
Posteriormente, se procede al entrenamiento del modelo MLP (Perceptrón Multicapa). La cantidad de capas ocultas y la configuración de nodos en cada capa se determinaron de manera automatizada mediante el método de validación cruzada de 5 pliegues. Asimismo, se eligió la función de activación como sigmoide, y el proceso de entrenamiento del modelo se ejecutó a lo largo de 20 iteraciones.
fitc <- mlp(train_confirmed, hd.auto.type="cv", reps=20, comb='median')
fitcMLP fit with 1 hidden node and 20 repetitions.
Univariate lags: (1,2,4)
Forecast combined using the median operator.
MSE: 57212.5871.plot(fitc)
Para llevar a cabo el pronóstico, se emplea el modelo de entrenamiento creado en la etapa anterior, manteniendo un horizonte de predicción que coincide en duración con los datos de prueba, tal como se hizo con la técnica anterior.
frcc <- forecast(fitc,h=length(test_confirmed))Las funciones aplicadas en esta sección son parte de la librería nnfor Kourentzes (2022) de R.
Con el propósito de llevar a cabo un análisis cuantitativo exhaustivo, se presenta a continuación una tabla comparativa de los resultados derivados de las dos técnicas implementadas y la base de datos de prueba. Posteriormente, se exhiben gráficas representativas de estos resultados. En la Figura 8.8 se muestra el pronóstico mediante Holt-Winters acompañado de su respectivo intervalo de confianza. En contraste, en la Figura 8.9, la gráfica punteada en color rojo representa el comportamiento real de los datos, mientras que en azul se representa el pronóstico obtenido a través de la red MLP.
Para evaluar la calidad o bondad de ajuste de los métodos utilizados en este estudio y seleccionar el modelo más apropiado, se aplican tres métricas de rendimiento, Error Cuadrático Medio (5.32), Error Absoluto Medio (5.30) y Error Porcentual Absoluto Medio (5.33) tanto en las fases de entrenamiento como en las de prueba. Los resultados correspondientes a éstas métricas se presentan en la Tabla 8.3 .
| Training | Testing | |||||
| RMSE | MAE | MAPE | RMSE | MAE | MAPE | |
| Holt-Winters | 262.9925 | 190.0482 | 20.8415 | 234.0094 | 165.8208 | 6.2967 |
| MLP | 239.3483 | 180.9937 | 14.6079 | 177.0605 | 136.4799 | 5.4441 |
Basándose en los resultados extraídos tanto de la tabla de pronósticos (Tabla 8.2) como de la tabla de errores (Tabla 8.3), se llega a la conclusión de que, para esta base de datos en particular, la técnica de redes neuronales MLP demuestra ser más efectiva en la predicción realizada. Esto se fundamenta en la evidencia de un menor error registrado en las tres métricas calculadas, tanto durante la fase de entrenamiento como en la fase de prueba.
Tras la identificación del modelo óptimo, se procedió a prever el comportamiento futuro de la serie temporal de casos confirmados para los próximos 30 días utilizando dicho modelo. Se elaboraron representaciones gráficas de la predicción de casos confirmados de COVID-19 a 30 días, realizando una comparación de la efectividad entre las implementaciones de redes neuronales en la paquetería de R y la paquetería nativa de Python, las cuales se encuentran en las figuras Figura 8.10 y Figura 8.11, respectivamente.
En la implementación de R, siguiendo el mismo procedimiento que en las fases de entrenamiento y prueba, se empleó un número específico de capas y nodos ocultos determinados automáticamente a través del método de validación cruzada de 5 pliegues. Esta configuración se llevó a cabo con una función de activación sigmoide, ejecutando 20 iteraciones para el entrenamiento de la red neuronal.
fit.mlp = mlp(ts(Confirmed_df$Confirmed), reps = 20, hd.auto.type = 'cv',
comb="median")
fore.mlp = forecast(fit.mlp, h = 30)
Los resultados del pronóstico indican que, el 14 de septiembre de 2020, se proyectan aproximadamente 2494 nuevos casos confirmados de COVID-19. Estos valores correspondientes al período de 30 días se detallan a continuación en la Tabla 8.4 .
En esta sección, se realizaron ajustes en el método para su implementación. Cada día proyectado se forma utilizando el dato del día anterior. En cada paso, se actualiza la secuencia de entrada eliminando el valor más antiguo e incorporando la predicción más reciente como el dato más reciente. Esta dinámica se representa esquemáticamente a continuación, donde \(n\) representa la extensión de la secuencia de entrada y \(T\) es la longitud de la serie temporal.
\[ \begin{split} y:\text{Observado}\quad &\quad \hat{y}:\text{Pronosticado}\\ y_{T-n+1}\quad y_{T-n+2}\quad y_{T-n+3}\quad&\cdots\quad y_{T-2}\quad y_{T-1}\quad y_T\quad \to\quad \hat{y}_{T+1}\\ y_{T-n+2}\quad y_{T-n+3}\quad y_{T-n+4}\quad&\cdots\quad y_{T-1}\quad y_{T}\quad \hat{y}_{T+1}\quad \to\quad \hat{y}_{T+2}\\ y_{T-n+3}\quad y_{T-n+4}\quad y_{T-n+5}\quad&\cdots\quad y_{T}\quad \hat{y}_{T+1}\quad \hat{y}_{T+2}\quad \to\quad \hat{y}_{T+3}\\ &\ddots\\ \end{split} \]
Se exhibe a continuación el código utilizado y el gráfico correspondiente al pronóstico generado por la red neuronal.
Los resultados del pronóstico en Figura 8.11 indican que, el 14 de Septiembre de 2020, se proyectan aproximadamente 2476 nuevos casos confirmados de COVID-19.