Prediction of confirmed active cases for next days

Conditionally to:

  1. the observed data up to 2020-04-29 (only the official confirmed cases)
  2. assuming that these data reflect covid-19 epidemic evolution;
  3. the model detailed below;

Then below is the prediction of confirmed cases for last and next days updated at the time of this report (see above).

It also estimate the probability of observing the peak (defined as less increment in cases than the previous day). Prediction limits are at 95% of probability.

Day Inf. Expected Sup. Prob of Peak Obs.
63 2020-04-26 1209 1281 1353 49 1280
64 2020-04-27 1218 1289 1362 51 1283
65 2020-04-28 1222 1295 1367 51 1285
66 2020-04-29 1226 1298 1370 52 1290
67 2020-04-30 1231 1302 1375 49 NA
68 2020-05-01 1229 1308 1390 49 NA
69 2020-05-02 1229 1315 1404 51 NA
70 2020-05-03 1229 1320 1417 50 NA
71 2020-05-04 1231 1326 1427 51 NA
72 2020-05-05 1232 1331 1438 50 NA
73 2020-05-06 1233 1336 1447 51 NA
74 2020-05-07 1234 1339 1457 50 NA
75 2020-05-08 1233 1343 1463 50 NA
76 2020-05-09 1233 1347 1470 50 NA
77 2020-05-10 1236 1350 1479 50 NA
78 2020-05-11 1236 1353 1486 50 NA
79 2020-05-12 1237 1356 1492 51 NA
80 2020-05-13 1238 1358 1500 50 NA
81 2020-05-14 1238 1360 1503 51 NA
82 2020-05-15 1238 1363 1507 50 NA
83 2020-05-16 1239 1364 1510 50 NA
84 2020-05-17 1241 1367 1517 50 NA
85 2020-05-18 1241 1369 1517 51 NA
86 2020-05-19 1242 1369 1523 51 NA
87 2020-05-20 1241 1371 1524 50 NA
88 2020-05-21 1243 1372 1530 50 NA
89 2020-05-22 1242 1374 1530 50 NA
90 2020-05-23 1243 1375 1534 50 NA
91 2020-05-24 1243 1376 1539 50 NA
92 2020-05-25 1243 1377 1540 50 NA
93 2020-05-26 1244 1378 1539 50 NA
94 2020-05-27 1243 1378 1543 50 NA
95 2020-05-28 1244 1379 1544 50 NA
96 2020-05-29 1244 1380 1548 50 NA

Evolution of increments

increment=NULL
for(i in 2:ncol(preds)) increment=cbind(increment,preds[,i]-preds[,i-1]) 
obsinc=cvirus$casos[-1]-cvirus$casos[-nrow(cvirus)]

dbinc=data.frame(day=dd[-1],yinf=apply(increment,2,quantile,p=0.025),
                 ymean=apply(increment,2,mean),ysup=apply(increment,2,quantile,p=0.975),
                 obs=c(obsinc,rep(NA,np)))

p=ggplot(dbinc, aes(x=day, y=ymean,ymin=yinf,ymax=ysup)) + 
  geom_point()+
 geom_ribbon(alpha = 0.5, colour = "yellow")+ geom_point(data=dbinc,aes(x=day,y=obs),color="red")+
  xlab("Day")+ylab("Cases")+
  ggtitle("Predicted and Observed increment of cases")

print(p)

Data

Here is the DataBase from protezione civile of only confirmed cases at the end of the Day.

These data can be smileading and the rest of analysis is subject to assuming that confirmed cases reflect evolution of covid-19 spread.

rm(list=ls())
url <- "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-regioni/dpc-covid19-ita-regioni.csv"
cvirus <- read.csv(text = getURL(url),sep=",")
cvirus <- cvirus[cvirus$denominazione_regione%in%c("Sardegna"),]
cvirus$fecha <- as.Date(as.character(cvirus$data))
cvirus$casos=as.numeric(cvirus$totale_casi)
cvirus=cvirus[order(cvirus$fecha),]
#cvirus=na.omit(cvirus)
n=nrow(cvirus)

Model for marginal counts cases

Let \(Y_t \in \mathcal{N}_0\) represents the number of cases at time (days) \(t\) where \(t=1\) is 2020-02-24.

The fitted model is an ARMA(1,1) on the Poisson mean:

\[ \begin{aligned} Y_t | \lambda_t & \sim Poisson(\lambda_t), \mbox{ for } t>0\\ \log(\lambda_t) & = \omega+\alpha\log(1+y_{t-1})+\beta\log(\lambda_{t-1}), \mbox{ for } t>1\\ \alpha,\beta,\omega & \sim N(0,10) (i.i.d.)\\ log(\lambda_1) & \sim N(-99,0.001) \end{aligned} \]

Interpetation of parameters:

  • \(\omega\) is the mean number (in log scale) of infected (actually the certified infected);
  • \(\alpha\) is the short term component (i.e. the proportion of new infected with respect to the day before);
  • \(\beta\) is the long term component that represents the evolution with respect to the mean (this is analogue to posing a GARCH on Poisson counts);

The non Bayesian and slighlty less complicated version of this model can be found here:

https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3551626

rstan_options(auto_write = FALSE)
#Sys.setenv(LOCAL_CPPFLAGS = '-march=native -mtune=native -axCORE-AVX2')
options(mc.cores = parallel::detectCores())

dotR <- file.path(Sys.getenv("HOME"), ".R")
if (!file.exists(dotR)) dir.create(dotR)
M <- file.path(dotR, ifelse(.Platform$OS.type == "windows", "Makevars.win", "Makevars"))
if (!file.exists(M)) file.create(M)
cat("\nCXX14FLAGS=-O3 -march=native -mtune=native",
    if( grepl("^darwin", R.version$os)) "CXX14FLAGS += -arch x86_64 -ftemplate-depth-256" else 
    if (.Platform$OS.type == "windows") "CXX11FLAGS=-O3 -march=corei7 -mtune=corei7" else
    "CXX14FLAGS += -fPIC",
    file = M, sep = "\n", append = TRUE)

mod.cov ="
data {
  int<lower=2> n;// number of observations
  int<lower=2> np;// number of predicted days
  int<lower=0> y[n]; // Cases
  
}

parameters {
  real alpha;
  real beta; 
  real omega;
}

transformed parameters {
  vector[n] llambda;
  llambda[1]=-99;
  for(t in 2:n) llambda[t]=omega+alpha*log(1+y[t-1])+beta*llambda[t-1];
}

model {
  // Priors
  alpha ~ normal(0,10);
  beta ~ normal(0,10);
  omega ~ normal(0,10);
  // Likelihood
  y ~ poisson_log(llambda);

}

generated quantities {
  int<lower=0> yp[n];
  int<lower=0> ypf[np];
  real llambdaf[np];
  yp[1]=0;
  for(t in 2:n) yp[t] = poisson_rng(exp(llambda[t])); //y values predicted by the model
  llambdaf[1]=omega+alpha*log(1+y[n])+beta*llambda[n];
  ypf[1]=poisson_rng(exp(llambdaf[1]));
  for(t in 2:np){
      llambdaf[t]=omega+alpha*log(1+ypf[t-1])+beta*llambdaf[t-1];
      ypf[t]=poisson_rng(exp(llambdaf[t]));
  }
}
"
ii=list(omega=0.57,alpha=0.38,beta=0.54)
init_f <- function () list(ii,ii,ii,ii)
m1 <- stan_model(model_code = mod.cov)

Hamiltonian MC is used for obtaining the posterior:

np=30
niter=10000
fit = sampling(m1,
           data=list(y=cvirus$casos,n=n,np=np), 
           iter=niter,chains = 4,
          init = init_f(),
#          control = list(adapt_delta = 0.99,max_treedepth=50),
           seed = 17,refresh=0)
save(fit,cvirus,np,n,file="stanmod-sardegna.RData")

Goodness of Fit

The model is reliable if is able to predict what observed, when taking into account prediction uncertainty. Here is the predicted mean and 95% posterior credible interval (i.e. small with respect to the mean).

load(file="stanmod-sardegna.RData")
library(bayesplot)
preds=extract(fit,pars="yp")$yp
ppc_intervals(
  y = apply(preds,2,mean),
  yrep = preds,
  x = cvirus$casos,
  prob = 0.95
)+labs(
    x = "Observed Cases",
    y = "Predicted cases",
    title = "95% posterior predictive intervals \nvs Observed",
    subtitle = "by day"
  ) +
  panel_bg(fill = "gray95", color = NA) +
  grid_lines(color = "white")+geom_abline(intercept = 0,slope=1)

Posterior Parameters

Posterior distributions of model paramters: mean and 95% credible intervals.

int.par=c("omega","alpha","beta")
plot(fit,pars = int.par,ci.level=0.95,point_est="mean")

print(fit,pars = int.par)
## Inference for Stan model: 6e55d19cd3f700e424d1da0a5fb69cf6.
## 4 chains, each with iter=10000; warmup=5000; thin=1; 
## post-warmup draws per chain=5000, total post-warmup draws=20000.
## 
##       mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
## omega 0.59    0.00 0.03 0.50 0.57 0.59 0.61  0.64   865    1
## alpha 0.39    0.01 0.16 0.25 0.29 0.32 0.41  0.84   457    1
## beta  0.53    0.01 0.15 0.09 0.51 0.59 0.63  0.66   459    1
## 
## Samples were drawn using NUTS(diag_e) at Thu Apr 30 15:10:03 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Iterative estimation of parameters and outof sample prediction.

This is interesting to monitor the evolution of covid-19.

If model is reliable, then also parameter evolution is reliable.

Further we check its out-of-sample prediction of cases since beginning of march. This is more reliable than the above goodness of fit.

ndays=sum(cvirus$fecha>"2020-03-10")
evpars=data.frame(NULL)
windows=list(1:(n-ndays+1))
for(i in (n-ndays+2):n) windows=c(windows,list(1:i))
nw=length(windows)
for(i in 1:nw){
  fit = sampling(m1,
           data=list(y=cvirus$casos[windows[[i]]],
                     n=length(windows[[i]]),np=2),
           iter=niter,chains = 4,
            init = init_f(),
#           control = list(adapt_delta = 0.99,max_treedepth=50),
           seed = 11,refresh=0)
  post.par=extract(fit,pars=c(int.par,"ypf"))
  for(j in 1:length(int.par)){
    evpars=rbind(evpars,data.frame(day=cvirus$fecha[max(windows[[i]])],
                              value=post.par[[j]],
                              param=int.par[j]))
  }
  evpars=rbind(evpars,data.frame(day=cvirus$fecha[max(windows[[i]])],
                              value=post.par[[j+1]][,1],
                              param="outpred"))

  cat("i=",i,"/",nw," - ")
}
save(evpars,file="evpars-sardegna.RData")

Evolution of parameters (using data since the beginning)

load(file="evpars-sardegna.RData")
ggplotly(ggplot(evpars[evpars$param!="outpred",], aes(x=day, y=value, colour=param)) + geom_smooth()+
  geom_vline(xintercept=as.Date("2020-03-10"),linetype="dashed",color = "red", size=2)+
  geom_hline(yintercept=0)+
  xlab("Days")+  ylab("Posterior mean and 95% C.I."))

Reaching a peak means \(\alpha_2<0\) and \(\beta_2<0\), while disappear of the covid means all parameters less than 0.

Out of sample prediction

outpred=evpars[evpars$param=="outpred",]
outpred=aggregate(value~day,outpred,
                             function(xx) c(quantile(xx,0.025),mean(xx),quantile(xx,0.975)))
outpred=data.frame(do.call(cbind,outpred))
colnames(outpred)=c("day","inf","mean","sup")
outpred$predday=sort(unique(evpars$day))+1
ggplot(outpred, aes(x=predday, y=mean)) + geom_line() +
geom_point(data=cvirus[cvirus$fecha%in%outpred$day,],
           aes(x=fecha,y=casos),color="red")+
  geom_errorbar(aes(ymin=inf, ymax=sup), width=.2,position=position_dodge(0.05))+
  xlab("Days")+ ylab("Posterior outofsample mean and 95% C.I.")