not only Stefano Cabras … also thanks to statisticians at uc3m.es

Final Prediction for next two days

This is the prediction for next two days updated at the time of this report (above)

Data

Here is the DataBase of only confirmed cases at OMS at the end of the Day.

rm(list=ls())
url <- "<https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv>"
cvirus <- read.table(url, sep = ",", header = T)
cvirus$Lat <- cvirus$Long <- NULL
cvirus$Province.State <- NULL
 
cvirus <- melt(cvirus, id.vars = "Country.Region")
 
colnames(cvirus) <- c("pais", "fecha", "casos")
cvirus <- cvirus[cvirus$pais%in%c("Spain"),]
cvirus$fecha <- as.Date(as.character(cvirus$fecha), format = "X%m.%d.%y")
cvirus$casos=as.numeric(cvirus$casos)
cvirus=cvirus[order(cvirus$fecha),]

y=cvirus$casos[-1]
ya=cvirus$casos[-length(cvirus)]
n=nrow(cvirus)-1
plot(1:n,y,type="l",xlab=paste("Days from",cvirus$fecha[1]),
     ylab="Confirmed Cases at the end of the day")

Model for marginal counts cases

Let represents the number of cases at time (days) where the day after 2020-01-22.

The fitted model is

Parameter interpretation:

• is the effect of under counted assuming that at least at the beginning all Infected cases 7 days before has been reported.

The non Bayesian version of this model is here:

rstan_options(auto_write = TRUE)
Sys.setenv(LOCAL_CPPFLAGS = '-march=corei7 -mtune=corei7')
options(mc.cores = parallel::detectCores())

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

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

transformed parameters {
  vector[n] llambda;
  llambda[1]=fakecero;
  for(t in 2:7) llambda[t]=omega+alpha*log(1+ya[t])+beta*llambda[t-1];
  for(t in 8:n) llambda[t]=omega+alpha*log(1+ya[t])+beta*llambda[t-1]+gamma*log(1+ya[t-7]);
}

model {
  // Priors
  fakecero  ~ normal(-99,0.001);
  alpha ~ normal(0,10);
  beta ~ normal(0,10);
  omega ~ normal(0,10);
  gamma ~ normal(0,10);
  // Likelihood
  y ~ poisson(exp(llambda));

}

generated quantities {
  int<lower=0> yp[n];
  int<lower=0> ypf1;
  real llambdaf1;
  int<lower=0> ypf2;
  real llambdaf2;
  yp[1]=0;
  for(t in 2:n) yp[t] = poisson_rng(exp(llambda[t])); //y values predicted by the model
  llambdaf1=omega+alpha*log(1+y[n])+beta*llambda[n];
  ypf1=poisson_rng(exp(llambdaf1));
  llambdaf2=omega+alpha*log(1+ypf1)+beta*llambdaf1;
  ypf2=poisson_rng(exp(llambdaf2));
}
"

Hamiltonian MC is used for obtaining the posterior:

fit = stan(model_code=mod.cov,
           data=list(y=y,
                     ya=ya,
                     n=n), 
           iter=1000,chains = 4,
#           control = list(adapt_delta = 0.99,max_treedepth=50),
           seed = 11)
save(fit,cvirus,file="stanmod-cov.RData")

Predictive checks

load(file="stanmod-cov.RData")
post.par=extract(fit,pars=c("alpha","beta","omega","gamma"))
post.pred=extract(fit,pars="yp")$yp
pred.mean=log(apply(post.pred,2,mean),10)
pred.ic.inf=log(apply(post.pred,2,quantile,p=0.025),10)
pred.ic.sup=log(apply(post.pred,2,quantile,p=0.975),10)
plot(1:n,log(y,10),xlab=paste("Days from",cvirus$fecha[1]))
lines(1:n,pred.mean)
segments(1:n,pred.ic.inf,1:n,pred.ic.sup,col=2)
points(1:n,pred.mean,pch=19)
legend("topleft",c("Observed","Posterior Predicted","95% C.I."),
       col=c(1,1,2),pch=c(21,19,22),bty="n")

Posterior Parameters