R/multi_stoch_shp.R

In qleclerc/efficientspatial: Efficient Spatial Modelling of Epidemics

#' @title Runs multiple stochastic simulations (Shapefile)
#'
#' @description Simulates multiple stochastic epidemics using the provided Shapefile, spatial kernel, contact
#'              matrix, and infection parameters.
#'
#' @param num_runs The number of stochastic runs to perform.
#' @param shp_data The dataframe object containing the population data extracted from the Shapefile object
#'                 (generated by the \code{\link{prep_simulation_shp}} function)
#' @param D The expanded kernel matrix to use for FOI calculation (generated by the \code{\link{calc_beta_shp}}
#'          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_shp}}
#'             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_shp}} function).
#' @param step Size of time step for stochastic simulation, in days (default is 1 day).
#' @param start_area Where to start the epidemic. You can provide the name of the area directly (you can check
#'                   which are possible by looking inside the shp_data dataframe). If using the District, County
#'                   or Region resolution, the default is London. Otherwise, no default is provided.
#' @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_shp}}
#'          function (e.g. if you want to simulate an epidemic using the Shapefile object "region_shp", you must
#'          prep_simulation_shp(region_shp) 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.
#'
#'
#' @export

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

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

  num_areas = length(shp_data$Population)   #derive number of areas
  num_ages = dim(contact_mat)[1]    #derive number of age categories from contact matrix


  for(z in 1:num_runs){

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

    results = run_simulation_shp(shp_data, 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, 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)
}