Introduction

This vignette demonstrates how to implement particle particle Markov-chain Monte Carlo (PMCMC) for disease transmission using the agentSIS package. The small MCMC iterations and particles are used to reduce computation time; in practice, larger values (say, 20000 MCMC iterations and 300 particles) should be used.

Statistcal Agent-based SIS Model

Agent-based model (ABM) is a simulation based modeling technique that aims to describe complex dynamic processes with a bottom-up approach. ABM provides considerable flexibility by explaining the complex dynamic process using simple rules that incorporate characteristics of individual entities (agents) and their interaction. For example, infectious disease is transmitted by the interaction between neighbors and by the attributes of each agent whose states fall into susceptible or infected states.

In the disease transmission, the observed prevalence at time \(t\) can be modeled as
\[\begin{equation}\label{model:bin} y_t \vert I_t,\rho \sim \operatorname{Binomial}(I_t,\rho) \end{equation}\] where the detection probability \(\rho\) accounts for the under-reported prevalence.

The observed prevalence at time \(t\) depends on unobserved latent agent states \(\mathbf{X}_t = (X_t^1, \ldots, X_t^N)\), which can be modeled as an independent Bernoulli distribution; \[\begin{equation}\label{model:ber} X_t^n \ \vert \ \xi_{t-1}^n \sim \operatorname{Beroulli}(\xi_{t-1}^n). \end{equation}\] where \(\xi_{t-1}^n\) corresponds to the transition probability between susceptible and infected state for agent \(n\). Define agent-specific initial infection rate \(\alpha_0^n\), infection rate \(\lambda^n\) and recovery rate \(\gamma^n\) as \[\begin{align*} \alpha_0^n = \left(1 + \exp(-\beta_{\alpha_0}^Tz^n)\right)^{-1}, \quad \lambda^n = \left(1 + \exp(-\beta_\lambda^Tz^n)\right)^{-1}, \quad \gamma^n = \left(1+\exp(-\beta_\gamma^Tz^n)\right)^{-1} \end{align*}\] where \(\beta_{\alpha_0},\beta_\lambda,\beta_\gamma\in\mathbb{R}^d\) are parameters and \(z^n\in\mathbb{R}^d\) are the covariates of agent \(n\).

Given the initial transition probability \(\xi_0^n=\alpha_0^n\), the latent state \(X_t\) for \(t=2,\ldots,T\) evolves according to a Markov process with transition probability \[\begin{align*} \xi^{n}_{t-1}= \begin{cases}\lambda^{n} \mathcal{D}(n)^{-1} \sum_{m \in \mathcal{N}_n} X_{t-1}^{m}, & \text { if } x_{t-1}^{n}=0 \\ 1-\gamma^{n}, & \text { if } x_{t-1}^{n}=1\end{cases} \end{align*}\] where \(\mathcal{N}_n\) is a neighborhood for agent \(n\) and \(\mathcal{D}(n)\) is the number of neighbors of agent \(n\). This transition probability depends only on the last value of the state of agents which is defined as a first-order Markovian. Since attributes of agents \(z^n\) account for the infection and recovery rate, the transition probability \(\xi_{t-1}^n\) is defined uniquely for each agent.

While the underlying dynamics of the disease transmission are governed by agent states \(\mathbf{X}_t\), these states are not observed directly. For example, public health agencies only report aggregate data at the population level rather than at the individual level to protect individuals’ privacy; the number of infected agents is reported on a daily basis in COVID-19 data. As a result, the hidden Markov model represents our best understanding of ABM dynamics where the individual agent’s states are treated as hidden and estimated using filtering approaches such as particle filter (PF).

This figure highlights that it is suitable to treat agent states as hidden states and the heterogeneous transition rates between agents cannot be ignored whereas aggregate counts are observed.

## Load library
library(tidyverse) # Load the tidyverse, mainly for ggplot

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(kableExtra) # For fancy tables

## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(agentSIS) # Our package
library(zeallot) 
library(reshape2)

## 
## Attaching package: 'reshape2'
## 
## The following object is masked from 'package:tidyr':
## 
##     smiths

options(knitr.table.format = 'markdown')

set.seed(2022)
N <- 30;
steps <- 30;
X <- matrix(nrow = 2, ncol = N);
X[1,] <- 1;
X[2,] <- rnorm(N);

b_lam0 <- -1; b_lam1 <- 2; 
b_gam0 <- -1; b_gam1 <- -1;
Lambda <- fun_rate(b_lam0, b_lam1);
Gamma <- fun_rate(b_gam0, b_gam1);
Alpha <- rep(0.1, N)
Rho <- 0.8

## simulate observations
true_config <- list(N = N, alpha0 = Alpha, lambda = Lambda,
                    gamma = Gamma, rho = Rho)


sis_simulate <- function (days, true_config) {
  agent_states <- matrix(NA, nrow = true_config$N, ncol = 1 + days)
  agent_states[,1] <- (runif(N) < true_config$alpha0)
  for (d in c(2:(days + 1))) {
    alpha_t <- rep(NA, N)
    alpha_t <- sis_get_alpha(agent_states[, d - 1], true_config)
    agent_states[, d] <- (runif(N) < alpha_t)
  }
  y <- apply(agent_states, 2, function(xx) rbinom(n = 1, size = sum(xx), 
                                                  prob = true_config$rho))
  return(list(agent_states = agent_states, y = y))
}

complete_data <- sis_simulate(days = steps, true_config = true_config);
y <- complete_data$y;

hypers <- Hypers()

opt <- Opt(num_burn = 100, num_save = 100)

## genrate figure relecting hidden states and observations
A <- complete_data$agent_states
longData <- melt(A)
ind <- longData$value == TRUE

longData$status <- "Susceptible"
longData$status[ind] <- "Infected"
longData$status <- factor(longData$status, levels=c("Susceptible", "Infected"))

df <- data.frame(x=c(1:length(y)), y = y)

ggplot(longData, aes(x = Var2, y=Var1)) + 
  geom_tile(aes(fill = status, width=0.8, height=0.8), size=2) + 
  scale_fill_manual(values=c("grey90","grey60")) + 
  geom_line(data=df, aes(x = x, y = y), size=1) + 
  xlab("time") + ylab("agents") + 
  theme_bw() + 
  theme(legend.position = "bottom", legend.title=element_blank(), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) 

fit_PMMH <- PMMH_SIS(X, y, opt = opt)

## parameter est :  -2.714 2.261 1.95 4.325 -0.582 1.403 0.673 ( iteration: 100 ) 
## parameter est :  -3.2 2.308 1.024 3.081 -0.514 2.645 0.805 ( iteration: 200 )

fit_PG <- PG_SIS(X, y, opt = opt)

## parameter est :  3.821 6.849 2.209 -0.043 1.104 -0.521 0.722 ( iteration: 100 ) 
## parameter est :  3.049 5.32 1.793 -0.218 0.946 -0.18 0.712 ( iteration: 200 )

fit_PG_AS <- PG_AS_SIS(X, y, opt = opt)

## parameter est :  0.393 3.197 0.636 -2.566 -1.06 0.118 0.274 ( iteration: 100 ) 
## parameter est :  -0.677 1.429 0.893 -2.998 -1.041 0.106 0.395 ( iteration: 200 )

est_mat <- data.frame(time = rep(c(0:steps), 4), 
                      est = c(y, fit_PMMH$obs_counts_mean, fit_PG$obs_counts_mean,
                              fit_PG_AS$obs_counts_mean), 
                      grp = rep(c("observation","PMMH", "PG", "PGAS"), 
                                each=length(fit_PG$obs_counts_mean)))

ggplot(est_mat) + 
  geom_line(aes(x = time, y = est, color = grp, linetype = grp)) + theme_bw() +
  scale_color_manual(values=c("#999999", "#E69F00", "#117733", "#CC6677")) + 
  scale_linetype_manual(values=c("dashed","solid","solid","solid")) + 
  theme(legend.title=element_blank(),legend.position="bottom") + ylab("estimated counts")

table_PMMH <- rbind(post_mean = colMeans(fit_PMMH$param_mat),
apply(fit_PMMH$param_mat, 2, function(x) quantile(x, c(0.025, 0.975))))
colnames(table_PMMH) <-  c("$\\alpha_0$", "$\\alpha_1$", "$\\lambda_0$", "$\\lambda_1$","$\\gamma_0$","$\\gamma_1$","$\\rho$")
kable(table_PMMH)

table_PG <- rbind(post_mean = colMeans(fit_PG$param_mat),
apply(fit_PG$param_mat, 2, function(x) quantile(x, c(0.025, 0.975))))
colnames(table_PG) <-  c("$\\alpha_0$", "$\\alpha_1$", "$\\lambda_0$", "$\\lambda_1$","$\\gamma_0$","$\\gamma_1$","$\\rho$")
kable(table_PG)

table_PG_AS <- rbind(post_mean = colMeans(fit_PG_AS$param_mat),
apply(fit_PG_AS$param_mat, 2, function(x) quantile(x, c(0.025, 0.975))))
colnames(table_PG_AS) <- c("$\\alpha_0$", "$\\alpha_1$", "$\\lambda_0$", "$\\lambda_1$","$\\gamma_0$","$\\gamma_1$","$\\rho$")
kable(table_PG_AS)