---
title: "Supplementary material to the paper: Spatio-temporal small area surveillance of the Covid-19 pandemics."
author: "Martinez-Beneito, Miguel A.; Mateu, Jorge and Botella-Rocamora, Paloma"
date: "October 25th, 2020"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

```{r, include=FALSE}
knitr::opts_chunk$set(cache = TRUE)

# This parameter controls if WinBUGS should be run in order to compute the fit of the spatio temporal models. These models are expensive to run so it could be convenient to save them once run and not running them again.
runWB <- FALSE
```

# Preliminaries

```{r initial setup, cache = FALSE, warning = FALSE, message = FALSE}
#library(mapview)
library(rgdal)
library(sf)
library(spdep)
library(corrplot)
library(R2WinBUGS)
library(RColorBrewer)
library(EpiEstim)
library(splines)

if (!require(pbugs)) {
  if (!require(devtools)) {
    install.packages("devtools")
    devtools::install_github("fisabio/pbugs")
  } else {
    install_github("fisabio/pbugs")
  }
}
```

# Model fitting, Castilla y León (CyL)
```{r load data for CyL}
load("DataCyL.RData")
```

## Temporal analysis for the whole CyL

### WinBUGS temporal model for the whole region
Some posterior estimates of this model will be used to initialize some parameters of the spatio-temporal model. This is done to speed up the convergence of that model.
```{r WinBUGS temporal model for the whole region}
#Splines temporal model for the whole region of study
TemporalModel.WB<-function(){
  for(j in 1:nDays){
    O[j]~dpois(mu[j])
    log(mu[j])<-log(pob.tot)+gamma[DoW[j]]+inprod2(design[j,],beta[])
  }
  for(k in 1:(nCoefs)){beta[k]~dflat()}
  gamma[1]<-0
  for(k in 2:7){gamma[k]~dflat()}
}
```

### WinBUGS call
```{r Data preparation for WinBUGS call CyL, temporal model}
nWeeks <- nrow(CasesCyL)%/%7

# design matrix for the B-splines linear basis
design.CyL <- ns(x = 1:nrow(CasesCyL), knots = 14*(1:((nWeeks)/2)), intercept = TRUE)

data <- list(O = as.vector(apply(CasesCyL,1,sum)), 
              nCoefs = ncol(design.CyL),
              nDays = nrow(CasesCyL),
              DoW = 1+(0:(nrow(CasesCyL)-1)%%7),
              pob.tot = sum(PopCyL),
              design = design.CyL)

inits <- function() {
  list(beta = rnorm(ncol(design.CyL), 0, 5), 
       gamma = c(NA, runif(6, -1, 1))
       )
}

# Parameters to be saved
param <- c("gamma","beta")

niters <- 5000
nburns <- 2000
```

```{r WinBUGS call CyL, temporal model}
set.seed(1)
resulTempCyL <- pbugs(data = data, 
                      inits = inits, 
                      parameters.to.save = param, 
                      model = TemporalModel.WB, 
                      n.iter = niters, 
                      n.burnin = nburns, 
                      n.chains = 3, 
                      DIC = F, 
                      debug = F, 
                      bugs.seed = 1)
```

## Spatio-temporal analysis for the whole CyL
### WinBUGS spatio-temporal model
```{r WinBUGS spatio-temporal model}
# Splines spatio-temporal model 
SpatioTemporalModel.WB <- function(){
  for(i in 1:nHD){
    for(j in 1:nDays){
      O[i,j]~dpois(mu[i,j])
      log(mu[i,j])<-log(pob[i])+gamma[DoW[j]]+inprod2(design[j,],beta[i,])+sd.epsi*epsi[i,j]
      epsi[i,j]~dnorm(0,1)
    }
  }
  
  for(i in 1:nHD){
    for(j in 1:(nBasis)){
      beta[i,j]<-beta.mean[j]+beta.var[i,j]
      # Leroux process for the coefficients of each spline node
      beta.var[i,j]~dnorm(mean.beta[i,j],prec.beta[i,j])
      prec.beta[i,j]<-(1-rho.beta+rho.beta*nneigh[i])/(sd.beta[j]*sd.beta[j])
    }
  }
  
  #keep on with the Leroux process
  for(j in 1:(nBasis)){
    mean.beta[1,j]<-(rho.beta/(1-rho.beta+rho.beta*nneigh[1]))*sum(beta.map[index[1]:index[2],j])
	  for(i in 2:(nHD)){
			mean.beta[i,j]<-(rho.beta/(1-rho.beta+rho.beta*nneigh[i]))*sum(beta.map[(index[i]+1):index[i+1],j])
	  }
    for(k in 1:nneigh.tot){
		  beta.map[k,j]<-beta.var[map[k],j]
    }
  }

  #Prior distributions
  for(j in 1:(nBasis)){
    beta.mean[j]~dflat()
  }
  
  gamma[1]<-0
  for(k in 2:7){gamma[k]~dflat()}
  
  for(j in 1:(nBasis)){
    sd.beta[j]~dunif(0,10)
  }
  rho.beta~dunif(0,1)
  sd.epsi~dunif(0,10)
}
```

### WinBUGS call
```{r Data preparation for WinBUGS call CyL, spatio-temporal model}
# design matrix for the B-splines linear basis
design.CyL <- ns(x = 1:nrow(CasesCyL), knots = 14*(1:(nWeeks/2)), intercept = TRUE)

data <- list(O = t(CasesCyL),
              nBasis = dim(design.CyL)[2],
              nDays = nrow(CasesCyL),
              DoW = 1+(0:(nrow(CasesCyL)-1)%%7),
              pob = PopCyL,
              nHD = ncol(CasesCyL),
              design = design.CyL,
              nneigh.tot = length(NeighsCyL$adj), 
              map = NeighsCyL$adj, 
              nneigh = NeighsCyL$num, 
              index = c(1,cumsum(NeighsCyL$num))
              )

inits <- function() {
  list(beta.mean = rnorm(dim(design.CyL)[2], resulTempCyL$mean$beta,1), 
       gamma = c(NA, runif(6, -1, 1)),
       sd.beta = runif(dim(design.CyL)[2], 0.5, 2),
       sd.epsi = runif(1, 0.5, 2))
}

# Parameters to save
param <- c("sd.epsi", "sd.beta", "rho.beta", "beta.mean", "gamma", "beta")

niters <- 5000
nburns <- 2000
```

```{r WinBUGS call CyL spatio-temporal model, eval = runWB}
set.seed(1)
resulSTCyL <- pbugs(data = data, 
                    inits = inits, 
                    parameters.to.save = param, 
                    model = SpatioTemporalModel.WB, 
                    n.iter = niters, 
                    n.burnin = nburns, 
                    n.chains = 5, 
                    DIC = F, 
                    debug = F, 
                    bugs.seed = 1, 
                    slice = TRUE)
#save(resulSTCyL,file = "ResulCyL.Rdata")
```

```{r, eval = !runWB, include = FALSE}
load(file = "ResulCyL.Rdata")
```

# Model fitting, Comunidad Valenciana (CV)
Data for this region are not provided, so you would possible rather to skip all the analysis code corresponding to this region.

```{r load data for CV}
load("DataCV.RData")
```

## Temporal analysis for the whole CV

### WinBUGS call
```{r Data preparation for WinBUGS call CV, temporal model}
nWeeks <- nrow(CasesCV)%/%7

# design matrix for the B-splines linear basis
design.CV <- ns(x = 1:nrow(CasesCV), knots = 14*(1:((nWeeks)/2)), intercept = TRUE)

data <- list(O = as.vector(apply(CasesCV,1,sum)), 
              nCoefs = ncol(design.CV),
              nDays = nrow(CasesCV),
              DoW = 1+(0:(nrow(CasesCV)-1)%%7),
              pob.tot = sum(PopCV),
              design = design.CV)

inits <- function() {
  list(beta = rnorm(ncol(design.CV), 0, 5), 
       gamma = c(NA, runif(6, -1, 1))
       )
}

# Parameters to save
param <- c("gamma", "beta")
niters <- 5000
nburns <- 2000
```

```{r WinBUGS call CV, temporal model}
set.seed(1)
resulTempCV <- pbugs(data = data, 
                     inits = inits, 
                     parameters.to.save = param, 
                     model = TemporalModel.WB, 
                     n.iter = niters, 
                     n.burnin = nburns, 
                     n.chains = 3, 
                     DIC = F, 
                     debug = F, 
                     bugs.seed = 1)
```

## Spatio-temporal analysis for the whole CV

### WinBUGS call
```{r Data preparation for WinBUGS call CV, spatio-temporal model}
# design matrix for the B-splines linear basis
design.CV <- ns(x = 1:nrow(CasesCV), knots = 14*(1:(nWeeks/2)), intercept = TRUE)

data <- list(O = t(CasesCV),
              nBasis = dim(design.CV)[2],
              nDays = nrow(CasesCV),
              DoW = 1+(0:(nrow(CasesCV)-1)%%7),
              pob = PopCV,
              nHD = ncol(CasesCV),
              design = design.CV,
              nneigh.tot = length(NeighsCV$adj), 
              map = NeighsCV$adj, 
              nneigh = NeighsCV$num, 
              index = c(1, cumsum(NeighsCV$num))
              )

inits <- function() {
  list(beta.mean = rnorm(dim(design.CV)[2], resulTempCV$mean$beta, 1), 
       gamma = c(NA, runif(6, -1, 1)),
       sd.beta = runif(dim(design.CV)[2], 0.5, 2),
       sd.epsi = runif(1), 0.5, 2)
}
param <- c("sd.epsi", "sd.beta", "rho.beta", "beta.mean", "gamma", "beta")

niters <- 5000
nburns <- 2000
```

```{r WinBUGS call CV, spatio-temporal model, eval = runWB}
set.seed(1)
resulSTCV <- pbugs(data = data, 
                   inits = inits, 
                   parameters.to.save = param, 
                   model = SpatioTemporalModel.WB, 
                   n.iter = niters, 
                   n.burnin = nburns, 
                   n.chains = 5, 
                   DIC = F, 
                   debug = F, 
                   bugs.seed = 1,
                   slice = TRUE)
#save(resulSTCV,file = "Paper/ResulCV.Rdata")
```

```{r, eval = !runWB, include = FALSE}
load(file = "ResulCV.Rdata")
```

# Results summaries

Calculation of rates and $R_t$s for the spatial units and whole regions to be represented later

## Summaries Castilla y León
```{r}
# rates
#######
rates.sim.CyL <- exp(apply(resulSTCyL$sims.list$beta, c(1, 2), function(x){design.CyL%*%matrix(x,ncol=1)}))*100000
# Effect of the unstructured random effect variability over the rates 
for(i in 1:(dim(rates.sim.CyL)[2])){
  rates.sim.CyL[,i , ] <- rates.sim.CyL[,i , ] * exp(((resulSTCyL$sims.list$sd.epsi[i])^2)/2)
}

# municipal daily rates (for each simulation) including the cyclic term
rates.sim.CyL.cyclic <- rates.sim.CyL
aux <- exp(cbind(rep(0,dim(resulSTCyL$sims.list$gamma)[1]),resulSTCV$sims.list$gamma))
for(i in 1:(dim(rates.sim.CyL.cyclic)[1])){
  rates.sim.CyL.cyclic[i,,] <- rates.sim.CyL[i,,] * aux[,1+(i-1)%%7]
}
# municipal daily rates (for each simulation) excluding the cyclic term
aux <- exp(apply(resulSTCyL$sims.list$gamma, 1, sum)/7)
rates.sim.CyL <- apply(rates.sim.CyL, c(1, 3), function(x){x*aux})
rates.sim.CyL <- aperm(a = rates.sim.CyL, perm = c(2, 1, 3))
# municipal daily rates summaries excluding the cyclic term
rates.CyL <- apply(rates.sim.CyL, c(1, 3), mean)
rates.CyL.CIL <- apply(rates.sim.CyL, c(1, 3), quantile, 0.025)
rates.CyL.CIU <- apply(rates.sim.CyL, c(1, 3), quantile, 0.975)
# daily rates for all the region and simulations, including the cyclic term 
rate.sim.CyL.cyclic <- apply(rates.sim.CyL.cyclic, 2, 
                             function(x){x%*%matrix(PopCyL, ncol=1)/sum(PopCyL)})
# daily rates for all the region and simulations, excluding the cyclic term
rate.sim.CyL <- apply(rates.sim.CyL, 2, 
                      function(x){x%*%matrix(PopCyL, ncol=1)/sum(PopCyL)})
# mean daily rates for all the region, including the cyclic term 
rate.CyL.cyclic <- apply(rate.sim.CyL.cyclic,1,mean)
# mean daily rates for all the region, excluding the cyclic term 
rate.CyL <- apply(rate.sim.CyL, 1, mean)
# ICs for daily rates for all the region, including the cyclic term 
rate.CyL.CIL.cyclic <- apply(rate.sim.CyL.cyclic, 1, quantile, 0.025)
rate.CyL.CIU.cyclic <- apply(rate.sim.CyL.cyclic, 1, quantile, 0.975)
# ICs for daily rates for all the region, excluding the cyclic term 
rate.CyL.CIL <- apply(rate.sim.CyL, 1, quantile, 0.025)
rate.CyL.CIU <- apply(rate.sim.CyL, 1, quantile, 0.975)

# Rt
#####
# We take into account "just" 25 days backwards for the Rts calcution since these cover 0.9999 of the probability of the serial interval
ws <- rev(discr_si(1:25, 4.7, 2.9))
# Expected cases for each day according to the observed cases in the previous days and Rt=1
rates.sim.CyL.E <- array(NA, dim = dim(rates.sim.CyL))
# We calculate Rt from the 5th day of analysis since for the first 4 days most of the period of the serial interval is unknown
for(i in 5:25){
    rates.sim.CyL.E[i,,] <- apply(rates.sim.CyL[1:(i-1), , ], c(2, 3),
                                  function(x){sum(x*ws[1:(i-1)])/sum(ws[1:(i-1)])})
}
for(i in 26:(dim(CasesCyL)[1])){
    rates.sim.CyL.E[i,,] <- apply(rates.sim.CyL[(i-25):(i-1), , ], c(2, 3),
                                  function(x){sum(x*ws)})
}
RTs.CyL <- apply(rates.sim.CyL/rates.sim.CyL.E, c(1, 3), mean)
RTs.CyL.CIL <- rbind(matrix(NA, ncol = (dim(CasesCyL)[2]), nrow = 4), apply((rates.sim.CyL/rates.sim.CyL.E)[-(1:4), , ], c(1, 3), quantile, 0.025))
RTs.CyL.CIU <- rbind(matrix(NA, ncol = (dim(CasesCyL)[2]), nrow = 4), apply((rates.sim.CyL/rates.sim.CyL.E)[-(1:4), , ], c(1, 3), quantile, 0.975))

# Expected cases for each day according to the observed cases in the previous days and Rt=1, whole region
rate.sim.CyL.E <- apply(rates.sim.CyL.E, 2, 
                        function(x){x%*%matrix(PopCyL, ncol = 1)/sum(PopCyL)})
RT.CyL <- apply(rate.sim.CyL/rate.sim.CyL.E, 1, mean)
RT.CyL.CIL <- c(rep(NA, 4), 
                apply((rate.sim.CyL/rate.sim.CyL.E)[-(1:4), ], 1, quantile, 0.025))
RT.CyL.CIU <- c(rep(NA, 4), 
                apply((rate.sim.CyL/rate.sim.CyL.E)[-(1:4), ], 1, quantile, 0.975))
```

## Summaries Comunidad Valenciana
```{r}
# rates
#######
rates.sim.CV <- exp(apply(resulSTCV$sims.list$beta, c(1, 2), 
                          function(x){design.CV%*%matrix(x, ncol = 1)}))*100000
# Effect of the unstructured random effect variability over the rates 
for(i in 1:(dim(rates.sim.CV)[2])){
  rates.sim.CV[,i , ] <- rates.sim.CV[,i , ] * exp(((resulSTCV$sims.list$sd.epsi[i])^2)/2)
}

# municipal daily rates (for each simulation) including the cyclic term
rates.sim.CV.cyclic <- rates.sim.CV
aux <- exp(cbind(rep(0, dim(resulSTCV$sims.list$gamma)[1]), resulSTCV$sims.list$gamma))
for(i in 1:(dim(rates.sim.CV.cyclic)[1])){
  rates.sim.CV.cyclic[i, , ] <- rates.sim.CV[i, , ] * aux[, 1+(i-1)%%7]
}
# municipal daily rates (for each simulation) excluding the cyclic term
aux <- exp(apply(resulSTCV$sims.list$gamma, 1, sum)/7)
rates.sim.CV <- apply(rates.sim.CV, c(1, 3), function(x){x*aux})
rates.sim.CV <- aperm(a = rates.sim.CV, perm = c(2, 1, 3))
# municipal daily rates summaries excluding the cyclic term
rates.CV <- apply(rates.sim.CV, c(1, 3), mean)
rates.CV.CIL <- apply(rates.sim.CV, c(1, 3), quantile, 0.025)
rates.CV.CIU <- apply(rates.sim.CV, c(1, 3), quantile, 0.975)
# daily rates for all the region and simulations, including the cyclic term 
rate.sim.CV.cyclic <- apply(rates.sim.CV.cyclic, 2, 
                            function(x){x%*%matrix(PopCV, ncol = 1)/sum(PopCV)})
# daily rates for all the region and simulations, excluding the cyclic term
rate.sim.CV <- apply(rates.sim.CV, 2, 
                     function(x){x%*%matrix(PopCV, ncol = 1)/sum(PopCV)})
# mean daily rates for all the region, including the cyclic term 
rate.CV.cyclic <- apply(rate.sim.CV.cyclic, 1, mean)
# mean daily rates for all the region, excluding the cyclic term 
rate.CV <- apply(rate.sim.CV, 1, mean)
# ICs for daily rates for all the region, including the cyclic term 
rate.CV.CIL.cyclic <- apply(rate.sim.CV.cyclic, 1, quantile, 0.025)
rate.CV.CIU.cyclic <- apply(rate.sim.CV.cyclic, 1, quantile, 0.975)
# ICs for daily rates for all the region, excluding the cyclic term 
rate.CV.CIL <- apply(rate.sim.CV, 1, quantile, 0.025)
rate.CV.CIU <- apply(rate.sim.CV, 1, quantile, 0.975)

# Rt
#####
# Expected cases for each day according to the observed cases in the previous days and Rt=1
rates.sim.CV.E <- array(NA, dim = dim(rates.sim.CV))
# We calculate Rt from the 5th day of analysis since for the first 4 days most of the period of the serial interval is unknown
for(i in 5:25){
    rates.sim.CV.E[i,,] <- apply(rates.sim.CV[1:(i-1), , ], c(2, 3),
                                 function(x){sum(x*ws[1:(i-1)])/sum(ws[1:(i-1)])})
}
for(i in 26:(dim(CasesCV)[1])){
    rates.sim.CV.E[i,,] <- apply(rates.sim.CV[(i-25):(i-1), , ], c(2, 3), 
                                 function(x){sum(x*ws)})
}
RTs.CV <- apply(rates.sim.CV/rates.sim.CV.E, c(1, 3), mean)
RTs.CV.CIL <- rbind(matrix(NA, ncol = (dim(CasesCV)[2]), nrow = 4), apply((rates.sim.CV/rates.sim.CV.E)[-(1:4), , ], c(1, 3), quantile, 0.025))
RTs.CV.CIU <- rbind(matrix(NA, ncol = (dim(CasesCV)[2]), nrow = 4), apply((rates.sim.CV/rates.sim.CV.E)[-(1:4), , ], c(1, 3), quantile, 0.975))

# Expected cases for each day according to the observed cases in the previous days and Rt=1, whole region
rate.sim.CV.E <- apply(rates.sim.CV.E, 2, 
                       function(x){x%*%matrix(PopCV, ncol = 1)/sum(PopCV)})
RT.CV <- apply(rate.sim.CV/rate.sim.CV.E, 1, mean)
RT.CV.CIL <- c(rep(NA, 4), 
               apply((rate.sim.CV/rate.sim.CV.E)[-(1:4), ], 1, quantile, 0.025))
RT.CV.CIU <- c(rep(NA, 4), 
               apply((rate.sim.CV/rate.sim.CV.E)[-(1:4), ], 1, quantile, 0.975))
```

## Figure 1
```{r}
FirstOfMonth <- which(names(apply(CasesCyL, 1, sum))%in%c("2020-04-01", "2020-05-01", "2020-06-01", "2020-07-01", "2020-08-01", "2020-09-01", "2020-10-01"))
# adjacency matrix for CV
Veci <- poly2nb(CartoCV)


#pdf(file = "Figure1.pdf")
par(mfrow = c(2, 2), bty = "l")
plot(sum(PopCV)*rate.CV/100000, 
     xlab = "Day", 
     ylab = "Observed cases", 
     main = "Time trend, whole CV", 
     ylim = range(apply(CasesCV, 1, sum)), 
     axes = F, 
     type = "n")
box()
axis(side = 2, at = 200*(0:4), labels = c("0", NA, "400", NA, "800"))
axis(side = 1, at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
abline(v = FirstOfMonth, col = "green", lty = 2)
points((1:223), apply(CasesCV, 1, sum), col = 2, cex= 0.5)
polygon(x = c(1:223, rev(1:223)),
        y = sum(PopCV)*c(rate.CV.CIL.cyclic, rev(rate.CV.CIU.cyclic))/100000, 
        col = rgb(0, 0, 1, 0.2), 
        border = FALSE)
lines(sum(PopCV)*rate.CV.cyclic/100000, col = rgb(0, 0, 1, 1))
polygon(x = c(1:223, rev(1:223)),
        y = sum(PopCV)*c(rate.CV.CIL, rev(rate.CV.CIU))/100000,
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
lines(sum(PopCV)*rate.CV/100000)

plot(sum(PopCyL)*rate.CyL/100000, 
     xlab = "Day", 
     ylab = "Observed cases", 
     main = "Time trend, whole CyL", 
     ylim = range(apply(CasesCyL, 1, sum)), 
     axes = F, 
     type = "n")
box()
axis(side = 2, at = 500*(0:6), labels = c("0", NA, "1000", NA, "2000", NA, "3000"))
axis(side = 1,at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
abline(v = FirstOfMonth, col = "green", lty = 2)
points((1:223), apply(CasesCyL, 1, sum), col = 2, cex = 0.5)
polygon(x = c(1:223, rev(1:223)),
        y = sum(PopCyL)*c(rate.CyL.CIL.cyclic, rev(rate.CyL.CIU.cyclic))/100000, 
        col = rgb(0, 0, 1, 0.2),
        border = FALSE)
lines(sum(PopCyL)*rate.CyL.cyclic/100000, col = rgb(0, 0, 1, 1))
polygon(x = c(1:223,rev(1:223)),
        y = sum(PopCyL)*c(rate.CyL.CIL, rev(rate.CyL.CIU))/100000, 
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
lines(sum(PopCyL)*rate.CyL/100000)

plot(RT.CV, xlab = "Day", 
     ylab = expression(R[t]), 
     main="Instantaneous reproduction number", 
     axes = F, 
     type = "n", 
     ylim =c(0.6, 2.4))
box()
axis(side = 2, at = 0.2*(4:12))
axis(side = 1, at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
abline(h = 1, lty = 2)
abline(v = FirstOfMonth, col = "green", lty = 2)
polygon(x = c(1:223, rev(1:223)),
        y = c(RT.CV.CIL,rev(RT.CV.CIU)), 
        col = rgb(0, 0, 1, 0.2), 
        border = FALSE)
lines(RT.CV, col = rgb(0, 0, 1, 1))
polygon(x=c(1:223, rev(1:223)),
        y = c(RT.CyL.CIL, rev(RT.CyL.CIU)), 
        col = rgb(1, 0, 0, 0.2), 
        border = FALSE)
lines(RT.CyL, col = rgb(1, 0, 0, 1))
legend(x="topright", 
       lty = 1, 
       col = c(rgb(0, 0, 1, 1), rgb(1, 0, 0, 1)), 
       legend = c("CV", "CyL"), 
       inset = 0.03, 
       bg = "white")

#Valencia
i <- which(names(PopCV) == "46250")
PosVal <- i
plot(rates.CV[, i], 
     type = "l", 
     axes = F, 
     xlab = "Day", 
     ylab = "Rate", 
     main = "Incidence rate")
box()
axis(side = 2)
axis(side = 1, at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
abline(v = FirstOfMonth, col = "green", lty = 2)
#Alicante
i <- which(names(PopCV) == "03014")
PosAli <- i
lines(rates.CV[,i],col = 2)
# Massalfassar
i <- which(names(PopCV) == names(which.min(PopCV[Veci[[PosVal]]])))
PosMas <- i
lines(rates.CV[, i], lty = 2)
#Aigues
i <- which(names(PopCV) == names(which.min(PopCV[Veci[[PosAli]]])))
PosAig <- i
lines(rates.CV[, i], lty = 2, col = 2)
legend(x = "topleft", 
       lty = c(1, 1, 2, 2), 
       col = c("black", "red", "black", "red"), 
       legend = c("Valencia", "Alicante", "Massalfassar", "Aigües"), 
       inset = 0.03, 
       bg = "white", 
       cex = 0.8)
#dev.off()
```

## Figure 2
```{r, cache=FALSE}
pdf(file = "Figure2.pdf")

par(mfrow = c(4, 2), bty = "l")
#Valencia
plot(PopCV[PosVal] * rates.CV[, PosVal]/100000, 
     type = "l", 
     axes = F, 
     xlab = "Day", 
     ylab = "Observed cases", 
     main = "Time trend, Valencia", 
     ylim = range(CasesCV[, PosVal]))
box()
axis(side = 2)
axis(side = 1,at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
abline(v = FirstOfMonth, col = "green", lty = 2)
polygon(x = c(1:223, rev(1:223)),
        y = PopCV[PosVal]*c(rates.CV.CIL[, PosVal], rev(rates.CV.CIU[, PosVal]))/100000,
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
points(CasesCV[, PosVal], col = 2, cex = 0.5)

plot(RTs.CV[, PosVal], 
     type = "l", 
     axes = F, 
     xlab = "Day", 
     ylab = expression(R[t]), 
     main = expression(paste(R[t], " Valencia")))
box()
axis(side = 2, at = 0.2 * (4:12))
axis(side = 1,at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
polygon(x = c(1:223, rev(1:223)),
        y = c(RTs.CV.CIL[, PosVal], rev(RTs.CV.CIU[, PosVal])), 
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
abline(h = 1, lty = 2)
abline(v = FirstOfMonth, col = "green", lty = 2)

# Alicante
plot(PopCV[PosAli] * rates.CV[, PosAli]/100000, 
     type = "l", 
     axes = F, 
     xlab = "Day",
     ylab = "Observed cases", 
     main = "Time trend, Alicante",
     ylim = range(CasesCV[, PosAli]))
box()
axis(side = 2)
axis(side = 1,at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
abline(v = FirstOfMonth, col = "green", lty = 2)
polygon(x = c(1:223, rev(1:223)),
        y = PopCV[PosAli] * c(rates.CV.CIL[, PosAli], rev(rates.CV.CIU[, PosAli]))/100000, 
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
points(CasesCV[, PosAli], col = 2, cex = 0.5)

plot(RTs.CV[, PosAli], 
     type = "l", 
     axes = F, 
     xlab = "Day", 
     ylab = expression(R[t]), 
     main=expression(paste(R[t], " Alicante")))
box()
axis(side = 2, at = 0.2*(4:12))
axis(side = 1,at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
polygon(x = c(1:223, rev(1:223)),
        y = c(RTs.CV.CIL[, PosAli], rev(RTs.CV.CIU[, PosAli])), 
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
abline(h = 1, lty = 2)
abline(v = FirstOfMonth, col = "green", lty = 2)

# Massalfassar
plot(PopCV[PosMas] * rates.CV[, PosMas]/100000, 
     type = "l", 
     axes = F, 
     xlab = "Day", 
     ylab = "Observed cases", 
     main = "Time trend, Massalfassar", 
     ylim = range(CasesCV[, PosMas]))
box()
axis(side = 2)
axis(side = 1, at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
abline(v = FirstOfMonth, col = "green", lty = 2)
polygon(x = c(1:223, rev(1:223)),
        y = PopCV[PosMas] * c(rates.CV.CIL[,PosMas], rev(rates.CV.CIU[, PosMas]))/100000, 
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
points(CasesCV[, PosMas], col = 2, cex = 0.5)

plot(RTs.CV[, PosMas], 
     type = "l", 
     axes = F, 
     xlab = "Day", 
     ylab = expression(R[t]), 
     main = expression(paste(R[t], " Massalfassar")))
box()
axis(side = 2, at = 0.2 * (4:12))
axis(side = 1,at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
polygon(x = c(1:223,rev(1:223)),
        y = c(RTs.CV.CIL[, PosMas], rev(RTs.CV.CIU[, PosMas])), 
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
abline(h = 1, lty = 2)
abline(v = FirstOfMonth, col = "green", lty = 2)

# Aigues
plot(PopCV[PosAig] * rates.CV[, PosAig]/100000, 
     type = "l", 
     axes = F, 
     xlab = "Day", 
     ylab = "Observed cases", 
     main = "Time trend, Aigües", 
     ylim = range(CasesCV[, PosAig]))
box()
axis(side = 2)
axis(side = 1, at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
abline(v = FirstOfMonth, col = "green", lty = 2)
polygon(x= c(1:223, rev(1:223)),
        y = PopCV[PosAig] * c(rates.CV.CIL[, PosAig], rev(rates.CV.CIU[, PosAig]))/100000, 
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
points(CasesCV[, PosAig], col = 2, cex = 0.5)

plot(RTs.CV[, PosAig], 
     type = "l", 
     axes = F, 
     xlab = "Day", 
     ylab = expression(R[t]), 
     main = expression(paste(R[t], " Aigües")))
box()
axis(side = 2, at = 0.2 * (4:12))
axis(side = 1, at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
polygon(x = c(1:223, rev(1:223)),
        y = c(RTs.CV.CIL[, PosAig], rev(RTs.CV.CIU[, PosAig])), 
        col = rgb(0, 0, 0, 0.2), 
        border = FALSE)
abline(h = 1, lty = 2)
abline(v = FirstOfMonth, col = "green", lty = 2)

dev.off()
knitr::include_graphics(path = "Figure2.pdf")
```

## Figure 3
```{r}
#pdf(file = "Figure3.pdf",width=13)
par(mfrow = c(1, 2), bty = "l")
betasCyL <- resulSTCyL$mean$beta[, c(16, 1:15, 17)]
betasCV <- resulSTCV$mean$beta[,c(16, 1:15, 17)]
colnames(betasCyL) <- c("3/6", "3/20", "4/2", "4/16", "4/30", "5/14", "5/28", "6/11", "6/24", "7/9", "7/23", "8/6", "8/20", "9/3", "9/17", "9/30", "10/14")
colnames(betasCV) <- colnames(betasCyL)
corrplot(cor(betasCV), title = "Comunidad Valenciana, correlations.", mar = c(0, 0, 3, 0))
corrplot(cor(betasCyL), title = "Castilla y León, correlations.", mar = c(0, 0, 3, 0))
#dev.off()
```

## Figure 4
```{r, warning = FALSE}
#pdf(file = "Figure4.pdf",width=10)
par(mfrow = c(1, 2))
pal.col <- brewer.pal(n = 9,name = "RdYlGn")

cutsRates <- c(10, 25, 75, 125)
cutsRTs <- c(1, 1.1, 1.5, 2)  

LastRates <- 7 * rates.CV[223, ]
LastRTs <- RTs.CV[223, ]

plot(LastRates, LastRTs, type = "n",
     xlim = c(0, max(7 * rates.CV[223,]) * 1.1), 
     ylim = c(0, max(RTs.CV[223, ]) * 1.1),
     xlab = "Incidence rate for the last 7 days",
     ylab = expression(R[t]), 
     main = "Risk levels per municipality")

auxCutsRates <- c(0, cutsRates, max(LastRates) * 1.1)
auxCutsRts <- c(0, cutsRTs, max(LastRTs) * 1.1)

for(i in 1:5){
  for(j in 1:5){
polygon(x = c(auxCutsRates[i], auxCutsRates[i+1], auxCutsRates[i+1], auxCutsRates[i]),
        y = c(auxCutsRts[j], auxCutsRts[j], auxCutsRts[j+1], auxCutsRts[j+1]), 
        col = pal.col[11-i-j], 
        border = FALSE)
polygon(x = c(auxCutsRates[i+1], auxCutsRates[i+2], auxCutsRates[i+2], auxCutsRates[i+1]),
        y = c(auxCutsRts[j], auxCutsRts[j], auxCutsRts[j+1], auxCutsRts[j+1]), 
        col = pal.col[10-i-j], 
        border = FALSE)
polygon(x = c(auxCutsRates[i+2], auxCutsRates[i+3], auxCutsRates[i+3], auxCutsRates[i+2]),
        y = c(auxCutsRts[j], auxCutsRts[j], auxCutsRts[j+1], auxCutsRts[j+1]), 
        col = pal.col[9-i-j], 
        border = FALSE)
polygon(x = c(auxCutsRates[i+3], auxCutsRates[i+4], auxCutsRates[i+4], auxCutsRates[i+3]),
        y = c(auxCutsRts[j], auxCutsRts[j], auxCutsRts[j+1], auxCutsRts[j+1]), 
        col = pal.col[8-i-j], 
        border = FALSE)
polygon(x = c(auxCutsRates[i+4], auxCutsRates[i+5], auxCutsRates[i+5], auxCutsRates[i+4]),
        y = c(auxCutsRts[j], auxCutsRts[j], auxCutsRts[j+1], auxCutsRts[j+1]),
        col=pal.col[7-i-j], 
        border = FALSE)
  }
}

points(LastRates, LastRTs, pch = 19, cex = 0.6)

LastRates.q <- sapply(LastRates, function(x){sum(x > c(10, 25, 50, 125))})
LastRTs.q <- sapply(LastRTs, function(x){sum(x > c(1, 1.1, 1.5, 2))})
grupo.q <- LastRates.q + LastRTs.q + 1
plot(CartoCV, col = pal.col[10-grupo.q], main = "Risk levels map")
#dev.off()
```

## Figure 5
```{r}
UrbanRural<-vector()
UrbanRural.ICL<-vector()
UrbanRural.ICU<-vector()
for(i in 1:(dim(rates.CyL)[1])){
  aux<-lm(log(rates.CyL[i,])~HDtype)
  UrbanRural[i]<-coef(aux)[2]
  UrbanRural.ICL[i]<-confint(aux)[2,1]
  UrbanRural.ICU[i]<-confint(aux)[2,2]
}

par(mfrow=c(1,1),bty="l")
plot(UrbanRural,xlab="Day",ylab="Urban risk excess (positive values)",main="(log-)risk excess in urban areas, CyL",type="l",ylim=c(-0.4,0.8),axes=F)
box()
axis(side = 2, at = 0.2 * (-2:4))
axis(side = 1, at = FirstOfMonth, labels = c("4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1"))
polygon(x = c(1:223,223:1),
        y = c(UrbanRural.ICL,rev(UrbanRural.ICU)),
        col=rgb(0, 0, 0, 0.2), 
        border = FALSE)
abline(h=0,col=2)
abline(v = FirstOfMonth, col = "green", lty = 2)
```



