Jointly estimate effective reproduction number and gravity model parameters • spatialbranchr

library(EpiEstim)
library(geosphere)
library(ggplot2)
library(MASS)
library(spatialbranchr)

#> Warning: replacing previous import 'Rcpp::LdFlags' by 'RcppParallel::LdFlags'
#> when loading 'spatialbranchr'

library(tidyr)

We will data from the West African Ebola epidemic to jointly estimate the effective reproduction number and two parameters of the gravity model capturing human population movement - (a) pstay, which characterises the probability that an infectious individual will not move out of a spatial unit (here country) during their infectious period, and (b) gamma, which modulates the effect of the distance between two spatial units on the flow of populations between them.

Data

We will use data from the West African Ebola epidemic bundled with the package.

data("promed_ebola")
data("metadata")
incid <- promed_ebola
## Get rid of dates column as we don't need it for fitting
incid_raw <- incid[, -1]
nloc <- ncol(incid_raw)
centroids <- metadata
si <- discr_si(k = 0:35, mu = 6, sigma = 4)
priors <- spatial_priors()

For illustration, we will focus on a few countries instead of including all countries in the estimation, which can take a very long time to run.

## countries of interest
cntrs <- c("SLE", "GIN", "LBR", "GHA", "NGA")

` The data set has more than 600 rows (each row is a day). We will estimate the model parameters using the first 35 days of data.

x <- incid_raw[1:35, cntrs]
pop <- centroids$Pop[centroids$ISO3 %in% cntrs]
coords <- centroids[centroids$ISO3 %in% cntrs, c("Centroid_Lon", "Centroid_Lat")]
## Generate distance matrix from centroids using geosphere package
dist <- distm(coords)
## This can take some time. 120 seconds on my computer with small-ish
## numver of iterations
out <- spatial_estimate(
  x, si, window = 7L, NULL, pop, dist, 1, 1, 1, "poisson", priors, iter = 2000, chains = 2
)

#> 
#> SAMPLING FOR MODEL 'estimate_both' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001845 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 18.45 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 30.3991 seconds (Warm-up)
#> Chain 1:                22.3668 seconds (Sampling)
#> Chain 1:                52.7659 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'estimate_both' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 0.001585 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 15.85 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 32.792 seconds (Warm-up)
#> Chain 2:                25.3166 seconds (Sampling)
#> Chain 2:                58.1086 seconds (Total)
#> Chain 2:

The object out now has the fitted model. You can check ggmcmc to check for convergence. Note that spatialbranchr does not carry out any convergence checks. As a final step, you might extract the variables from out into a more friendly format.

processed <- process_spatial_estimate(out, x, window = 7L)

The first component of this list is three-dimensional array where the dimensions are time, space, and samples from the posterior ditribution of Rt. We can summarise these:

rt <- processed$Rt
rt_summary <- apply(rt, c(1, 2), quantile, probs = c(0.025, 0.5, 0.975))

Since we have estimated a single Rt over a window for each location, the same samples have been “filled in” for each day over the window.

We can similarly summarise the parameters gamma and ptsay. We can also use them to obtain a flow matrix which can then be used as input in spatial_project.