R/run_multi_stoch.R

In qleclerc/epicspatial: Flexible Spatial Modelling of Epidemics

#' @title Runs multiple stochastic simulations
#'
#' @description Simulates multiple stochastic epidemics using the provided spatial dataset, spatial kernel, contact
#'              matrix, and infection parameters.
#'
#' @param num_runs The number of stochastic runs to perform.
#' @param spatial_data The spatial dataset 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. Default: area with highest population density. If the spatial
#'                   dataset is a SpatialPolygonsDataFrame, you can specify a string indicating the name of your
#'                   desired starting area. Alternatively, no matter which type of spatial dataset you have, you can
#'                   indicate the index of you desired starting area.
#' @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 spatial dataset "testdata", you must
#'          prep_simulation(test_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 spatial dataset:
#' 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 = run_multi_stoch(100, test_data, expanded_D, contact_mat, beta = beta)
#' View(results)
#'
#'
#' @export


run_multi_stoch = function(num_runs=100, spatial_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_ages = dim(contact_mat)[1]    #derive number of age categories from contact matrix
  num_areas = dim(D)[1]/num_ages    #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(spatial_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, 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)
}