R/multi_stoch.R

In qleclerc/efficientspatial: Efficient Spatial Modelling of Epidemics

#' @title Runs multiple stochastic simulations (Raster)
#'
#' @description Simulates multiple stochastic epidemics using the provided Rasterlayer, spatial kernel, contact
#'              matrix, and infection parameters.
#'
#' @param num_runs The number of stochastic runs to perform.
#' @param rasterl The RasterLayer object containing the population data.
#' @param D The expanded kernel matrix to use for FOI calculation (generated by the \code{\link{calc_beta}}
#'          function).
#' @param contact_mat The contact matrix for mixing between age groups.
#' @param beta The beta value for the epidemic (calculated from a given R0 using the \code{\link{calc_beta}}
#'             function).
#' @param sigma The recovery rate for the epidemic (must match the one used to calculate beta from R0
#'              using the \code{\link{calc_beta}} function).
#' @param step Size of time step for stochastic simulation, in days (default is 1 day).
#' @param start_area Where to start the epidemic. 1: Most highly populated area (default), 2: A
#'                   random area in the middle of the country (typically medium population density), 3: A
#'                   random area in the north of the country (typically low population density). NOTE: CURRENTLY
#'                   ONLY SUPPORTS OPTION 1
#' @param start_num Number of infected individuals to start the epidemic.
#' @param t_max How many days to run the simulation for.
#'
#' @details This functions requires specific objects to run. These can be generated using the \code{\link{prep_simulation}}
#'          function (e.g. if you want to simulate an epidemic using the RasterLayer object "toy_data", you must
#'          prep_simulation(toy_data) first). The model used is an SIR model, where individuals can be either
#'          Susceptible, Infected or Recovered with regards to the disease. This assumes that Infected individuals
#'          are infectious, and that Recovered individuals are immune and cannot be reinfected.
#'
#' @return Returns one dataframe object containing the median and 95% range of all of the simulation results.
#'
#' @examples
#'
#' #Create a RasterLayer object:
#' test_data = raster(nrow=10, ncol=10, xmn=1, xmx=100000, ymn=1, ymx=100000)
#' values(test_data) = runif(100, 1, 1000)
#'
#' #Calculate the parameters for the simulation:
#' prep_simulation(test_data)
#'
#' results = multi_stoch(100, test_data, expanded_D, contact_mat, beta = beta)
#' View(results)
#'
#'
#' @export


multi_stoch = function(num_runs=100, rasterl, D, contact_mat, beta, sigma=1/2.6, step=1, start_area=1, start_num=10, t_max=200){

  #create matrix to store the results:
  total_results = matrix(0, nrow = (t_max+1), ncol = num_runs)

  num_ages = dim(contact_mat)[1]    #derive number of age categories from contact matrix
  num_areas = dim(D)[1]/4    #derive number of areas



  for(z in 1:num_runs){

    #counter to display progress during model execution:
    print(paste0("Run ", z, " out of ", num_runs))

    results = run_simulation(rasterl, D, contact_mat, beta, sigma=sigma, t_max=t_max, stoch=T, step=step, start_area = start_area, start_num = start_num)

    #calculate incidence:
    results_mem = rowSums(results[seq(1,max(results[,1])+1, 1/step),2:(num_areas*num_ages+1)])
    results_mem = c(0,abs(diff(results_mem)))

    #save incidence values:
    total_results[,z] = results_mem

  }


  #calculate median and 95% range:
  median_results = apply(total_results, 1, median)
  top_results = apply(total_results, 1, quantile, 0.975)
  bot_results = apply(total_results, 1, quantile, 0.025)

  #save results:
  final_results = data.frame(Time=seq(0,t_max,1), Median=median_results, Upper=top_results, Lower=bot_results)

  return(final_results)
}