Package {primarycensored}

Title:	Primary Event Censored Distributions
Version:	1.5.2
Description:	Provides functions for working with primary event censored distributions and 'Stan' implementations for use in Bayesian modeling. Primary event censored distributions are useful for modeling delayed reporting scenarios in epidemiology and other fields (Charniga et al. (2024) <doi:10.48550/arXiv.2405.08841>). It also provides support for arbitrary delay distributions, a range of common primary distributions, and allows for truncation and secondary event censoring to be accounted for (Park et al. (2024) <doi:10.1101/2024.01.12.24301247>). A subset of common distributions also have analytical solutions implemented, allowing for faster computation. In addition, it provides multiple methods for fitting primary event censored distributions to data via optional dependencies.
License:	MIT + file LICENSE
URL:	https://primarycensored.epinowcast.org, https://github.com/epinowcast/primarycensored
BugReports:	https://github.com/epinowcast/primarycensored/issues
Depends:	R (≥ 4.0.0)
Imports:	lifecycle
Suggests:	bookdown, cmdstanr, dplyr, fitdistrplus, flexsurv, knitr, ggplot2, rmarkdown, spelling, testthat (≥ 3.1.9), usethis, withr
Additional_repositories:	https://stan-dev.r-universe.dev
Config/Needs/hexsticker:	hexSticker, sysfonts, ggplot2
Config/Needs/website:	r-lib/pkgdown, epinowcast/enwtheme
Config/testthat/edition:	3
Encoding:	UTF-8
Language:	en-GB
LazyData:	true
VignetteBuilder:	knitr
RoxygenNote:	7.3.3
NeedsCompilation:	no
Packaged:	2026-09-10 21:05:08 UTC; seabbs
Author:	Sam Abbott [aut, cre, cph], Sam Brand [aut], Adam Howes [ctb], James Mba Azam [aut], Carl Pearson [aut], Sebastian Funk [aut], Kelly Charniga [aut]
Maintainer:	Sam Abbott <contact@samabbott.co.uk>
Repository:	CRAN
Date/Publication:	2026-09-11 03:00:02 UTC

primarycensored: Primary Event Censored Distributions

Description

Provides functions for working with primary event censored distributions and 'Stan' implementations for use in Bayesian modeling. Primary event censored distributions are useful for modeling delayed reporting scenarios in epidemiology and other fields (Charniga et al. (2024) doi:10.48550/arXiv.2405.08841). It also provides support for arbitrary delay distributions, a range of common primary distributions, and allows for truncation and secondary event censoring to be accounted for (Park et al. (2024) doi:10.1101/2024.01.12.24301247). A subset of common distributions also have analytical solutions implemented, allowing for faster computation. In addition, it provides multiple methods for fitting primary event censored distributions to data via optional dependencies.

Author(s)

Maintainer: Sam Abbott contact@samabbott.co.uk (ORCID) [copyright holder]

Authors:

• Sam Brand usi1@cdc.gov (ORCID)

• James Mba Azam james.azam@lshtm.ac.uk (ORCID)

• Carl Pearson carl.ab.pearson@gmail.com (ORCID)

• Sebastian Funk sebastian.funk@lshtm.ac.uk (ORCID)

• Kelly Charniga kelly.charniga@gmail.com (ORCID)

Other contributors:

• Adam Howes adamthowes@gmail.com (ORCID) [contributor]

See Also

Useful links:

• https://primarycensored.epinowcast.org

• https://github.com/epinowcast/primarycensored

• Report bugs at https://github.com/epinowcast/primarycensored/issues

Build the (dpcens_dist, ppcens_dist) closure pair for fitdistdoublecens

Description

Single path: scalar named parameters are gathered from the closure's environment, optionally folded into the vector argument named by vector_param via the distribution's param_transform, and dispatched through .dpcens/.ppcens. Formals on the closures are derived from the supplied start list (parametric) or from a vector_param-aware naming convention (non-parametric).

Usage

.build_pcens_closures(
  pdist,
  ddist,
  params,
  dprimary,
  primary_args,
  pprimary = NULL,
  vector_param,
  param_transform = NULL,
  fit_penalty,
  prior,
  N,
  start,
  pdist_extras = list(),
  check_once = function() FALSE
)

Arguments

Validate truncation bounds L and D

Description

Internal function to validate that L (lower truncation) and D (upper truncation) parameters are valid. L must be less than D. L may be negative or -Inf for delay distributions with support below zero.

Usage

.check_truncation_bounds(L, D)

Arguments

Value

Invisible NULL if valid, otherwise stops with an error message.

Validate truncation bounds in a data frame

Description

Internal function to validate that L (lower truncation) is less than D (upper truncation) for all rows in a data frame.

Usage

.check_truncation_bounds_df(data, L_col, D_col)

Arguments

Value

Invisible NULL if valid, otherwise stops with an error message.

Define a fitdistrplus compatible wrapper around dprimarycensored

Description

Define a fitdistrplus compatible wrapper around dprimarycensored

Usage

.dpcens(
  x,
  params,
  pdist,
  dprimary,
  primary_args,
  pprimary = NULL,
  check = TRUE,
  ...
)

Arguments

Extract base function name

Description

This helper function extracts the base name of a function, removing namespace prefixes. Base R distribution functions are identified from the C routine they call. Functions exported by another package (for example flexsurv::pgengamma.orig()) are identified by their exported name.

Usage

.extract_function_name(func)

Arguments

Value

Character string representing the base name of the function, or "unknown" if it cannot be determined.

Extract function names or content from Stan code

Description

Extract function names or content from Stan code

Usage

.extract_stan_functions(content, names_only = FALSE, functions = NULL)

Arguments

Value

Character vector of function names or content

Extract and Combine Distribution Names

Description

This helper function attempts to determine distribution names and uses those to establish a class hierarchy for potential analytical solutions.

Usage

.format_class(pdist, dprimary)

Arguments

Value

A character vector of class names: specific (delay + primary), delay-only, and base class.

Shared MAP penalty for the logit-hazard families

Description

Both RW and RE variants share priors on alpha, log_sigma, and eps_*; only how the hazards are built from those scalars differs (see .make_hazard_transform()).

Usage

.hazard_fit_penalty(par_named, N, prior_settings = NULL)

Test whether a numeric vector lies on the simplex

Description

Returns TRUE when the vector contains no missing values, no negative entries, and sums to 1 within 10^{-8}. Used by the discretestep family to apply a soft penalty inside fitting closures rather than erroring.

Usage

.is_valid_simplex(p, tol = 1e-08)

Look up the primary event CDF from the registry

Description

Given a primary event density function dprimary, looks up the corresponding CDF function from pcd_primary_distributions using the "name" attribute. Returns NULL silently when no match is found so that callers can fall back to numerical integration.

Usage

.lookup_pprimary(dprimary)

Arguments

Value

A function (the primary CDF) or NULL.

Build a hazard-vector transform for a given model

Description

Returns a function that maps a named list of scalar parameters (alpha, log_sigma, eps_1, ..., eps_{K-1}) to a length-K hazard vector with the final entry pinned to 1. model = "rw" uses a Gaussian random walk on the logit hazard; model = "re" treats the innovations as IID logit random effects around the intercept.

Usage

.make_hazard_transform(model = c("rw", "re"))

Normalise a primary event censored CDF

Description

Internal function to normalise a primary event censored CDF when truncation is applied. The CDF is normalised using (F(q) - F(L)) / (F(D) - F(L)) and values outside [L, D] are clamped to 0 or 1.

Usage

.normalise_cdf(result, q, L, D, pcens_obj, pwindow)

Arguments

Value

Normalised CDF values as a numeric vector

Analytical primary event censored CDF for the generalised gamma

Description

Shared implementation for the Stacy parameterisation used by both generalised gamma pcens_cdf() methods.

Usage

.pcens_cdf_gengamma_unif(q, pwindow, shape, scale, k)

Arguments

Value

Vector of computed primary event censored CDFs

Define a fitdistrplus compatible wrapper around pprimarycensored

Description

Define a fitdistrplus compatible wrapper around pprimarycensored

Usage

.ppcens(
  q,
  params,
  pdist,
  dprimary,
  primary_args,
  pprimary = NULL,
  check = TRUE,
  ...
)

Arguments

Resolve the prior settings for the hazard fit_penalty

Description

Defaults: alpha ~ N(0, 5), log_sigma ~ N(log(0.5), 1). User overrides are merged in: each component of prior_settings may itself be a list with mean and sd entries.

Usage

.resolve_hazard_prior(prior_settings = NULL)

Resolve a delay distribution function from a name or function

Description

Accepts either a function (returned as-is, with its existing "name" attribute preserved) or a character string that is looked up against pcd_distributions. When a string is supplied, the corresponding base R p<name> function is returned with the "name" attribute attached so analytical solutions can dispatch.

Usage

.resolve_pdist(pdist, type = c("p", "d"))

Arguments

Value

A function with a "name" attribute.

Resolve the primary CDF, validating against dprimary if both supplied

Description

Returns the primary CDF to use. If the user supplies pprimary explicitly (either a function or a string name), it is returned (after resolving the string via pcd_dist_name). When both dprimary and pprimary carry a "name" attribute, the names must agree on everything except the leading d/p; otherwise we error to catch typos like dunif + pexpgrowth. If pprimary is not supplied, falls back to a registry lookup against dprimary via .lookup_pprimary, which may return NULL.

Usage

.resolve_pprimary(dprimary, pprimary = NULL)

Arguments

Value

A primary CDF function, or NULL if no match was found.

Resolve primary_args / dprimary_args with a deprecation

Description

The resolver distinguishes "user supplied nothing" from "user supplied an empty list" so the deprecated dprimary_args path can be detected. The returned list is never NULL.

Usage

.resolve_primary_args(
  primary_args,
  dprimary_args,
  fn,
  env = parent.frame(),
  user_env = parent.frame(2)
)

Arguments

Details

The deprecation is soft (lifecycle::deprecate_soft()): a warning is shown when the exported function is called from the global environment or from the package under test, and calls from other packages stay silent.

Value

A list (possibly empty) of primary distribution arguments.

Sort eps_* parameter names by trailing numeric suffix

Description

Lexicographic sorting of c("eps_1", ..., "eps_10") returns eps_10 before eps_2, which would scramble the random-walk innovations when K > 10. This helper extracts the trailing integer and orders by it.

Usage

.sort_eps_names(eps_names)

Validate boundaries and pmf for step distribution functions

Description

Validate boundaries and pmf for step distribution functions

Usage

.validate_discretestep_args(boundaries, pmf)

Helper method for custom distributions

Description

pprimarycensored() and related functions can identify which distributions are provided via the pdist and dprimary arguments when those are base R functions (e.g. punif, dexp) via the name attribute.

Usage

add_name_attribute(func, name)

Arguments

Details

If you need to use a non-base R implementation, but know the distribution name, you can use this helper function to set it in a way that will be detected by pprimarycensored() and related functions.

This is useful as it enables the automatic use of analytical solutions for distributions where they exist. You can check which analytical solutions are available using methods(pcens_cdf) and check distribution names using pcd_dist_name().

Value

Function, with a "name" attribute added

See Also

Utility functions for working with distributions pcd_dist_name(), pcd_distributions, pcd_primary_distributions

Examples

dist <- add_name_attribute(pnorm, "hello")
attr(dist, "name")

Check if a function is a valid bounded probability density function (PDF)

Description

This function tests whether a given function behaves like a valid PDF by checking if it integrates to approximately 1 over the specified range and if it takes the arguments min and max.

Usage

check_dprimary(dprimary, pwindow, dprimary_args = list(), tolerance = 0.001)

Arguments

Value

NULL. The function will stop execution with an error message if dprimary is not a valid PDF.

See Also

Distribution checking functions check_pdist(), check_truncation()

Examples

check_dprimary(dunif, pwindow = 1)

Check if a function is a valid cumulative distribution function (CDF)

Description

This function tests whether a given function behaves like a valid CDF by checking if it's monotonically increasing and bounded between 0 and 1.

Usage

check_pdist(pdist, D = Inf, ...)

Arguments

Value

NULL. The function will stop execution with an error message if pdist is not a valid CDF.

See Also

Distribution checking functions check_dprimary(), check_truncation()

Examples

check_pdist(pnorm, D = 10)

Check if truncation time is appropriate relative to the maximum delay

Description

This function checks if the truncation time D is appropriate relative to the maximum delay. If D is much larger than necessary, it suggests considering setting it to Inf for better efficiency with minimal accuracy cost.

Usage

check_truncation(delays, D, multiplier = 2)

Arguments

Value

Invisible NULL. Prints a message if the condition is met.

See Also

Distribution checking functions check_dprimary(), check_pdist()

Examples

check_truncation(delays = c(1, 2, 3, 4), D = 10, multiplier = 2)

Hazard-parameterised piecewise-constant PMF

Description

Returns the probability mass for each value in x. Converts hazards to a PMF via hazards_to_pmf() then delegates to ddiscretestep().

Usage

ddiscretehazard(x, boundaries = NULL, hazards)

Arguments

Details

Outside fitting it is a deterministic wrapper around ddiscretestep(). It earns its keep as a fitting parameterisation in fitdistdoublecens() because the random walk on the logit hazard smooths the recovered PMF.

Value

Numeric vector of PMF values, the same length as x.

See Also

ddiscretestep(), hazards_to_pmf(), pmf_to_hazards(), fitdistdoublecens()

Other pdiscretehazard: discretehazard_start(), pdiscretehazard(), rdiscretehazard()

Examples

hazards <- c(0.3, 0.5, 1)
ddiscretehazard(c(1, 2, 3), boundaries = 0:3, hazards = hazards)

Step (piecewise-constant) PMF

Description

Returns the probability mass for each value in x. Mass pmf[i] is located at boundaries[i+1] (the right edge of bin i); all other values return 0.

Usage

ddiscretestep(x, boundaries = NULL, pmf)

Arguments

Details

Like pdiscretestep, this function applies a soft simplex penalty: if pmf contains negative entries or fails to sum to 1 (within 10^{-8}), the function returns near-zero density (.Machine$double.eps) rather than erroring. This makes it safe to call from inside fitdistrplus::fitdist() closures driven by fitdistdoublecens.

Value

Numeric vector of PMF values, the same length as x.

See Also

Other pdiscretestep: hazards_to_pmf(), pdiscretestep(), pmf_to_hazards(), rdiscretestep()

Examples

ddiscretestep(c(0, 1, 2, 3), boundaries = 0:3, pmf = c(0.2, 0.5, 0.3))

Start values for the logit-hazard parameterisation

Description

Builds a named list of starting values for fitdistdoublecens() with distr = "discretehazard". The free parameters are the logit intercept alpha, the log random-walk (or random-effect) scale log_sigma, and the K - 1 innovations eps_1, ..., eps_{K-1}. The final bin hazard is pinned to 1 inside the parameterisation so the implied PMF sums to 1.

Usage

discretehazard_start(K, alpha = -2, log_sigma = log(1), eps = 0)

Arguments

Value

Named list suitable for the start argument of fitdistdoublecens().

See Also

Other pdiscretehazard: ddiscretehazard(), pdiscretehazard(), rdiscretehazard()

Examples

discretehazard_start(K = 5)

Compute the primary event censored PMF for delays

Description

This function computes the primary event censored probability mass function (PMF) for a given set of quantiles. It adjusts the PMF of the primary event distribution by accounting for the delay distribution and potential truncation at a maximum delay (D) and minimum delay (L). The function allows for custom primary event distributions and delay distributions.

Usage

dprimarycensored(
  x,
  pdist,
  pwindow = 1,
  swindow = 1,
  L = -Inf,
  D = Inf,
  dprimary = dunif,
  primary_args = NULL,
  pprimary = NULL,
  dprimary_args = NULL,
  log = FALSE,
  ...,
  check = TRUE
)

dpcens(
  x,
  pdist,
  pwindow = 1,
  swindow = 1,
  L = -Inf,
  D = Inf,
  dprimary = dunif,
  primary_args = NULL,
  pprimary = NULL,
  dprimary_args = NULL,
  log = FALSE,
  ...,
  check = TRUE
)

Arguments

Details

The primary event censored PMF is computed by taking the difference of the primary event censored cumulative distribution function (CDF) at two points, d + \text{swindow} and d. The primary event censored PMF, f_{\text{cens}}(d), is given by:

f_{\text{cens}}(d) = F_{\text{cens}}(d + \text{swindow}) - F_{\text{cens}}(d)

where F_{\text{cens}} is the primary event censored CDF.

The function first computes the CDFs for all unique points (including both d and d + \text{swindow}) using pprimarycensored(). It then creates a lookup table for these CDFs to efficiently calculate the PMF for each input value. For delays less than L, the function returns 0.

When the secondary censoring interval extends past the upper truncation point (d + \text{swindow} > D) but the lower endpoint satisfies d < D, the upper endpoint is internally clipped to D before evaluating the CDF. The likelihood for such an observation is P(X \in [d, \min(d + \text{swindow}, D)] \mid L \le X \le D), which equals the usual interval probability when d + \text{swindow} \le D. This avoids erroring when an observation's secondary window straddles the truncation point (relevant for non-parametric delays such as pdiscretestep()).

Observations with d \ge D are rejected with an error: under the truncation X \le D, no event with latent value d \ge D is observable, and accepting such inputs would otherwise yield a 0/0 likelihood.

The PMF is normalised to ensure it sums to 1 over the range [L, D\). This normalization uses:

f_{\text{cens,norm}}(d) = \frac{f_{\text{cens}}(d)}{ F_{\text{cens}}(D) - F_{\text{cens}}(L)}

where f_{\text{cens,norm}}(d) is the normalized PMF. For the explanation and mathematical details of the CDF, refer to the documentation of pprimarycensored().

Value

Vector of primary event censored PMFs, normalized over [L, D] if truncation is applied

See Also

Primary event censored distribution functions pprimarycensored(), qprimarycensored(), rprimarycensored()

Examples

# Example: Weibull distribution with uniform primary events
dprimarycensored(c(0.1, 0.5, 1), pweibull, shape = 1.5, scale = 2.0)

# Example: Weibull distribution with exponential growth primary events
dprimarycensored(
  c(0.1, 0.5, 1), pweibull,
  dprimary = dexpgrowth,
  primary_args = list(r = 0.2), shape = 1.5, scale = 2.0
)

# Example: Left-truncated distribution (e.g., for generation intervals)
dprimarycensored(1:9, pweibull, L = 1, D = 10, shape = 1.5, scale = 2.0)

Exponential growth distribution functions

Description

Density, distribution function, and random generation for the exponential growth distribution.

Usage

dexpgrowth(x, min = 0, max = 1, r, log = FALSE)

pexpgrowth(q, min = 0, max = 1, r, lower.tail = TRUE, log.p = FALSE)

rexpgrowth(n, min = 0, max = 1, r)

Arguments

Details

The exponential growth distribution is defined on the interval [min, max] with rate parameter (r). Its probability density function (PDF) is:

f(x) = \frac{r \cdot \exp(r \cdot x)}{\exp(r \cdot max) - \exp(r \cdot min)}

The cumulative distribution function (CDF) is:

F(x) = \frac{\exp(r \cdot x) - \exp(r \cdot min)}{ \exp(r \cdot max) - \exp(r \cdot min)}

For random number generation, we use the inverse transform sampling method:

1. Generate u \sim \text{Uniform}(0,1)

2. Set F(x) = u and solve for x:

   x = \frac{1}{r} \cdot \log(u \cdot \exp(r \cdot max) + (1 - u) \cdot \exp(r \cdot min))

This method works because of the probability integral transform theorem, which states that if X is a continuous random variable with CDF F(x), then Y = F(X) follows a \text{Uniform}(0,1) distribution. Conversely, if U is a \text{Uniform}(0,1) random variable, then F^{-1}(U) has the same distribution as X, where F^{-1} is the inverse of the CDF.

In our case, we generate u from \text{Uniform}(0,1), then solve F(x) = u for x to get a sample from our exponential growth distribution. The formula for x is derived by algebraically solving the equation:

u = \frac{\exp(r \cdot x) - \exp(r \cdot min)}{\exp(r \cdot max) - \exp(r \cdot min)}

When r is very close to 0 (|r| < 1e-10), the distribution approximates a uniform distribution on [min, max], and we use a simpler method to generate samples directly from this uniform distribution.

Value

dexpgrowth gives the density, pexpgrowth gives the distribution function, and rexpgrowth generates random deviates.

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

Examples

x <- seq(0, 1, by = 0.1)
probs <- dexpgrowth(x, r = 0.2)
cumprobs <- pexpgrowth(x, r = 0.2)
samples <- rexpgrowth(100, r = 0.2)

Fit a distribution to doubly censored data

Description

This function wraps the custom approach for fitting distributions to doubly censored data using fitdistrplus and primarycensored. It handles primary censoring (when the primary event time is not known exactly), secondary censoring (when the secondary event time is interval-censored), and truncation (when events are only observed within a delay range [L, D]).

Usage

fitdistdoublecens(
  censdata,
  distr,
  left = "left",
  right = "right",
  pwindow = "pwindow",
  L = "L",
  D = "D",
  dprimary = dunif,
  primary_args = NULL,
  pprimary = NULL,
  dprimary_args = NULL,
  truncation_check_multiplier = 2,
  prior = NULL,
  hazard_model = c("rw", "re"),
  check = TRUE,
  ...
)

Arguments

Details

How distribution functions are resolved

The distr argument names a distribution. The function looks up the density and CDF functions by prepending d and p to the name (e.g. distr = "gamma" resolves to dgamma() and pgamma()). Custom distributions can be used as long as the corresponding ⁠d<distr>()⁠ and ⁠p<distr>()⁠ functions are defined.

Non-parametric distributions

Two non-parametric distributions are supported. They share a common fitting machinery: the dist function carries a vector_param attribute ("pmf" for pdiscretestep()/ddiscretestep(), "hazards" for pdiscretehazard()/ddiscretehazard()) that drives this function to build a closure mapping flat scalar parameters into the underlying vector argument.

• distr = "discretestep": free parameters ⁠p1, ..., p_{K-1}⁠ (in ⁠[0, 1]⁠); the last bin probability is ⁠1 - sum(p1, ..., p_{K-1})⁠. See pdiscretestep() for parameterisation details and the soft simplex penalty applied when probabilities are infeasible.

• distr = "discretehazard": free parameters alpha, log_sigma, ⁠eps_1, ..., eps_{K-1}⁠. The hazard form parameterises the same family of step distributions as "discretestep", but its free parameters drive either a Gaussian random walk on the logit hazard (hazard_model = "rw", the default, logit(h_i) = alpha + sigma * cumsum(eps)) or an IID logit random-effect transform (hazard_model = "re", logit(h_i) = alpha + sigma * eps_i with eps_i ~ N(0, 1)). The smoothing of the random walk regularises the recovered PMF against over-fitting in sparse data and replaces the simplex constraint with an unconstrained optimisation; the random-effect variant models hazards as independent draws around alpha rather than a smoothed trajectory. See pdiscretehazard() for full parameterisation details and the MAP-equivalent prior penalty applied during fitting; pass prior to override the default prior settings.

For non-parametric distributions K is implied by length(start): K = length(start) + 1 for "discretestep" and K = length(start) - 1 for "discretehazard". start is therefore required.

Value

An object of class "fitdist" as returned by fitdistrplus::fitdist.

See Also

pdiscretestep() pdiscretehazard()

Modelling wrappers for external fitting packages pcd_as_stan_data(), pcd_cmdstan_model()

Examples

# Example with normal distribution
set.seed(123)
n <- 1000
true_mean <- 5
true_sd <- 2
pwindow <- 2
swindow <- 2
D <- 10
samples <- rprimarycensored(
  n, rnorm,
  mean = true_mean, sd = true_sd,
  pwindow = pwindow, swindow = swindow, D = D
)

delay_data <- data.frame(
  left = samples,
  right = samples + swindow,
  pwindow = rep(pwindow, n),
  D = rep(D, n)
)

fit_norm <- fitdistdoublecens(
  delay_data,
  distr = "norm",
  start = list(mean = 0, sd = 1)
)

summary(fit_norm)


# Example with discretestep (non-parametric PMF) distribution
set.seed(42)
true_pmf <- c(0.1, 0.3, 0.4, 0.15, 0.05)
step_samples <- rprimarycensored(
  500, rdiscretestep,
  boundaries = 0:5, pmf = true_pmf,
  pwindow = 1, swindow = 1, D = 6
)
step_data <- data.frame(
  left = step_samples,
  right = step_samples + 1,
  pwindow = rep(1, 500),
  D = rep(6, 500)
)
fit_step <- fitdistdoublecens(
  step_data,
  distr = "discretestep",
  boundaries = 0:5,
  start = as.list(setNames(rep(0.2, 4), paste0("p", 1:4)))
)

# Example with discretehazard (logit-hazard random walk) distribution
fit_haz <- fitdistdoublecens(
  step_data,
  distr = "discretehazard",
  boundaries = 0:5,
  start = c(
    list(alpha = -2, log_sigma = log(0.5)),
    as.list(setNames(rep(0, 4), paste0("eps_", 1:4)))
  )
)

Convert discrete-time hazards to a PMF

Description

Given a vector of discrete-time hazards h_1, \ldots, h_K, computes the corresponding PMF:

pmf_i = h_i \prod_{j < i} (1 - h_j)

Usage

hazards_to_pmf(hazards)

Arguments

Details

The last hazard must equal 1 (to ensure the PMF sums to 1). If a vector of length K-1 is supplied (all hazards except the final exit hazard), 1 is appended automatically.

Value

Numeric vector of PMF values.

See Also

Other pdiscretestep: ddiscretestep(), pdiscretestep(), pmf_to_hazards(), rdiscretestep()

Examples

hazards_to_pmf(c(0.2, 0.3, 1))

S3 class for primary event censored distribution computation

Description

S3 class for primary event censored distribution computation

Usage

new_pcens(
  pdist,
  dprimary,
  primary_args = NULL,
  pprimary = NULL,
  dprimary_args = NULL,
  ...
)

Arguments

Value

An object with class hierarchy c("pcens_{pdist_name}_{dprimary_name}", "pcens_{pdist_name}", "pcens"). This contains the primary event distribution, the delay distribution, the delay distribution arguments, the primary event CDF (if available), and any additional arguments. It can be used with the pcens_cdf() function to compute the primary event censored CDF.

See Also

Low level primary event censored distribution objects and methods pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile(), pcens_quantile.default()

Examples

new_pcens(
  pdist = pgamma, dprimary = dunif,
  primary_args = list(min = 0, max = 1),
  shape = 1, scale = 1
)

Prepare data for primarycensored Stan model

Description

This function takes in delay data and prepares it for use with the primarycensored Stan model.

Usage

pcd_as_stan_data(
  data,
  delay = "delay",
  delay_upper = "delay_upper",
  n = "n",
  pwindow = "pwindow",
  start_relative_obs_time = "start_relative_obs_time",
  relative_obs_time = "relative_obs_time",
  dist_id,
  primary_id,
  param_bounds,
  primary_param_bounds,
  priors,
  primary_priors,
  compute_log_lik = FALSE,
  use_reduce_sum = FALSE,
  truncation_check_multiplier = 2,
  dist_options = NULL
)

Arguments

Value

A list containing the data formatted for use with pcd_cmdstan_model()

See Also

Modelling wrappers for external fitting packages fitdistdoublecens(), pcd_cmdstan_model()

Examples

data <- data.frame(
  delay = c(1, 2, 3),
  delay_upper = c(2, 3, 4),
  n = c(10, 20, 15),
  pwindow = c(1, 1, 2),
  relative_obs_time = c(10, 10, 10)
)
stan_data <- pcd_as_stan_data(
  data,
  dist_id = 1,
  primary_id = 1,
  param_bounds = list(lower = c(0, 0), upper = c(10, 10)),
  primary_param_bounds = list(lower = numeric(0), upper = numeric(0)),
  priors = list(location = c(1, 1), scale = c(1, 1)),
  primary_priors = list(location = numeric(0), scale = numeric(0))
)

Create a CmdStanModel with primarycensored Stan functions

Description

This function creates a CmdStanModel object using the Stan model and functions from primarycensored and optionally includes additional user-specified Stan files.

Usage

pcd_cmdstan_model(include_paths = primarycensored::pcd_stan_path(), ...)

Arguments

Details

The underlying Stan model (pcens_model.stan) supports various features:

• Multiple probability distributions for modeling delays

• Primary and secondary censoring

• Truncation

• Optional use of reduce_sum for improved performance (via within chain parallelism).

• Flexible prior specifications

• Optional computation of log-likelihood for model comparison

Value

A CmdStanModel object.

See Also

Modelling wrappers for external fitting packages fitdistdoublecens(), pcd_as_stan_data()

Examples

if (!is.null(cmdstanr::cmdstan_version(error_on_NA = FALSE))) {
  model <- pcd_cmdstan_model(compile = FALSE)
  model
}

Get distribution function cdf or pdf name

Description

Get distribution function cdf or pdf name

Usage

pcd_dist_name(name, type = c("delay", "primary"))

Arguments

Value

String distribution function name or NA if no base R implementation

See Also

Utility functions for working with distributions add_name_attribute(), pcd_distributions, pcd_primary_distributions

Examples

pcd_dist_name("lnorm")
pcd_dist_name("lognormal")
pcd_dist_name("gamma")
pcd_dist_name("weibull")
pcd_dist_name("exp")
pcd_dist_name("unif", type = "primary")
pcd_dist_name("expgrowth", type = "primary")

Supported delay distributions

Description

A dataset containing information about the supported delay distributions in primarycensored. Includes both distributions with base R implementations and those only available in Stan. Distributions beyond these are not supported in the stan code but any user functions can be used in the R code.

Usage

pcd_distributions

Format

A data.frame with 28 rows and 4 columns:

name

Distribution name

pdist

R distribution function name (e.g. plnorm), NA if there is no base R implementation

aliases

Alternative names/identifiers

stan_id

Stan distribution ID used in the stan code

See Also

Utility functions for working with distributions add_name_attribute(), pcd_dist_name(), pcd_primary_distributions

Load Stan functions as a string

Description

Load Stan functions as a string

Usage

pcd_load_stan_functions(
  functions = NULL,
  stan_path = primarycensored::pcd_stan_path(),
  wrap_in_block = FALSE,
  write_to_file = FALSE,
  output_file = "pcd_functions.stan",
  dependencies = FALSE
)

Arguments

Value

A character string containing the requested Stan functions

See Also

Tools for working with package Stan functions pcd_stan_dist_id(), pcd_stan_files(), pcd_stan_function_deps(), pcd_stan_functions(), pcd_stan_path()

Supported primary event distributions

Description

A dataset containing information about the supported primary event distributions in primarycensored. Distributions beyond these are not supported in the stan code but any user functions can be used in the R code.

Usage

pcd_primary_distributions

Format

A data.frame with 2 rows and 4 columns:

name

Distribution name

dprimary

R density function name

aliases

Alternative names/identifiers

stan_id

Stan distribution ID used in the stan code

See Also

Utility functions for working with distributions add_name_attribute(), pcd_dist_name(), pcd_distributions

Get distribution stan ID by name

Description

Get distribution stan ID by name

Usage

pcd_stan_dist_id(name, type = c("delay", "primary"))

Arguments

Value

Numeric distribution ID

See Also

Tools for working with package Stan functions pcd_load_stan_functions(), pcd_stan_files(), pcd_stan_function_deps(), pcd_stan_functions(), pcd_stan_path()

Examples

pcd_stan_dist_id("lnorm")
pcd_stan_dist_id("lognormal")
pcd_stan_dist_id("gamma")
pcd_stan_dist_id("weibull")
pcd_stan_dist_id("exp")
pcd_stan_dist_id("unif", type = "primary")

Get Stan files containing specified functions

Description

This function retrieves Stan files from a specified directory, optionally filtering for files that contain specific functions.

Usage

pcd_stan_files(functions = NULL, stan_path = primarycensored::pcd_stan_path())

Arguments

Value

A character vector of file paths to Stan files.

See Also

Tools for working with package Stan functions pcd_load_stan_functions(), pcd_stan_dist_id(), pcd_stan_function_deps(), pcd_stan_functions(), pcd_stan_path()

Get dependencies for a Stan function

Description

Returns all Stan functions that the specified function depends on, in topological order (dependencies before the functions that use them).

Usage

pcd_stan_function_deps(
  function_name,
  stan_path = primarycensored::pcd_stan_path()
)

Arguments

Value

A character vector of function names that the specified function depends on, ordered so that dependencies come before functions that use them. The requested function itself is included as the last element.

See Also

Tools for working with package Stan functions pcd_load_stan_functions(), pcd_stan_dist_id(), pcd_stan_files(), pcd_stan_functions(), pcd_stan_path()

Examples

# See what primarycensored_lpmf depends on
pcd_stan_function_deps("primarycensored_lpmf")

# A function with no dependencies
pcd_stan_function_deps("expgrowth_pdf")

Get Stan function names from Stan files

Description

This function reads all Stan files in the specified directory and extracts the names of all functions defined in those files.

Usage

pcd_stan_functions(stan_path = primarycensored::pcd_stan_path())

Arguments

Value

A character vector containing unique names of all functions found in the Stan files.

See Also

Tools for working with package Stan functions pcd_load_stan_functions(), pcd_stan_dist_id(), pcd_stan_files(), pcd_stan_function_deps(), pcd_stan_path()

Get the path to the Stan code

Description

Get the path to the Stan code

Usage

pcd_stan_path()

Value

A character string with the path to the Stan code

See Also

Tools for working with package Stan functions pcd_load_stan_functions(), pcd_stan_dist_id(), pcd_stan_files(), pcd_stan_function_deps(), pcd_stan_functions()

Compute primary event censored CDF

Description

This function dispatches to either analytical solutions (if available) or numerical integration via the default method. To see which combinations have analytical solutions implemented, use methods(pcens_cdf). For example, pcens_cdf.gamma_unif indicates an analytical solution exists for gamma delay with uniform primary event distributions.

Usage

pcens_cdf(object, q, pwindow, use_numeric = FALSE)

Arguments

Value

Vector of computed primary event censored CDFs

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile(), pcens_quantile.default()

Default method for computing primary event censored CDF

Description

This method serves as a fallback for combinations of delay and primary event distributions that don't have specific implementations. It uses a numeric integration method.

Usage

## Default S3 method:
pcens_cdf(object, q, pwindow, use_numeric = FALSE)

Arguments

Details

This method implements the numerical integration approach for computing the primary event censored CDF. It uses the same mathematical formulation as described in the details section of pprimarycensored(), but applies numerical integration instead of analytical solutions.

Value

Vector of computed primary event censored CDFs

See Also

pprimarycensored() for the mathematical details of the primary event censored CDF computation.

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile(), pcens_quantile.default()

Examples

# Create a primarycensored object with gamma delay and uniform primary
pcens_obj <- new_pcens(
  pdist = pgamma,
  dprimary = dunif,
  primary_args = list(min = 0, max = 1),
  shape = 3,
  scale = 2
)

# Compute CDF for a single value
pcens_cdf(pcens_obj, q = 9, pwindow = 1)

# Compute CDF for multiple values
pcens_cdf(pcens_obj, q = c(4, 6, 8), pwindow = 1)

Method for hazard CDF delay with general primary event distribution

Description

Computes the analytic primary event censored CDF for a hazard-parameterised piecewise-constant delay distribution. Converts hazards to a PMF via hazards_to_pmf() then dispatches back through pcens_cdf() using a freshly constructed step-distribution object. The same primary event distribution and arguments are preserved.

Usage

## S3 method for class 'pcens_pdiscretehazard'
pcens_cdf(object, q, pwindow, use_numeric = FALSE)

Arguments

Value

Vector of computed primary event censored CDFs

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile(), pcens_quantile.default()

Method for step CDF delay with general primary event distribution

Description

Computes the analytic primary event censored CDF for a piecewise-constant (step) delay distribution and an arbitrary primary event distribution whose CDF F_{primary} is available via object$pprimary.

Usage

## S3 method for class 'pcens_pdiscretestep'
pcens_cdf(object, q, pwindow, use_numeric = FALSE)

Arguments

Details

The observation CDF is

F_{obs}(q) = \int_0^{pwindow} F_{step}(q-p)\,dF_{primary}(p)

Because F_{step} is piecewise constant, the integral reduces to

F_{obs}(q) = \sum_k c_k \,[F_{primary}(p^{end}_k) - F_{primary}(p^{start}_k)]

where c_k is the constant value of F_{step} on the k-th sub-interval of the primary event window induced by the step-function knots. The partition is exact for any bin widths, so bins may be wider or narrower than pwindow, and for boundaries that start below zero.

Falls back to pcens_cdf.default when use_numeric = TRUE or when no primary CDF is available on the object.

Value

Vector of computed primary event censored CDFs

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile(), pcens_quantile.default()

Method for Gamma delay with uniform primary

Description

Method for Gamma delay with uniform primary

Usage

## S3 method for class 'pcens_pgamma_dunif'
pcens_cdf(object, q, pwindow, use_numeric = FALSE)

Arguments

Value

Vector of computed primary event censored CDFs

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile(), pcens_quantile.default()

Method for generalised gamma delay with uniform primary

Description

Analytical solution for the generalised gamma distribution in the Stacy parameterisation used by flexsurv::pgengamma.orig(), with parameters shape, scale and k. The delay CDF is F_T(t) = P(k, (t / \theta)^a) with P the regularised lower incomplete gamma function, a the shape and \theta the scale. The mean is E[T] = \theta \Gamma(k + 1/a) / \Gamma(k) and the partial expectation distribution is \tilde F_T(t) = P(k + 1/a, (t / \theta)^a), so the solution generalises the gamma (shape = 1) and Weibull (k = 1) cases. See vignette("analytic-solutions") for the derivation.

Usage

## S3 method for class 'pcens_pgengamma.orig_dunif'
pcens_cdf(object, q, pwindow, use_numeric = FALSE)

Arguments

Value

Vector of computed primary event censored CDFs

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile(), pcens_quantile.default()

Examples

pcens_obj <- new_pcens(
  pdist = flexsurv::pgengamma.orig,
  dprimary = dunif,
  dprimary_args = list(min = 0, max = 1),
  shape = 1.5,
  scale = 2,
  k = 0.8
)
pcens_cdf(pcens_obj, q = c(1, 4, 8), pwindow = 1)

Method for generalised gamma (Prentice parameterisation) delay with uniform primary

Description

Analytical solution for the generalised gamma distribution in the Prentice parameterisation used by flexsurv::pgengamma(), with parameters mu, sigma and Q. For Q > 0 this is mapped to the Stacy parameterisation of pcens_cdf.pcens_pgengamma.orig_dunif() via shape = Q / sigma, scale = exp(mu) * Q^(2 * sigma / Q) and k = 1 / Q^2. For Q <= 0 (the lognormal and reflected cases) the numerical pcens_cdf.default() method is used.

Usage

## S3 method for class 'pcens_pgengamma_dunif'
pcens_cdf(object, q, pwindow, use_numeric = FALSE)

Arguments

Value

Vector of computed primary event censored CDFs

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile(), pcens_quantile.default()

Method for Log-Normal delay with uniform primary

Description

Method for Log-Normal delay with uniform primary

Usage

## S3 method for class 'pcens_plnorm_dunif'
pcens_cdf(object, q, pwindow, use_numeric = FALSE)

Arguments

Value

Vector of computed primary event censored CDFs

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile(), pcens_quantile.default()

Method for Weibull delay with uniform primary

Description

Method for Weibull delay with uniform primary

Usage

## S3 method for class 'pcens_pweibull_dunif'
pcens_cdf(object, q, pwindow, use_numeric = FALSE)

Arguments

Value

Vector of computed primary event censored CDFs

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_quantile(), pcens_quantile.default()

Compute primary event censored quantiles

Description

This function inverts the primary event censored CDF to compute quantiles. It uses numerical optimisation via optim to find the value q such that pcens_cdf() is close to the specified probability. Currently, only the default numerical inversion method is implemented. Future analytical solutions may be added.

Usage

pcens_quantile(object, p, pwindow, L = -Inf, D = Inf, use_numeric = FALSE, ...)

Arguments

Value

Vector of primary event censored quantiles.

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile.default()

Default method for computing primary event censored quantiles

Description

This method inverts the primary event censored CDF by root-finding via stats::uniroot() with extendInt = "upX". The censored CDF is monotone in q, so stats::uniroot() extends its starting bracket outward as needed and handles infinite L or D without special casing.

Usage

## Default S3 method:
pcens_quantile(
  object,
  p,
  pwindow,
  L = -Inf,
  D = Inf,
  use_numeric = FALSE,
  init = 5,
  tol = 1e-08,
  max_iter = 10000,
  ...
)

Arguments

Value

A numeric vector containing the computed primary event censored quantiles.

See Also

Low level primary event censored distribution objects and methods new_pcens(), pcens_cdf(), pcens_cdf.default(), pcens_cdf.pcens_pdiscretehazard(), pcens_cdf.pcens_pdiscretestep(), pcens_cdf.pcens_pgamma_dunif(), pcens_cdf.pcens_pgengamma.orig_dunif(), pcens_cdf.pcens_pgengamma_dunif(), pcens_cdf.pcens_plnorm_dunif(), pcens_cdf.pcens_pweibull_dunif(), pcens_quantile()

Examples

# Create a primarycensored object with gamma delay and uniform primary
pcens_obj <- new_pcens(
  pdist = pgamma,
  dprimary = dunif,
  primary_args = list(min = 0, max = 1),
  shape = 3,
  scale = 2
)

# Compute quantile for a single probability
pcens_quantile(pcens_obj, p = 0.8, pwindow = 1)

# Compute quantiles for multiple probabilities
pcens_quantile(pcens_obj, p = c(0.25, 0.5, 0.75), pwindow = 1)

# Compute quantiles for multiple probabilities with truncation
pcens_quantile(pcens_obj, p = c(0.25, 0.5, 0.75), pwindow = 1, D = 10)

# Compute quantiles with left truncation
pcens_quantile(pcens_obj, p = c(0.25, 0.5, 0.75), pwindow = 1, L = 1, D = 10)

Hazard-parameterised piecewise-constant CDF

Description

Returns the CDF of a discrete distribution specified by its bin-wise discrete-time hazards. Converts hazards to a PMF via hazards_to_pmf() then delegates to pdiscretestep().

Usage

pdiscretehazard(q, boundaries = NULL, hazards)

Arguments

Details

Outside fitting it is a deterministic wrapper around pdiscretestep(). It earns its keep as a fitting parameterisation in fitdistdoublecens() because the random walk on the logit hazard smooths the recovered PMF.

Value

Numeric vector of CDF values, the same length as q.

See Also

pdiscretestep(), hazards_to_pmf(), pmf_to_hazards(), fitdistdoublecens()

Other pdiscretehazard: ddiscretehazard(), discretehazard_start(), rdiscretehazard()

Examples

hazards <- c(0.3, 0.5, 1)
pdiscretehazard(c(0.5, 1, 2, 3), boundaries = 0:3, hazards = hazards)

Step (piecewise-constant) CDF

Description

Returns the CDF of a discrete distribution whose mass is concentrated at the right boundary of each bin. The CDF is right-continuous and piecewise constant with jumps at boundaries[2], ..., boundaries[K+1].

Usage

pdiscretestep(q, boundaries = NULL, pmf)

Arguments

Details

Below boundaries[1] the function returns 0. At boundaries[i+1] (the right edge of bin i), F jumps by pmf[i], so F(boundaries[i+1]) = \sum_{j=1}^{i} pmf_j. For q in [boundaries[i], boundaries[i+1]), F equals \sum_{j=1}^{i-1} pmf_j. At or above boundaries[K+1] the function returns 1.

Use with fitdistdoublecens()

This function carries the attribute vector_param = "pmf" so that fitdistdoublecens can drive it from a flat list of scalar parameters p1, ..., p_{K-1}. The free parameters are the first K-1 bin probabilities; the last is set to 1 - sum(p1, ..., p_{K-1}). When the implied probabilities violate the simplex (any negative entry, or sum departing from 1 by more than 10^{-8}), the function returns 0 (or near-zero density in ddiscretestep) rather than erroring; this drives the optimiser back to the feasible region.

Value

Numeric vector of CDF values, the same length as q.

See Also

Other pdiscretestep: ddiscretestep(), hazards_to_pmf(), pmf_to_hazards(), rdiscretestep()

Examples

# Two-bin PMF: mass 0.3 at x=1, mass 0.7 at x=2
pdiscretestep(c(0.5, 1, 1.5, 2), boundaries = 0:2, pmf = c(0.3, 0.7))

Convert a PMF to discrete-time hazards

Description

Inverts hazards_to_pmf(): given a PMF, computes the discrete-time conditional hazard at each time point.

h_i = pmf_i / (1 - \sum_{j < i} pmf_j)

Usage

pmf_to_hazards(pmf)

Arguments

Details

The returned vector has the same length as pmf, with the last entry equal to 1.

Value

Numeric vector of hazards in [0, 1].

See Also

Other pdiscretestep: ddiscretestep(), hazards_to_pmf(), pdiscretestep(), rdiscretestep()

Examples

pmf_to_hazards(c(0.2, 0.3, 0.5))

Compute the primary event censored CDF for delays

Description

This function computes the primary event censored cumulative distribution function (CDF) for a given set of quantiles. It adjusts the CDF of the primary event distribution by accounting for the delay distribution and potential truncation at a maximum delay (D) and minimum delay (L). The function allows for custom primary event distributions and delay distributions.

Usage

pprimarycensored(
  q,
  pdist,
  pwindow = 1,
  L = -Inf,
  D = Inf,
  dprimary = dunif,
  primary_args = NULL,
  pprimary = NULL,
  dprimary_args = NULL,
  ...,
  check = TRUE
)

ppcens(
  q,
  pdist,
  pwindow = 1,
  L = -Inf,
  D = Inf,
  dprimary = dunif,
  primary_args = NULL,
  pprimary = NULL,
  dprimary_args = NULL,
  ...,
  check = TRUE
)

Arguments

Details

The primary event censored CDF is computed by integrating the product of the delay distribution function (CDF) and the primary event distribution function (PDF) over the primary event window. The integration is adjusted for truncation if specified.

The primary event censored CDF, F_{\text{cens}}(q), is given by:

F_{\text{cens}}(q) = \int_{0}^{pwindow} F(q - p) \cdot f_{\text{primary}}(p) \, dp

where F is the CDF of the delay distribution, f_{\text{primary}} is the PDF of the primary event times, and pwindow is the primary event window.

If truncation is applied (finite D or finite L), the CDF is normalized:

F_{\text{cens,norm}}(q) = \frac{F_{\text{cens}}(q) - F_{\text{cens}}(L)}{ F_{\text{cens}}(D) - F_{\text{cens}}(L)}

where F_{\text{cens,norm}}(q) is the normalized CDF. For values q \leq L, the function returns 0; for values q \geq D, it returns 1.

This function creates a primarycensored object using new_pcens() and then computes the primary event censored CDF using pcens_cdf(). This abstraction allows for automatic use of analytical solutions when available, while seamlessly falling back to numerical integration when necessary.

See methods(pcens_cdf) for which combinations have analytical solutions implemented.

Value

Vector of primary event censored CDFs, normalized over [L, D] if truncation is applied

See Also

new_pcens() and pcens_cdf()

Primary event censored distribution functions dprimarycensored(), qprimarycensored(), rprimarycensored()

Examples

# Example: Lognormal distribution with uniform primary events
pprimarycensored(c(0.1, 0.5, 1), plnorm, meanlog = 0, sdlog = 1)

# Example: Lognormal distribution with exponential growth primary events
pprimarycensored(
  c(0.1, 0.5, 1), plnorm,
  dprimary = dexpgrowth,
  primary_args = list(r = 0.2), meanlog = 0, sdlog = 1
)

# Example: Left-truncated distribution (e.g., for generation intervals)
pprimarycensored(
  c(1, 2, 3), plnorm,
  L = 1, D = 10,
  meanlog = 0, sdlog = 1
)

Compute quantiles corresponding to target probabilities for primary event censored delays

Description

This function computes the quantiles (delay values) that correspond to specified probabilities in the primary event censored distribution. For a given probability p, it computes the delay value q such that the cumulative probability up to q equals p in the primary event censored distribution. The distribution accounts for both the delay distribution and the primary event timing distribution.

Usage

qprimarycensored(
  p,
  pdist,
  pwindow = 1,
  L = -Inf,
  D = Inf,
  dprimary = dunif,
  primary_args = NULL,
  pprimary = NULL,
  dprimary_args = NULL,
  ...,
  check = TRUE
)

qpcens(
  p,
  pdist,
  pwindow = 1,
  L = -Inf,
  D = Inf,
  dprimary = dunif,
  primary_args = NULL,
  pprimary = NULL,
  dprimary_args = NULL,
  ...,
  check = TRUE
)

Arguments

Details

For each probability, the function finds the delay value where that proportion of events have occurred by that time in the primary event censored distribution. This is done by inverting the cumulative distribution function.

The function creates a primarycensored object using new_pcens() and then computes the quantiles using pcens_quantile(). This approach allows for analytical solutions when available, falling back to numerical methods when necessary.

For example, if p = 0.5, the function returns the median delay (truncated over [L, D] if specified) where 50% of censored events occur by this time and 50% occur after.

See methods(pcens_quantile) for which combinations have analytical solutions implemented.

Value

Vector of delay values (quantiles) corresponding to the input probabilities

See Also

new_pcens() and pcens_quantile()

Primary event censored distribution functions dprimarycensored(), pprimarycensored(), rprimarycensored()

Examples

# Compute delays where 25%, 50%, and 75% of events occur by (quartiles)
# Using lognormal delays with uniform primary events
qprimarycensored(c(0.25, 0.5, 0.75), plnorm, meanlog = 0, sdlog = 1)

# Same quartiles but with exponential growth in primary events
qprimarycensored(
  c(0.25, 0.5, 0.75), plnorm,
  dprimary = dexpgrowth,
  primary_args = list(r = 0.2), meanlog = 0, sdlog = 1
)

# Same quartiles but with truncation at 10
qprimarycensored(
  c(0.25, 0.5, 0.75), plnorm,
  dprimary = dexpgrowth,
  primary_args = list(r = 0.2), meanlog = 0, sdlog = 1, D = 10
)

# Left-truncated distribution (e.g., for generation intervals)
qprimarycensored(
  c(0.25, 0.5, 0.75), plnorm,
  L = 1, D = 10, meanlog = 0, sdlog = 1
)

Sample from a hazard-parameterised step distribution

Description

Draws n independent samples from the discrete distribution defined by hazards. Converts hazards to a PMF via hazards_to_pmf() then delegates to rdiscretestep().

Usage

rdiscretehazard(n, boundaries = NULL, hazards)

Arguments

Details

Outside fitting it is a deterministic wrapper around rdiscretestep(). It earns its keep as a fitting parameterisation in fitdistdoublecens() because the random walk on the logit hazard smooths the recovered PMF.

Value

Numeric vector of length n.

See Also

rdiscretestep(), hazards_to_pmf(), pmf_to_hazards(), fitdistdoublecens()

Other pdiscretehazard: ddiscretehazard(), discretehazard_start(), pdiscretehazard()

Examples

set.seed(42)
rdiscretehazard(10, boundaries = 0:3, hazards = c(0.3, 0.5, 1))

Sample from a step distribution

Description

Draws n independent samples from the discrete distribution with probability mass pmf[i] at boundaries[i+1].

Usage

rdiscretestep(n, boundaries = NULL, pmf)

Arguments

Value

Numeric vector of length n.

See Also

Other pdiscretestep: ddiscretestep(), hazards_to_pmf(), pdiscretestep(), pmf_to_hazards()

Examples

set.seed(42)
rdiscretestep(10, boundaries = 0:3, pmf = c(0.2, 0.5, 0.3))

Generate random samples from a primary event censored distribution

Description

This function generates random samples from a primary event censored distribution. It adjusts the distribution by accounting for the primary event distribution and potential truncation at a maximum delay (D) and minimum delay (L). The function allows for custom primary event distributions and delay distributions.

Usage

rprimarycensored(
  n,
  rdist,
  pwindow = 1,
  swindow = 1,
  L = -Inf,
  D = Inf,
  rprimary = stats::runif,
  rprimary_args = list(),
  oversampling_factor = 1.2,
  ...
)

rpcens(
  n,
  rdist,
  pwindow = 1,
  swindow = 1,
  L = -Inf,
  D = Inf,
  rprimary = stats::runif,
  rprimary_args = list(),
  oversampling_factor = 1.2,
  ...
)

Arguments

Details

The mathematical formulation for generating random samples from a primary event censored distribution is as follows:

1. Generate primary event times (p) from the specified primary event distribution (f_p) with parameters phi, defined between 0 and the primary event window (pwindow):

   p \sim f_p(\phi), \quad p \in [0, pwindow]

2. Generate delays (d) from the specified delay distribution (f_d) with parameters theta:

   d \sim f_d(\theta)

3. Calculate the total delays (t) by adding the primary event times and the delays:

   t = p + d

4. Apply upper truncation to remove delays >= D:

   t_{upper} = \{t \mid t < D\}

5. Round the delays to the nearest secondary event window (swindow):

   t_{rounded} = \lfloor \frac{t_{upper}}{swindow} \rfloor \times swindow

6. Apply lower truncation on the rounded values to ensure observed delays are >= L:

   t_{valid} = \{t_{rounded} \mid t_{rounded} \geq L\}

The function oversamples to account for potential truncation and generates additional samples if needed to reach the desired number of valid samples.

Value

Vector of random samples from the primary event censored distribution censored by the secondary event window.

See Also

Primary event censored distribution functions dprimarycensored(), pprimarycensored(), qprimarycensored()

Examples

# Example: Lognormal distribution with uniform primary events
rprimarycensored(10, rlnorm, meanlog = 0, sdlog = 1)

# Example: Lognormal distribution with exponential growth primary events
rprimarycensored(
  10, rlnorm,
  rprimary = rexpgrowth, rprimary_args = list(r = 0.2),
  meanlog = 0, sdlog = 1
)

# Example: Left-truncated distribution (e.g., for generation intervals)
rprimarycensored(10, rlnorm, L = 1, D = 10, meanlog = 0, sdlog = 1)