# Simple MCMC under SIR In MultiBD: Multivariate Birth-Death Processes

We describe how to set-up and run a simple Metropolis-Hastings-based Markov chain Monte Carlo (MCMC) sampler under the susceptible-infected-removed (SIR) model.

library(MultiBD)


This example uses the Eyam data that consist the population counts of susceptible, infected and removed individuals across several time points.

data(Eyam)
Eyam


The log likelihood function is the sum of the log of the transition probabilities between two consecutive observations. Note that, we will use $(\log \alpha, \log \beta)$ as parameters instead of $(\alpha, \beta)$. The rows and columns of the transition probability matrix returned by dbd_prob() correspond to possible values of $S$ (from $a$ to $a0$) and $I$ (from $0$ to $B$) respectively.

loglik_sir <- function(param, data) {
alpha <- exp(param[1]) # Rates must be non-negative
beta  <- exp(param[2])

# Set-up SIR model
drates1 <- function(a, b) { 0 }
brates2 <- function(a, b) { 0 }
drates2 <- function(a, b) { alpha * b     }
trans12 <- function(a, b) { beta  * a * b }

sum(sapply(1:(nrow(data) - 1), # Sum across all time steps k
function(k) {
log(
dbd_prob(  # Compute the transition probability matrix
t  = data$time[k + 1] - data$time[k], # Time increment
a0 = data$S[k], b0 = data$I[k],       # From: S(t_k), I(t_k)
drates1, brates2, drates2, trans12,
a = data$S[k + 1], B = data$S[k] + data$I[k] - data$S[k + 1],
computeMode = 4, nblocks = 80         # Compute using 4 threads
)[1, data$I[k + 1] + 1] # To: S(t_(k+1)), I(t_(k+1)) ) })) }  Here, we choose$\text{Normal}(0, 100^2)$as the prior for both$\log \alpha$and$\log \beta$. logprior <- function(param) { log_alpha <- param[1] log_beta <- param[2] dnorm(log_alpha, mean = 0, sd = 100, log = TRUE) + dnorm(log_beta, mean = 0, sd = 100, log = TRUE) }  We will use the random walk Metropolis algorithm implemented in the function MCMCmetrop1R() (MCMCpack package) to explore the posterior distribution. So, we first need to install the package and its dependencies. source("http://bioconductor.org/biocLite.R") biocLite("graph") biocLite("Rgraphviz") install.packages("MCMCpack", repos = 'http://cran.us.r-project.org') library(MCMCpack)  # Provide manual caching because knitr's caching # is not working in my hands file <- system.file("vignetteCache", "post_sample.RData", package="MultiBD") if (!file.exists(file)) { <<loadStuff>> }  The starting point of our Markov chain is the estimated value of$(\alpha, \beta)$from Raggett (1982). alpha0 <- 3.39 beta0 <- 0.0212  We discard the first$200$iterations and keep the next$1000$iterations of the chain. post_sample <- MCMCmetrop1R(fun = function(param) { loglik_sir(param, Eyam) + logprior(param) }, theta.init = log(c(alpha0, beta0)), mcmc = 1000, burnin = 200)  # Provide manual caching because knitr's caching # is not working in my hands if (file.exists(file)) { load(file) } else { <<longMCMCRun>> # dir.create("cache", showWarnings = FALSE) save(post_sample, file = "../inst/vignetteCache/post_sample.RData") }  The trace plots of both$\log \alpha$and$\log \beta$look good. plot(as.vector(post_sample[,1]), type = "l", xlab = "Iteration", ylab = expression(log(alpha))) plot(as.vector(post_sample[,2]), type = "l", xlab = "Iteration", ylab = expression(log(beta)))  We can visualize the joint posterior distribution of$\log \alpha$and$\log \beta$using the ggplot2 package. library(ggplot2) x = as.vector(post_sample[,1]) y = as.vector(post_sample[,2]) df <- data.frame(x, y) ggplot(df,aes(x = x,y = y)) + stat_density2d(aes(fill = ..level..), geom = "polygon", h = 0.26) + scale_fill_gradient(low = "grey85", high = "grey35", guide = FALSE) + xlab(expression(log(alpha))) + ylab(expression(log(beta)))  We can also construct the$95\%$Bayesian credible intervals for$\alpha$and$\beta\$.

quantile(exp(post_sample[,1]), probs = c(0.025,0.975))
quantile(exp(post_sample[,2]), probs = c(0.025,0.975))


## Try the MultiBD package in your browser

Any scripts or data that you put into this service are public.

MultiBD documentation built on May 2, 2019, 11:50 a.m.