Package {BayesPIM}

BayesPIM: Bayesian Prevalence-Incidence Mixture Model

Description

Models time-to-event data from interval-censored screening studies. It accounts for latent prevalence at baseline and incorporates misclassification due to imperfect test sensitivity. For usage details, see the package vignette "BayesPIM_intro". Further details can be found in Klausch, Lissenberg-Witte and Coupé (2026) doi:10.1002/sim.70433.

Author(s)

Maintainer: Thomas Klausch t.klausch@amsterdamumc.nl

Authors:

• Thomas Klausch t.klausch@amsterdamumc.nl

See Also

Useful links:

• https://github.com/thomasklausch2/BayesPIM

• Report bugs at https://github.com/thomasklausch2/BayesPIM/issues

Fitting Bayesian Prevalence-Incidence Mixture Model

Description

Estimates the prevalence-incidence mixture model of Klausch et al. (2026) using a Bayesian Gibbs sampler. The model is formulated as an interval-censored survival model over successive intervals, with the possibility of missed events due to imperfect test sensitivity. In addition, baseline tests at time zero may fail to detect pre-study events (prevalence).

Usage

bayespim(
  v_obs,
  x_t = NULL,
  x_g = NULL,
  r = NULL,
  dist = "weibull",
  kappa = NULL,
  update_kappa = FALSE,
  kappa_prior = NULL,
  ndraws = 1000,
  warmup = floor(ndraws/2),
  warmup_updated = FALSE,
  prop_sd = NULL,
  slice_width = 1,
  chains,
  seed_chains = NULL,
  save_every = 1,
  update_till_converge = FALSE,
  maxit = Inf,
  max_rhat = 1.01,
  min_effss = chains * 100,
  sampler = "slice_collapsed",
  log_prior_fun = log_aft_prior,
  beta_prior = "norm",
  tau_t = 1,
  sig_prior = 1,
  tau_g = 1,
  fix_sigma = FALSE,
  q_prior_sd = 1,
  fix_q = FALSE,
  prev_run = NULL,
  ndraws_update = NULL,
  prev = TRUE,
  par_exp = FALSE,
  rescale_times = TRUE,
  standardize_covariates = TRUE,
  ini_spread = 0.5,
  silent = FALSE
)

Arguments

Details

This Bayesian prevalence-incidence mixture model (PIM) characterizes time to incidence through accelerated failure time (AFT) scaling. With \eta_i = \bm{x}_{ti}'\bm{\beta}_t, covariates multiply event times by \exp(\eta_i) relative to the corresponding baseline distribution. For the Weibull, lognormal, and log-logistic families this is the familiar log-location-scale representation

\log(t_i) = \eta_i + \sigma_t\epsilon_i.

For dist = "gamma", BayesPIM instead uses

t_i \mid \bm{x}_{ti} \sim \mathrm{Gamma}\{\sigma_t^{-2}, \sigma_t^{-2}\exp(-\eta_i)\},

so E(t_i \mid \bm{x}_{ti}) = \exp(\eta_i) and \sigma_t is the conditional coefficient of variation. For dist = "gengamma", BayesPIM uses the Prentice generalized gamma with location \eta_i, scale \sigma_t, and signed shape Q. Generalized-gamma fitting is available only with sampler = "slice_collapsed". The covariate matrix for the latent-time model is supplied through x_t.

Baseline prevalence is modeled using a probit formulation Pr(g_i=1 | \bm{x}_{gi}) = Pr(w_i > 0 | \bm{x}_{gi}) with

w_i = \bm{x}_{gi}' \bm{\beta}_g + \psi_i

where \psi_i follows a standard normal distribution, and the covariate vector \bm{x}_{gi} is given in the x_g matrix. The latent continuous probit variable w_i determines the modeled binary prevalence variable: g_i = 1 if w_i > 0 and g_i = 0 otherwise.

The argument v_obs provides the observed testing times for all individuals. It is a list of numeric vectors, where each vector starts with 0 (representing the baseline time) and is followed by one or more screening times. The final entry is Inf in the case of right censoring or indicates the time of a positive test if an event is observed. Specifically:

• If the baseline test is positive, the vector consists solely of c(0).

• If the baseline test is negative and right censoring occurs before the first regular screening, the vector is c(0, Inf).

• Otherwise, the vector ends with Inf in the case of right censoring (e.g., c(0, 1, 3, 6, Inf)) or ends at the event time (e.g., c(0, 1, 3, 6) for an event detected at time 6).

By convention, every vector in v_obs starts with 0. However, the binary vector r of length n indicates whether the baseline test was conducted (r[i] = 1) or missing (r[i] = 0) for each individual i in v_obs. For further details on coding, see Section 2 of the main paper.

Test sensitivity can be fixed to a value kappa by setting update_kappa = FALSE, or it can be estimated if update_kappa = TRUE. When estimated with kappa_prior = c(m, s), a Beta prior with mean m and standard deviation s is used. Its shape parameters are calculated analytically as

\nu = \frac{m(1-m)}{s^2} - 1,\qquad a = m\nu,\qquad b = (1-m)\nu.

The requested mean must lie strictly between zero and one, and the standard deviation must satisfy 0 < s < \sqrt{m(1-m)}. Malformed or infeasible specifications stop with an informative error. The calculated shape parameters must be finite and positive. If kappa_prior = NULL, the function warns and uses an uninformative \mathrm{Beta}(1,1) prior. In general, we advise against using an uninformative prior, but this default avoids favoring any specific informative prior.

The Gibbs sampler runs for ndraws iterations for each of chains total chains. Incidence parameters can be updated by Metropolis-Hastings, by univariate slice sampling conditional on augmented exact incidence times, or by collapsed univariate slice sampling conditional only on sampled screening intervals. The collapsed slice sampler (sampler = "slice_collapsed") is the default because it has shown lower autocorrelation and faster convergence than the two exact-time samplers. This differs from the original implementation described by Klausch et al., which used a Metropolis-Hastings sampler. The collapsed sampler also supports the generalized-gamma model, for which convergence with the other samplers is typically slow; consequently, dist = "gengamma" is available only with the collapsed sampler.

When sampler = "mh", the Metropolis step applies a normal proposal distribution with standard deviation prop_sd, which must be selected by trial and error. An optimal acceptance rate is approximately 23%, which can be computed per MCMC run from the model output. The function search_prop_sd provides a heuristic for selecting an effective proposal standard deviation.

We recommend running at least two chains to facilitate standard MCMC diagnostics such as the Gelman-Rubin statistic. For larger analyses, users may run more chains, but should ensure that their computing environment permits the requested parallel workers. CRAN examples use small sequential runs. Additionally, we suggest first running the sampler for a moderate number of iterations to assess its behavior before using the updating functionality in prev_run to extend sampling (see below).

The option update_till_converge = TRUE allows bayespim to run until convergence. Convergence is achieved when the rank-normalized split R-hat calculated by posterior::rhat() is at most max_rhat for every sampled parameter and the minimum effective sample size calculated by posterior::ess_mean() reaches min_effss. The diagnostics use every stored post-warm-up draw; no additional diagnostic thinning is applied. The current diagnostics are always calculated, stored, and printed. If automated updating is enabled, the sampler continues updating and printing diagnostics until convergence is attained or maxit is reached.

Setting silent = TRUE suppresses this printing for both forms of updating, the automated loop just described and a manual update requested through prev_run. Only the progress output is withheld: the diagnostics are still calculated and stored in the convergence element of the returned object, and warnings and errors are still signaled.

The priors for the regression coefficients and distributional parameters can be controlled using beta_prior, tau_t, sig_prior, tau_g, and q_prior_sd. Specifically:

• beta_prior determines the prior type for \beta_{tj} (either normal or Student-t t).

• tau_t specifies either the standard deviation (for normal priors) or degrees of freedom (for Student-t priors). The default is a standard normal prior.

• A half-normal prior is used for the positive family scale/dispersion parameter \sigma, with sig_prior controlling the standard deviation. In the gamma model, \sigma is the conditional coefficient of variation.

• A zero-centered normal prior is assigned to \beta_{gj}, with tau_g controlling its standard deviation (default: standard normal).

• In the generalized-gamma model, a zero-centered normal prior is assigned to the signed shape parameter Q, with standard deviation q_prior_sd.

Sometimes model fitting can be improved by fixing the \sigma parameter to a value, which is achieved through setting fix_sigma = TRUE. Then, the value specified as sig_prior is regarded as the fixed value of \sigma. For the gamma model this fixes the conditional coefficient of variation. The functionality can also be used to obtain the exponential distribution, akin to a Markov model. For this choose dist = "weibull", sig_prior = 1, and fix_sigma = TRUE.

The prev_run argument allows updating a previous run with additional MCMC draws. The MCMC chain resumes from the last draws, continues, and merges with the original run. If an initial model was fit using mod <- bayespim(...), it can be updated using mod_update <- bayespim(prev_run = mod). By default, ndraws additional iterations are added unless otherwise specified via ndraws_update. If warmup_updated = TRUE, the warmup increases with each update; otherwise the initial warmup is retained.

The Gibbs sampler requires starting values. Incidence-model coefficients are anchored at estimates from an interval-censored survival model, while prevalence-model coefficients are anchored at a penalized probit fit. With multiple chains, the two coefficient vectors are multiplied by equally spaced factors from ini_spread to 1; a single chain uses the unscaled fitted values. Each complete candidate state, including the latent prevalence indicators and incidence intervals or times, must have a finite posterior. Each chain makes up to ten initialization attempts. The first eight redraw the latent states at that chain's own spread factor; the ninth uses the fitted estimates unscaled; and the tenth uses zero-centered coefficients. If none of the ten yields a finite posterior, initialization stops with an informative error.

Value

A list containing the following elements:

The returned list has class "bayespim". Additionally, most input arguments are returned as part of the output for reference.

References

T. Klausch, B. I. Lissenberg-Witte, and V. M. H. Coupé (2026). "A Bayesian prevalence-incidence mixture model for screening outcomes with misclassification.", Statistics in Medicine, 45(8-9), e70433. doi:10.1002/sim.70433

J. S. Liu and Y. N. Wu, “Parameter Expansion for Data Augmentation,” Journal of the American Statistical Association, vol. 94, no. 448, pp. 1264–1274, 1999, https://doi.org/10.2307/2669940.

See Also

summary.bayespim() and plot.bayespim()

Examples

# A deliberately short fit for illustrating the interface.
# Use substantially more draws for scientific inference.
set.seed(2025)
dat <- gen_data(kappa = 0.7, n = 100, theta = 0.2,
               p = 1, p_discrete = 1,
               beta_t = c(0.2, 0.2), beta_g = c(0.2, 0.2),
               v_min = 20, v_max = 30, mean_rc = 80,
               sigma_t = 0.2, mu_t = 5, dist = "weibull",
               prob_r = 1)

fit <- bayespim(
  v_obs = dat$v_obs,
  x_t = dat$x,
  x_g = dat$x,
  r = dat$r,
  kappa = 0.7,
  update_kappa = FALSE,
  ndraws = 100,
  warmup = 10,
  chains = 2,
  seed_chains = c(202501, 202502),
  update_till_converge = FALSE,
  sampler = "slice_collapsed",
  dist = "weibull",
  silent = TRUE
)
fit$runtime

gen_data: Simulate Screening Data for a Prevalence-Incidence Mixture Model

Description

Generates synthetic data according to the Bayesian prevalence-incidence mixture (PIM) framework of Klausch et al. (2026) with interval-censored screening outcomes. The function simulates continuous or discrete baseline covariates, event times from one of several parametric families, and irregular screening schedules, yielding interval-censored observations suitable for testing or demonstrating PIM-based or other interval-censored survival methods.

Usage

gen_data(
  kappa = 0.7,
  n = 1000,
  p = 2,
  p_discrete = 0,
  rho = 0,
  s = 1,
  sigma_t = 1/2,
  mu_t = 4,
  beta_t = NULL,
  beta_g = NULL,
  theta = 0.15,
  v_min = 1,
  v_max = 6,
  mean_rc = 40,
  dist = "weibull",
  q = 1,
  sel_mod = "probit",
  prob_r = 0
)

Arguments

Details

The data-generating process includes:

1. Covariates X: Continuous multivariate normal distributed covariates are simulated using a correlation structure specified by rho and a common standard deviation s. If p_discrete = 1, a single discrete covariate is added, drawn from \mathrm{Bernoulli}(0.5).

2. Latent Times t: An accelerated failure time (AFT) model is used, with linear predictor

   \eta_i = \beta_{t0} + \beta_t^\top x_{ti},

   where \beta_{t0} is set by mu_t. Covariates multiply event times by \exp(\eta_i) relative to the corresponding baseline family.

       For \code{"weibull"}, \code{"lognormal"}, and \code{"loglog"}
       (log-logistic), latent times are generated from
       \deqn{\log(t_i) = \eta_i + \sigma_t \epsilon_i.}{
                 log(t_i) = eta_i + sigma_t * epsilon_i.}
       For \code{"gamma"},
       \deqn{t_i \mid x_{ti} \sim
             \mathrm{Gamma}{\sigma_t^{-2},
             \sigma_t^{-2}\exp(-\eta_i)},}{
             t_i | x_ti ~ Gamma(shape = sigma_t^-2,
             rate = sigma_t^-2 * exp(-eta_i)),}
       so that \eqn{E(t_i \mid x_{ti}) = \exp(\eta_i)} and
       \eqn{\sigma_t} is the conditional coefficient of variation.
       For \code{"gengamma"}, the Prentice generalized gamma is generated
       with location \eqn{\eta_i}, scale \eqn{\sigma_t}, and signed shape
       parameter \code{q}.
   

3. Irregular Screening Schedules V_i: Each individual has multiple screening times generated randomly between v_min and v_max, ending in right censoring or the time of detection. These screening times (including a 0 for baseline and Inf for censoring) are returned in v_obs.

4. Prevalence Indicator g_i: Baseline prevalence is modeled via either a probit or logit link, consistent with:

   w_i = \beta_{g0} + \beta_g^\top x_{gi} + \psi_i,

   where \beta_{g0} is determined by theta, and \beta_g by beta_g. Specifically:

   • If sel_mod = "probit", then \beta_{g0} = \mathrm{qnorm}(\theta).

   • If sel_mod = "logit", then \beta_{g0} = \log(\theta / (1-\theta)).

   We set g_i = 1 if w_i > 0, and g_i = 0 otherwise.

5. Baseline Test Missingness r_i: A baseline test indicator r_i \in \{0,1\} is drawn from a Bernoulli distribution with success probability prob_r, so r_i = 1 means the baseline test is performed and r_i = 0 means it is missing.

6. Test Sensitivity \kappa: A misclassification parameter \kappa (test sensitivity) can be specified via kappa. If \kappa < 1, some truly positive cases are missed.

Value

A list with the following elements:

v_obs

A list of length n, each entry containing screening times. The first element is 0 (baseline), and Inf may indicate right censoring. The right-censoring time is drawn as the first post-baseline screening time plus an exponential increment, so it always exceeds that screening time and every individual receives at least one screening after baseline. Consequently the coding c(0, Inf), a negative or missing baseline test followed by right censoring before the first regular screening, is never generated here, even though bayespim accepts and models that case.

times_true

Numeric vector of length n giving the true latent times t_i.

x

Numeric matrix of dimension n \times p (plus an extra column if p_discrete = 1) containing the covariates.

g

Binary vector of length n, indicating whether an individual is truly positive at baseline (g_i = 1).

r

Binary vector of length n, indicating whether the baseline test was performed (r_i = 1) or missing (r_i = 0). This is the vector passed to the r argument of bayespim.

prob_g

Numeric vector of length n giving the true prevalence probabilities, P(g_i = 1).

References

T. Klausch, B. I. Lissenberg-Witte, and V. M. H. Coupé (2026). "A Bayesian prevalence-incidence mixture model for screening outcomes with misclassification.", Statistics in Medicine, 45(8-9), e70433. doi:10.1002/sim.70433

Examples

# Generate a small dataset for testing
set.seed(2025)
sim_data <- gen_data(n = 20, p = 1, p_discrete = 1,
                    sigma_t = 0.5, mu_t = 2,
                    beta_t = c(0.2, 0.2), beta_g = c(0.5, -0.2),
                    theta = 0.2,
                    dist = "weibull", sel_mod = "probit")
names(sim_data)

Compute Information Criteria for a Bayesian Prevalence-Incidence Mixture Model

Description

Computes and returns information criteria for a fitted Bayesian prevalence-incidence mixture model, including the Widely Applicable Information Criterion 1 (WAIC-1), WAIC-2, and the Deviance Information Criterion (DIC). These criteria are commonly used for model comparison and evaluation in Bayesian analysis. See Gelman et al. (2014) for further details on these criteria.

Usage

get_ic(mod, samples = NULL, cores = NULL)

Arguments

Details

This function calculates information criteria for a fitted Bayesian prevalence-incidence mixture model (bayespim). The information criteria include:

• WAIC-1: Based on the sum of posterior variances of log-likelihood contributions.

• WAIC-2: Similar to WAIC-1 but incorporates an alternative variance estimate.

• DIC: Measures model fit by penalizing complexity via the effective number of parameters.

The computation is performed by evaluating log-likelihood values for MCMC samples. By default, all MCMC samples after burn-in are used, though a subset can be specified via the samples argument. All incidence-time distributions supported by bayespim() (Weibull, log-logistic, log-normal, gamma, and generalized gamma) are handled. For a model fitted with covariate standardization, likelihood calculations reconstruct the standardized parameterization and operate on standardized covariates to reduce cancellation on large original scales; this is algebraically equivalent to using the returned original-scale coefficients.

Parallelization is available via the foreach package, utilizing multiple cores if cores is set accordingly. If cores = NULL, all available cores will be used.

Value

A matrix containing WAIC-1, WAIC-2, and DIC values for the model.

References

Gelman, A., Hwang, J., & Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Stat Comput, 24(6), 997–1016.

Examples

data(mod)
set.seed(2025)
get_ic(mod, samples = 20, cores = 1)

Log prior for accelerated failure time models

Description

Evaluates the log-prior density of the incidence-model parameters in an accelerated failure time (AFT) model. For the Weibull, lognormal, log-logistic, and gamma distributions, the parameter vector is

Usage

log_aft_prior(eta, dist, beta_prior, tau_t, sig_prior, q_prior_sd)

Arguments

Details

\eta = (\beta_0, \beta_1, \ldots, \beta_p, \log(\sigma)).

For the Prentice generalized-gamma distribution, the signed shape parameter Q is appended without transformation:

\eta = (\beta_0, \beta_1, \ldots, \beta_p, \log(\sigma), Q).

The same prior-function contract is used for every supported distribution and incidence sampler. A custom function supplied through bayespim() must accept the named arguments eta, dist, beta_prior, tau_t, sig_prior, and q_prior_sd, and must return one numeric log-density value. Arguments that are irrelevant to a custom prior may be ignored. When bayespim() is called with standardize_covariates = TRUE, eta contains the internally standardized incidence coefficients; returned posterior draws are transformed back to the original covariate scale only after sampling.

The default prior places either independent Student-t priors or independent zero-centered normal priors on the regression coefficients. A zero-centered half-normal prior is placed on the positive family scale/dispersion parameter \sigma. For the gamma model, \sigma is the conditional coefficient of variation; for the other families it is their AFT scale parameter. The log density includes the Jacobian for the transformation from \sigma to \log(\sigma). For the generalized-gamma model, Q additionally receives a zero-centered normal prior with standard deviation q_prior_sd.

Value

A single numeric value giving the log-prior density.

Examples

# Use the default prior but give the incidence intercept five times its
# default prior scale. The function retains the complete prior contract and
# therefore works for all supported AFT distributions.
log_aft_prior_relaxed_intercept <- function(
    eta, dist, beta_prior, tau_t, sig_prior, q_prior_sd
) {
  log_prior <- log_aft_prior(
    eta = eta,
    dist = dist,
    beta_prior = beta_prior,
    tau_t = tau_t,
    sig_prior = sig_prior,
    q_prior_sd = q_prior_sd
  )

  if (beta_prior == "t") {
    log_prior -
      dt(eta[1], df = tau_t, log = TRUE) +
      dt(eta[1] / 5, df = tau_t, log = TRUE) -
      log(5)
  } else {
    log_prior -
      dnorm(eta[1], sd = tau_t, log = TRUE) +
      dnorm(eta[1], sd = 5 * tau_t, log = TRUE)
  }
}

Converged example BayesPIM fit

Description

A converged Weibull prevalence-incidence mixture model fitted to simulated screening data. The object is provided so examples of post-estimation functions can run without refitting the model.

Usage

data(mod)

Format

A "bayespim" object containing four MCMC chains with 2,000 generated and stored draws each (save_every = 1). The first 1,000 draws per chain are marked as warm-up. Its continuous covariate was standardized internally and its binary covariate was left unchanged; stored coefficients are on the original covariate scale. The fit contains the simulated screening observations, incidence and prevalence covariates, fitted parameter chains, convergence diagnostics, and the settings needed by BayesPIM post-estimation methods.

Details

The data were generated with set.seed(2025) using the data-generating settings shown in the main bayespim() documentation: 1,000 individuals, Weibull incidence times, two covariates, prevalence probability 0.2 at zero-valued covariates, and fixed test sensitivity 0.7.

The collapsed slice sampler was run once with automatic convergence updating enabled, although no update was required. All monitored parameters met R-hat at most 1.01 and effective sample size at least 400 after the initial 2,000 draws per chain.

Source

Simulated and fitted with gen_data() and bayespim().

Examples

data(mod)
inherits(mod, "bayespim")
mod$convergence$converged

Plot method for bayespim objects

Description

Plot method for bayespim objects

Usage

## S3 method for class 'bayespim'
plot(x, warmup = x$warmup, thinning = 1L, ...)

Arguments

Value

Invisibly returns x. The method is called for its side effect of producing trace and density plots for the requested parameter blocks.

Examples

data(mod)
plot(mod, thinning = 20)

Plot posterior predictive cumulative incidence functions

Description

Plot the mixture CIF, the non-prevalent CIF, or both CIFs from an object returned by ppCIF.

Usage

## S3 method for class 'ppCIF'
plot(
  x,
  y = NULL,
  type = c("mixture", "nonprevalent", "both"),
  ci = TRUE,
  main = NULL,
  xlab = "Time",
  ylab = "Cumulative incidence",
  xlim = NULL,
  ylim = c(0, 1),
  col = "#0072B2",
  ci_col = NULL,
  lwd = 2,
  ...
)

Arguments

Details

The plotting contract is independent of the fitted incidence distribution and therefore also applies to generalized-gamma posterior predictions.

For ppd_type = "percentiles", time is shown on the x-axis and uncertainty is represented vertically. For ppd_type = "quantiles", the stored inverse-CDF representation is drawn with uncertainty in time. Pointwise regions can be omitted with ci = FALSE.

Value

Invisibly returns the supplied "ppCIF" object. The method is called for its base-graphics plotting side effect.

Examples

data(mod)
set.seed(2025)
cif <- ppCIF(
  mod,
  pst_samples = 50,
  quant = seq(0, 300, length.out = 51)
)
plot(cif, type = "both")

Posterior predictive cumulative incidence functions

Description

Compute posterior predictive cumulative incidence functions (CIFs) from a fitted bayespim model.

Usage

ppCIF(
  mod,
  fix_x_t = NULL,
  fix_x_g = NULL,
  pst_samples = 1000,
  perc = seq(0, 1, 0.01),
  ppd_type = c("percentiles", "quantiles"),
  quant = NULL
)

Arguments

Details

A prevalence-incidence mixture model defines two related CIFs. One ppCIF() call always computes and stores both; the type argument of plot.ppCIF() controls which one is displayed. The mixture CIF includes prevalent cases as a point mass at time zero. The nonprevalent CIF describes event times among individuals who are non-prevalent at baseline. Both are computed from the same posterior draws and posterior predictive replicates and retained in the returned object. The candidate draws comprise every stored post-warm-up draw; no additional thinning is applied before the explicit random subsample of pst_samples draws. Fixed covariate values are supplied on their original scale. For a model fitted with covariate standardization, linear predictors are evaluated on the standardized computational scale to reduce numerical cancellation.

For a model fitted with prev = FALSE, the two CIFs are identical. Posterior prediction supports every incidence distribution available in bayespim(), including the Prentice generalized gamma. For dist = "gengamma", event times are drawn using the fitted location \mu = \mathbf{x}_t^\top\boldsymbol{\beta}_t, positive scale \sigma_t, and signed shape Q.

Value

An object of class "ppCIF" with:

mixture, nonprevalent

Lists containing med_cdf, the pointwise posterior predictive median, and med_cdf_ci, a two-row matrix with the pointwise 2.5 and 97.5 percent posterior predictive quantiles.

ppd_type

The selected representation.

quant, perc

The applicable time or probability grid; the unused grid is NULL.

pst_samples

The number of posterior draws used.

distribution

The fitted incidence distribution.

prevalence_model

Whether prevalence was fitted.

call

The matched call.

See Also

plot.ppCIF, bayespim

Examples

data(mod)
set.seed(2025)
cif <- ppCIF(
  mod,
  pst_samples = 50,
  ppd_type = "percentiles",
  quant = seq(0, 300, length.out = 51)
)
plot(cif)
plot(cif, type = "nonprevalent")
plot(cif, type = "both")

Automated Heuristic Search of a Proposal Standard Deviation for bayespim

Description

When bayespim is fitted with sampler = "mh", it uses a Metropolis step for sampling the incidence-model parameters and requires a standard deviation for the normal proposal distribution. This function uses a heuristic algorithm to find a proposal standard deviation such that the Metropolis sampler accepts proposed draws at a rate within the user-defined interval (by default around 20–25%). The default sampler in bayespim() is "slice_collapsed", so the Metropolis sampler must be requested explicitly before using this function.

Usage

search_prop_sd(m, ndraws = 1000, succ_min = 3, acc_bounds = c(0.2, 0.25))

Arguments

Details

Starting from an initial bayespim model object m, the function attempts to calibrate the standard deviation of the proposal distribution. Specifically, it:

1. Runs an initial update of ndraws iterations and computes an acceptance rate.

2. If the acceptance rate lies within acc_bounds, the number of MCMC draws ndraws is doubled, and the process repeats.

3. Otherwise, the proposal standard deviation \sigma is adjusted based on whether the acceptance rate p is below the lower bound a or above the upper bound b of acc_bounds.

4. The formula for adjustment is:

   \sigma \leftarrow \sigma \times (1 - \frac{ (a-p)}{a}) \quad\text{if } p < a, \quad \sigma \leftarrow \sigma \times (1 + \frac{ (p-b)}{b}) \quad\text{if } p > b.

By default, if the acceptance rate falls within [0.2, 0.25], that \sigma is considered acceptable, and the process continues until succ_min consecutive successes (doubles) are achieved.

Value

A list with the following elements:

prop_sd

The final (adjusted) proposal standard deviation.

ac

The acceptance rate in the last iteration.

Examples

## Not run: 
# search_prop_sd() requires an MH fit. This deliberately short initial fit
# illustrates the interface; use more draws for an actual calibration.
set.seed(2025)
dat <- gen_data(
  kappa = 0.7, n = 100, theta = 0.2,
  p = 1, p_discrete = 1,
  beta_t = c(0.2, 0.2), beta_g = c(0.2, 0.2),
  v_min = 20, v_max = 30, mean_rc = 80,
  sigma_t = 0.2, mu_t = 5, dist = "weibull",
  prob_r = 1
)

fit_mh <- bayespim(
  v_obs = dat$v_obs, x_t = dat$x, x_g = dat$x, r = dat$r,
  kappa = 0.7, update_kappa = FALSE,
  ndraws = 100, warmup = 10, chains = 2,
  seed_chains = c(202501, 202502),
  update_till_converge = FALSE,
  sampler = "mh", prop_sd = 0.005, dist = "weibull",
  silent = TRUE
)

search_sd <- search_prop_sd(
  m = fit_mh,
  ndraws = 100,
  succ_min = 1
)
print(search_sd)

## End(Not run)

Summary method for bayespim objects

Description

Summary method for bayespim objects

Usage

## S3 method for class 'bayespim'
summary(object, warmup = object$warmup, probs = c(0.025, 0.5, 0.975), ...)

Arguments

Value

Invisibly returns a list of class "summary.bayespim" containing the fitted latent-time distribution and incidence sampler, posterior and convergence tables for the latent-time model and prevalence model, and kappa when estimated, together with the convergence criteria, formatted draw information, and the underlying draw counts. The same information is printed to the console.

Examples

data(mod)
summary(mod)

Trim and thin an mcmc.list

Description

Convenience function for trimming burn-in iterations and applying thinning to an object of class mcmc.list.

Usage

trim_mcmc(obj, burnin = 0L, end = NULL, thinning = 1L)

Arguments

Details

The function selects stored rows seq(burnin + 1, end, by = thinning) from each chain and reconstructs the result as an mcmc.list, preserving its original iteration numbering and multiplying its existing storage interval by thinning.

Value

An object of class mcmc.list containing the trimmed and thinned chains.

Examples

data(mod)
trimmed <- trim_mcmc(
  mod$par,
  burnin = mod$warmup,
  thinning = 20
)
nrow(as.matrix(trimmed[[1]]))