Package 'ggoutbreak'

Description

Time periods are just a zero based numeric representation of dates with a time unit baked in. This allows variable length periods (e.g. days or weeks), and fractional days to be represented in a consistent(ish) way between things that want to deal in dates (like ggplot) and things that want to deal in numbers (like model fitting)

Usage

as.time_period(x, ...)

## S3 method for class 'time_period'
as.time_period(x, unit = NULL, start_date = NULL, ...)

## S3 method for class 'Date'
as.time_period(x, unit = NULL, anchor = NULL, ...)

## S3 method for class 'numeric'
as.time_period(x, unit = NULL, start_date = NULL, ...)

## S3 method for class 'grates_epiweek'
as.time_period(x, ...)

## S3 method for class 'grates_isoweek'
as.time_period(x, ...)

## S3 method for class 'grates_period'
as.time_period(x, ...)

## S3 method for class 'time_period'
seq(from, to = from, ...)

is.time_period(x)

date_to_time(dates, unit = NULL, start_date = NULL)

time_to_date(
  timepoints,
  unit = attr(timepoints, "unit"),
  start_date = attr(timepoints, "start_date")
)

Arguments

Value

a vector of class time_period

Functions

• seq(time_period): Create a sequence using time_periods

• is.time_period(): Check is a time_period

• date_to_time(): Convert a set of dates to numeric timepoints

• time_to_date(): Convert a set of time points to dates

Unit tests

tmp = as.time_period(grates::as_epiweek("2019-W12", format = "yearweek")+0:2)
testthat::expect_equal(
  lubridate::epiweek(as.Date(tmp)),
  c(12, 13, 14)
)

tmp = as.time_period(grates::as_isoweek("2019-W12", format = "yearweek")+0:2)
testthat::expect_equal(
  lubridate::isoweek(as.Date(tmp)),
  c(12, 13, 14)
)

x = as.time_period(as.Date("2025-01-01")+0:2, anchor="start", unit="1 day")
y = as.time_period(as.Date("2025-01-07")+0:2, anchor="start", unit="1 day")
testthat::expect_equal(as.numeric(c(x, y)), c(0, 1, 2, 6, 7, 8))
testthat::expect_equal(as.numeric(c(y, x)), c(0, 1, 2, -6, -5, -4))

z = as.time_period(as.Date("2025-01-01")+0:2*7, anchor="start", unit="1 week")
testthat::expect_equal(as.numeric(c(x, z)), c(0, 1, 2, 0, 7, 14))
testthat::expect_equal(
  as.numeric(c(z, x)),
  c(0, 1, 2, 0, 0.142857142857143, 0.285714285714286)
)

testthat::expect_equal(
  as.Date(c(y, z)),
  structure(c(20095, 20096, 20097, 20089, 20096, 20103), class = "Date")
)

seq(x[[1]], y[[1]])

Examples

#' # 100 weeks from 2020-01-01

tmp = as.time_period(0:100, 7, "2020-01-01")
as.Date(tmp)

range(tmp)
min(tmp)
tmp2 = as.integer(as.Date(tmp))
# testthat::expect_true(all(na.omit(tmp2-lag(tmp2)) == 7))

tmp2 = as.time_period(0:23, 1/24, "2020-01-01")
as.POSIXct(tmp2)

# convert timeseries to new "unit"
tmp = as.time_period(0:100, 7, "2020-01-01")
tmp2 = as.time_period(tmp,1)
testthat::expect_equal(as.numeric(tmp2), 0:100*7)

# 100 weeks from 2020-01-01

tmp = as.time_period(0:100, 7, "2020-01-01")
as.Date(tmp)

range(tmp)
min(tmp)
tmp2 = as.integer(as.Date(tmp))
# testthat::expect_true(all(na.omit(tmp2-lag(tmp2)) == 7))

tmp2 = as.time_period(0:23, 1/24, "2020-01-01")
as.POSIXct(tmp2)

# convert timeseries to new "unit"
tmp = as.time_period(0:100, 7, "2020-01-01")
tmp2 = as.time_period(tmp,1)
testthat::expect_equal(as.numeric(tmp2), 0:100*7)

# Time to date
times = date_to_time(as.Date("2019-12-29")+0:100, "1 week")
dates = time_to_date(times)

# Date to time
times = date_to_time(as.Date("2019-12-29")+0:100, "1 week")
dates = time_to_date(times)

Description

A scales breaks generator for log1p scales

Usage

breaks_log1p(n = 5, base = 10)

Arguments

Value

a function for ggplot scale breaks

Examples

ggplot2::ggplot(ggplot2::diamonds, ggplot2::aes(x=price))+
  ggplot2::geom_density()+
  ggplot2::scale_x_continuous(trans="log1p", breaks=breaks_log1p())

Description

Generate a random probability based on features of the simulation

Usage

cfg_beta_prob_rng(probability_fn = ~0.8, kappa_fn = ~0.1)

Arguments

Value

a time dependent function that inputs a time (as a time_period) and returns an probability for each day defined by the probability_fn and kappa_fn

See Also

rbeta2()

Examples

fn = cfg_beta_prob_rng(~ ifelse(.x<=5,0.1,0.9))
fn(1:10)

Description

Get a IP generating function from time varying mean and SD of a gamma function

Usage

cfg_gamma_ip_fn(mean_fn = ~2, sd_fn = function(mean) sqrt(mean))

Arguments

Value

a time dependent function that inputs a time (as a time_period) and returns an ip delay distribution for each day defined by the mean_fn and sd_fn

Examples

fn = cfg_gamma_ip_fn(mean_fn = function(t) ifelse(t < 5, 4, 2))
# a gamma function that changes mean at time 5
fn(4)
fn(7)

Description

This is used for random sampling from the infectivity profile for times to infection, for example. There is nothing to stop you putting in a delay distribution with negative times but strange things may happen in your simulation.

Usage

cfg_ip_sampler_rng(ip = i_empirical_ip)

Arguments

Value

a function which accepts n parameter which produces random samples from the ip distribution

Examples

tmp = cfg_ip_sampler_rng(example_ganyani_ip())(10000)

# This discretised ganyani distribution is based on these figures:
# mean: 5.2 (3.78-6.78) and sd: 1.72 (0.91-3.93)
format_ip(example_ganyani_ip())

mean(tmp) # Should be about 5.2
stats::sd(tmp) # Should be about 1.72

Description

Linear function from dataframe

Usage

cfg_linear_fn(changes, ..., col_name = setdiff(colnames(changes), "t"))

Arguments

Value

a function that inputs a vector t and returns a linearly interpolated value from ⁠<col_name>⁠

Description

Step function from dataframe

Usage

cfg_step_fn(changes, ..., col_name = setdiff(colnames(changes), "t"))

Arguments

Value

a function that inputs a vector t and returns the next smallest corresponding value in ⁠<col_name>⁠ (or the first one)

Description

This is particularly designed for use within the fn_list_next_gen parameter of sim_branching_process() to allow age group contract matrices to be applied for example (assuming a categorical age).

Usage

cfg_transition_fn(transition)

Arguments

Value

a function that given an input will return samples of the output class according to the probability distributions.

Examples

age = rep(c("child","adult","elderly"),100)

fn = cfg_transition_fn(dplyr::tribble(
  ~input, ~output, ~probability,
  "child", "child", 0.5,
  "child", "adult", 0.5,
  "adult", "child", 0.5,
  "adult", "adult", 0.3,
  "adult", "elderly", 0.2,
  "elderly","elderly", 0.5,
  "elderly","adult", 0.5,
))

table(fn(age),age)

Description

This function returns a random number generator from a gamma distribution that has a weekly period and configurable degree of additional variability around weekly pattern. It is suited to delays of things like testing which may depend on the day of week.

Usage

cfg_weekly_gamma_rng(
  mean = c(1, 1, 1, 1, 4, 3, 2),
  sd = sqrt(mean),
  week_starts = weekdays(as.Date("2024-10-14"))
)

Arguments

Value

a random number generator taking t time parameter and returning a duration the that time t depending on the day of the week.

Examples

fn = cfg_weekly_gamma_rng(c(1,1,1,1,4,3,2))
matrix(fn(1:42),ncol=7,byrow=TRUE)

Description

Weekly convolution distribution function

Usage

cfg_weekly_ip_fn(
  mean = c(1, 1, 1, 1, 4, 3, 2),
  sd = sqrt(mean),
  week_starts = weekdays(as.Date("2024-10-14"))
)

Arguments

Value

a time dependent function that inputs a time (as a time_period) and generates an IP delay distribution for each day varying by day of week

Examples

cat(sapply(cfg_weekly_ip_fn()(1:7),format_ip),sep = "\n")

Description

This function returns a random probability generator that has a weekly period and configurable degree of additional variability around weekly pattern. It is suited to probabilities of things like testing which may depend on the day of week.

Usage

cfg_weekly_proportion_rng(
  prob = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.5, 0.5),
  kappa = 0.1,
  week_starts = weekdays(as.Date("2024-10-14"))
)

Arguments

Value

a random number generator function taking t time parameter and returning a probability of ascertainment for that time, depending on day of week etc.

Examples

fn = cfg_weekly_proportion_rng(c(0.9,0.9,0.9,0.9,0.9,0.1,0.1))
matrix(fn(1:42),ncol=7,byrow=TRUE)

Description

The counterpart to date_seq_dates(). Take an original set of data and place it within a regular time series where the periodicity of the time series may be expressed as numbers of days, weeks, months quarters, or years, and the periods are defined by an anchoring date, day of the week or by reference to the start or end of the input dates. This can either return the periods as dates or factors (e.g. for plotting) or as a time_period for analysis that relies on a numeric representation of the date or duration from the anchor.

Usage

cut_date(
  dates,
  unit,
  anchor = "start",
  output = c("date", "factor", "time_period"),
  dfmt = "%d/%b/%y",
  ifmt = "{start} — {end}",
  ...
)

Arguments

Value

a set of dates, times or a factor level, representing the start of the period the date falls into, where the period is defined by the duration and the anchor

Examples

dates = as.Date(c("2020-01-01","2020-02-01","2020-01-15","2020-02-03",NA))
fs = date_seq(dates, "2 days")
dates - cut_date(dates, "2 days")
cut_date(dates,unit="2 days", output="time_period")

# A weekly set of dates:
dates2 = Sys.Date() + floor(stats::runif(50,max=10))*7

# in this specific situation the final date is not truncated because the
# input data is seen as an exact match for the whole output period.
cut_date(dates2, "1 week", "sun", output="factor")
cut_date(dates2, dfmt = "%d/%b", output="factor", unit = "2 weeks", anchor="sun")

Description

This is useful if you want to fill in missing values that should have been observed but weren't. For example, date_seq(c(1, 2, 4, 6), 1) will return 1:6.

Usage

date_seq(x, period, ...)

Arguments

Value

a vector of the same type as the input

Examples

date_seq(c(1, 2, 4, 5, 10), 1)

Description

Derive from a vector of observation dates, a complete ordered sequence of periods in a regular time series, where the length of the periods is specified, as a number of days, weeks, years etcetera. E.g. this can convert a random set of dates to a ordered complete list of 1 week intervals (or 2 month intervals) spanning the same range as the dates. This has some interesting problems regarding where to put breaks within a month or week. Often this is either based on a specific date (e.g. yearly periods starting at 2020-01-01) or a day of week (e.g. 2 weekly periods staring on a Sunday) or maybe relative to the input time series (weekly ending on the last date of the data). There is also a problem when we consider data that may have incomplete starting and end periods, which may not be comparable to other periods, and we may need to exclude these from the result.

Usage

## S3 method for class 'Date'
date_seq(x, period = .day_interval(x), anchor = "start", complete = FALSE, ...)

Arguments

Value

a vector of dates for regular periods between the minimum and maximum of dates, with the boundaries defined by the anchor.

Examples

date_seq(as.Date(c("2020-01-01","2020-02-01","2020-01-15","2020-02-01",NA)), "2 days")

Description

This is useful if you want to fill in missing values that should have been observed but weren't. For example, date_seq(c(1, 2, 4, 6), 1) will return 1:6.

Usage

## S3 method for class 'numeric'
date_seq(x, period = 1, tol = 1e-06, ...)

Arguments

Value

a vector of the same type as the input

Examples

date_seq(c(1, 2, 4, 5, 10), 1)

Description

Derive from a vector of observation time_periods, a complete ordered sequence of periods in a regular time series, where the length of the periods is specified, as a number of days, weeks, years etc. E.g. this can convert a random set of times to a ordered complete list of 1 week intervals (or 2 month intervals) spanning the same range as the dates. This has some interesting problems regarding where to put breaks within a month or week. Often this is either based on a specific date (e.g. yearly periods starting at 2020-01-01) or a day of week (e.g. 2 weekly periods staring on a sunday) or maybe relative to the input time series (weekly ending on the last date of the data). There is also a problem when we consider data that may have incomplete starting and end periods, which may not be comparable to other periods, and we may need to exclude these from the result.

Usage

## S3 method for class 'time_period'
date_seq(x, period = attributes(x)$unit, complete = FALSE, ...)

Arguments

Value

a vector of time_periods for regular periods between the minimum and maximum of dates, with the boundaries defined by the anchor.

Examples

tmp = as.time_period(c(0,10,100), 7, "2020-01-01")
date_seq(tmp, "7 days")
date_seq(tmp, "1 day")

Description

Density, distribution function, quantile function and random generation for the Beta distribution with parameters shape1 and shape2 (and optional non-centrality parameter ncp).

Usage

dbeta2(x, prob, kappa, log = FALSE)

Arguments

Value

dbeta gives the density, pbeta the distribution function, qbeta the quantile function, and rbeta generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rbeta, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::dbeta()

Examples

dbeta2(c(0.25,0.5,0.75), 0.5, 0.25)

Description

The following conversion describes the parameters mean and kappa

Usage

dcgamma(x, mean, kappa = 1/mean, log = FALSE)

Arguments

Details

$\text{shape: } \alpha = \frac{1}{\kappa} \\ \text{rate: } \beta = \frac{1}{\mu \times \kappa} \\ \text{scale: } \sigma = \mu \times \kappa \\$

Value

dgamma gives the density, pgamma gives the distribution function, qgamma gives the quantile function, and rgamma generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgamma, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::dgamma()

Examples

dcgamma(seq(0,4,0.25), 2, 0.5)

Description

Density, distribution function, quantile function and random generation for the Gamma distribution with parameters shape and scale.

Usage

dgamma2(x, mean, sd = sqrt(mean), log = FALSE)

Arguments

Value

dgamma gives the density, pgamma gives the distribution function, qgamma gives the quantile function, and rgamma generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgamma, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::dgamma()

Examples

dgamma2(seq(0,4,0.25), 2, 1)

Description

Density, distribution function, quantile function and random generation for the log normal distribution whose logarithm has mean equal to meanlog and standard deviation equal to sdlog.

Usage

dlnorm2(x, mean = 1, sd = sqrt(exp(1) - 1), log = FALSE)

Arguments

Details

The log normal distribution has density

$f(x) = \frac{1}{\sqrt{2\pi}\sigma x} e^{-(\log(x) - \mu)^2/2 \sigma^2}$

where $\mu$ and $\sigma$ are the mean and standard deviation of the logarithm. The mean is $E(X) = exp(\mu + 1/2 \sigma^2)$, the median is $med(X) = exp(\mu)$, and the variance $Var(X) = exp(2\mu + \sigma^2)(exp(\sigma^2) - 1)$ and hence the coefficient of variation is $\sqrt{exp(\sigma^2) - 1}$ which is approximately $\sigma$ when that is small (e.g., $\sigma < 1/2$).

Value

dlnorm gives the density, plnorm gives the distribution function, qlnorm gives the quantile function, and rlnorm generates random deviates.

The length of the result is determined by n for rlnorm, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

Note

The cumulative hazard $H(t) = - \log(1 - F(t))$ is -plnorm(t, r, lower = FALSE, log = TRUE).

Source

dlnorm is calculated from the definition (in ‘Details’). [pqr]lnorm are based on the relationship to the normal.

Consequently, they model a single point mass at exp(meanlog) for the boundary case sdlog = 0.

References

Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole.

Johnson, N. L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions, volume 1, chapter 14. Wiley, New York.

See Also

Distributions for other standard distributions, including dnorm for the normal distribution.

Examples

dlnorm2(seq(0,4,0.25), 2, 1)

Description

The logit-normal distribution has a support of 0 to 1.

Usage

dlogitnorm(x, meanlogit = 0, sdlogit = 1, log = FALSE)

Arguments

Value

a vector of probabilities, quantiles, densities or samples.

Examples

dlogitnorm(seq(0.1,0.9,0.1), 0, 1)

Description

The logit-normal distribution has a support of 0 to 1.

Usage

dlogitnorm2(x, prob.0.5 = 0.5, kappa = 1 - exp(-1), log = FALSE)

Arguments

Value

a vector of probabilities, quantiles, densities or samples.

Examples

q = seq(0.1,0.9,0.1)
eps = sqrt(.Machine$double.eps)
dlogitnorm2(q, 0.75, 0.2)

(plogitnorm2(q+eps, 0.75, 0.2) -
  plogitnorm2(q-eps, 0.75, 0.2)) / (2*eps)

Description

Density, distribution function, quantile function and random generation for the negative binomial distribution with parameters size and prob.

Usage

dnbinom2(x, mean, sd = sqrt(mean), log = FALSE)

Arguments

Details

The negative binomial distribution with size $= n$ and prob $= p$ has density

$p(x) = \frac{\Gamma(x+n)}{\Gamma(n) x!} p^n (1-p)^x$

for $x = 0, 1, 2, \ldots$, $n > 0$ and $0 < p \le 1$.

This represents the number of failures which occur in a sequence of Bernoulli trials before a target number of successes is reached. The mean is $\mu = n(1-p)/p$ and variance $n(1-p)/p^2$.

A negative binomial distribution can also arise as a mixture of Poisson distributions with mean distributed as a gamma distribution (see pgamma) with scale parameter (1 - prob)/prob and shape parameter size. (This definition allows non-integer values of size.)

An alternative parametrization (often used in ecology) is by the mean mu (see above), and size, the dispersion parameter, where prob = size/(size+mu). The variance is mu + mu^2/size in this parametrization.

If an element of x is not integer, the result of dnbinom is zero, with a warning.

The case size == 0 is the distribution concentrated at zero. This is the limiting distribution for size approaching zero, even if mu rather than prob is held constant. Notice though, that the mean of the limit distribution is 0, whatever the value of mu.

The quantile is defined as the smallest value $x$ such that $F(x) \ge p$, where $F$ is the distribution function.

Value

dnbinom gives the density, pnbinom gives the distribution function, qnbinom gives the quantile function, and rnbinom generates random deviates.

Invalid size or prob will result in return value NaN, with a warning.

The length of the result is determined by n for rnbinom, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

rnbinom returns a vector of type integer unless generated values exceed the maximum representable integer when double values are returned.

Source

dnbinom computes via binomial probabilities, using code contributed by Catherine Loader (see dbinom).

pnbinom uses pbeta.

qnbinom uses the Cornish–Fisher Expansion to include a skewness correction to a normal approximation, followed by a search.

rnbinom uses the derivation as a gamma mixture of Poisson distributions, see

Devroye, L. (1986) Non-Uniform Random Variate Generation. Springer-Verlag, New York. Page 480.

See Also

Distributions for standard distributions, including dbinom for the binomial, dpois for the Poisson and dgeom for the geometric distribution, which is a special case of the negative binomial.

Examples

dnbinom2(0:5, 2, sqrt(2))

Description

Null distributions always returns NA

Usage

dnull(x, ...)

Arguments

Examples

dnull(c(0.25,0.5,0.75), 0.5, 0.25)

Description

The unit of doubling times is always days.

Usage

doubling_time(x, ...)

Arguments

Value

the same dataframe with additional columns for doubling time or relative doubling time plus confidence intervals.

Examples

example_poisson_rt_smooth() %>%
  poisson_locfit_model(window=21) %>%
  doubling_time() %>%
  dplyr::glimpse()

Description

The wedge distribution has a support of 0 to 1 and has a linear probability density function over that support.

Usage

dwedge(x, a, lower.tail = TRUE, log = FALSE)

Arguments

Details

The rwedge can be combined with quantile functions to skew standard distributions, or introduce correlation or down weight certain parts of the distribution.

Value

a vector of probabilities, quantiles, densities or samples.

Examples

pwedge(seq(0,1,0.1), a=1)
dwedge(seq(0,1,0.1), a=1)
qwedge(c(0.25,0.5,0.75), a=-1)

stats::cor(
  stats::qnorm(rwedge(1000, a=2)),
  stats::qnorm(rwedge(1000, a=-2))
)

Description

Format date as dmy

Usage

fdmy(date)

Arguments

Value

the formatted date

Examples

fdmy(Sys.Date())

Description

Print a summary of an infectivity profile

Usage

format_ip(ip = i_empirical_ip)

Arguments

Value

an infectivity profile description

Examples

format_ip(example_ganyani_ip())

Description

Delayed GAM reporting model function generator

Usage

gam_delayed_reporting(
  window,
  max_delay = 40,
  ...,
  knots_fn = ~gam_knots(.x, window, ...)
)

Arguments

Details

This function is used to configure a delayed reporting GAM model. The model is of the form:

count ~ s(time, bs = "cr", k = length(kts)) + s(log(tau), k = 4, pc = max_delay)

where tau is the difference between time series observation time and the time of the data point in the time series, and we have multiple observations of the same time series. This function helps specify the knots of the GAM and the maximum expected delay

Value

a list with 2 entries - model_fn and predict suitable as the input for poisson_gam_model(model_fn = ..., predict=...).

Examples

data = example_delayed_observation() %>% dplyr::group_by(obs_time)
cfg = gam_delayed_reporting(14,40)
fit = cfg$model_fn(data)
summary(fit)

Description

At the moment this does nothing sophisticated. An mostly equally spaced grid of knots with gaps at start and end to prevent over-fitting there. In the future this could look at the number of observations or areas where there is a lot of change to add in more knots.

Usage

gam_knots(data = i_incidence_data, window, ..., k = NULL)

Arguments

Value

a vector of times (as a numeric)

Examples

gam_knots(example_poisson_rt(), 14)
gam_knots(example_poisson_rt(), k=10)

Description

This function configures a very simple negative binomial count model and using:

Usage

gam_nb_model_fn(window, ..., knots_fn = ~gam_knots(.x, window, ...))

Arguments

Details

count ~ s(time, bs = "cr", k = length(kts))

The control of knot positions is defined by the knots_fn and window parameters.

Value

a function suitable as the input for poisson_gam_model(model_fn = ...).

Examples

model_fn = gam_nb_model_fn(14)
fit = model_fn(example_poisson_rt() %>% dplyr::ungroup())
summary(fit)

Description

This function configures a very simple poisson count model and using:

Usage

gam_poisson_model_fn(window, ..., knots_fn = ~gam_knots(.x, window, ...))

Arguments

Details

count ~ s(time, bs = "cr", k = length(kts))

The control of knot positions is defined by the knots_fn and window parameters.

Value

a function suitable as the input for poisson_gam_model(model_fn = ...).

Examples

model_fn = gam_poisson_model_fn(14)
fit = model_fn(example_poisson_rt() %>% dplyr::ungroup())
summary(fit)

Description

The x axis must be a date.

Usage

geom_events(
  events = i_events,
  event_label_size = 7,
  event_label_colour = "black",
  event_label_angle = -30,
  event_line_colour = "grey50",
  event_fill_colour = "grey50",
  hide_labels = FALSE,
  guide_axis = ggplot2::derive(),
  x_axis_style = c("date", "time_period"),
  ...
)

Arguments

Value

a set of geoms for a time series.

Description

This assumes a modelled incidence estimate that is log-normal. The exponential growth rate is the first derivative of the mu parameters of this log-normal. On the link scale these are normally distributed. This function assumes that the time series incidence estimates are uncorrelated to estimate the error in the growth rate, which is a conservative approach resulting in more uncertainty in growth rate than might be possible through other methods. This is all based on Savitsky-Golay filters applied to the normally distributed log-incidence estimates.

Usage

growth_rate_from_incidence(d = i_incidence_model, window = 7, deg = 2)

Arguments

Value

the timeseries with growth rate columns: A dataframe containing the following columns:

• time (ggoutbreak::time_period + group_unique) - A (usually complete) set of singular observations per unit time as a time_period

• growth.fit (double) - an estimate of the growth rate

• growth.se.fit (positive_double) - the standard error the growth rate

• growth.0.025 (double) - lower confidence limit of the growth rate

• growth.0.5 (double) - median estimate of the growth rate

• growth.0.975 (double) - upper confidence limit of the growth rate

Any grouping allowed.

Examples

data = example_poisson_growth_rate()
tmp2 = data %>%
  poisson_glm_model(window=7,deg=1) %>%
  growth_rate_from_incidence(window = 13, deg=1)

if(interactive()) {
  plot_incidence(tmp2,
      date_labels="%b %y")

  plot_growth_rate(
      tmp2,
      date_labels="%b %y"
    )+
    sim_geom_function(data,colour="red")
}

Description

This assumes a prevalence that is logit-normally distributed. The exponential growth rate is the first derivative of the mu parameters of this logit-normal. On the link scale these are normally distributed. This function assumes that the time series incidence estimates are uncorrelated to estimate the error in the growth rate, which is a conservative approach resulting in more uncertainty in growth rate than might be possible through other methods. This is all based on Savitsky-Golay filters applied to the normally distributed logit-proportion estimates.

Usage

growth_rate_from_prevalence(d = i_prevalence_model, window = 7, deg = 2)

Arguments

Value

the timeseries with growth rate columns: A dataframe containing the following columns:

• time (ggoutbreak::time_period + group_unique) - A (usually complete) set of singular observations per unit time as a time_period

• proportion.fit (double) - an estimate of the proportion on a logit scale

• proportion.se.fit (positive_double) - the standard error of proportion estimate on a logit scale

• proportion.0.025 (proportion) - lower confidence limit of proportion (true scale)

• proportion.0.5 (proportion) - median estimate of proportion (true scale)

• proportion.0.975 (proportion) - upper confidence limit of proportion (true scale)

• relative.growth.fit (double) - an estimate of the relative growth rate

• relative.growth.se.fit (positive_double) - the standard error the relative growth rate

• relative.growth.0.025 (double) - lower confidence limit of the relative growth rate

• relative.growth.0.5 (double) - median estimate of the relative growth rate

• relative.growth.0.975 (double) - upper confidence limit of the relative growth rate

Any grouping allowed.

Examples

data = example_poisson_rt_2class()

tmp = data %>%
  proportion_glm_model(window=7,deg=2) %>%
  dplyr::select(time, dplyr::starts_with("proportion")) %>%
  dplyr::rename_with(~ gsub("proportion","prevalence",.x)) %>%
  dplyr::select(-prevalence.fit, -prevalence.se.fit)

tmp = tmp %>%
  growth_rate_from_prevalence(window = 25, deg=1)

plot_growth_rate(
      tmp,
      date_labels="%b %y"
  )

data1 = ukc19::ons_infection_survey %>%
  dplyr::filter(name=="England") %>%
  timeseries(count=FALSE)

tmp2 = data1 %>%
  growth_rate_from_prevalence(window = 25, deg=1)

if(interactive()) {
  plot_growth_rate(
      tmp2,
      date_labels="%b %y",
      mapping = ggplot2::aes(colour=name),
      events = ukc19::timeline
  )
}

Description

This assumes a modelled proportion that is logit-normally distributed. The exponential growth rate is the first derivative of the mu parameters of this logit-normal. On the link scale these are normally distributed. This function assumes that the time series incidence estimates are uncorrelated to estimate the error in the growth rate, which is a conservative approach resulting in more uncertainty in growth rate than might be possible through other methods. This is all based on Savitsky-Golay filters applied to the normally distributed logit-proportion estimates.

Usage

growth_rate_from_proportion(d = i_proportion_model, window = 7, deg = 2)

Arguments

Value

the timeseries with growth rate columns: A dataframe containing the following columns:

• time (ggoutbreak::time_period + group_unique) - A (usually complete) set of singular observations per unit time as a time_period

• proportion.fit (double) - an estimate of the proportion on a logit scale

• proportion.se.fit (positive_double) - the standard error of proportion estimate on a logit scale

• proportion.0.025 (proportion) - lower confidence limit of proportion (true scale)

• proportion.0.5 (proportion) - median estimate of proportion (true scale)

• proportion.0.975 (proportion) - upper confidence limit of proportion (true scale)

• relative.growth.fit (double) - an estimate of the relative growth rate

• relative.growth.se.fit (positive_double) - the standard error the relative growth rate

• relative.growth.0.025 (double) - lower confidence limit of the relative growth rate

• relative.growth.0.5 (double) - median estimate of the relative growth rate

• relative.growth.0.975 (double) - upper confidence limit of the relative growth rate

Any grouping allowed.

Examples

data = example_poisson_rt_2class()

tmp2 = data %>%
  proportion_glm_model(window=7,deg=2) %>%
  growth_rate_from_proportion(window = 13, deg=1)

if(interactive()) {
  plot_proportion(tmp2,
      date_labels="%b %y")

  plot_growth_rate(
      tmp2,
      date_labels="%b %y"
    )
}

Description

This function augments any timeseries with a population denominator. The population data may be static estimates, or a set of estimates at time points. The population data may be grouped in which case the grouping might be geographical area or age group or gender for example. The two inputs must have compatible grouping (i.e. all the groups in the population data must be present in the timeseries).

Usage

infer_population(df = i_timeseries, pop = i_population_data)

Arguments

Value

the df timeseries with additional population column

Examples

# The COVID data already has a population column so we are just double checking:
data = example_england_covid_by_age() %>%
  dplyr::rename(pop_old = population)

demog = ukc19::uk_population_2019_by_5yr_age %>% dplyr::group_by(code, name, class)

data %>%
  infer_population(demog) %>%
  dplyr::glimpse()

Description

Log-scaled incidence estimates are used to generate samples of incidence. These are convolved with the duration of infection (as the probability still infected) on any given day (as a set of discrete distributions) and the result is evaluated as a fraction of population. The result is fitted to a logit-normal (using moments) and quantiles produced from samples. If test sensitivity is used instead of duration of infection then the prevalence estimate will be the expected proportion of test positives by screening test rather than a true prevalence. Differences between this and the proportion of test positives observed may be due to ascertainment bias or test sensitivity misspecification.

Usage

infer_prevalence(
  modelled = i_incidence_model,
  pop = i_population_data,
  ip = i_discrete_ip,
  bootstraps = 1000,
  seed = Sys.time(),
  ...
)

Arguments

Value

the modelled input with additional proportion columns

Examples

tmp = example_poisson_rt_smooth() %>%
  poisson_locfit_model(window=14) %>%
  infer_prevalence(
     pop = 10000,
     ip = example_ip()
  )

if(interactive()) {
  plot_prevalence(tmp)
}

Description

This enables incidence rates are able to be compared to a baseline figure for incidence. The baseline could come for example from a population average or average incidence over time. The output is an incidence rate ratio. The incidence_baseline column is a rate of events per unit time. The time unit is expected to be the same as that of the date in modelled and this is not checked.

Usage

infer_rate_ratio(
  modelled = i_incidence_model,
  base = i_baseline_incidence_data,
  ...
)

Arguments

Value

a dataframe with incidence rate ratios for each of the classes in modelled. A dataframe containing the following columns:

• time (ggoutbreak::time_period + group_unique) - A (usually complete) set of singular observations per unit time as a time_period

• rate_ratio.0.025 (positive_double) - lower confidence limit of the rate ratio for a population group

• rate_ratio.0.5 (positive_double) - median estimate of the rate ratio for a population group

• rate_ratio.0.975 (positive_double) - upper confidence limit of the rate ratio for a population group

Any grouping allowed.

Examples

# not age stratified
baseline = example_poisson_locfit() %>%
  dplyr::mutate(baseline_incidence = incidence.0.5)

# age stratified rate ratios:
tmp = example_poisson_age_stratified() %>%
  infer_rate_ratio(baseline) %>%
  dplyr::glimpse()

Description

This assumes that for example, case distribution proportions are stratified by a population grouping, e.g. geography or age, and we have estimates of the size of that population during that time period. Normalising by population proportion allows us to compare the relative risk of the outcome in groups, compared to the expected population risk if the outcome is evenly distributed across the population. There may be other proportions other than population fractions that may be useful to compare. At the moment this does not handle any uncertainty.

Usage

infer_risk_ratio(
  modelled = i_proportion_model,
  base = i_baseline_proportion_data,
  ...
)

Arguments

Value

a dataframe with relative risk / risk ratio columns. A dataframe containing the following columns:

• time (ggoutbreak::time_period + group_unique) - A (usually complete) set of singular observations per unit time as a time_period

• risk_ratio.0.025 (positive_double) - lower confidence limit of the excess risk ratio for a population group

• risk_ratio.0.5 (positive_double) - median estimate of the excess risk ratio for a population group

• risk_ratio.0.975 (positive_double) - upper confidence limit of the excess risk ratio for a population group

Any grouping allowed.

Examples

demog = ukc19::uk_population_2019_by_5yr_age %>%
  dplyr::filter(name == "England")


tmp = example_proportion_age_stratified() %>%
  infer_risk_ratio(demog) %>%
  dplyr::glimpse()

if(interactive()) {
  plot_growth_phase(tmp,duration = 14*4)+
  ggplot2::scale_colour_viridis_d()+
  ggplot2::coord_cartesian(xlim=c(-0.25,0.25),ylim=c(0.02,20))
}

Description

Strictly integer breaks for continuous scale

Usage

integer_breaks(n = 5, ...)

Arguments

Value

a ggplot breaks function

Description

This solves the relationship between $R_t$ and growth rates as described by Wallinga and Lipsitch, to get a growth rate from $R_t$ and infectivity profile.

Usage

inv_wallinga_lipsitch(
  Rt,
  y = i_empirical_ip,
  a1 = seq(0.5, length.out = length(y)),
  a0 = dplyr::lag(a1, default = 0)
)

Arguments

Details

This function uses a single empirical distribution for the infectivity profile / generation time. If multiple are provided then the average central value is chosen (i.e. this does not propagate uncertainty in infectivity profile)

Value

an vector of growth rates

Examples

inv_wallinga_lipsitch(Rt=seq(0.5,2.5,length.out=9), y=example_ip())

Description

Check whether vector is a date

Usage

is.Date(x)

Arguments

Value

TRUE if dates, FALSE otherwise

Examples

is.Date(Sys.Date())

Description

Extract the weekday, month or quarter, or the Julian time (days since some origin). These are generic functions: the methods for the internal date-time classes are documented here.

Usage

## S3 method for class 'time_period'
julian(x, ...)

Arguments

Value

weekdays and months return a character vector of names in the locale in use, i.e., Sys.getlocale("LC_TIME").

quarters returns a character vector of "Q1" to "Q4".

julian returns the number of days (possibly fractional) since the origin, with the origin as a "origin" attribute. All time calculations in R are done ignoring leap-seconds.

Note

Other components such as the day of the month or the year are very easy to compute: just use as.POSIXlt and extract the relevant component. Alternatively (especially if the components are desired as character strings), use strftime.

See Also

DateTimeClasses, Date; Sys.getlocale("LC_TIME") crucially for months() and weekdays().

Examples

## first two are locale dependent:
weekdays(.leap.seconds)
months  (.leap.seconds)
quarters(.leap.seconds)

## Show how easily you get month, day, year, day (of {month, week, yr}), ... :
## (remember to count from 0 (!): mon = 0..11, wday = 0..6,  etc !!)

##' Transform (Time-)Date vector  to  convenient data frame :
dt2df <- function(dt, dName = deparse(substitute(dt))) {
    DF <- as.data.frame(unclass(as.POSIXlt( dt )))
    `names<-`(cbind(dt, DF, deparse.level=0L), c(dName, names(DF)))
}
## e.g.,
dt2df(.leap.seconds)    # date+time
dt2df(Sys.Date() + 0:9) # date

##' Even simpler:  Date -> Matrix - dropping time info {sec,min,hour, isdst}
d2mat <- function(x) simplify2array(unclass(as.POSIXlt(x))[4:7])
## e.g.,
d2mat(seq(as.Date("2000-02-02"), by=1, length.out=30)) # has R 1.0.0's release date


## Julian Day Number (JDN, https://en.wikipedia.org/wiki/Julian_day)
## is the number of days since noon UTC on the first day of 4317 BCE.
## in the proleptic Julian calendar.  To more recently, in
## 'Terrestrial Time' which differs from UTC by a few seconds
## See https://en.wikipedia.org/wiki/Terrestrial_Time
julian(Sys.Date(), -2440588) # from a day
floor(as.numeric(julian(Sys.time())) + 2440587.5) # from a date-time

Description

Create a set of labels for a time period based on the start and duration of the period. The format is configurable using the start and end dates and the dfmt and ifmt parameters, however if the time period has names then these are used in preference.

Usage

## S3 method for class 'time_period'
labels(
  object,
  ...,
  dfmt = "%d/%b",
  ifmt = "{start} — {end}",
  na.value = "Unknown"
)

Arguments

Value

a set of character labels for the time

Examples

eg = as.time_period(Sys.Date()+0:10*7, unit="1 week", anchor="start")

labels(eg)
labels(eg, ifmt="{start}", dfmt="%d/%b/%y")
labels(eg, ifmt="until {end}", dfmt="%d %b %Y")

# labels retained in constructor:
eg2 = Sys.Date()+0:10*7
names(eg2) = paste0("week ",0:10)
labels(eg2)
labels(as.time_period(eg2, anchor="start"))

Description

Coerce an object to a ggoutbreak compatible case linelist.

Usage

linelist(x, ...)

Arguments

Value

minimally a time stamped linelist dataframe A dataframe containing the following columns:

• time (ggoutbreak::time_period) - A set of events with a timestamp as a time_period

Any grouping allowed.

Examples

(Sys.Date()+stats::runif(100)*7) %>% linelist()

Description

Perform logit scaling with correct axis formatting. To not be used directly but with ggplot (e.g. ggplot2::scale_y_continuous(trans = "logit") )

Usage

logit_trans(n = 5, ...)

Arguments

Value

A scales object

Examples

dplyr::tibble(pvalue = c(0.001, 0.05, 0.1), fold_change = 1:3) %>%
 ggplot2::ggplot(ggplot2::aes(fold_change , pvalue)) +
 ggplot2::geom_point() +
 ggplot2::scale_y_continuous(trans = "logit")

Description

Recover a long format infectivity profile from an EpiEstim style matrix

Usage

make_empirical_ip(omega, normalise = TRUE)

Arguments

Value

a long format ip delay distribution

Examples

format_ip(make_empirical_ip(c(0,0,1,1,1,2,2,2,1,1)))

Description

Generate a simple discrete infectivity profile from a gamma distribution

Usage

make_fixed_ip(mean, sd = sqrt(mean), epiestim_compat = FALSE)

Arguments

Value

a discrete infectivity profile with 1 bootstrap

Examples

tmp = make_fixed_ip(5,2)

if(interactive()) {
  plot_ip(tmp, alpha=1)
}

Description

The infectivity profile is typically fitted to data by MCMC and reported as median and 95% credible intervals, of the mean, and the SD of (usually) a gamma distribution. This function generates a discrete infectivity probability distribution representing the chance that an infectee was infected on any specific day after the infector was infected (given that the infectee was infected).

Usage

make_gamma_ip(
  median_of_mean,
  lower_ci_of_mean = median_of_mean,
  upper_ci_of_mean = median_of_mean,
  median_of_sd = sqrt(median_of_mean),
  lower_ci_of_sd = median_of_sd,
  upper_ci_of_sd = median_of_sd,
  correlation = NA,
  n_boots = 100,
  epiestim_compat = FALSE,
  epiestim_sampler = epiestim_compat,
  z_crit = 0.95,
  seed = Sys.time()
)

Arguments

Details

EpiEstim generates these distributions by sampling from a truncated normal distribution for both mean and sd. The means and sds thus produced are discretised using a gamma distribution offset by 1 day, to enforce that the probability of infection on day zero is zero.

This constraint changes the shape of the distribution somewhat and may cause a small bias (although there is no ground truth to evaluate). In this function two different sampling and discretisation strategy are provided. The sampler uses log-normal distributions for both mean and SD with a degree of correlation. The discretizer assigns probabilities direct from the CDF of the gamma distribution without offset. This results in non zero values for the probability at time zero and can only be used with Rt estimation methods that can handle zero/negative serial intervals (e.g. rt_from_incidence or rt_from_renewal, or rt_from_growth_rate). The alternative follows EpiEstims algorithm.

Value

a long format infectivity profile data frame, or a list of dataframes if input is a vector.

Examples

# COVID-19 estimates from Ganyani et al 2020.
tmp = make_gamma_ip(5.2, 3.78, 6.78, 1.72, 0.91, 3.93,
  epiestim_sampler=FALSE, epiestim_compat=FALSE)

tmp %>%
  dplyr::group_by(boot) %>%
  dplyr::summarise(
    mean = sum(tau*probability),
    sd = sqrt(sum((tau-sum(tau*probability))^2*probability))
  ) %>%
  dplyr::summarise(
    mean = sprintf("%1.2f [%1.2f-%1.2f]",
      stats::quantile(mean,0.5),
      stats::quantile(mean,0.025),
      stats::quantile(mean,0.975)),
    sd = sprintf("%1.2f [%1.2f-%1.2f]",
      stats::quantile(sd,0.5),
      stats::quantile(sd,0.025),
      stats::quantile(sd,0.975))
  )

if(interactive()) {
  plot_ip(tmp, alpha=0.1) +
    ggplot2::coord_cartesian(xlim=c(0,15))
}

means = c(3,4,5)
ips = make_gamma_ip(means)

if (interactive()) {
  purrr::map2(ips,means, ~ .x %>% dplyr::mutate(label = sprintf("Mean: %1.2f",.y))) %>%
    purrr::map( ~ plot_ip(.x,alpha=0.1)+ggplot2::facet_wrap(~label))
}

Description

The infectivity profile is typically fitted to data by MCMC as a gamma distribution. This function generates a discrete infectivity probability distribution representing the chance that an infectee was infected on any specific day after the infector was infected (given that the infectee was infected), from posterior samples.

Usage

make_posterior_ip(
  ...,
  mean,
  sd,
  shape,
  rate,
  scale,
  epiestim_compat = FALSE,
  n_boots = 100
)

Arguments

Details

If using EpiEstim and coarseDataTools::dic.fit.mcmc the output of the MCMC will be a S4 object with a samples slot, containing a dataframe of shape=var1 and scale=var2 columns. to use this output with make_posterior_ip invoke it like this:

do.call(make_posterior_ip, SI_fit_clever@samples %>% dplyr::rename(shape=var1, scale=var2))

N.b. only one combination of mean and sd, shape and rate, or shape and scale, are required.

Value

a long format ip delay distribution

Examples

tmp = make_posterior_ip(
  mean = stats::rnorm(100,5,0.1),
  sd = stats::rnorm(100,1.5,0.1)
)
tmp %>% dplyr::glimpse()
if (interactive()) plot_ip(tmp)

Description

Suits larger contact tracing data sets where there is a delay between 2 events which may or may not be precisely known.

Usage

make_resampled_ip(
  tau,
  min_tau = pmax(tau - 0.5, truncate),
  max_tau = tau + 0.5,
  add_noise = TRUE,
  truncate = 0,
  n_boots = 100,
  seed = Sys.time()
)

Arguments

Value

a long format ip delay distribution

Examples

tau = rgamma2(100, 5,2)
ip = make_resampled_ip(min_tau = tau-1, max_tau = tau+1, seed = 100)
if(interactive()) {
  plot_ip(ip,alpha=0.1)
}

Description

max.Date returns an integer and -Inf for a set of NA dates. This is usually inconvenient.

Usage

max_date(x, ...)

Arguments

Value

a date. '0001-01-01“ if there is no well defined minimum.

Examples

max_date(NA)

Description

min.Date returns an integer and Inf for a set of NA dates. This is usually inconvenient.

Usage

min_date(x, ...)

Arguments

Value

a date. 9999-12-31 if there is no well defined minimum.

Examples

min_date(NA)

Description

Extract the weekday, month or quarter, or the Julian time (days since some origin). These are generic functions: the methods for the internal date-time classes are documented here.

Usage

## S3 method for class 'time_period'
months(x, abbreviate = FALSE)

Arguments

Value

weekdays and months return a character vector of names in the locale in use, i.e., Sys.getlocale("LC_TIME").

quarters returns a character vector of "Q1" to "Q4".

julian returns the number of days (possibly fractional) since the origin, with the origin as a "origin" attribute. All time calculations in R are done ignoring leap-seconds.

Note

Other components such as the day of the month or the year are very easy to compute: just use as.POSIXlt and extract the relevant component. Alternatively (especially if the components are desired as character strings), use strftime.

See Also

DateTimeClasses, Date; Sys.getlocale("LC_TIME") crucially for months() and weekdays().

Examples

## first two are locale dependent:
weekdays(.leap.seconds)
months  (.leap.seconds)
quarters(.leap.seconds)

## Show how easily you get month, day, year, day (of {month, week, yr}), ... :
## (remember to count from 0 (!): mon = 0..11, wday = 0..6,  etc !!)

##' Transform (Time-)Date vector  to  convenient data frame :
dt2df <- function(dt, dName = deparse(substitute(dt))) {
    DF <- as.data.frame(unclass(as.POSIXlt( dt )))
    `names<-`(cbind(dt, DF, deparse.level=0L), c(dName, names(DF)))
}
## e.g.,
dt2df(.leap.seconds)    # date+time
dt2df(Sys.Date() + 0:9) # date

##' Even simpler:  Date -> Matrix - dropping time info {sec,min,hour, isdst}
d2mat <- function(x) simplify2array(unclass(as.POSIXlt(x))[4:7])
## e.g.,
d2mat(seq(as.Date("2000-02-02"), by=1, length.out=30)) # has R 1.0.0's release date


## Julian Day Number (JDN, https://en.wikipedia.org/wiki/Julian_day)
## is the number of days since noon UTC on the first day of 4317 BCE.
## in the proleptic Julian calendar.  To more recently, in
## 'Terrestrial Time' which differs from UTC by a few seconds
## See https://en.wikipedia.org/wiki/Terrestrial_Time
julian(Sys.Date(), -2440588) # from a day
floor(as.numeric(julian(Sys.time())) + 2440587.5) # from a date-time

Description

Takes a list of times, classes and counts, e.g. a COGUK variant like data set with time, (multinomial) class (e.g. variant) and count being the count in that time period. Fits a quadratic B-spline on time to the proportion of the data using nnet::multinom, with approx one degree of freedom per class and per window units of the time series.

Usage

multinomial_nnet_model(
  d = i_multinomial_input,
  ...,
  window = 14,
  frequency = "1 day",
  predict = TRUE,
  .progress = interactive()
)

Arguments

Details

Additional groupings are treated as distinct proportions models.

Value

a new dataframe with time (as a time period), class, and proportion.0.5, or a model object

Examples

data = example_poisson_rt_2class()
tmp = data %>% multinomial_nnet_model(window=14)

if (interactive()) {
  plot_multinomial(tmp, date_labels="%b %y")
}

Description

This assumes positive disease counts are stratified by a population grouping, e.g. geography or age, and we have estimates of the size of that population during that time period. Normalising by population size allows us to compare groups.

Usage

normalise_count(
  raw = i_incidence_data,
  pop = i_population_data,
  ...,
  population_unit = 1e+05,
  normalise_time = FALSE
)

Arguments

Value

a dataframe with incidence rates per unit capita. A dataframe containing the following columns:

• population (positive_integer) - Size of population

• count (positive_integer) - Positive case counts associated with the specified time frame

• time (ggoutbreak::time_period + group_unique) - A (usually complete) set of singular observations per unit time as a time_period

Any grouping allowed.

Examples

data = example_england_covid_by_age()
demog = ukc19::uk_population_2019_by_5yr_age

data %>%
  normalise_count(demog) %>%
  dplyr::glimpse()

Description

This assumes positive disease counts are stratified by a population grouping, e.g. geography or age, and we have estimates of the size of that population during that time period. Normalising by population size allows us to compare groups.

Usage

normalise_incidence(
  modelled = i_incidence_model,
  pop = i_population_data,
  ...,
  population_unit = 1e+05,
  normalise_time = FALSE
)

Arguments

Value

a dataframe with incidence rates per unit capita. A dataframe containing the following columns:

• time (ggoutbreak::time_period + group_unique) - A (usually complete) set of singular observations per unit time as a time_period

• incidence.per_capita.fit (double) - an estimate of the incidence per capita rate on a log scale

• incidence.per_capita.se.fit (positive_double) - the standard error of the incidence per capita rate estimate on a log scale

• incidence.per_capita.0.025 (positive_double) - lower confidence limit of the incidence per capita rate (true scale)

• incidence.per_capita.0.5 (positive_double) - median estimate of the incidence per capita rate (true scale)

• incidence.per_capita.0.975 (positive_double) - upper confidence limit of the incidence per capita rate (true scale)

• population_unit (double) - The population unit on which the per capita incidence rate is calculated

• time_unit (lubridate::as.period) - The time period over which the per capita incidence rate is calculated

Any grouping allowed.

Examples

model = example_poisson_age_stratified()
demog = ukc19::uk_population_2019_by_5yr_age

model %>%
  normalise_incidence(demog) %>%
  dplyr::glimpse()

Description

This converts a long format infectivity profile to something that will work with EpiEstim. In general ggoutbreak will do the conversion internally. This is provided as a utility for examples mostly.

Usage

omega_matrix(ip = i_discrete_ip, epiestim_compat = TRUE)

Arguments

Value

a matrix with 0:max(tau) rows and 1:max(boot) columns.

Examples

omega_matrix(example_ip())

Description

Density, distribution function, quantile function and random generation for the Beta distribution with parameters shape1 and shape2 (and optional non-centrality parameter ncp).

Usage

pbeta2(q, prob, kappa, lower.tail = TRUE, log.p = FALSE)

Arguments

Value

dbeta gives the density, pbeta the distribution function, qbeta the quantile function, and rbeta generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rbeta, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::pbeta()

Examples

pbeta2(c(0.25,0.5,0.75), 0.5, 0.25)

Description

The following conversion describes the parameters mean and kappa

Usage

pcgamma(q, mean, kappa = 1/mean, lower.tail = TRUE, log.p = FALSE)

Arguments

Details

$\text{shape:} \alpha = \frac{1}{\kappa} \\ \text{rate:} \beta = \frac{1}{\mu \times \kappa} \\ \text{scale:} \sigma = \mu \times \kappa \\$

Value

dgamma gives the density, pgamma gives the distribution function, qgamma gives the quantile function, and rgamma generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgamma, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::pgamma()

Examples

pcgamma(seq(0,4,0.25), 2, 0.5)

Description

Density, distribution function, quantile function and random generation for the Gamma distribution with parameters shape and scale.

Usage

pgamma2(q, mean, sd = sqrt(mean), lower.tail = TRUE, log.p = FALSE)

Arguments

Value

dgamma gives the density, pgamma gives the distribution function, qgamma gives the quantile function, and rgamma generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgamma, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::pgamma()

Examples

pgamma2(seq(0,4,0.25), 2, 1)

Description

Density, distribution function, quantile function and random generation for the log normal distribution whose logarithm has mean equal to meanlog and standard deviation equal to sdlog.

Usage

plnorm2(q, mean = 1, sd = sqrt(exp(1) - 1), lower.tail = TRUE, log.p = FALSE)

Arguments

Details

The log normal distribution has density

$f(x) = \frac{1}{\sqrt{2\pi}\sigma x} e^{-(\log(x) - \mu)^2/2 \sigma^2}$

where $\mu$ and $\sigma$ are the mean and standard deviation of the logarithm. The mean is $E(X) = exp(\mu + 1/2 \sigma^2)$, the median is $med(X) = exp(\mu)$, and the variance $Var(X) = exp(2\mu + \sigma^2)(exp(\sigma^2) - 1)$ and hence the coefficient of variation is $\sqrt{exp(\sigma^2) - 1}$ which is approximately $\sigma$ when that is small (e.g., $\sigma < 1/2$).

Value

dlnorm gives the density, plnorm gives the distribution function, qlnorm gives the quantile function, and rlnorm generates random deviates.

The length of the result is determined by n for rlnorm, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

Note

The cumulative hazard $H(t) = - \log(1 - F(t))$ is -plnorm(t, r, lower = FALSE, log = TRUE).

Source

dlnorm is calculated from the definition (in ‘Details’). [pqr]lnorm are based on the relationship to the normal.

Consequently, they model a single point mass at exp(meanlog) for the boundary case sdlog = 0.

References

Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole.

Johnson, N. L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions, volume 1, chapter 14. Wiley, New York.

See Also

Distributions for other standard distributions, including dnorm for the normal distribution.

Examples

plnorm2(seq(0,4,0.25), 2, 1)

Description

The logit-normal distribution has a support of 0 to 1.

Usage

plogitnorm(q, meanlogit = 0, sdlogit = 1, lower.tail = TRUE, log.p = FALSE)

Arguments

Value

a vector of probabilities, quantiles, densities or samples.

Examples

plogitnorm(seq(0.1,0.9,0.1), 0, 1)

Description

The logit-normal distribution has a support of 0 to 1.

Usage

plogitnorm2(
  q,
  prob.0.5 = 0.5,
  kappa = 1 - exp(-1),
  lower.tail = TRUE,
  log.p = FALSE
)

Arguments

Value

a vector of probabilities, quantiles, densities or samples.

Examples

plogitnorm2(seq(0.1,0.9,0.1), 0.75, 0.2)
plogitnorm2(qlogitnorm2(seq(0.1,0.9,0.1), 0.75, 0.2), 0.75, 0.2)

Description

Plot a line-list of cases as a histogram

Usage

plot_cases(
  raw,
  ...,
  mapping = .check_for_aes(raw, ..., class_aes = "fill"),
  events = i_events
)

Arguments

Value

a ggplot object

Examples

with_defaults("2025-01-01" ,"1 day", {
  tmp = dplyr::tibble(
    time = as.time_period( rexpgrowth(100,0.05,40)),
    class = rep(c("one","two","three"), length.out=100)
  )
})

if(interactive()) {
  plot_cases(tmp, mapping=ggplot2::aes(fill = class),linewidth=0.1,colour="white")
}

Description

Plot a raw case count timeseries

Usage

plot_counts(raw, ..., mapping = .check_for_aes(raw, ...), events = i_events)

Arguments

Value

a ggplot object

Examples

# example code

tmp = example_england_covid_by_age() %>%
  time_aggregate(count=sum(count)) %>%
  normalise_count(pop=56489700, population_unit=1000, normalise_time=TRUE)

# normalised by England population (56489700 people)

if(interactive()) {
  plot_counts(tmp, colour="blue",size=0.25)
}

Description

Plot an incidence or proportion versus growth phase diagram

Usage

plot_growth_phase(
  modelled,
  timepoints = NULL,
  duration = max(dplyr::count(modelled)$n),
  interval = 7,
  mapping = if (interfacer::is_col_present(modelled, class)) {
     ggplot2::aes(colour =
    class)
 } else {
     ggplot2::aes()
 },
  cis = TRUE,
  ...
)

Arguments

Value

a ggplot

Examples

data = example_poisson_rt_2class()
tmp2 = data %>% poisson_locfit_model()

timepoints = as.Date(tmp2$time[c(40,80,120,160)])

if(interactive()) {
  plot_growth_phase(tmp2, timepoints, duration=108)
}