ci.curve {REffectivePred}

Confidence bands

Description

Computes the pointwise confidence interval of the epidemic curve.

Usage

ci.curve(
  fit = NULL,
  H.E = NULL,
  H.W = NULL,
  scenario = NULL,
  cases = NULL,
  cfg = NULL,
  restrictions = NULL,
  restriction.starts = NULL,
  ranges = NULL,
  rt_func = 1,
  fit.t.pred = NULL,
  predict.beyond = 0,
  lt = NULL,
  adj.period = NULL,
  population = NULL,
  rho = NULL,
  serial_mean = NULL,
  serial_var = NULL,
  window_size = NULL,
  eps = .Machine$double.eps^(1/2)
)

Arguments

fit

Output from function estimate.mle.

H.E

Mobility metrics for category Retail & Entertainment. Currently unsupported.

H.W

Mobility metrics for category Workplaces. Currently unsupported.

scenario

A character string describing options to deal with restrictions. Currently unsupported.

cases

vector of case counts.

cfg

The object that contains all variables from the configuration file. fit, H.E, H.W, scenario, and cases are also required for the method to execute. All other parameters will not be used if cfg is passed to the method.

restrictions

A numeric integer vector giving the severity of restrictions. Zero means no restriction, and higher numbers means greater severity/disruption. The ordered unique values should be consecutive integers starting from zero. Each number (other than 0) adds a new parameter to the fit.

restriction.starts

A vector of same length as restrictions, of times when restrictions came into effect. Note: the first index time should be 1.

ranges

A vector of time ranges for the different waves. The wave ranges should be contiguous, with at least one unit of time between consecutive waves.

rt_func

The parametric form of function c(). Options are listed under function c_helper.

fit.t.pred

Time from which prediction is done. If use.actual.not.predicted is TRUE, values of S_t before this time will be computed using actual counts.

predict.beyond

Number of days to predict beyond the end of cases. See Details for usage notes.

lt

The length of cases.

adj.period

Adjustment period following a change in severity level. Restriction level (psi) is linearly interpolated from the old to the new value over this period.

population

Total population size.

rho

A vector of under-reporting rates of the same length as cases. If a scalar is supplied, the vector will be constant with this value.

serial_mean

Mean of the serial interval on the log scale.

serial_var

Variance of the serial interval on the log scale.

window_size

The maximum value for the serial interval.

eps

The epsilon value for computing finite differences.

Value

Returns a matrix with two rows containing Wald-style confidence bounds:

• ci_lower - lower bound of a 95% pointwise CI of the best fit curve.

• ci_upper - upper bound of a 95% pointwise CI of the best fit curve.

Examples

library(REffectivePred)
## Read in the data
path_to_data <- system.file("extdata/NY_OCT_4_2022.csv", package = "REffectivePred")
data <- read.csv(path_to_data)
head(data)
cases <- diff(c(0, data$cases)) # Convert cumulative cases into daily cases
lt <- length(cases)             # Length of cases
Time <- as.Date(data$date, tryFormats = c("%d-%m-%Y", "%d/%m/%Y"))

navigate_to_config() # Open the config file, make any necessary changes here.
path_to_config <- system.file("config.yml", package = "REffectivePred")  # Read config file
cfg <- load_config()    # Build the cfg object

# Estimate parameters
est <- estimate.mle(
    cases = cases,
    cfg = cfg,
    hessian = TRUE
    )
a1 <- est$a1
a2 <- est$a2
a3 <- est$a3
a4 <- est$a4
nu <- est$nu
vt <- c(1, est$vt_params_est)
psi <- est$Psi
betas <- est$betas

# Predict curve
r1 <- pred.curve(
a1 = a1,
a2 = a2,
a3 = a3,
a4 = a4,
nu = nu,
variant.transm = vt,
Psi = psi,
betas = betas,
cases = cases,
cfg = cfg
)

plot_outputs(Time = Time,
cases = cases,
cfg = cfg,
curve = r1,
option = 2
)

bounds <- ci.curve(fit = est,
                   cases = cases,
                   cfg = cfg)

# Adding CI bands
# lines(c(Time, Time[length(Time)]+(1:predict.beyond)), bounds[2,], lty = 2)
# lines(c(Time, Time[length(Time)]+(1:predict.beyond)), bounds[1,], lty = 2)