#' @include generics.R
#'
NULL
#' An S4 object representing the SEIRDV.
#'
#' This class represents the SEIRV model, showing how populations of susceptible,
#' exposed, infectious and recovered individuals evolve over time. Vaccinated
#' individuals are considered in their own compartment.
#'
#' @slot output_names list of compartments name which are used by the model and
#' incidence.
#' @slot initial_condition_names list of names of initial conditions
#' (characters). Default is list("S0", "E0", "I0", R0", "V0", "VR0").
#' @slot transmission_parameter_names list of names of transmission parameters
#' (characters). Default is list("beta", "kappa", "gamma", "mu", "nu",
#' "delta_V", "delta_R", "delta_VR").
#' @slot intervention_parameter_names list of names of intervention parameters
#' (characters). Default is list ("starts", "stops", "coverages")
#' @slot initial_conditions list of values for initial conditions (double).
#' @slot transmission_parameters list of values for transmission parameters
#' (double).
#' @slot intervention_parameters list of values for intervention parameters
#' (double).
#'
#' @import deSolve
#' @import tidyverse
#' @import glue
#' @import reshape2
#'
#' @export SEIRDV
SEIRDV <- setClass("SEIRDV",
# slots
slots = c(
output_names = "list",
initial_condition_names = "list",
transmission_parameter_names = "list",
intervention_parameter_names = "list",
initial_conditions = "list",
transmission_parameters = "list",
intervention_parameters = "list"
),
# prototypes for the slots, automatically set parameter names and
# its data type
prototype = list(
output_names = list("S", "E", "I", "R", "V", "VR", "D", "Incidence", "Deaths"),
initial_condition_names = list("S0", "E0", "I0", "R0", "V0", "VR0"),
transmission_parameter_names = list("beta", "kappa", "gamma",
"mu", "nu", "delta_V",
"delta_R", "delta_VR"),
intervention_parameter_names = list("starts", "stops", "coverages"),
initial_conditions = vector(mode = "list", length = 6),
transmission_parameters = vector(mode = "list", length = 8),
intervention_parameters = vector(mode = "list", length = 3)
)
)
#' @describeIn SEIRDV Retrieves initial conditions of SEIRV model.
#'
#' @param object An object of the class SEIRDV.
#'
#' @export
setMethod("initial_conditions", "SEIRDV",
function(object) object@initial_conditions)
#' @describeIn SEIRDV Retrieves transmission parameters of SEIRV model.
#'
#' @param object An object of the class SEIRDV.
#'
#' @export
setMethod("transmission_parameters", "SEIRDV",
function(object) object@transmission_parameters)
#' @describeIn SEIRDV Setter method for initial conditions (S0, E0, I0, R0, V0, VR0)
#' of the SEIRV model.
#'
#' All initial conditions must sum up to 1.
#' If the initial conditions provided to do not sum to 1, an error is thrown.
#'
#' @param object an object of the class SEIRDV
#' @param value (list) list of initial conditions S0, E0, I0, R0, V0, VR0.
#'
#' @return object of class SEIRDV with initial conditions assigned.
#'
#' @export
setMethod(
"initial_conditions<-", "SEIRDV",
function(object, value) {
if (mean(names(value) %in% object@initial_condition_names) != 1)
stop(paste0("Initial conditions must contain: ",
object@initial_condition_names))
init_cond <- value
# raise errors if age category dimensions do not match initial state vectors
# also raise errors if initial state and parameter values are not doubles
for (p in list("S0", "E0", "I0", "R0", "V0", "VR0")) {
if (!is.numeric(init_cond[[p]])) {
stop(glue("{p} format must be numeric"))
}
}
# check that the initial conditions are properly normalized
if (init_cond$S0 + init_cond$E0 + init_cond$I0 +
init_cond$R0 + init_cond$V0 + init_cond$VR0 != 1) {
stop("Invalid initial conditions. Must add up to 1.")
}
object@initial_conditions <- init_cond
object
})
#' @describeIn SEIRDV Set transmission parameters (beta, kappa, gamma, mu, nu,
#' delta_V, delta_R) of the SEIRV model.
#'
#' If the transmission parameters provided to are not 1-dimensional an error is
#' thrown.
#'
#' @param object (SEIRDV model)
#' @param value (list) list of values for beta, kappa, gamma, mu, nu, delta_V and
#' delta_R respectively.
#'
#' @return object of class SEIRDV with transmission parameter values
#' assigned.
#'
#' @export
setMethod(
"transmission_parameters<-", "SEIRDV",
function(object, value) {
# create list of parameter values
if (mean(names(value) %in% object@transmission_parameter_names) != 1)
stop(paste0("Transmission parameters must contain: ",
object@transmission_parameter_names))
trans_params <- value
# check format of parameters
if (length(trans_params$beta) != 1
| length(trans_params$kappa) != 1
| length(trans_params$gamma) != 1
| length(trans_params$mu) != 1
| length(trans_params$nu) != 1
| length(trans_params$delta_V) != 1
| length(trans_params$delta_R) != 1
| length(trans_params$delta_VR) != 1) {
stop("The parameter values should be 1-dimensional.")
}
# if all above tests are passed, assign the trans_params namelist to the
# object
object@transmission_parameters <- trans_params
object
})
# SEIRDV class specific functions
#' @describeIn SEIRDV Retrieves intervention parameters of SEIRDV model.
#'
#' @param object An object of the class SEIRDV.
#'
#' @export
setMethod("intervention_parameters", "SEIRDV",
function(object) object@intervention_parameters)
#' @describeIn SEIRDV Setter method for intervention parameters of the SEIRV model.
#'
#' Intervention parameters have same size. A tanh function is used to smooth interventions during simulation. This class is designed for interventions
#' which last several days at least and have several days between them; interventions involving rapid fluctuations may be distorted.
#'
#' @param object an object of the class SEIRDV
#' @param value (list) list of intervention parameters: starts, stops and
#' coverages.
#'
#' @return object of class SEIRDV with intervention parameters assigned.
#'
#' @export
setMethod(
"intervention_parameters<-", "SEIRDV",
function(object, value) {
if (mean(names(value) %in% object@intervention_parameter_names) != 1)
stop(paste0("Intervention parameters must contain: ",
object@intervention_parameter_names))
interv_par <- value
# raise errors if intervention parameters are not doubles
for (p in list("starts", "stops", "coverages")) {
if (!is.numeric(interv_par[[p]])) {
stop(glue("{p} format must be numeric"))
}
}
# check that the intervention parameters are all of the same size
if (length(interv_par$starts) != length(interv_par$stops)|
length(interv_par$starts) != length(interv_par$coverages)|
length(interv_par$coverages) != length(interv_par$stops)) {
stop("Invalid intervention parameters. Must have same size.")
}
object@intervention_parameters <- interv_par
object
})
#' @describeIn SEIRDV Solves ODEs of the SEIRDV specified in object
#' for the time points specified in times and integration method specified in
#' solve_method.
#'
#' \deqn{\frac{dS(t)}{dt} = - beta S(t) I(t) - nu Inter(t) S(t) + delta_V V(t) + delta_R R(t) + delta_VR VR(t)}
#' \deqn{\frac{dE(t)}{dt} = beta S(t) I(t) - kappa E(t)}
#' \deqn{\frac{dI(t)}{dt} = kappa E(t) - (gamma + mu) I(t)}
#' \deqn{\frac{dR(t)}{dt} = gamma I(t) - delta_R R(t) - nu Inter(t) R(t)}
#' \deqn{\frac{dV(t)}{dt} = nu Inter(t) S(t) - delta_V V(t)}
#' \deqn{\frac{dVR(t)}{dt} = nu Inter(t) R(t) - delta_VR VR(t)}
#' \deqn{\frac{dC(t)}{dt} = beta S(t) I(t)}
#' \deqn{\frac{dD(t)}{dt} = mu I(t)}
#'
#' where Inter(t) is the value at time t of the intervention protocol defined by
#' the intervention parameters. This function relies on the package deSolve.
#'
#' @param object an object of the class SEIRDV
#' @param times (double) a sequence of time points at which the solution to
#' the system of ODEs should be returned. Must be of the form
#' seq(t_start, t_end, by=t_step). Default time series is seq(0, 100, by = 1).
#' @param solve_method (string) a string of chosen numerical integration method
#' for solving the ode system. Default is lsoda which is also the default for
#' the ode function in the deSolve package used in this function.
#'
#' @return two dataframes: one with the time steps, age range, time series of S,
#' E, I, R and V population fractions, and one with the time steps, age range,
#' time series of incidences and deaths population fraction.
#'
#' @export
#'
setMethod(
"run", "SEIRDV",
function(object, times, solve_method = "lsoda") {
if (!is.double(times)) {
stop("Evaluation times of the model storage format must be a vector.")
}
if (is.null(unlist(object@transmission_parameters)))
stop("Transmission parameters must be set before running.")
if (is.null(unlist(object@initial_conditions)))
stop("Initial conditions must be set before running.")
if (is.null(unlist(object@intervention_parameters)))
stop("Intervention parameters must be set before running.")
# set initial state vector
state <- c(S = initial_conditions(object)$S0,
E = initial_conditions(object)$E0,
I = initial_conditions(object)$I0,
R = initial_conditions(object)$R0,
V = initial_conditions(object)$V0,
VR = initial_conditions(object)$VR0,
C = 0,
D = 0)
# set transmission parameters vector
parameters <- c(b = transmission_parameters(object)$beta,
k = transmission_parameters(object)$kappa,
g = transmission_parameters(object)$gamma,
m = transmission_parameters(object)$mu,
n = transmission_parameters(object)$nu,
d_v = transmission_parameters(object)$delta_V,
d_r = transmission_parameters(object)$delta_R,
d_vr = transmission_parameters(object)$delta_VR)
# set intervention parameters vector
int_parms <-
InterventionParameters(
start=intervention_parameters(object)$starts,
stop=intervention_parameters(object)$stops,
coverage= intervention_parameters(object)$coverages)
# use tstep=0.1 and tanh_slope=1 for good nice step-function-like shape of
# the intervention wave
sim_parms <- SimulationParameters(start =0, stop = tail(times, n=1),
tstep = 0.1)
inter_prot <- intervention_protocol(int_parms, sim_parms, 1)
intervention <- approxfun(inter_prot$time, inter_prot$coverage, rule=2)
# function for RHS of ode system
right_hand_side <- function(t, state, parameters, input) {
with(
as.list(c(state, parameters)), {
inter <- input(t)
s <- state[1]
e <- state[2]
i <- state[3]
r <- state[4]
v <- state[5]
vr <- state[6]
c <- state[7]
d <- state[8]
# rate of change
ds <- -b * s * i - n * inter * s + d_v * v + d_r * r + d_vr * vr
de <- b * s * i - k * e
di <- k * e - (g + m ) * i
dr <- g * i - d_r * r - n * inter * r
dv <- n * inter * s - d_v * v
dvr <- n * inter * r - d_vr * vr
dc <- b * s * i
d_death <- m * i
# return the rate of change
list(c(ds, de, di, dr, dv, dvr, dc, d_death))
})
}
# call ode solver
out <- ode(
y = state, times = times, func = right_hand_side,
parms = parameters, method = solve_method, input = intervention)
output <- as.data.frame.array(out)
# Compute incidences and deaths
cases <- c(0, diff(output$C))
deaths <- c(0, diff(output$D))
output$Incidence <- cases
output$Deaths <- deaths
output <- output[, c("time", unlist(object@output_names))]
# Create long format of output
output <- melt(output, id.vars = "time")
output <- output[, c("time", "value", "variable")]
names(output) <- c("time", "value", "compartment")
# Added for consistency of output format across models
output$age_range <- rep("0-150", length(output$time))
# Split output into 2 dataframes: one with S,E,I,R and V and one with C and D
states <- subset(output, !output$compartment %in% c("Incidence", "Deaths"))
states <- droplevels(states)
changes <- subset(output, output$compartment %in% c("Incidence", "Deaths"))
changes <- droplevels(changes)
list("states" = states, "changes" = changes)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.