SVEIRD compartmental model with optional Bayesian data assimilation
Source:R/SVEIRD.BayesianDataAssimilation.R
SVEIRD.BayesianDataAssimilation.RdRun a SVEIRD compartmental model of an epidemic, optionally using Bayesian data assimilation.
Usage
SVEIRD.BayesianDataAssimilation(
alpha,
beta,
gamma,
sigma,
delta,
n.days,
seedData,
neighbourhood.order = 0,
lambda,
populationSpatRaster,
subregionsSpatVector,
aggregationFactor,
startDate,
countryCodeISO3C,
incidenceData,
deathData,
simulationIsDeterministic = TRUE,
dataAssimilationEnabled = FALSE,
healthZoneCoordinates,
variableCovarianceFunction,
forecastError.cov.sdBackground,
forecastError.cor.length,
neighbourhood.Bayes,
psi.diagonal,
callback
)Arguments
- alpha
The rate of vaccination (per day)
- beta
The rate of exposure (per day)
- gamma
The rate of becoming infectious (per day)
- sigma
The rate of recovery (per day)
- delta
The fatality rate (per day)
- n.days
The number of days the simulation will run for, beginning from the startDate.
- seedData
a dataframe, as described in initialInfections.fourCities.
- neighbourhood.order
The order of the queen's neighbourhood over which to distribute the Exposed and Infected intial state variable; zero is no neighbourhood, and corresponds to only the same cell where the queen already is. Higher orders follow the simple formula
(2 * x + 1)^2to determine the number of cells over which to evenly disperse the exposed and infected state variables. All other state variables are placed directly in the grid cell corresponding to the health zone. At the moment orders higher than one are prohibited, so the only valid values for this argument are zero and one.- lambda
the average dialy movement distance of an individual (in kilometers).
- populationSpatRaster
a SpatRaster of population count data; it must be pre-aggregated if any aggregation is to be used during the simulation; none is performed by this function or downstream functions in the overall implementation of SVEIRD.BayesianDataAssimilation simulations.
- subregionsSpatVector
a SpatVector object of subregions used to crop the population SpatRaster before creating the other layers in the raster object to represent the other compartments in a SVEIRD epidemic model. If it is missing no cropping is performed.
- aggregationFactor
The number of adjacent cells in any one direction to aggregate into a single cell. The aggregation factor must be the same as that used to generate the SpatRaster for layers.
- startDate
The date (in YYYY-MM-DD format) the simulation begins.
- countryCodeISO3C
The ISO three character code for a recognized country.
- incidenceData
A "situation report" dataframe. The first column provides the date of the officially reported, observed incidence of the disease, in YYYY-MM-DD format.
Date Beni Biena Butembo Goma Kalunguta 05-Aug 3 0 2 0 0 12-Aug 2 0 0 0 0 20-Aug 1 0 0 0 0 26-Aug 5 0 0 0 0 02-Sep 8 0 0 0 1 09-Sep 5 0 2 0 0 16-Sep 5 0 3 0 0 23-Sep 4 0 1 0 0 02-Oct 10 0 1 0 0 07-Oct 14 0 4 0 0 15-Oct 28 0 2 0 1 21-Oct 19 0 3 0 0 28-Oct 28 0 8 0 0 04-Nov 16 0 6 0 1 11-Nov 14 0 0 0 21- deathData
Data of the same format as incidenceData, but observations represent deaths, not infections.
- simulationIsDeterministic
Whether stochasticity is enabled or not; if the simulation is deterministic then no stochastic processes are used and the simulation is entirely deterministic. Either
TRUEorFALSE.- dataAssimilationEnabled
Whether Bayesian data assimilation will be used for state reporting data.
- healthZoneCoordinates
The coordinates of health zones in the country of interest, TODO: describe the use of the data so users understand why they must include it; e.g. (which will be used to group and summarize the compartmental model data at the end of the simmulation.)
HealthZone Latitude Longitude Alimbongo -0.365515 29.1911818 Beni 0.49113 29.47306 Biena 0.57923 29.115633 Butembo 0.140692 29.335014 Goma -1.658271 29.220132 Kalunguta 0.323085 29.354677 Katwa 0.116985 29.411838 Kayna -0.603936 29.174066 Kyondo -0.005622 29.408813 Lubero -0.15633 29.24057 Mabalako 0.461257 29.210687 Manguredjipa 0.353433 28.765479 Masereka -0.133333 29.333333 Musienene 0.04022 29.26246 Mutwanga 0.32893 29.74576 Nyiragongo -1.516667 29.25 Oicha 0.698681 29.518834 Pinga -0.9830504 28.687911 Vuhovi 0.1416 29.4075 Ariwara 3.136873 30.706615 Bunia 1.566667 30.25 Komanda 1.367496 29.774322 Lolwa 1.352969 29.498455 Mambasa 1.359731 29.029226 Mandima 1.35551 29.08173 Nyakunde 1.431271 30.029283 Rwampara 1.4053 30.3449 Tchomia 1.4412 30.4845- variableCovarianceFunction
Passed directly to
forecastError.cov()to generate a forecastErrorCovariance matrix.- forecastError.cov.sdBackground
the "background" or default amount of error, in standard deviations.
- forecastError.cor.length
TODO
- neighbourhood.Bayes
The neighbourhood used in Bayesian data assimilation; this is different from that used in the casting of seed data about a neighbourhood of cells.
- psi.diagonal
A replacement value for elements of the Psi matrix' diagonal which are zero.
- callback
a callback function to run, with no arguments, which will be called every time the main loop of the simulation iterates, or a list of callback functions which run before, during, and after the loop, and which have the following structure. For each component, if the component is a list it must have members fun and args, where fun is a function symbol and args is a list of named arguments to the function; if it is not a list, the component of the list (before, during, or after) must be a function. See the examples.
Value
a list with components table, a tibble, and timeseries, a SpatRasterDataset representing timeseries of different variables.
Examples
IturiNordKivu <- terra::vect(
system.file(
"extdata",
## COD: Nord-Kivu and Ituri (Democratic Republic of Congo)
"subregionsSpatVector",
package = "spatialEpisim.foundation",
mustWork = TRUE
)
)
populationDemocraticRepublicCongo <- terra::rast(
system.file(
"extdata",
"susceptibleSpatRaster.tif", # Congo population
package = "spatialEpisim.foundation",
mustWork = TRUE
)
)
data("healthZonesCongo", package = "spatialEpisim.foundation")
data("initialInfections.fourCities", package = "spatialEpisim.foundation")
data("Congo.EbolaIncidence", package = "spatialEpisim.foundation")
rasterAggregationFactor = 9
simulation.days <- 31
if (requireNamespace("cli", quietly = TRUE)) {
callback <- list(before = list(fun = cli::cli_progress_bar,
args = list(name = "Simulating epidemic (SEI-type)",
total = simulation.days)),
during = cli::cli_progress_update,
after = cli::cli_progress_done)
} else {
callback <- "{" # does nothing
}
SVEIRD.BayesianDataAssimilation(
## Parameters
alpha = 3.5e-5,
beta = 7e-3,
gamma = 1/7,
sigma = 1/36,
delta = 2/36,
lambda = 18,
## Model runtime
n.days = simulation.days,
## Model data
seedData = initialInfections.fourCities,
neighbourhood.order = 1,
populationSpatRaster = populationDemocraticRepublicCongo,
subregionsSpatVector = IturiNordKivu,
aggregationFactor = rasterAggregationFactor,
startDate = "2018-08-01",
countryCodeISO3C = "COD",
## Model options
simulationIsDeterministic = TRUE,
dataAssimilationEnabled = FALSE,
callback = callback
)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'global': object 'rasterAggregationFactor' not found
## A large, lengthy simulation with Bayesian data assimilation (runtime
## approximately four minutes).
simulation.days <- 28
SVEIRD.BayesianDataAssimilation(
## Parameters
alpha = 3.5e-5,
beta = 7e-3,
gamma = 1/7,
sigma = 1/36,
delta = 2/36,
lambda = 18,
## Model runtime
n.days = simulation.days,
## Model data
seedData = initialInfections.fourCities,
neighbourhood.order = 1,
populationSpatRaster = populationDemocraticRepublicCongo,
aggregationFactor = rasterAggregationFactor,
subregionsSpatVector = IturiNordKivu,
startDate = "2018-08-01",
countryCodeISO3C = "COD",
incidenceData = head(Congo.EbolaIncidence, n = 4),
## Model options
simulationIsDeterministic = TRUE,
dataAssimilationEnabled = TRUE,
healthZoneCoordinates = healthZonesCongo,
variableCovarianceFunction = "DBD",
## Special parameters
forecastError.cov.sdBackground = 0.55,
forecastError.cor.length = 6.75e-1,
neighbourhood.Bayes = 3,
psi.diagonal = 1e-3
)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'global': object 'rasterAggregationFactor' not found