Todo modelo predictivo debe probarse con datos reales, para determinar su validez. Una forma de medir el grado de error es el error medio cuadrático, que se calcula

\[ EMC= \sum_i^k \frac{(y_i-\hat{y_i})^2}{n} \]

Donde \(y_i\) es el dato observado en el momento \(i\), entanto que \(\hat{y_i}\) es valor predicho. En nuestros ejemplos, como todo lo hacemos sobre el logaritmo de las curvas, se calcular el EMC sobre los valores al logaritmo.

En este día, usaremos un modelo lineal cuadrático sobre los datos en bruto, con una estructura de autocorrelación AR(2). ## Modelos de casos totales

Probaremos la adecuación de los tres modelos, que usamos hasta el momento - ARIMA(1,1,0) con trend, casos nuevos + AR(1), casos nuevos + AR(4), usando el modelo calculado el día 8 de Abril.

Se debe considerar que como los reportes del MINSAL se realizan con un día de retraso, realmente estamos usando los datos del día 7 para predecir hasta el 14.

data.frame(inicio=convertToDate(f.inicial),
           final=convertToDate(f.final))

Se puede ver que la predicción cumplió en todo momento el intervalo de confianza para ARIMA(1,1,0) y T+AR(1) y T+AR(4), pero siempre con valores menores a los valores observados. El modelo cuadrático se acerca más a los datos observados, pero no se cumple su intervalo de confianza; además, tiende a tener valores mayores a los observados.

calcular.mse<-function(x,y) {
  mean((log(x$casos)-log(y$casos))^2)
}



dato.a.inicial<-datos.casos$total[datos.casos$total$fecha <= f.inicial,]
dato.a.final<-datos.casos$total[datos.casos$total$fecha   <= f.final,]

obtener.prediccion<-function(x.comp,min.casos, n.ahead=7) {
  n.comp<-nrow(x.comp)
  x.parc<-x.comp[1:(n.comp-n.ahead), ]
  
  dato.obs<-data.frame(
    dia=x.comp$dia,
    casos=x.comp$casos,
    li=x.comp$casos, 
    ls=x.comp$casos,
    tipo="Observado",
    pais="Chile"
)
  prediccion<-prediccion.casos(list(Chile=x.parc),min.casos = min.casos, modelos = c("tar2", "tar4", "arima","cuad"),n.ahead = n.ahead)
  out<-list()
  out$dato.obs<-dato.obs
  out$dato.obs.parc<-tail(x.comp,n.ahead)
  out$p.tar21<-tail(prediccion$tar2$Chile,n.ahead)
  out$p.tar4 <-tail(prediccion$tar4$Chile,n.ahead)
  out$p.arima1<-tail(prediccion$arima$Chile,n.ahead)
  out$p.cuad<-tail(prediccion$cuad$Chile,n.ahead)
  out$para.ver<-rbind(out$dato.obs, out$p.tar21,out$p.arima1,out$p.tar4,out$p.cuad)
  
  out
}

pred.50<-obtener.prediccion(dato.a.final, 50)

ggplot(pred.50$para.ver, aes(x=dia, y=casos, color=tipo, ymin=li,ymax=ls))+geom_point()+geom_line()+geom_ribbon(alpha=0.1)+geom_hline(yintercept = 50)+scale_y_continuous(trans="log10",limits = c(3000,60000))+xlim(c(20,45))
## Warning: Removed 28 rows containing missing values (geom_point).
## Warning: Removed 28 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_hline).

El MSE del ARIMA sigue están sobre 0.20, al igual que antes, lo que es un error muy grande. AR(4) presenta menor MSE que AR(1), pero es un tanto peor que el error de la semana anterior, que estaba en 0.003 para T+AR(1). El modelo cuadrático presenta un menor MSE.

data.frame("T+AR(1)"=calcular.mse(pred.50$dato.obs.parc, pred.50$p.tar21),
           "ARIMA(1,0,0)"=calcular.mse(pred.50$dato.obs.parc, pred.50$p.arima1),
           "T+AR(4)"=calcular.mse(pred.50$dato.obs.parc, pred.50$p.tar4),
           "Cuadrático"=calcular.mse(pred.50$dato.obs.parc, pred.50$p.cuad)
           )

¿Qué pasaría si hubiesemos considerado solo desde el caso 250 en adelante, como hace dos semanas? Se puede ver que habríamos caído fuera del intervalo de confianza para las predicciones de T+AR(1) y T+AR(4). Además, el modelo exponencial también hubiera fallado en su intervalo de confianza, también sobrestimando.

pred.250<-obtener.prediccion(dato.a.final, 250)

ggplot(pred.250$para.ver, aes(x=dia, y=casos, color=tipo, ymin=li,ymax=ls))+geom_point()+geom_line()+geom_ribbon(alpha=0.1)+geom_hline(yintercept = 250)+scale_y_continuous(trans="log10",limits = c(1000,60000))+xlim(c(20,45))
## Warning: Removed 21 rows containing missing values (geom_point).
## Warning: Removed 21 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_hline).

El MSE es peor para T+AR(1) y para T+AR(4), aunque un poco mejor para ARIMA. El modelo cuadrático sigue siendo mejor.

data.frame("T+AR(1)"=calcular.mse(pred.250$dato.obs.parc, pred.250$p.tar21),
           "ARIMA(1,0,0)"=calcular.mse(pred.250$dato.obs.parc, pred.250$p.arima1),
           "T+AR(4)"=calcular.mse(pred.250$dato.obs.parc, pred.250$p.tar4),
           "Cuadrático"=calcular.mse(pred.250$dato.obs.parc, pred.50$p.cuad)
           )

Si utilizamos el corte del periodo anterior, con 20 casos, los modelos T+AR siguen subestimando, y el modelo cuadrático claramente sobreestimando.

pred.20<-obtener.prediccion(dato.a.final, 20)
ggplot(pred.20$para.ver, aes(x=dia, y=casos, color=tipo, ymin=li,ymax=ls))+geom_point()+geom_line()+geom_ribbon(alpha=0.1)+scale_y_continuous(trans="log10")+geom_hline(yintercept = 20)+scale_y_continuous(trans="log10",limits = c(1000,60000))+xlim(c(20,45))
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 21 rows containing missing values (geom_point).
## Warning: Removed 21 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_hline).

Los errores siguen siendo superiores en un orden de magnitud. Lo interesante es que el modelo cuadrático se mantiene bastante estable

data.frame("T+AR(1)"=calcular.mse(pred.20$dato.obs.parc, pred.20$p.tar21),
           "ARIMA(1,0,0)"=calcular.mse(pred.20$dato.obs.parc, pred.20$p.arima1),
           "T+AR(4)"=calcular.mse(pred.20$dato.obs.parc, pred.20$p.tar4),
           "Cuadrático"=calcular.mse(pred.20$dato.obs.parc, pred.20$p.cuad)
           )

Tanto para AR(1) como para AR(4), podemos ver que la predicción subestima sistemáticamente el valor total, en tanto que el cuadrático sobreestima, siendo sus intervalos de confianza erróneos, pero siempre con menor error cuadrático.

Predicción por regiones

Durante las últimas semanas hemos probado un modelo de predicción que suma los resultados parciales por región. Para cada región se elige un punto de corte distinto, el cual se puede consultar en el código a continuación

min.corte<-list(
    "Arica.y.Parinacota"=10,
    "Tarapacá"=10,
    "Antofagasta"=10,
    "Atacama"=2,
    "Coquimbo"=10,
    "Valparaíso"=10,
    "Metropolitana"=50,
    "O’Higgins"=10,
    "Maule"=10,
    "Ñuble"=20,
    "Biobío"=10,
    "Araucanía"=10,
    "Los.Ríos"=20,
    "Los.Lagos"=20,
    "Aysén"=2,
    "Magallanes"=10)

dc.inicio<-lapply(datos.casos,function(x) {x[x$fecha<=f.inicial,]})
dc.final<-lapply(datos.casos,function(x) {x[x$fecha<=f.final,]})

ppz.inic<-prediccion.por.zonas(min.corte,dc.inicio,n.ahead = 7, modelos = c("tar2","cuad"))
obs<-dato.a.final[dato.a.final$fecha>=f.inicial+1 & dato.a.final$fecha<=f.final,]

Se observamos la serie posterior, con detalle, podemos ver que AR(1) subestima el resultado final, en tanto que el modelo exponencial lo sobrestima.

ggplot(rbind(ppz.inic, data.frame(tipo="observado",obs[,c("dia","fecha","casos")])), aes(x=dia,y=casos, color=tipo))+geom_line()+scale_y_continuous(trans="log10",limit=c(3000,10000))+xlim(25,45)
## Warning: Removed 28 rows containing missing values (geom_path).

El MSE del TAR regional es un poco menor que el de la serie completa con T+AR(4), en tanto que el lineal cuadrático es aun mejor.

data.frame("Suma regional T+AR(1)"=calcular.mse( obs, ppz.inic[ppz.inic$tipo=="Casos nuevo: Tendencia + AR(1)",]),
           "Suma regional Lineal cuadrático"=calcular.mse(obs, ppz.inic[ppz.inic$tipo=="General: Lineal cuadrático",])
           )

Conclusión

Durante la próxima semana, la predicción se realizará usando un modelo que parta del día 50. Se eliminará el modelo ARIMA, ya que no muestra ser útil. Se agregará el modelo lineal cuadrático que si bien presenta malos intervalos de confianza, muestra un EMC bastante menor.