ABC-SMC algorithm with two populations

Here is a possible solution for the ABC-SMC with two populations.

library(fitR)
library(MASS)
source(here::here("scripts", "snippets", "sumstat-examples.r"))
source(here::here("scripts", "snippets", "distance-examples.r"))

initState <- c(S = 250, E = 0, I = 4, T = 0, L = 30, Inc = 0)


# use the ABC rejection algorithm to find population 1 in the ABC-SMC algorithm

pop_1 <- my_abcAlgorithm(
  N = 1000, epsilon = 50,
  sumStats = list(ssMax, ssSize),
  distanceAbc = ssMeanRelDistance,
  fitmodel = seitlStoch,
  initState = initState,
  data = fluTdc1971
)


N <- dim(pop_1)[1]

# specify a smaller tolerance for the second population
epsilon_2 <- 5

# Sigma is the covariance matrix for the multivarite normal distribution pertubation
Sigma <- matrix(c(0.5, 0, 0, 0.5), 2, 2)

# set up empty matrix to store results
results <- matrix(nrow = 0, ncol = 6)

# initialise with i=0
i <- 0

# while the length of the accepted values (result) is less than the desired length (N)
while (i < N) {
  # - draw a random number between 1 and 1000
  rowNo <- sample(1000, 1)
  # extract corresponding row from pop_1 using this number
  pars <- pop_1[rowNo, c("D_lat", "D_inf")]

  # perturb these parameters using multivariate Gaussian distribution
  parsPerturb <- mvrnorm(n = 1, pars, Sigma)

  theta <- c(R_0 = 2, D_lat = parsPerturb[["D_lat"]], D_inf = parsPerturb[["D_inf"]], alpha = 0.9, D_imm = 13, rho = 0.85)

  # if any of the proposed parameters are less than 0, reject them (set distance to infinity)
  if (any(theta < 0)) {
    dist <- Inf
  } else {
    # use computeDistanceAbc to calculate a distance between the model
    # and data
    dist <- computeDistanceAbc(
      sumStats = list(ssMax, ssSize),
      distanceAbc = ssMeanRelDistance,
      fitmodel = seitlStoch,
      theta = theta,
      initState = initState,
      data = fluTdc1971
    )
  }

  ## if the model distance is within the epsilon window
  # store the accepeted parameter values
  if (dist <= epsilon_2) {
    results <- rbind(results, theta)
  }

  # update i (dimension of results store)
  i <- dim(results)[1]
}

# return the accepted values
head(results)

Return to the ABC session.

 

---

This web site and the material contained in it were originally created in support of an annual short course on Model Fitting and Inference for Infectious Disease Dynamics at the London School of Hygiene & Tropical Medicine. All material is under a MIT license. Please report any issues or suggestions for improvement on the corresponding GitHub issue tracker. We are always keen to hear about any uses of the material here, so please do get in touch using the Discussion board if you have any questions or ideas, or if you find the material here useful or use it in your own teaching.