BayesPIM: user guide

Gibbs sampler

The function bayespim estimates the parameters of the model in equations @ref(eq:aft) and @ref(eq:probit) with a Gibbs sampler. The observed data for a non-prevalent individual is a series of negative screening tests followed either by a positive test or by right censoring. The exact event time is therefore never observed; it is only known to lie in one of the screening intervals. Moreover, because a test can miss a truly positive individual with probability \(1-\kappa\), even the interval containing the event is uncertain. Bayesian data augmentation handles this missing information by imputing plausible latent quantities and then updating the model parameters given these quantities.

Klausch et al. (2026) originally augmented the exact latent event time \(t_i\) for each individual and updated the incidence parameters with a Metropolis-Hastings step on the resulting complete-data likelihood. The current version of BayesPIM instead uses a collapsed slice sampler, which is the default (sampler = "slice_collapsed"). Rather than imputing an exact time, the collapsed sampler augments only the latent screening interval \((v_{i,k-1}, v_{ik}]\) in which the (possibly missed) event occurred, and updates the incidence parameters from the interval-censored likelihood, marginalizing over the exact position of \(t_i\) within the interval. Marginalizing over \(t_i\) in this way is known as collapsing (Liu, 1994). A collapsed Gibbs sampler preserves the target posterior but reduces the autocorrelation that exact-time augmentation induces, which yields faster convergence. This is why the collapsed sampler replaced the original Metropolis-Hastings sampler as the default; the earlier samplers remain available (Section @ref(sec-samplers)).

The main steps of the Gibbs sampler are, after initialization:

1. Augment the latent prevalence indicators \(g_i\) for individuals whose baseline status is unknown, collapsing over the latent time.
   
2. Augment the latent screening interval containing the event for each non-prevalent individual, collapsing over the exact time within the interval.
   
3. Draw the incidence parameters \((\boldsymbol{\beta}_t, \sigma_t)\) from the interval-censored complete-data posterior by univariate slice sampling.
   
4. Draw the prevalence parameters \(\boldsymbol{\beta}_g\) using the latent-normal representation of the probit model.
   
5. If test sensitivity is estimated, update \(\kappa\).
   

These steps are run repeatedly. After a warm-up period, the sampler yields auto-correlated draws from the target posterior distribution of the parameters. Conditional on the augmented intervals and prevalence indicators, the incidence parameters are drawn from the complete-data posterior

\[\begin{align} q(\boldsymbol{\beta}_t, \sigma_t \mid \cdot) \propto L(\boldsymbol{\beta}_t, \sigma_t \mid \cdot)\,\pi(\boldsymbol{\beta}_t, \sigma_t), \end{align}\]

where \(L\) is the interval-censored complete-data likelihood and \(\pi\) the prior density (Section @ref(sec-prior-assumptions)). For an individual whose event has been augmented into interval \((v_{k-1}, v_k]\), the likelihood contribution is the probability that the event falls in that interval, \(F_t(v_k) - F_t(v_{k-1})\), where \(F_t\) is the cumulative distribution function of the incidence distribution. The full derivation is given in Klausch et al. (2026).

bayespim input data structure

bayespim accepts four data inputs: v_obs, x_t, x_g, and r. These have to be available for all individuals in the data; missing information has to be handled outside of BayesPIM.

Argument v_obs provides the observed screening times as a list of length \(n\), one numeric vector per individual. Every vector starts with 0, the baseline time. The remaining entries encode the screening series and its outcome:

• If the baseline test is positive, the vector is 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 for right censoring (e.g. c(0, 22, 48, 77, Inf)) or ends at the time of the positive test if an event is detected (e.g. c(0, 23, 52, 78, 105, 132)).

Because a longer screening series carries no extra information under the model’s assumptions, only the times up to the event or to the last follow-up need to be supplied. Argument r is a binary vector of length \(n\) indicating whether the baseline test was actually carried out (r[i] = 1) or is missing (r[i] = 0). Together, v_obs and r encode which cases are known prevalent (i.e. positive test at baseline), incident, or right censored, and whether their baseline status is known.

Arguments x_t and x_g take the covariates for the incidence and prevalence models, respectively, each an \(n \times p\) matrix. The two sets of covariates may differ, although in practice the same covariates are often used for both. If no covariates are supplied, x_t = NULL and/or x_g = NULL give intercept-only models. Categorical variables have to be dummy-coded by the user; stats::model.matrix is useful for this.

By default, standardize_covariates = TRUE centers each non-binary covariate at its sample mean and divides it by its sample standard deviation before fitting. A column with exactly two observed values is detected as binary/dummy and is left unchanged. This standardization generally improves sampling when continuous covariates have very different scales. The returned chains and summaries are transformed back to the original covariate scale, and the original matrices remain in $x_t and $x_g. Set standardize_covariates = FALSE to fit directly on the supplied scale. Numerically coded categorical variables with more than two levels must still be dummy-coded; otherwise they are treated as continuous.

BayesPIM has a built-in data-generating function gen_data that simulates data under the model of Klausch et al. (2026); see ?gen_data for all options. We use it to illustrate the input structure and model estimation.

library(BayesPIM)

# Generate data under the PIM of Klausch et al. (2026)
set.seed(2025)
dat <- gen_data(
  kappa   = 0.7,          # Test sensitivity
  n       = 1e3,          # Sample size
  theta   = 0.2,          # Baseline prevalence when all covariates are zero
  p       = 1,            # Number of continuous covariates
  p_discrete = 1,         # Add one Bernoulli(0.5) covariate
  beta_t  = c(0.2, 0.2),  # True incidence slopes
  beta_g  = c(0.2, 0.2),  # True prevalence slopes
  mu_t    = 5,            # True incidence intercept
  sigma_t = 0.2,          # True incidence AFT scale
  dist    = "weibull",    # Incidence distribution
  v_min   = 20,           # Minimum time between screening moments
  v_max   = 30,           # Maximum time between screening moments
  mean_rc = 80,           # Mean time to right censoring (exponential)
  prob_r  = 1             # Probability that a baseline test is done
)

gen_data simulates the incidence times from the AFT model in equation @ref(eq:aft), with intercept mu_t and slopes beta_t, and generates prevalence from the probit model in equation @ref(eq:probit), with intercept qnorm(theta) (so theta is the prevalence probability when all covariates are zero) and slopes beta_g. A screening series is then superimposed: successive screening times are drawn uniformly between v_min and v_max apart, and the right-censoring time is exponential with mean mean_rc. At each screening time after the true event, a test detects the event with probability kappa. Besides weibull, the distributions lognormal, loglog (log-logistic), gamma, and gengamma (Prentice generalized gamma) are available.

The screening times are returned in dat$v_obs, in exactly the format bayespim expects. Looking at the first individuals illustrates the coding.

head(dat$v_obs)

[[1]]
[1] 0

[[2]]
[1]   0.00000  22.53705  51.89452  78.35944 104.52213 131.74624

[[3]]
[1]  0.00000 25.37066

[[4]]
[1] 0

[[5]]
[1]   0.00000  27.37410  53.60168  80.54386 108.43628       Inf

[[6]]
[1]  0.00000 22.31108 48.06398 77.10746      Inf

Individuals 1 and 4 have c(0) and are known prevalent (due to a positive baseline test). Individual 2 has an event detected at time 131.7, and individual 3 at time 25.4. Individuals 5 and 6 end in Inf and are right-censored without a detected event. In addition, gen_data returns the covariate matrix dat$x, the baseline-test indicator dat$r, and the latent quantities dat$times_true (true event times), dat$g (true prevalence status), and dat$prob_g (true prevalence probabilities), which are unobserved in practice and are not passed to bayespim.

Before running the sampler, it is useful to inspect the case mix.

v <- dat$v_obs
prevalent <- vapply(v, function(x) length(x) == 1L, logical(1))
censored  <- vapply(v, function(x) length(x) > 1L && is.infinite(x[length(x)]), logical(1))
c(prevalent = sum(prevalent),
  incident  = sum(!prevalent & !censored),
  censored  = sum(censored))

prevalent  incident  censored
      181       209       610

About 18% of individuals are prevalent and 21% have a detected incident event, which is ample for bayespim to run successfully. Care in model specification (number of covariates and choice of distribution) should be taken when advanced-state events are infrequent. BayesPIM is designed to cope with sparse events and usually converges if run long enough, but the precision of estimation may then be low, especially with many covariates.

Prior assumptions

Klausch et al. (2026) showed that weakly informative priors help to obtain stable estimates even when events are infrequent and the likelihood is nearly flat. bayespim uses such priors by default. For the incidence regression coefficients \(\boldsymbol{\beta}_t\), a normal prior (beta_prior = 'norm') is used by default, with standard deviation tau_t (default 1); a Student-\(t\) prior is available with beta_prior = 't', in which case tau_t sets the degrees of freedom. A half-normal prior with standard deviation sig_prior (default 1) is placed on the positive scale/dispersion parameter \(\sigma_t\). For the prevalence coefficients \(\boldsymbol{\beta}_g\), a zero-centered normal prior with standard deviation tau_g (default 1) is used. In the generalized-gamma model, a zero-centered normal prior with standard deviation q_prior_sd is placed on the signed shape parameter \(Q\). With the default covariate standardization, tau_t and tau_g apply to effects per one-standard-deviation change in a continuous covariate, while binary effects remain contrasts between their two supplied values. The internal intercept concerns an individual at the mean continuous-covariate values; returned intercepts are transformed back to covariate value zero.

The posterior can be sensitive to these choices when data are sparse, so a prior sensitivity analysis can be useful; see Klausch et al. (2026) for an example. Users can also supply a custom log-prior through log_prior_fun; see Section @ref(sec-user-defined-prior) and ?log_aft_prior.

Basic bayespim run

We now fit a Weibull model to the simulated data. The incidence covariates are passed to x_t, the prevalence covariates to x_g, and the baseline-test indicator to r. Here the test sensitivity is known, so we set kappa = 0.7 and update_kappa = FALSE; estimating \(\kappa\) is discussed in Section @ref(sec-kappa). The sampler runs chains = 4 MCMC chains in parallel, which is needed for reliable convergence diagnostics and to generate draws efficiently. The user’s machine should have at least chains free CPUs available.

mod_slice <- bayespim(
  v_obs = dat$v_obs,
  x_t = dat$x,
  x_g = dat$x,
  r = dat$r,
  kappa = 0.7,
  update_kappa = FALSE,
  ndraws = 1e3,
  warmup = 5e2,
  save_every = 1,
  standardize_covariates = TRUE,
  chains = 4,
  seed_chains = 1:4,
  min_effss = 800,
  update_till_converge = FALSE,
  sampler = "slice_collapsed",
  dist = "weibull"
)

The desired incidence distribution is passed via dist. The basic behavior of the Gibbs sampler is controlled by ndraws, warmup, chains, save_every, and seed_chains. The call above runs ndraws = 1000 iterations per chain and then stops (update_till_converge = FALSE). Each chain is randomly initialized, but the initialization and sampling are reproducible through seed_chains: one integer seed per chain, here seed_chains = 1:4 for chains = 4. Repeating the run with the same seeds reproduces the posterior chains exactly; seed_chains = NULL initializes the chains randomly.

A good warmup value is not known a priori, so it is advisable to run the sampler for a moderate number of iterations first and inspect the trace plots (Section @ref(sec-parameter-summaries)) to choose it.

By default, save_every = 1 stores every generated parameter draw in the returned $par chains. Convergence diagnostics and summary discard the first warmup generated iterations and then use every stored post-warm-up draw. warmup is always specified on the generated-iteration scale, independently of save_every: for example, warmup = 500 and save_every = 10 omit the 50 stored draws at iterations 10 through 500 in each chain, and the first retained draw is iteration 510. For values larger than one only every save_every-th draw is retained in memory and returned as $par; for example save_every = 10 saves every 10th posterior draw in each chain in memory. Such ‘thinning’ reduces memory and chain auto-correlation but decreases precision of estimation of posterior statistics and convergence diagnostics. If storing all draws would exceed available memory, especially when update_till_converge = TRUE may append several large updates, set save_every to a value greater than one. Storage selection does not alter the Markov-chain trajectory, but the intervening parameter draws are permanently discarded and diagnostics and summaries can then use only the stored draws. The setting is inherited by manual and automatic updates. Because retaining all draws generally gives the most informative diagnostics and posterior estimates, values above one are recommended only when memory is a practical constraint. Plotting has a separate thinning argument that can reduce rendering cost without discarding draws from the fitted object.

After the requested draws, bayespim evaluates and prints convergence diagnostics: the rank-normalized split R-hat of Vehtari et al. (2021) and the effective sample size (ESS), computed with the posterior package (Bürkner et al., 2026). These are checked against max_rhat (default 1.01) and min_effss (default chains * 100). To make the update workflow visible in this example, we use the stricter requirement min_effss = 800 rather than the four-chain default of 400.

Convergence diagnostics after 1000 iterations per chain (500 stored post-warm-up draws used).
Convergence criteria: R-hat <= 1.010 and ESS >= 800.0
 block        parameter R_hat   ESS
     t beta_t_intercept 1.005 446.5
     t       beta_t_x_1 1.002 778.8
     t       beta_t_x_2 1.005 576.8
     t          sigma_t 1.010 554.4
     g beta_g_intercept 1.004 549.7
     g       beta_g_x_1 1.006 678.0
     g       beta_g_x_2 1.005 631.1

The table reports diagnostics for the incidence (t) and prevalence (g) parameter blocks. After 1000 draws all R-hat values meet their threshold, but several effective sample sizes are below the deliberately strict requirement of 800. We therefore need more draws under this example’s convergence rule, which can be added by updating the run (Section @ref(sec-updating-runs)) or by letting bayespim update automatically (Section @ref(sec-auto-convergence)). The runtime of the sampler is stored in the fitted object.

mod_slice$runtime

Time difference of 4.020705 secs

Updating previous bayespim runs

If an initial run has not converged, it can be continued by passing the fitted object to the prev_run argument of a new bayespim call. The chains resume from their last draws and the new draws are merged with the previous run. The number of added draws is controlled by ndraws_update (default: the original ndraws).

mod_slice_update <- bayespim(
  prev_run = mod_slice,
  ndraws_update = 2e3,
  min_effss = 800
)

Updating previous MCMC run.

Convergence diagnostics after 3000 iterations per chain (2500 stored post-warm-up draws used).
Convergence criteria: R-hat <= 1.010 and ESS >= 800.0
 block        parameter R_hat    ESS
     t beta_t_intercept 1.001 2186.3
     t       beta_t_x_1 1.002 4351.5
     t       beta_t_x_2 1.001 2685.5
     t          sigma_t 1.001 2672.2
     g beta_g_intercept 1.001 2951.2
     g       beta_g_x_1 1.001 3355.5
     g       beta_g_x_2 1.001 3287.6

After a further 2000 iterations per chain (3000 in total), all R-hat values are at most 1.002 and all effective sample sizes exceed 2100, so the sampler has converged and we have ample draws for posterior inference. Data, sampler, priors, covariate-scaling constants, storage interval, and the end-of-chain random-number state are inherited from prev_run; we only repeat the stricter ESS criterion because convergence thresholds are call-specific. By default the original warm-up is retained; warmup_updated = TRUE instead increases the warm-up with each update.

Automatic updating till convergence

Instead of adding draws manually, bayespim can update automatically until convergence with update_till_converge = TRUE. If the initial ndraws draws are not sufficient, another ndraws_update draws are added and convergence is re-evaluated, and so on, until the criteria (max_rhat, min_effss) are met or the maximum number of draws maxit (default Inf) is reached. This can be requested directly, without a manual initial run.

mod_weibull <- bayespim(
  v_obs = dat$v_obs,
  x_t = dat$x,
  x_g = dat$x,
  r = dat$r,
  kappa = 0.7,
  update_kappa = FALSE,
  ndraws = 1e3,
  warmup = 5e2,
  save_every = 1,
  standardize_covariates = TRUE,
  chains = 4,
  seed_chains = 1:4,
  min_effss = 800,
  update_till_converge = TRUE,
  ndraws_update = 1e3,
  sampler = "slice_collapsed",
  dist = "weibull"
)

Not converged after 1000 iterations per chain; updating with 1000 iterations.
Convergence criteria: R-hat <= 1.010 and ESS >= 800.0
 block        parameter R_hat   ESS
     t beta_t_intercept 1.005 446.5
     t       beta_t_x_1 1.002 778.8
     t       beta_t_x_2 1.005 576.8
     t          sigma_t 1.010 554.4
     g beta_g_intercept 1.004 549.7
     g       beta_g_x_1 1.006 678.0
     g       beta_g_x_2 1.005 631.1
Updating previous MCMC run.

Converged after 2000 iterations per chain (1500 stored post-warm-up draws used).
Convergence criteria: R-hat <= 1.010 and ESS >= 800.0
 block        parameter R_hat    ESS
     t beta_t_intercept 1.001 1275.2
     t       beta_t_x_1 1.002 2508.5
     t       beta_t_x_2 1.002 1651.9
     t          sigma_t 1.001 1660.3
     g beta_g_intercept 1.002 1769.4
     g       beta_g_x_1 1.001 2044.3
     g       beta_g_x_2 1.002 1904.0

Here a single automatic update of 1000 iterations was enough to reach convergence at 2000 iterations per chain. The whole fit, including the automatic update, took about seven seconds. This example uses save_every = 1, so all 2000 draws per chain are retained. For a substantially longer automatic run where the growing $par object might exceed available R memory, rerun the model with, for example, save_every = 5; the fit will store every fifth generated state, and convergence checks and summaries will use all of those stored states without further thinning.

mod_weibull$runtime

Time difference of 7.404977 secs

We use mod_weibull as the fitted Weibull model in the remainder of this guide. If silent = TRUE is passed, the progress output shown above is suppressed; the diagnostics are still computed and stored in the convergence element of the returned object.

Posterior summaries of the model parameters

The posterior parameter draws are stored in the fitted object as a coda::mcmc.list named $par (Plummer et al., 2006). With the default save_every = 1, it contains every generated draw, including warm-up. With save_every > 1, it contains only the selected storage iterations. The plotting method discards stored draws belonging to the warm-up and may apply additional plot-only thinning for rendering speed; it does not modify the fitted object.

plot(mod_weibull, thinning = 5)

It can be helpful to thin the chains for faster plotting and to discard some additional warm-up. Here we show the trace and density plots for the incidence (latent-time) parameters after a thinning interval of five; the four chains are overlaid.

Trace and density plots for the incidence (latent-time) model parameters of the Weibull fit.

The chains mix well and are indistinguishable across the four colors, consistent with the R-hat and ESS diagnostics. Posterior summaries are obtained efficiently through the summary method, which reports posterior medians, 95% credible intervals, and the convergence diagnostics for each parameter block.

summary(mod_weibull)

Latent-time distribution: weibull
Incidence sampler: slice_collapsed

Parameters of the latent-time model
                  2.5%   50% 97.5% R_hat    ESS
beta_t_intercept 4.922 4.973 5.029 1.001 1275.2
beta_t_x_1       0.156 0.197 0.239 1.002 2508.5
beta_t_x_2       0.124 0.204 0.286 1.002 1651.9
sigma_t          0.177 0.207 0.241 1.001 1660.3

Parameters of the prevalence model
                   2.5%    50%  97.5% R_hat    ESS
beta_g_intercept -0.986 -0.857 -0.728 1.002 1769.4
beta_g_x_1        0.058  0.145  0.236 1.001 2044.3
beta_g_x_2        0.093  0.273  0.452 1.002 1904.0

Convergence criteria: R-hat <= 1.010 and ESS >= 800.0

MCMC iterations generated: 8000 (2000 per chain)
Parameter draws stored: 8000 (2000 per chain; save_every = 1)
Warm-up cutoff: 500 generated iterations per chain
Stored warm-up draws omitted: 2000 (500 per chain)
Posterior draws used: 6000 (1500 per chain)

The estimates are close to the values used in the data generation: the incidence intercept is near mu_t = 5, both incidence slopes near beta_t = 0.2, and the AFT scale near sigma_t = 0.2; the credible intervals cover the true values. The prevalence slopes are likewise near beta_g = 0.2, and the prevalence intercept near qnorm(0.2) = -0.84. By default summary discards the warm-up stored in the object; a different warm-up can be supplied through the warmup argument for this summary only.

For posterior means, additional quantiles, and Monte Carlo standard errors, coda’s own summary method can be used on the parameter chains. The helper trim_mcmc conveniently discards warm-up and thins an mcmc.list.

summary(trim_mcmc(mod_weibull$par, burnin = mod_weibull$warmup))

Iterations = 501:2000
Thinning interval = 1
Number of chains = 4
Sample size per chain = 1500

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                    Mean      SD  Naive SE Time-series SE
beta_t_intercept  4.9735 0.02684 0.0003465      0.0007390
beta_t_x_1        0.1970 0.02129 0.0002748      0.0004249
beta_t_x_2        0.2037 0.04127 0.0005328      0.0010258
sigma_t           0.2078 0.01662 0.0002146      0.0004079
beta_g_intercept -0.8572 0.06593 0.0008511      0.0015724
beta_g_x_1        0.1451 0.04586 0.0005921      0.0010301
beta_g_x_2        0.2728 0.09145 0.0011806      0.0021003

2. Quantiles for each variable:

                     2.5%     25%     50%     75%   97.5%
beta_t_intercept  4.92248  4.9560  4.9728  4.9904  5.0290
beta_t_x_1        0.15619  0.1826  0.1967  0.2110  0.2394
beta_t_x_2        0.12380  0.1756  0.2038  0.2308  0.2857
sigma_t           0.17677  0.1963  0.2072  0.2187  0.2414
beta_g_intercept -0.98616 -0.9014 -0.8574 -0.8131 -0.7280
beta_g_x_1        0.05766  0.1137  0.1447  0.1771  0.2362
beta_g_x_2        0.09282  0.2127  0.2732  0.3315  0.4522

The generalized-gamma fit is summarized in the same way and additionally reports the signed shape parameter q_t, for which q_t = 1 recovers the Weibull and q_t = 0 the log-normal. Here q_t is estimated with considerable posterior uncertainty, reflecting that interval-censored screening data carry limited information about the precise distributional shape; this is consistent with the small difference in information criteria between the Weibull and generalized-gamma fits found in Section @ref(sec-information-criteria).

Posterior cumulative incidence functions

A central question in screening research is the cumulative probability of progression as a function of time. BayesPIM estimates posterior predictive cumulative incidence functions (CIFs) with the function ppCIF. As motivated in the Overview, BayesPIM distinguishes two CIFs that handle prevalence differently, and one ppCIF call computes both.

The non-prevalent CIF is the incidence function among individuals who are disease-free at baseline. It disentangles incidence from prevalence and equals the cumulative distribution function of the latent event time,

\[\begin{align} F_t(t_0 \mid \boldsymbol{\beta}_t, \sigma_t) = \Pr(t \le t_0 \mid \boldsymbol{\beta}_t, \sigma_t). \end{align}\]

Like an ordinary CIF, it starts at zero at time zero. The mixture CIF additionally represents prevalence as a point mass at time zero. Writing \(\theta\) for the prevalence probability,

\[\begin{align} F_{\mathrm{mix}}(t_0 \mid \cdot) = \theta + (1 - \theta)\,F_t(t_0 \mid \boldsymbol{\beta}_t, \sigma_t), \end{align}\]

so the mixture CIF starts at \(\theta > 0\) at time zero and then accumulates incident events among the non-prevalent fraction. The mixture CIF is the joint representation, useful when the total baseline-to-disease burden (prevalent plus incident) is of interest; the non-prevalent CIF isolates incidence.

Both CIFs come in a marginal and a conditional variant. A marginal CIF integrates over the empirical covariate distribution and describes a randomly selected individual from the population. A conditional CIF fixes one or more covariates and describes a specific subgroup (for example, individuals with a given biomarker value). Conditionally, the prevalence probability is \(\theta(\mathbf{x}_g) = \Phi(\mathbf{x}_{g}'\boldsymbol{\beta}_g)\) and the incidence CIF is \(F_t(t_0 \mid \mathbf{x}_t, \boldsymbol{\beta}_t, \sigma_t)\); marginal CIFs average these over the covariates. Because closed forms are not available for all quantities, ppCIF obtains both CIFs by posterior predictive simulation. It considers every stored post-warm-up draw and then randomly selects the explicitly requested number pst_samples; for each selected draw it simulates prevalence indicators and event times and summarizes the results by the pointwise posterior median and 95% credible band.

We distinguish two use cases: obtaining predictive probabilities at a few time points, and obtaining a full curve for plotting.

set.seed(2025)
cif_pts <- ppCIF(mod_weibull, ppd_type = "percentiles", quant = c(0, 100, 200))
cif_pts$mixture$med_cdf         # posterior median mixture CIF
cif_pts$mixture$med_cdf_ci      # 2.5% and 97.5% posterior band
cif_pts$nonprevalent$med_cdf    # posterior median non-prevalent CIF

> cif_pts$mixture$med_cdf
[1] 0.240 0.362 0.897
> cif_pts$mixture$med_cdf_ci
         [,1]  [,2]  [,3]
[1,] 0.204000 0.320 0.853
[2,] 0.276025 0.407 0.936
> cif_pts$nonprevalent$med_cdf
[1] 0.0000000 0.1596530 0.8646518

ppCIF(mod_weibull, ppd_type = "percentiles", quant = c(0, 100, 200),
      fix_x_t = c(NA, 1))

Plotting CIFs

For plotting, ppCIF is evaluated on a fine grid of times and the resulting curve is drawn with plot.ppCIF. By default the plot method shows the mixture CIF; type = "nonprevalent" shows the non-prevalent CIF, and type = "both" shows them side by side on a common time scale.

set.seed(2026)
cif <- ppCIF(mod_weibull, pst_samples = 1e3, ppd_type = "percentiles",
             quant = seq(0, 300, length.out = 601))
plot(cif, type = "both", xlim = c(0, 300))

Posterior predictive mixture CIF (left) and non-prevalent CIF (right) with 95% credible bands.

The two panels make the role of prevalence explicit. The mixture CIF starts at about 0.24, the prevalence mass, and rises to one, whereas the non-prevalent CIF starts at zero and describes incidence in the initially healthy sub-population. If quant is left NULL, a default grid from zero to the maximum finite follow-up time is used. Setting ci = FALSE omits the credible bands.

Conditional CIFs are obtained in the same way, by fixing covariates through fix_x_t and fix_x_g. For a compact comparison, the median curves returned in $nonprevalent$med_cdf (and the corresponding grid in $quant) can be drawn together. The figure below overlays the conditional non-prevalent CIFs for the two levels of the discrete covariate, marginalizing over the continuous covariate.

set.seed(2027)
cif_x2_0 <- ppCIF(mod_weibull, fix_x_t = c(NA, 0), pst_samples = 1e3,
                  ppd_type = "percentiles", quant = seq(0, 300, length.out = 601))
set.seed(2028)
cif_x2_1 <- ppCIF(mod_weibull, fix_x_t = c(NA, 1), pst_samples = 1e3,
                  ppd_type = "percentiles", quant = seq(0, 300, length.out = 601))

plot(cif_x2_0, type = "nonprevalent", ci = FALSE, xlim = c(0, 300),
     main = "Conditional non-prevalent CIF by baseline covariate")
lines(cif_x2_1$quant, cif_x2_1$nonprevalent$med_cdf, col = "#D55E00", lwd = 2)
legend("bottomright", bty = "n", lwd = 2, col = c("#0072B2", "#D55E00"),
       legend = c("x2 = 0", "x2 = 1"))

Conditional non-prevalent CIFs for the two levels of the discrete covariate, marginalizing over the continuous covariate.

Because the incidence slope of the discrete covariate is positive, individuals with x2 = 1 reach a given cumulative incidence later than individuals with x2 = 0, so their curve lies to the right.

ppCIF also supports the inverse representation with ppd_type = "quantiles", which returns event times at cumulative probabilities supplied in perc (default a fine grid on \([0,1]\)). This is convenient when the quantities of interest are times-to-probability rather than probabilities-at-time. The plotting method handles both representations automatically.

Estimating test sensitivity

So far the test sensitivity \(\kappa\) was fixed at a known value with update_kappa = FALSE. When \(\kappa\) is unknown, it can be estimated by setting update_kappa = TRUE and supplying an informative Beta prior through kappa_prior = c(mean, sd). The prior mean must lie strictly between zero and one, and the standard deviation must satisfy \(0 < \mathrm{sd} < \sqrt{\mathrm{mean}(1-\mathrm{mean})}\); the Beta shape parameters are computed from these two moments. An informative prior is advisable because sensitivity and incidence can be weakly identified from screening data alone.

mod_kappa <- bayespim(
  v_obs = dat$v_obs, x_t = dat$x, x_g = dat$x, r = dat$r,
  update_kappa = TRUE,
  kappa_prior = c(0.7, 0.1),   # Beta prior with mean 0.7 and sd 0.1
  ndraws = 1e3, warmup = 5e2, chains = 4, seed_chains = 1:4,
  update_till_converge = TRUE, dist = "weibull"
)

When \(\kappa\) is estimated, a kappa block is added to the convergence diagnostics and to summary. If kappa_prior = NULL, a warning is issued and an uninformative \(\mathrm{Beta}(1,1)\) prior is used; we generally advise against this default.

Slice, collapsed, and Metropolis samplers

bayespim offers three samplers for the incidence parameters, selected with sampler. The default "slice_collapsed" augments only the latent screening interval and updates the parameters from the interval-censored likelihood, as described in Section @ref(sec-gibbs). The alternative "slice" augments the exact latent event times and then updates the parameters by univariate slice sampling, and "mh" augments the exact times and updates the parameters with a random-walk Metropolis step, as in the original implementation of Klausch et al. (2026). The two exact-time samplers are retained mainly for comparison. The collapsed sampler is the default because collapsing over the exact times reduces MCMC autocorrelation and speeds up convergence, and because the generalized-gamma model is supported only with the collapsed sampler.

Metropolis sampler and the proposal standard deviation

Unlike the slice samplers, the Metropolis sampler (sampler = "mh") requires a tuning parameter: the standard deviation prop_sd of the normal random-walk proposal (\(\sigma_t\) is proposed on the log scale). A good prop_sd yields an acceptance rate near 23%, which is often quoted as efficient for random-walk Metropolis under approximate normality (Roberts et al., 1997). A suitable value is not known a priori and has to be tuned; in our experience it usually lies between 0.001 and 0.1.

The helper search_prop_sd automates this search. It takes a short initial Metropolis fit and adjusts prop_sd with a heuristic rule until the acceptance rate falls within acc_bounds (default c(0.2, 0.25)), doubling the number of draws succ_min times (default 3) to confirm stability.

mod_mh_ini <- bayespim(
  v_obs = dat$v_obs, x_t = dat$x, x_g = dat$x, r = dat$r,
  kappa = 0.7, update_kappa = FALSE,
  ndraws = 1e3, warmup = 5e2, chains = 4, seed_chains = 1:4,
  sampler = "mh", prop_sd = 0.005, dist = "weibull"
)

search_sd <- search_prop_sd(m = mod_mh_ini)
search_sd$prop_sd

Iteration 1
Acceptance rate was: 0.746
prop_sd is set to 0.015
...
Success. Doubling number of MCMC draws: 4000
Finished calibrating proposal variance.

> search_sd$prop_sd
[1] 0.01590472

The calibrated prop_sd can then be passed to a full Metropolis run. In this example the acceptance rate is brought close to the target, but the incidence parameters still mix slowly under Metropolis sampling, needing many more draws than the collapsed slice sampler to reach comparable effective sample sizes. This illustrates why the collapsed slice sampler is the default; the Metropolis sampler is mainly of interest for reproducing the original method.

Slice sampler step size

The slice samplers have their own tuning parameter, the initial bracket width used in the step-out procedure (Neal, 2003), controlled by slice_width (default 1). It does not affect the target posterior, only computational efficiency. If the width is too small, the sampler may need many step-out steps to bracket the slice; if too large, it may spend more time shrinking the bracket. The default works well in most applications; if the slice sampler mixes slowly or is unexpectedly slow, adjusting slice_width can help.

User-defined prior function

By default bayespim uses the weakly informative priors of Section @ref(sec-prior-assumptions) through its internal log-prior function log_aft_prior. Users can supply a custom log-prior for the incidence parameters through log_prior_fun. The function must accept the named arguments eta, dist, beta_prior, tau_t, sig_prior, and q_prior_sd, and return a single numeric log-density. Here eta is the vector of incidence coefficients with log(sigma_t) appended in the last position (and the signed shape Q after it for dist = "gengamma"). When standardize_covariates = TRUE, this function receives the internally standardized incidence coefficients, so a custom prior is also defined on that scale. See ?log_aft_prior for the exact contract.

Internal scaling of times and covariates

The scale of the screening times passed to bayespim is often arbitrary (days, weeks, months, or years). Changing the time scale mainly shifts the AFT intercept \(\beta_{t0}\), which is regularized by its prior, so the amount of regularization would otherwise depend on the chosen units. By default (rescale_times = TRUE), bayespim rescales the screening times internally by the median finite observation time before fitting, and restores the returned times and incidence intercept to the original scale afterward. Dividing times by a constant \(c\) shifts the intercept by \(-\log(c)\), which is added back after fitting; the slopes and scale parameter are unaffected. This keeps the default prior on the intercept meaningful regardless of the time unit. Rescaling can be turned off with rescale_times = FALSE, in which case care should be taken with the intercept prior.

Covariate standardization is another affine reparameterization. For a continuous column \(x_j\), let \(z_j=(x_j-m_j)/s_j\), where \(m_j\) and \(s_j\) are its sample mean and standard deviation. If \(\widetilde\beta_j\) is sampled as the coefficient of \(z_j\), the coefficient per original covariate unit is

\[ \beta_j=\frac{\widetilde\beta_j}{s_j}. \]

For both the AFT and probit prevalence components, the original-scale intercept is obtained by subtracting \(\sum_j m_j\beta_j\) from the internally sampled intercept. The AFT intercept additionally receives the \(+\log(c)\) time-scale correction described above. Thus

\[ \beta_{t0}=\widetilde\beta_{t0}+\log(c)-\sum_jm_{tj}\beta_{tj}, \qquad \beta_{g0}=\widetilde\beta_{g0}-\sum_jm_{gj}\beta_{gj}. \]

Every posterior draw is transformed before it is returned, so credible intervals and parameter dependence are preserved. Scale/shape parameters such as \(\sigma_t\) and generalized-gamma \(Q\), and test sensitivity \(\kappa\), are unaffected. Information-criterion and posterior-predictive calculations reconstruct the centered parameterization internally and use standardized covariates for numerical stability, while fixed covariate values supplied to ppCIF remain on the original scale. The fitted centers and standard deviations are stored in $covariate_scaling and inherited unchanged by model updates.