R/SEIRfansy.R
In Ritwik-Bhaduri/SEIR-fansy: Extended Susceptible-Exposed-Infected-Recovery Model
Defines functions
SEIRfansy
Documented in
SEIRfansy
#' Estimate SEIRfansy Model Parameters
#'
#' This function is used to estimate the different parameters of interest like
#' the time varying transmission rates, proportion of reported cases, and the
#' basic reproduction rates.
#'
#' @param N (mandatory): The population size.
#' @param data (mandatory): If the model is Multinomial, then the data matrix
#' should contain the 3 columns Confirmed, Recovered, and Death. If the model
#' is Poisson or Binomial, then the data should contain only the column Confirmed.
#' Please ensure that the names of the columns are exactly as stated above.
#' @param data_init (mandatory): These are the initial data values
#' provided by the user as a numeric vector of length six. The entries should be
#' the Total Confirmed, Total Recovered, Total Death, Daily Confirmed,
#' Daily Recovered, and Daily Death for the Starting Date. Note: If the starting total
#' confirmed is 0, please replace it by 1.
#' @param init_pars NULL(default): If not equal to NULL, then the user can give a user input initial
#' parameters which should consist of the initial values of the time varying beta,
#' proportion of testing for the different periods.
#' @param niter 1e5 (default): Number of iterations for the MCMC Metropolis Hastings
#' algorithm.
#' @param BurnIn 5e4 (default): This is the number of burn-in iterations for the MCMC algorithm
#' @param model "Multinomial" (default): This is the likelihood function that will
#' be used. There are three options available: "Multinomial", "Poisson", or "Binomial".
#' @param plot TRUE (default): This will give the box plot for the basic reproduction
#' number for the different periods.
#' @param auto.initialize TRUE (default): This is the option for using a mle
#' based initial parameter.
#' @param period_start The total time period is divided into small periods depending
#' on the lock down measures imposed by the government. So this is a numeric vector
#' consisting of the start dates for the different time periods.
#' @param ... arguments passed to the function model_initializeR and model_plotR
#' which is used for initializing the parameters. The parameters are
#' described below:
#' \itemize{
#' \item{step_pars: }{init_pars/500 (default): It is the variance
#' of the proposal distribution for the Metropolis Hastings Algorithm
#' which is assumed to be a Random Walk.}
#' \item{alpha_p: }{0.5 (default): It is defined as the ratio
#' of rate of spread of infection by tested positive patients to
#' that by false negatives. We have taken $\alpha_p < 1$ as patients
#' who are tested positive are subjected to quarantine where the
#' chance of spreading the disease is less than that of false
#' negative patients who are mostly unaware of their infectious
#' nature. So, false negative individuals continue to spread the
#' disease freely at a higher rate than tested positive individuals.}
#' \item{alpha_u: }{0.7 (default): It is defined as the scaling
#' factor for the rate of spread of infection by untested individuals.
#' $\alpha_u$ is assumed to be < 1 as U mostly consists of
#' asymptomatic or mildly symptomatic cases who are known to spread
#' the disease at a much lower rate than those with higher levels
#' of symptoms.}
#' \item{beta_1: }{0.6 (default): It is the scaling factor for
#' rate of recovery for untested individuals. $\beta_1$ is
#' assumed to be less than 1. The condition of Untested individuals
#' is not as severe as they consist of mostly asymptomatic people.
#' So, they are assumed to recover faster than the Current
#' Positive ones.}
#' \item{beta_2: }{0.7 (default): It is the inverse of the
#' scaling factor for rate of recovery for false negative individuals.
#' $\beta_2$ is assumed to be less than 1. It is assumed that the
#' recovery rate is slower than the detected ones for the False
#' Negative ones because they are not getting any hospital treatments.}
#' \item{delta_1: }{0.3 (default): It is the scaling factor for
#' death rate for undetected individuals. $\delta_1$ is assumed
#' to be less than 1. Similarly for the Untested ones, the death
#' rate is taken to be lesser because they are mostly asymptomatic.
#' So, their probability of dying is much lower.}
#' \item{delta_2: }{0.7 (default): It is the inverse of the
#' scaling factor for death rate for false negative individuals.
#' $\delta_2$ is assumed to be less than 1. Same as before, the
#' death rate for False Negative ones are assumed to be higher than
#' the Current detected Positive as they are not receiving proper
#' treatment.}
#' \item{lambda: }{1 / (66.26 * 365) (default): Natural birth
#' rate. The value given here as the default value is the world's
#' common birth rate.}
#' \item{mu: }{1 / (66.26 * 365) (default): Natural death rates.
#' This is assumed to be equal with natural birth rate for the sake
#' of simplicity.}
#' \item{D_d: }{17.8 (default): Mean days until death for positive
#' individuals.}
#' \item{D_e: }{5.2 (default): Incubation period.}
#' \item{Dr: }{17.8 (default): Mean days until recovery for positive
#' individuals.}
#' \item{f: }{0.15 (default): False negative probability of
#' RT-PCR test.}
#' \item{mCFR: }{NULL (It is calculated from the data by default)
#' It is defined as the ratio of the total reported deceased cases
#' and the total removed cases until that day.}
#' \item{init.exposed.ratio: }{3 (default): This is the scaling
#' factor for the calculation of the initial number of exposed people
#' from the sum of the initial number of unreported, reported people.}
#' \item{init.confirmed.ratio: }{0.15 (default): This is the initial
#' value of the probability of being tested.}
#' \item{opt_num: }{100 (default): The number of times an user
#' wants to run the mle optimization before deciding on the best
#' initial parameter.}
#' \item{trace_plot.common_axis: }{FALSE (default): This will give
#' the trace plot for the convergence of the MCMC estimated time
#' varying parameters.}
#' \item{save plot: }{TRUE (default): It is the option for saving
#' the plots in the directory folder.}
#' }
#'
#'
#'
#' @return A list with class "SEIRfansy", which contains the items described
#' below:
#' \itemize{
#' \item{mcmc_pars: }{a matrix of the mcmc draws for the parameters}
#' \item{plots: }{a list of ggplot objects}
#' }
#'
#' @export
#'
#' @examples
#' library(dplyr)
#' train = covid19[which(covid19$Date == "01 April "):which(covid19$Date == "30 June "),]
#' test = covid19[which(covid19$Date == "01 July "):which(covid19$Date == "31 July "),]
#'
#' train_multinom =
#' train %>%
#' rename(Confirmed = Daily.Confirmed,
#' Recovered = Daily.Recovered,
#' Deceased = Daily.Deceased) %>%
#' dplyr::select(Confirmed, Recovered, Deceased)
#'
#' test_multinom =
#' test %>%
#' rename(Confirmed = Daily.Confirmed,
#' Recovered = Daily.Recovered,
#' Deceased = Daily.Deceased) %>%
#' dplyr::select(Confirmed, Recovered, Deceased)
#'
#' N = 1341e6 # population size of India
#' data_initial = c(2059, 169, 58, 424, 9, 11)
#' pars_start = c(c(1,0.8,0.6,0.4,0.2), c(0.2,0.2,0.2,0.25,0.2))
#' phases = c(1,15,34,48,62)
#'
#' cov19est = SEIRfansy(data = train_multinom, init_pars = pars_start,
#' data_init = data_initial, niter = 1e3, BurnIn = 1e2,
#' model = "Multinomial", N = N, lambda = 1/(69.416 * 365),
#' mu = 1/(69.416 * 365), period_start = phases, opt_num = 1,
#' auto.initialize = TRUE, f = 0.15)
#'
#' names(cov19est)
#' class(cov19est$mcmc_pars)
#' names(cov19est$plots)
#'
#' plot(cov19est, type = "trace")
#' plot(cov19est, type = "boxplot")
#'
#'
SEIRfansy <- function(data, data_init,N, init_pars = NULL, niter = 1e5, BurnIn = 5e3,
model = "Multinomial", plot = TRUE, period_start, auto.initialize = TRUE,
... ){
#library(dplyr)
# library(arm)
# library(pbapply)
# library(DescTools)
# library("patchwork")
if(model == "Multinomial"){
if(! all(c("Confirmed", "Recovered", "Deceased" ) %in% colnames(data))){
stop("Data does not contain required columns")
}
} else if (model == "Poisson" || model == "Binomial"){
if(! all(c("Confirmed") %in% colnames(data))){
stop("Data does not contain required columns")
}
} else
stop("Incorrect model specified")
if(length(data_init) != 6) stop("Incorrect initial values: length must be 6")
if(any(data < 0) || data_init < 0) stop("Negative values in data")
var_init = model_initializeR(data = data, data_init = data_init, init_pars = init_pars, model = model,
period_start = period_start, auto.initialize = auto.initialize, ...)
fix_pars = var_init$fix_pars
init_state_num = var_init$init_state_num
init_pars = var_init$init_pars
period_start = var_init$period_start
pars_estimate = mcmc_performR(data = data, init_state_num = init_state_num, init_pars = init_pars, model = model,
niter = niter, BurnIn = BurnIn, fix_pars = fix_pars, period_start = period_start, ...)
## Plot
if(plot == TRUE){
plots <- model_plotR(mcmc_pars = pars_estimate, n_period = length(period_start), data = data,
data_initial = data_initial, call = "estimateR", model = model, ...)
rlist = list(mcmc_pars = pars_estimate, "plots" = plots)
class(rlist) = "SEIRfansy"
return(rlist)
}
rlist = list(mcmc_pars = pars_estimate)
class(rlist) = "SEIRfansy"
return(rlist)
}
Ritwik-Bhaduri/SEIR-fansy documentation built on Oct. 2, 2020, 12:03 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions SEIRfansy
Documented in SEIRfansy
#' Estimate SEIRfansy Model Parameters
#'
#' This function is used to estimate the different parameters of interest like
#' the time varying transmission rates, proportion of reported cases, and the
#' basic reproduction rates.
#'
#' @param N (mandatory): The population size.
#' @param data (mandatory): If the model is Multinomial, then the data matrix
#' should contain the 3 columns Confirmed, Recovered, and Death. If the model
#' is Poisson or Binomial, then the data should contain only the column Confirmed.
#' Please ensure that the names of the columns are exactly as stated above.
#' @param data_init (mandatory): These are the initial data values
#' provided by the user as a numeric vector of length six. The entries should be
#' the Total Confirmed, Total Recovered, Total Death, Daily Confirmed,
#' Daily Recovered, and Daily Death for the Starting Date. Note: If the starting total
#' confirmed is 0, please replace it by 1.
#' @param init_pars NULL(default): If not equal to NULL, then the user can give a user input initial
#' parameters which should consist of the initial values of the time varying beta,
#' proportion of testing for the different periods.
#' @param niter 1e5 (default): Number of iterations for the MCMC Metropolis Hastings
#' algorithm.
#' @param BurnIn 5e4 (default): This is the number of burn-in iterations for the MCMC algorithm
#' @param model "Multinomial" (default): This is the likelihood function that will
#' be used. There are three options available: "Multinomial", "Poisson", or "Binomial".
#' @param plot TRUE (default): This will give the box plot for the basic reproduction
#' number for the different periods.
#' @param auto.initialize TRUE (default): This is the option for using a mle
#' based initial parameter.
#' @param period_start The total time period is divided into small periods depending
#' on the lock down measures imposed by the government. So this is a numeric vector
#' consisting of the start dates for the different time periods.
#' @param ... arguments passed to the function model_initializeR and model_plotR
#' which is used for initializing the parameters. The parameters are
#' described below:
#' \itemize{
#' \item{step_pars: }{init_pars/500 (default): It is the variance
#' of the proposal distribution for the Metropolis Hastings Algorithm
#' which is assumed to be a Random Walk.}
#' \item{alpha_p: }{0.5 (default): It is defined as the ratio
#' of rate of spread of infection by tested positive patients to
#' that by false negatives. We have taken $\alpha_p < 1$ as patients
#' who are tested positive are subjected to quarantine where the
#' chance of spreading the disease is less than that of false
#' negative patients who are mostly unaware of their infectious
#' nature. So, false negative individuals continue to spread the
#' disease freely at a higher rate than tested positive individuals.}
#' \item{alpha_u: }{0.7 (default): It is defined as the scaling
#' factor for the rate of spread of infection by untested individuals.
#' $\alpha_u$ is assumed to be < 1 as U mostly consists of
#' asymptomatic or mildly symptomatic cases who are known to spread
#' the disease at a much lower rate than those with higher levels
#' of symptoms.}
#' \item{beta_1: }{0.6 (default): It is the scaling factor for
#' rate of recovery for untested individuals. $\beta_1$ is
#' assumed to be less than 1. The condition of Untested individuals
#' is not as severe as they consist of mostly asymptomatic people.
#' So, they are assumed to recover faster than the Current
#' Positive ones.}
#' \item{beta_2: }{0.7 (default): It is the inverse of the
#' scaling factor for rate of recovery for false negative individuals.
#' $\beta_2$ is assumed to be less than 1. It is assumed that the
#' recovery rate is slower than the detected ones for the False
#' Negative ones because they are not getting any hospital treatments.}
#' \item{delta_1: }{0.3 (default): It is the scaling factor for
#' death rate for undetected individuals. $\delta_1$ is assumed
#' to be less than 1. Similarly for the Untested ones, the death
#' rate is taken to be lesser because they are mostly asymptomatic.
#' So, their probability of dying is much lower.}
#' \item{delta_2: }{0.7 (default): It is the inverse of the
#' scaling factor for death rate for false negative individuals.
#' $\delta_2$ is assumed to be less than 1. Same as before, the
#' death rate for False Negative ones are assumed to be higher than
#' the Current detected Positive as they are not receiving proper
#' treatment.}
#' \item{lambda: }{1 / (66.26 * 365) (default): Natural birth
#' rate. The value given here as the default value is the world's
#' common birth rate.}
#' \item{mu: }{1 / (66.26 * 365) (default): Natural death rates.
#' This is assumed to be equal with natural birth rate for the sake
#' of simplicity.}
#' \item{D_d: }{17.8 (default): Mean days until death for positive
#' individuals.}
#' \item{D_e: }{5.2 (default): Incubation period.}
#' \item{Dr: }{17.8 (default): Mean days until recovery for positive
#' individuals.}
#' \item{f: }{0.15 (default): False negative probability of
#' RT-PCR test.}
#' \item{mCFR: }{NULL (It is calculated from the data by default)
#' It is defined as the ratio of the total reported deceased cases
#' and the total removed cases until that day.}
#' \item{init.exposed.ratio: }{3 (default): This is the scaling
#' factor for the calculation of the initial number of exposed people
#' from the sum of the initial number of unreported, reported people.}
#' \item{init.confirmed.ratio: }{0.15 (default): This is the initial
#' value of the probability of being tested.}
#' \item{opt_num: }{100 (default): The number of times an user
#' wants to run the mle optimization before deciding on the best
#' initial parameter.}
#' \item{trace_plot.common_axis: }{FALSE (default): This will give
#' the trace plot for the convergence of the MCMC estimated time
#' varying parameters.}
#' \item{save plot: }{TRUE (default): It is the option for saving
#' the plots in the directory folder.}
#' }
#'
#'
#'
#' @return A list with class "SEIRfansy", which contains the items described
#' below:
#' \itemize{
#' \item{mcmc_pars: }{a matrix of the mcmc draws for the parameters}
#' \item{plots: }{a list of ggplot objects}
#' }
#'
#' @export
#'
#' @examples
#' library(dplyr)
#' train = covid19[which(covid19$Date == "01 April "):which(covid19$Date == "30 June "),]
#' test = covid19[which(covid19$Date == "01 July "):which(covid19$Date == "31 July "),]
#'
#' train_multinom =
#' train %>%
#' rename(Confirmed = Daily.Confirmed,
#' Recovered = Daily.Recovered,
#' Deceased = Daily.Deceased) %>%
#' dplyr::select(Confirmed, Recovered, Deceased)
#'
#' test_multinom =
#' test %>%
#' rename(Confirmed = Daily.Confirmed,
#' Recovered = Daily.Recovered,
#' Deceased = Daily.Deceased) %>%
#' dplyr::select(Confirmed, Recovered, Deceased)
#'
#' N = 1341e6 # population size of India
#' data_initial = c(2059, 169, 58, 424, 9, 11)
#' pars_start = c(c(1,0.8,0.6,0.4,0.2), c(0.2,0.2,0.2,0.25,0.2))
#' phases = c(1,15,34,48,62)
#'
#' cov19est = SEIRfansy(data = train_multinom, init_pars = pars_start,
#' data_init = data_initial, niter = 1e3, BurnIn = 1e2,
#' model = "Multinomial", N = N, lambda = 1/(69.416 * 365),
#' mu = 1/(69.416 * 365), period_start = phases, opt_num = 1,
#' auto.initialize = TRUE, f = 0.15)
#'
#' names(cov19est)
#' class(cov19est$mcmc_pars)
#' names(cov19est$plots)
#'
#' plot(cov19est, type = "trace")
#' plot(cov19est, type = "boxplot")
#'
#'
SEIRfansy <- function(data, data_init,N, init_pars = NULL, niter = 1e5, BurnIn = 5e3,
model = "Multinomial", plot = TRUE, period_start, auto.initialize = TRUE,
... ){
#library(dplyr)
# library(arm)
# library(pbapply)
# library(DescTools)
# library("patchwork")
if(model == "Multinomial"){
if(! all(c("Confirmed", "Recovered", "Deceased" ) %in% colnames(data))){
stop("Data does not contain required columns")
}
} else if (model == "Poisson" || model == "Binomial"){
if(! all(c("Confirmed") %in% colnames(data))){
stop("Data does not contain required columns")
}
} else
stop("Incorrect model specified")
if(length(data_init) != 6) stop("Incorrect initial values: length must be 6")
if(any(data < 0) || data_init < 0) stop("Negative values in data")
var_init = model_initializeR(data = data, data_init = data_init, init_pars = init_pars, model = model,
period_start = period_start, auto.initialize = auto.initialize, ...)
fix_pars = var_init$fix_pars
init_state_num = var_init$init_state_num
init_pars = var_init$init_pars
period_start = var_init$period_start
pars_estimate = mcmc_performR(data = data, init_state_num = init_state_num, init_pars = init_pars, model = model,
niter = niter, BurnIn = BurnIn, fix_pars = fix_pars, period_start = period_start, ...)
## Plot
if(plot == TRUE){
plots <- model_plotR(mcmc_pars = pars_estimate, n_period = length(period_start), data = data,
data_initial = data_initial, call = "estimateR", model = model, ...)
rlist = list(mcmc_pars = pars_estimate, "plots" = plots)
class(rlist) = "SEIRfansy"
return(rlist)
}
rlist = list(mcmc_pars = pars_estimate)
class(rlist) = "SEIRfansy"
return(rlist)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.