Version:	0.15.4
Date:	2025-02-14
License:	GPL-3
Encoding:	UTF-8
URL:	https://github.com/davidearn/epigrowthfit
BugReports:	https://github.com/davidearn/epigrowthfit/issues
Title:	Nonlinear Mixed Effects Models of Epidemic Growth
Description:	Maximum likelihood estimation of nonlinear mixed effects models of epidemic growth using Template Model Builder ('TMB'). Enables joint estimation for collections of disease incidence time series, including time series that describe multiple epidemic waves. Supports a set of widely used phenomenological models: exponential, logistic, Richards (generalized logistic), subexponential, and Gompertz. Provides methods for interrogating model objects and several auxiliary functions, including one for computing basic reproduction numbers from fitted values of the initial exponential growth rate. Preliminary versions of this software were applied in Ma et al. (2014) <doi:10.1007/s11538-013-9918-2> and in Earn et al. (2020) <doi:10.1073/pnas.2004904117>.
Depends:	R (≥ 4.3)
LinkingTo:	RcppEigen (≥ 0.3.4.0.0), TMB
Imports:	Matrix (≥ 1.6-2), TMB, grDevices, graphics, methods, nlme, stats, utils
LazyData:	true
NeedsCompilation:	yes
Packaged:	2025-02-18 13:50:46 UTC; mikael
Author:	Mikael Jagan [aut, cre], Ben Bolker [aut], Jonathan Dushoff [ctb], David Earn [ctb], Junling Ma [ctb]
Maintainer:	Mikael Jagan <jaganmn@mcmaster.ca>
Repository:	CRAN
Date/Publication:	2025-02-19 14:50:02 UTC

R Package epigrowthfit

Description

An R package for estimating nonlinear mixed effects models of epidemic growth.

Details

The “main” model estimating function is egf.

To render a list of available help topics, use help(package = "epigrowthfit"). Many of these document methods for the class of objects returned by egf.

To report a bug or request a change, use bug.report(package = "epigrowthfit").

Author(s)

Mikael Jagan jaganmn@mcmaster.ca

Compute the Basic Reproduction Number

Description

Computes the basic reproduction number as a function of the initial exponential growth rate, conditional on a binned generation interval distribution.

Usage

R0(r, breaks, probs)

Arguments

Value

A numeric vector listing basic reproduction numbers.

Computation

For an initial exponential growth rate r, the basic reproduction number is computed as

r / sum(probs * (exp(-r * breaks[-n]) - exp(-r * breaks[-1L])) / (breaks[-1L] - breaks[-n])) ,

where n = length(breaks).

References

Wallinga, J. & Lipsitch M. (2007). How generation intervals shape the relationship between growth rates and reproductive numbers. Proceedings of the Royal Society B: Biological Sciences, 274(1609), 599-604. doi:10.1098/rspb.2006.3754

See Also

timescale, finalsize.

Examples

r <- seq(0, 1, 0.02)
breaks <- seq(0, 20, 1)
probs <- diff(pgamma(breaks, shape = 1, scale = 2.5))

plot(r, R0(r, breaks, probs), las = 1,
     xlab = "initial exponential growth rate",
     ylab = "basic reproduction number")

Extract Coefficients and Random Effect Covariance Parameters

Description

Extracts the bottom level parameter vector c(beta, theta, b) or a subset. Segments beta, theta, and b contain (respectively) fixed effect coefficients, random effect covariance parameters, and random effect coefficients.

Usage

## S3 method for class 'egf'
coef(object, random = FALSE, full = FALSE, ...)

Arguments

Value

A numeric vector concatenating beta, theta, and b, without b if random = FALSE and without mapped elements if full = FALSE.

Attribute len is a named integer vector partitioning the result by segment.

Attribute map is a named list of integer vectors i such that that a full segment y and its condensed counterpart x are related by y = x[i], with the exception that i[j] is NA if y[j] is mapped to an initial value. NULL is used in place of an integer vector where x and y are identical.

The result inherits from class coef.egf, which has methods for print, as.list, and labels.

See Also

The generic function coef.

Examples

example("egf", package = "epigrowthfit")

for (random in c(FALSE, TRUE)) {
    for (full in c(FALSE, TRUE)) {
        cat(sprintf("random = %s, full = %s :\n\n", random, full))
        str(coef(m1, random = random, full = full))
        cat("\n")
    }
}

Confidence Intervals

Description

Computes confidence intervals on fixed effect coefficients, random effect covariance parameters, and linear combinations thereof, including population fitted values. Intervals on individual fitted values accounting for random effects are supported, but only by method = "wald".

Usage

## S3 method for class 'egf'
confint(object, parm, level = 0.95,
        A = seq_along(par), method = c("wald", "profile", "uniroot"), scale = 7,
        parallel = egf_parallel(), trace = FALSE,
        top = egf_top(object), subset = NULL, select = NULL,
        class = FALSE, link = TRUE, random = FALSE, ...)

## S3 method for class 'confint.egf'
plot(x, by = 12L,
     subset = NULL, order = NULL, label = NULL, main = NULL, ...)

Arguments

Details

Three methods for computing confidence intervals are available:

"wald"

confidence limits are calculated as

value + c(-1, 1) * sqrt(q) * se

where q = qchisq(level, df = 1).

"profile", "uniroot"

confidence limits are calculated as approximate solutions of the equation

2 * (f(x) - f(value)) = q

where q = qchisq(level, df = 1) and f is the negative log marginal likelihood function expressed as a function of the parameter x in question. Solutions are approximated by interpolating a likelihood profile ("profile") or by rootfinding ("uniroot").

"wald" assumes asymptotic normality of the maximum likelihood estimator. "profile" and "uniroot" avoid this contraint but are typically expensive, requiring estimation of many restricted models. They are parallelized at the C++ level when there is OpenMP support and object[["control"]][["omp_num_threads"]] is set to an integer greater than 1. When there is no OpenMP support, they can still be parallelized at the R level with appropriate setting of argument parallel.

Value

A numeric array in 2 or 3 dimensions containing the lower and upper confidence limits in the last dimension.

When confidence intervals on fitted values are desired, the user will set A = NULL and in that case have the option of passing class = TRUE to obtain an augmented result. Thus, alternatively:

A data frame inheriting from class confint.egf, with variables:

The confidence level level is preserved as an attribute.

See Also

The generic function confint.

Examples

example("egf", package = "epigrowthfit")

zz1 <- confint(m1, A = NULL, method =    "wald", class = TRUE,
               random = TRUE)
str(zz1)

op <- par(mar = c(4.5, 4, 2, 1), oma = c(0, 0, 0, 0))
plot(zz1)
par(op)

zz2 <- confint(m1, A = NULL, method = "profile", class = TRUE,
               top = "log(r)", subset = quote(country == "A" & wave == 1))
zz3 <- confint(m1, A = NULL, method = "uniroot", class = TRUE,
               top = "log(r)", subset = quote(country == "A" & wave == 1))

Compute a Packed Representation of a Covariance Matrix

Description

Transform covariances matrices to a “packed” representation or compute the inverse transformation.

Usage

cov2theta(Sigma)
theta2cov(theta)

Arguments

Details

An n-by-n real, symmetric, positive definite matrix \Sigma can be factorized as

\Sigma = R' R\,.

The upper triangular Cholesky factor R can be written as

R = R_{1} D^{-1/2} D_{\sigma}^{1/2}\,,

where R_{1} is a unit upper triangular matrix and D = \mathrm{diag}(\mathrm{diag}(R_{1}' R_{1})) and D_{\sigma} = \mathrm{diag}(\mathrm{diag}(\Sigma)) are diagonal matrices.

cov2theta takes \Sigma and returns the vector \theta of length n(n+1)/2 containing the log diagonal entries of D_{\sigma} followed by (in column-major order) the strictly upper triangular entries of R_{1}. theta2cov computes the inverse transformation.

Value

A vector like theta (cov2theta) or a matrix like Sigma (theta2cov); see ‘Details’.

COVID-19 in Ontario, Canada

Description

Time series of COVID-19 cases and tests in Ontario, Canada, daily from February 8, 2020 to May 1, 2022.

Usage

data(covid19.ontario, package = "epigrowthfit")

Format

A data frame with 814 rows and 3 variables:

date

a Date vector.

cases

an integer vector. cases[i] is the number of cases confirmed by Ontario public health units prior to date[i]. This number includes resolved and fatal cases as well as reinfections.

tests

an integer vector. tests[i] is the number of tests completed prior to date[i]. This number includes repeated tests by individuals except prior to April 15, 2020, when individuals were counted at most once.

Source

This data set is a processed subset of a larger data set downloaded on 2024-01-10 from the It is updated using an installed script:

\link{system.file}("scripts", "covid19.ontario.R", package = "epigrowthfit")

Examples

data(covid19.ontario, package = "epigrowthfit")
plot(1 + diff(c(NA, cases)) ~ date, data = covid19.ontario, log = "y")

Extract the Residual Degrees of Freedom

Description

Extracts from a model object the number of observations (see nobs) minus the number of estimated parameters (fixed effect coefficients and random effect covariance parameters).

Usage

## S3 method for class 'egf'
df.residual(object, ...)

Arguments

Value

An integer.

See Also

The generic function df.residual.

Fit Nonlinear Mixed Effects Models of Epidemic Growth

Description

Fits nonlinear mixed effects models of epidemic growth to collections of one or more disease incidence time series.

Usage

egf(model, ...)

## S3 method for class 'egf_model'
egf(model,
    formula_ts,
    formula_windows,
    formula_parameters = list(),
    formula_priors = list(),
    data_ts,
    data_windows,
    subset_ts = NULL,
    subset_windows = NULL,
    select_windows = NULL,
    na_action_ts = c("fail", "pass"),
    na_action_windows = c("fail", "omit"),
    control = egf_control(),
    init = list(),
    map = list(),
    fit = TRUE,
    se = FALSE,
    ...)

Arguments

Details

Users attempting to set arguments formula_priors, init, and map should know the structure of the bottom level parameter vector. It is described under topic egf-class.

If

formula_ts = cbind(time, x) ~ ts1
formula_windows = cbind(start, end) ~ ts2

then it is expected that time, start, and end (after coercion to numeric) measure time on the same scale. To be precise, numeric times should have a common unit of measure and, at least within time series, represent displacements from a common reference time. These conditions will always hold if time, start, and end all evaluate to Date, POSIXct, or POSIXlt vectors.

When day of week effects are estimated, numeric times are interpreted as numbers of days since midnight on January 1, 1970, so that time points can be mapped unambiguously to days of week. Furthermore, in this case, time (after coercion to numeric) is required to be integer-valued with one day spacing in all time series. This means that

isTRUE(all.equal(time, round(time))) &&
    all(range(diff(round(time))) == 1)

must be TRUE in each level of ts1. These conditions ensure that intervals between successive time points represent exactly one day of week.

Value

A list inheriting from class egf. See topic egf-class for class documentation.

See Also

The many methods for class egf, listed by methods(class = "egf").

Examples

## Simulate 'N' incidence time series exhibiting exponential growth
set.seed(180149L)
N <- 10L
f <- function(time, r, c0) {
    lambda <- diff(exp(log(c0) + r * time))
    c(NA, rpois(lambda, lambda))
}
time <- seq.int(0, 40, 1)
r <- rlnorm(N, -3.2, 0.2)
c0 <- rlnorm(N, 6, 0.2)
data_ts <-
    data.frame(country = gl(N, length(time), labels = LETTERS[1:N]),
               time = rep.int(time, N),
               x = unlist(Map(f, time = list(time), r = r, c0 = c0)))
rm(f, time)

## Define fitting windows (here, two per time series)
data_windows <-
    data.frame(country = gl(N, 1L, 2L * N, labels = LETTERS[1:N]),
               wave = gl(2L, 10L),
               start = c(sample(seq.int(0, 5, 1), N, TRUE),
                         sample(seq.int(20, 25, 1), N, TRUE)),
               end = c(sample(seq.int(15, 20, 1), N, TRUE),
                       sample(seq.int(35, 40, 1), N, TRUE)))

## Estimate the generative model
m1 <-
    egf(model = egf_model(curve = "exponential", family = "pois"),
        formula_ts = cbind(time, x) ~ country,
        formula_windows = cbind(start, end) ~ country,
        formula_parameters = ~(1 | country:wave),
        data_ts = data_ts,
        data_windows = data_windows,
        se = TRUE)

## Re-estimate the generative model with:
## * Gaussian prior on beta[1L]
## * LKJ prior on all random effect covariance matrices
##   (here there happens to be just one)
## * initial value of 'theta' set explicitly
## * theta[3L] fixed at initial value
m2 <-
    update(m1,
           formula_priors = list(beta[1L] ~ Normal(mu = -3, sigma = 1),
                                 Sigma ~ LKJ(eta = 2)),
           init = list(theta = c(log(0.5), log(0.5), 0)),
           map = list(theta = 3L))

Description of Objects of Class egf

Description

Class egf designates models estimated by function egf. Objects of this class hold information about an estimated model. Components can be accessed directly. However, as the components are subject to change without notice, portable code will rely on exported methods for interrogation.

Details

Currently, a legitimate egf object is a list with elements:

model

a copy of the so-named argument of egf.

frame

a list of the form list(ts, windows, parameters, extra). ts and windows are data frames preserving time series and fitting window endpoints. parameters is a list of mixed effects model frames, with one element for each top level nonlinear model parameter. extra is a data frame preserving additional variables specified in call[["select_windows"]]. windows, the model frames listed in parameters, and extra all correspond rowwise.

priors

a list of the form list(top, bottom = list(beta, theta, Sigma)), where top, beta, theta, and Sigma are all lists of egf_prior objects.

control

a copy of the so-named argument of egf.

tmb_out

the list output of MakeADFun.

optimizer_out

the list output of the optimizer specified by control[["optimizer"]].

init, best

numeric vectors giving the values of the condensed bottom level parameter vector used in the first and maximal likelihood evaluations.

random

a logical vector indexing the elements of the condensed bottom level parameter vector that are not arguments of the negative log marginal likelihood function. It indexes all elements of segment b (random effect coefficients) and (but only if control[["profile"]] = TRUE) all elements of segment beta (fixed effect coefficients).

value, gradient

numeric vectors giving the value and gradient of the negative log marginal likelihood function at best[!random].

hessian

a logical flag indicating whether the Hessian matrix of the negative log marginal likelihood function is positive definite at best[!random]. NA means that the matrix has not been computed.

coefficients

a list of the form list(fixed, random), where fixed and random are data frames preserving interpretive information about fixed and random effect coefficients.

contrasts

a list of the form list(fixed, random), where fixed and random are lists preserving contrasts used to construct the fixed and random effects design matrices.

call

the call to egf, enabling updates to the object by the default method of generic function update.

Bottom Level Parameter Vector

An estimated model is specified by a bottom level parameter vector that is the concatenation of three segments:

beta

the result of unlist(lbeta), where lbeta is a list of numeric vectors of fixed effect coefficients, with one vector for each top level nonlinear model parameter. The order of top level parameters is specified by egf_top(model).

theta

the result of unlist(ltheta), where ltheta is a list of numeric vectors of random effect covariance parameters, with one vector for each distinct random effect term in formula_parameters. Each vector parametrizes a random effect covariance matrix via theta2cov and its inverse cov2theta.

The list Sigma mentioned in the description of egf argument formula_priors is precisely lapply(ltheta, theta2cov).

b

the result of unlist(lb), where lb is a list of numeric matrices of scaled random effect coefficients, corresponding elementwise to ltheta. The columns of lb[[i]] (one per level of the grouping variable) are interpreted as samples from a zero mean, unit variance multivariate normal distribution with covariance matrix cov2cor(theta2cov(ltheta[[i]])).

When elements of this vector are “mapped” via egf argument map, likelihood is defined as a function of the condensed vector that excludes mapped elements.

Methods are defined for generic functions coef, fixef, and ranef to allow users to interrogate the structure of the vector.

Examples

methods(class = "egf")
help.search("\\.egf$", fields = "alias", package = "epigrowthfit")
## less verbosely: alias??`\\.egf$`

Define Control Parameters

Description

Set parameters controlling the behaviour of egf.

Usage

egf_control(outer_optimizer = egf_optimizer(nlminb),
            inner_optimizer = egf_optimizer(newton),
            trace = FALSE, profile = FALSE, sparse_X = FALSE,
            omp_num_threads = getOption("egf.cores", 1L))

Arguments

Details

trace affects the amount of information printed during likelihood evaluations:

0

likelihood evaluations are always silent.

1

a message is printed whenever a negative log marginal likelihood term is NaN or exceeds 1e+09.

2

all negative log marginal likelihood terms are printed.

egf passes silent = trace == 0L to MakeADFun. A corollary is that nonzero values of trace have a number of additional side effects:

• error messages are printed during function and gradient evaluations;

• the maximum absolute gradient element is printed with each gradient evaluation; and

• trace flags set by config are turned on.

Value

A list inheriting from class egf_control containing the validated arguments.

Warning

Setting trace > 0L and omp_num_threads > 0L simultaneously should be avoided, because tracing messages are printed using R API functions that are not thread-safe.

Examples

control <- egf_control()
str(control)

Define Control Parameters for Plotting

Description

Sets parameters controlling the graphical output of plot for objects of class egf. Supplied values override package defaults (retrievable as defaults <- egf_control_plot()), which in turn override global defaults set via par.

Below, x, type, time_as, and delta refer to the so-named arguments of plot.egf.

Usage

egf_control_plot(window, data, predict, asymptote, box, axis, title, doubling)

Arguments

Details

Setting an argument (or an element thereof in the case of nested lists) to NULL has the effect of suppressing the corresponding plot element.

Value

A named list containing the package defaults modified according to the arguments in the call.

Test for Convergence

Description

Performs simple diagnostic tests to assess whether the optimizer that produced an estimated model actually converged to a local minimum point of the negative log marginal likelihood function.

Usage

egf_has_converged(object, check = TRUE, tol = 1)

Arguments

Value

TRUE if all tests pass. FALSE if any test fails. NA if no test fails, but the test for a positive definite Hessian matrix is indeterminate because the matrix has not been computed.

Test for Random Effects

Description

Tests whether an object specifies a model with random effects.

Usage

egf_has_random(object, check = TRUE)

Arguments

Value

TRUE or FALSE.

Define a Top Level Nonlinear Model

Description

Sets flags defining a top level nonlinear model of epidemic growth to be estimated by egf.

Usage

egf_model(curve = c("logistic", "richards", "exponential",
                    "subexponential", "gompertz"),
          excess = FALSE,
          family = c("nbinom", "pois"),
          day_of_week = FALSE)

Arguments

Value

A list inheriting from class egf_model containing the validated arguments.

See Also

simulate.egf_model.

Examples

model <- egf_model()
str(model)

Define an Optimization Method

Description

Utilities for linking an optimizer with optional arguments and control parameters to define an optimization method for use by egf.

Usage

egf_optimizer(f = nlminb, args = list(), control = list())

Arguments

Details

An optim-like function is a function f such that:

• the first three arguments of f specify an initial parameter vector, an objective function, and a gradient function, respectively;

• f accepts control as a fourth (or later) argument; and

• f returns a named list with elements par, value, convergence, and message.

Value

A list inheriting from class egf_optimizer containing the validated arguments, wherein f may be a new function wrapping the supplied one to make it optim-like.

Examples

optimizer <- egf_optimizer(nlminb)
str(optimizer)

Define a Parallelization Method

Description

Defines instructions for parallelization by linking a method with options.

Usage

egf_parallel(method = c("serial", "multicore", "snow"),
             outfile = "", cores = getOption("egf.cores", 1L),
             args = list(), cl = NULL)

Arguments

Value

A list inheriting from class "egf_parallel" containing the arguments (after possible matching and coercion).

See Also

vignette("parallel", "parallel").

Examples

parallel <- egf_parallel()
str(parallel)

Prior Distributions

Description

Functions used by egf to specify prior distributions of bottom level mixed effects model parameters.

Usage

Normal(mu = 0, sigma = 1)
LKJ(eta = 1)
Wishart(df, scale, tol = 1e-06)
InverseWishart(df, scale, tol = 1e-06)

Arguments

Value

A list inheriting from class egf_prior, with elements:

Examples

Normal(mu = 0, sigma = 1)
Normal(mu = -5:5, sigma = c(0.1, 1))

LKJ(eta = 2)

u <- matrix(rnorm(9L), 3L, 3L)
utu <- crossprod(u)
uut <- tcrossprod(u)
Wishart(df = 6, scale = utu)
InverseWishart(df = 6, scale = list(utu, uut))

Top Level Nonlinear Model Parameter Names

Description

Retrieves the names used internally for top level nonlinear model parameters.

Usage

egf_top(object, ...)

## S3 method for class 'egf_model'
egf_top(object, link = TRUE, ...)
## S3 method for class 'egf'
egf_top(object, link = TRUE, ...)

Arguments

Value

A character vector listing names relevant to object.

Defunct Functions in Package epigrowthfit

Description

The functions and other objects listed here are no longer part of epigrowthfit as they are no longer needed.

Usage

## Nothing yet!

Details

These either are stubs reporting that they are defunct or have been removed completely (apart from being documented here).

See Also

Deprecated, base-deprecated, epigrowthfit-deprecated, epigrowthfit-notyet.

Deprecated Functions in Package epigrowthfit

Description

The functions and other objects listed here are provided only for compatibility with older versions of epigrowthfit and may become defunct as soon as the next release.

Usage

## Nothing yet!

See Also

Defunct, base-defunct, epigrowthfit-defunct, epigrowthfit-notyet.

Not Yet Implemented Functions in Package epigrowthfit

Description

The functions listed here are defined but not yet implemented. Use bug.report(package = "epigrowthfit") to request an implementation.

Usage

## Nothing yet!

See Also

NotYetImplemented, epigrowthfit-deprecated, epigrowthfit-defunct.

Extract the (Generalized) AIC

Description

Extracts from a model object the generalized Akaike Information Criterion (AIC).

Usage

## S3 method for class 'egf'
extractAIC(fit, scale, k = 2, ...)

Arguments

Value

An numeric vector of length 2 giving the equivalent degrees of freedom and criterion value.

See Also

The generic function extractAIC.

Compute the Expected Epidemic Final Size

Description

Computes the proportion of a population expected to be infected over the course of an epidemic, as a function of the basic reproduction number.

Usage

finalsize(R0, S0, I0)

Arguments

Details

At least one of S0 and I0 must be supplied. If S0 (I0) is supplied but not I0 (S0), then the latter is assigned the value of one minus the former.

R0, S0, and I0 are recycled to a common length (the maximum of their lengths).

Value

A numeric vector listing values in the interval [0,1] for the expected epidemic final size.

Computation

The basic reproduction number R0 defines the expected epidemic final size Z through an implicit equation,

Z = S0 * (1 - exp(-R0 * (Z + I0))) ,

which admits an explicit solution

Z = S0 + (1/R0) * W(-R0 * S0 * exp(-R0 * (S0 + I0))) .

Here, W denotes the Lambert W function. finalsize computes this solution, relying on function lambertW from package emdbook.

References

Ma, J. & Earn, D. J. D. (2006). Generality of the final size formula for an epidemic of a newly invading infectious disease. Bulletin of Mathetmatical Biology, 68(3), 679-702. doi:10.1007/s11538-005-9047-7

See Also

timescale, R0.

Examples

R0 <- 10^seq(-3, 1, length.out = 151L)
plot(R0, finalsize(R0, S0 = 1, I0 = 0), type = "l", las = 1,
     xlab = "basic reproduction number",
     ylab = "final size")

Fitted Values

Description

Retrieves fitted values of top level nonlinear model parameters. The fitted value of a given parameter for a given fitting window is obtained by adding (i) the population fitted value computed as a linear combination of fixed effect coefficients and (ii) all applicable random effects, with random effects set equal to their conditional modes.

Usage

## S3 method for class 'egf'
fitted(object,
       top = egf_top(object), subset = NULL, select = NULL,
       class = FALSE, se = FALSE, ...)

## S3 method for class 'fitted.egf'
confint(object, parm = seq_len(nrow(object)), level = 0.95,
        class = FALSE, ...)

Arguments

Value

A numeric matrix containing fitted values.

Users can pass class = TRUE to obtain an augmented result. Thus, alternatively:

A data frame inheriting from class fitted.egf, with variables:

See Also

The generic function fitted.

Examples

example("egf", package = "epigrowthfit")

zz <- fitted(m1, class = TRUE, se = TRUE)
str(zz)

confint(zz, class = TRUE)

Details about Fixed Effect Coefficients

Description

Extracts from a model object details about the fixed effect coefficients, namely segment beta of the bottom level parameter vector.

Usage

## S3 method for class 'egf'
fixef(object, ...)

Arguments

Value

A data frame with one row per coefficient and variables:

See Also

The generic function fixef.

Extract Model Formulae

Description

Extracts from a model object the mixed effects model formula corresponding to a top level nonlinear model parameter.

Usage

## S3 method for class 'egf'
formula(x, top = egf_top(x), split = FALSE, ...)

Arguments

Value

By default, the mixed effects model formula corresponding to top. If split = TRUE, then the same formula without random effect terms. The deleted terms are stored in an expression vector and preserved as attribute random of the result.

See Also

The generic function formula.

Extract Model Calls

Description

Extracts from a model object the call to egf that produced it. This method exists mainly to enable compatibility with the default method of generic function update.

Usage

## S3 method for class 'egf'
getCall(x, ...)

Arguments

Value

A call to egf.

See Also

The generic function getCall.

Generation Interval Distribution

Description

Generation interval density function (dgi), distribution function (pgi), quantile function (qgi), and sampling (rgi). Results are conditional on supplied latent and infectious period distributions. It is assumed

• that the latent period and infectious waiting time are independent,

• that infectiousness is constant over the infectious period, and

• that the latent and infectious periods are positive and integer-valued (in arbitrary but common units of time).

Usage

dgi(x, latent, infectious)
pgi(q, latent, infectious)
qgi(p, latent, infectious)
rgi(n, latent, infectious)

Arguments

Value

A numeric vector with length equal to the that of the first argument or length n in the case of rgi.

References

Svensson, Å. (2007). A note on generation times in epidemic models. Mathematical Biosciences, 208(1), 300-311. doi:10.1016/j.mbs.2006.10.010

Examples

latent <- c(0.026, 0.104, 0.182, 0.246, 0.318, 0.104,
            0.013, 0.004, 0.003)
m <- length(latent)

infectious <- c(0.138, 0.462, 0.256, 0.078, 0.041, 0.007,
                0.004, 0.004, 0.006, 0.004)
n <- length(infectious)

## Histogram of samples
y <- rgi(1e06, latent, infectious)
hist(y, breaks = seq(0, m + n + 1), freq = FALSE, las = 1,
     ylab = "relative frequency",
     main = "")

## Density and distribution functions
x <- seq(0, m + n + 1, by = 0.02)
fx <- dgi(x, latent, infectious)
Fx <- pgi(x, latent, infectious)
plot(x, fx, type = "l", las = 1, # consistent with histogram
     xlab = "generation interval",
     ylab = "density function")
plot(x, Fx, type = "l", las = 1,
     xlab = "generation interval",
     ylab = "distribution function")

## Quantile function
p <- seq(0, 1, by = 0.001)
qp <- qgi(p, latent, infectious)
plot(p, qp, type = "l", las = 1,
     xlab = "probability",
     ylab = "quantile function")

Extract the Log (Marginal) Likelihood

Description

Extracts from a model object the value of the log marginal likelihood. Whether the result represents a local maximum depends on, among other things, convergence of the optimizer.

Usage

## S3 method for class 'egf'
logLik(object, ...)

Arguments

Value

A numeric vector of length 1 inheriting from class logLik. Attribute df is the number of estimated parameters (fixed effect coefficients and random effect covariance parameters). Attribute nobs is the number of observations of disease incidence used in estimation.

See Also

The generic function logLik.

Extract Model Frames

Description

Extracts from a model object any of several data frames used to specify the model, including the mixed effects model frames.

Usage

## S3 method for class 'egf'
model.frame(formula,
            which = c("ts", "windows", "parameters", "extra", "combined"),
            top = egf_top(formula), full = FALSE, ...)

Arguments

Value

A data frame.

See Also

The generic function model.frame.

Extract Design Matrices

Description

Extracts from a model object fixed and random effects design matrices.

Usage

## S3 method for class 'egf'
model.matrix(object, which = c("fixed", "random"),
             top = NULL, random = NULL, ...)

Arguments

Details

model.matrix(which = "fixed", top = NULL) returns the result of combining (in the sense of cbind) all fixed effects design matrices.

model.matrix(which = "random", top = "<name>", random = NULL) returns the result of combining all random effects design matrices associated with parameter top.

model.matrix(which = "random", top = NULL, random = NULL) returns the result of combining all random effects design matrices and permuting the columns to obtain a convenient ordering of random effect coefficients. (Coefficients are sorted by relation to a common random vector. Random vectors are sorted by relation to a common covariance matrix.)

None of these “combined” design matrices possesses attributes assign and contrasts.

Value

A (sparse) dgCMatrix or a traditional (dense) matrix, with attributes assign and contrasts except in special cases; see ‘Details’.

See Also

The generic function model.matrix.

Extract the Number of Observations

Description

Returns the number of observations of disease incidence that were used in estimation of a model. This number excludes missing values and observations not belonging to a fitting window, which, despite being preserved in model objects, do not affect estimation.

Usage

## S3 method for class 'egf'
nobs(object, ...)

Arguments

Value

An integer.

See Also

The generic function nobs.

Plot Nonlinear Mixed Effects Models of Epidemic Growth

Description

A method for printing objects of class egf.

Usage

## S3 method for class 'egf'
plot(x, type = c("interval", "cumulative", "rt"),
     time_as = c("Date", "numeric"), delta = 1, log = TRUE, zero = NA,
     show_predict = TRUE, show_doubling = FALSE, level = 0.95,
     control = egf_control_plot(), cache = NULL, plot = TRUE,
     subset = NULL, order = NULL, xlim = NULL, ylim = NULL,
     main = NULL, sub = NULL, xlab = NULL, ylab = NULL, ...)

Arguments

Details

Computation of fitted and predicted values and standard errors is performed before any plots are created. To avoid waste of computation time, cached computations are returned even if an error is thrown during plotting. To ensure that the cache is preserved, assign the result of the function call to a name:

cache <- plot(x, \dots)

. Caching functionality must be used with care, as mismatch between x and cache will not be detected. Constructions such as plot(y, cache = plot(x, ...), ...), where x and y are different egf objects, should not be expected to produce correct results.

Value

A data frame inheriting from class plot.egf.

If argument cache was supplied in the function call, then this data frame is the result of augmenting cache with any new computations.

See Also

The generic function plot.

Examples

example("egf", package = "epigrowthfit")

l <- list(legend = list(cex = 0.8),
          value = list(cex = 0.8, font = 2),
          ci = list(cex = 0.8))
control <- egf_control_plot(doubling = l)

op <- par(mar = c(3.5, 5, 5, 1))
plot(m1,
     type = "interval",
     show_predict = 2L,
     show_doubling = 2L,
     control = control)
plot(m1,
     type = "cumulative",
     main = "Fitted exponential model",
     sub = quote(paste("Country", country)))
par(op)

op <- par(mar = c(3.5, 9.5, 5, 1))
plot(m1, type = "rt", subset = quote(country %in% LETTERS[4:6]))
par(op)

Predicted Values

Description

Computes predicted values of top level nonlinear model parameters. These are conditional on an estimated nonlinear mixed effects model and, optionally, new data.

Usage

## S3 method for class 'egf'
predict(object, newdata = NULL, class = FALSE, se = FALSE, ...)

Arguments

Value

A numeric matrix containing predicted values.

Users can pass class = TRUE to obtain an augmented result. Thus, alternatively:

A data frame inheriting from class predict.egf, with variables:

See Also

The generic function predict.

Printing Model Objects

Description

A method for printing objects of class egf.

Usage

## S3 method for class 'egf'
print(x, width = 0.9 * getOption("width"), indent = 2L, ...)

Arguments

Value

The argument x, unchanged but invisible.

See Also

The generic function print.

Univariate Likelihood Profiles

Description

Computes univariate likelihood profiles of fixed effect coefficients, random effect covariance parameters, and linear combinations thereof, including population fitted values.

Usage

## S3 method for class 'egf'
profile(fitted, level = 0.95,
        A = seq_along(par), grid = 12L,
        parallel = egf_parallel(), trace = FALSE,
        top = egf_top(fitted), subset = NULL, select = NULL, ...)

## S3 method for class 'profile.egf'
confint(object, parm = seq_along(object), level = attr(object, "level"),
        class = FALSE, ...)

## S3 method for class 'profile.egf'
plot(x, parm = seq_along(x), level = attr(x, "level"),
     type = c("z", "abs(z)", "z^2"), ...)

Arguments

Details

Computation of likelihood profiles is typically expensive, requiring estimation of many restricted models. It is parallelized at the C++ level when there is OpenMP support and fitted[["control"]][["omp_num_threads"]] is set to an integer greater than 1. When there is no OpenMP support, it can still be parallelized at the R level with appropriate setting of argument parallel.

Value

A list of length nrow(A) inheriting from classes profile.egf and profile. Each element is a data frame specifying a profile, with two variables:

The confidence level level is preserved as an attribute.

See Also

The generic function profile. The more basic “next” method for generic function plot, namely plot.profile.

Examples

example("egf", package = "epigrowthfit")

zz <- profile(m1, A = NULL,
              top = "log(r)", subset = quote(country == "A" & wave == 1))
str(zz)

confint(zz, class = TRUE)

pty <- c("z", "abs(z)", "z^2")
bty <- c("l", "u", "u")
for (i in 1:3)
    plot(zz, type = pty[i], bty = bty[i], las = 1)

Details about Random Effect Coefficients

Description

Extracts from a model object details about the random effect coefficients, namely segment b of the bottom level parameter vector.

Usage

## S3 method for class 'egf'
ranef(object, makeSigma = FALSE, ...)

Arguments

Value

A data frame with one row per coefficient and variables:

If makeSigma = TRUE, then the result has attribute Sigma, a list of covariance matrices corresponding to the levels of variable cov.

See Also

The generic function ranef.

Simulation and Parametric Bootstrapping

Description

Simulates incidence data conditional on a fitted nonlinear mixed effects model of epidemic growth. Optionally re-estimates the model given the simulated data, thus generating samples from the conditional distribution of the bottom level parameter vector.

Usage

## S3 method for class 'egf'
simulate(object, nsim = 1, seed = NULL,
         bootstrap = FALSE,
         control = list(), parallel = egf_parallel(),
         trace = FALSE, ...)

Arguments

Details

Bootstrap optimizations are typically expensive for nontrivial models. They are parallelized at the C++ level when there is OpenMP support and object[["control"]][["omp_num_threads"]] is set to an integer greater than 1. When there is no OpenMP support, they can still be parallelized at the R level with appropriate setting of argument parallel.

Arguments control, parallel, and trace are unused when bootstrap = FALSE.

Value

A list inheriting from class simulate.egf, with elements:

Attribute RNGstate preserves the RNG state prior to simulation, making the result reproducible.

See Also

The generic function simulate.

Examples

example("egf", package = "epigrowthfit")

zz <- simulate(m2, nsim = 6L, seed = 181952L, bootstrap = TRUE)
str(zz)

matplot(t(zz[["bootstrap"]][!m2[["random"]], ]),
        type = "o", las = 1, xlab = "simulation", ylab = "value")

Simulating Incidence Time Series

Description

Simulates equally spaced incidence time series according to a specified nonlinear model. Top level nonlinear model parameters vary between time series according to a fixed intercept model ~ts or random intercept model ~(1 | ts).

Usage

## S3 method for class 'egf_model'
simulate(object, nsim = 1, seed = NULL,
         mu, Sigma = NULL, tol = 1e-06, cstart = 0, tend = 100, ...)

Arguments

Value

A list inheriting from class simulate.egf_model, with elements:

Attribute RNGstate preserves the RNG state prior to simulation, making the result reproducible.

See Also

The generic function simulate.

Examples

r <- 0.04
c0 <- 400
s <- 0.2

mu <- log(c(r, c0))
Sigma <- diag(rep.int(s^2, length(mu)))

zz <- simulate(object = egf_model(curve = "exponential", family = "pois"),
               nsim = 20L,
               seed = 202737L,
               mu = mu,
               Sigma = Sigma,
               cstart = 10)
str(zz)

mm <- egf(zz)
(pp <- cbind(actual = coef(zz), fitted = coef(mm)))

Model Summaries

Description

Summarizes fitted values of top level nonlinear model parameters and gathers diagnostic information that can be used to quickly assess convergence of the optimizer.

Usage

## S3 method for class 'egf'
summary(object, ...)

Arguments

Value

A list inheriting from class summary.egf, with elements:

See Also

The generic function summary.

Examples

example("egf", package = "epigrowthfit")

zz <- summary(m1)
str(zz)

Model Terms

Description

Extracts the terms object corresponding to a top level nonlinear model parameter.

Usage

## S3 method for class 'egf'
terms(x, top = egf_top(x), ...)

Arguments

Value

A terms object.

See Also

The generic function terms.

Compute the Characteristic Time Scale

Description

Computes characteristic time scales corresponding to exponential growth rates.

Usage

timescale(r, units)

Arguments

Value

1/r, as a difftime if units is not missing.

See Also

R0, finalsize.

Examples

r <- 10^(-2:0)
units <- "days"
stopifnot(all.equal(timescale(r), 1 / r),
          all.equal(timescale(r, units), .difftime(1 / r, units)))

Model Covariance Matrices

Description

Extracts (or, if necessary, computes) the covariance matrix of bottom level parameters beta and theta, corresponding to the output of coef(object).

Usage

## S3 method for class 'egf'
vcov(object, ...)

Arguments

Details

If the resulting matrix is not finite or not positive definite, then the fit specified by object should be investigated, as the optimizer that produced the fit may have failed to converge.

Value

A real, symmetric matrix.

See Also

The generic function vcov.