In mrc-ide/spatialbranchr: Spatially explicit branching process model

knitr::opts_chunk$set(ollapse = TRUE, comment = "#>")
options(rmarkdown.html_vignette.check_title = FALSE)

library(EpiEstim)
library(geosphere)
library(ggplot2)
library(MASS)
library(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
)

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.

mrc-ide/spatialbranchr documentation built on Oct. 2, 2022, 5:56 p.m.