bayes.2S {BayesPIM} | R Documentation |
Fitting Bayesian Prevalence-Incidence Mixture Model
Description
Estimates the Pattern Mixture model of Klausch et al. (2025) 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
bayes.2S(
Vobs,
Z.X = NULL,
Z.W = NULL,
r = NULL,
dist.X = "weibull",
kappa = 0.5,
update.kappa = FALSE,
kappa.prior = NULL,
ndraws = 1000,
prop.sd.X = NULL,
chains = 3,
thining = 1,
parallel = TRUE,
update.till.converge = FALSE,
maxit = Inf,
conv.crit = "upper",
min_effss = chains * 10,
beta.prior = "norm",
beta.prior.X = 1,
sig.prior.X = 1,
tau.w = 1,
fix.sigma.X = FALSE,
prev.run = NULL,
update.burnin = TRUE,
ndraws.update = NULL,
prev = TRUE,
vanilla = FALSE,
ndraws.naive = 10000,
naive.run.prop.sd.X = prop.sd.X,
par.exp = FALSE,
collapsed.g = TRUE,
k.prior = 1,
fix.k = FALSE
)
Arguments
Vobs |
A list of length |
Z.X |
A numeric matrix of dimension |
Z.W |
A numeric matrix of dimension |
r |
A binary vector of length |
dist.X |
Character. Specifies the distribution for the time-to-incidence variable; choices are |
kappa |
Numeric. The test sensitivity value to be used if |
update.kappa |
Logical. If |
kappa.prior |
A numeric vector of length 2. When specified, a Beta distribution prior is used for |
ndraws |
Integer. The total number of MCMC draws for the main Gibbs sampler. |
prop.sd.X |
Numeric. The standard deviation for the proposal (jumping) distribution in the Metropolis sampler used for updating |
chains |
Integer. The number of MCMC chains to run. |
thining |
Integer. The thinning interval for the MCMC sampler. |
parallel |
Logical. If |
update.till.converge |
Logical. If |
maxit |
Numeric. The maximum number of MCMC draws allowed before interrupting the update process when |
conv.crit |
Character. Specifies whether the convergence check uses the point estimate ( |
min_effss |
Integer. The minimum effective sample size required for each parameter before convergence is accepted during iterative updating. |
beta.prior |
Character. Specifies the type of prior for the incidence regression coefficients ( |
beta.prior.X |
Numeric. The hyperparameter for the prior distribution of the regression coefficients ( |
sig.prior.X |
Numeric. The hyperparameter (standard deviation) for a half-normal prior on the scale parameter ( |
tau.w |
Numeric. The hyperparameter (standard deviation) for the normal prior distribution of the regression coefficients ( |
fix.sigma.X |
Logical. If |
prev.run |
Optional. An object of class |
update.burnin |
Logical. If |
ndraws.update |
Integer. The number of MCMC draws for updating a previous run or for convergence updates. If unspecified, |
prev |
Logical. If |
vanilla |
Logical. If |
ndraws.naive |
Integer. The number of MCMC draws for a preliminary vanilla run used to obtain starting values. Increase if initial values lead to issues (e.g., an infinite posterior). |
naive.run.prop.sd.X |
Numeric. The standard deviation for the proposal distribution used in the vanilla run. Adjust only if the acceptance rate is significantly off target, as indicated by an interruption message. |
par.exp |
Logical. If |
collapsed.g |
Logical. If |
k.prior |
Experimental prior parameter for generalized gamma; currently not used. |
fix.k |
Experimental fixing of prior parameter for generalized gamma; currently not used. |
Details
This Bayesian prevalence-incidence mixture model (PIM) characterizes the time to incidence using an accelerated failure time (AFT) model of the form:
\log(x_i) = \bm{z}_{xi}' \bm{\beta}_x + \sigma \epsilon_i
where \epsilon_i
is chosen such that x_i
follows a weibull
, lognormal
, or loglog
(log-logistic) distribution, as specified by the dist
argument. The covariate vector \bm{z}_{xi}
for individual i
is provided in the Z.X
matrix.
Baseline prevalence is modeled using a probit formulation Pr(g_i=1 | \bm{z}_{wi}) = Pr(w_i > 0 | \bm{z}_{wi})
with
w_i = \bm{z}_{wi}' \bm{\beta}_w + \psi_i
where \psi_i
follows a standard normal distribution, and the covariate vector \bm{z}_{wi}
is given in the Z.W
matrix. The latent variable w_i
determines prevalence status, such that g_i = 1
if w_i > 0
and g_i = 0
otherwise.
The argument Vobs
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 time6
).
By convention, every vector in Vobs
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 Vobs
. 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, a Beta prior is used, centered on the first element of kappa.prior
, with a standard deviation equal to its second element. An internal optimization process finds the Beta prior hyperparameters that best match this choice. If the chosen prior is not feasible, unexpected behavior may occur. If kappa.prior
is not specified (the default), an uninformative uniform(0,1) prior is used. 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. The Metropolis step used for sampling the parameters of the incidence model applies a normal proposal (jumping) distribution with a standard deviation prop.sd.X
, 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. Alternatively, the function search.prop.sd provides a heuristic for selecting an effective proposal distribution standard deviation.
If parallel = TRUE
, the Gibbs sampler runs in parallel with one chain per CPU (if possible), using the foreach
package. If this package causes issues on some operating systems, set parallel = FALSE
or use the bayes.2S_seq function, which iterates over 1:chains
using a for
loop. This sequential function may also be useful in Monte Carlo simulations that parallelize experimental replications using foreach
.
We recommend running at least two chains in parallel, and preferably more, to facilitate standard MCMC diagnostics such as the Gelman-Rubin R
statistic. 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.convergence = TRUE
allows bayes.2S
to run until convergence. Convergence is achieved when R < 1.1
for all parameters and the minimum effective sample size min_effs
is reached. The sampler continues updating until convergence is attained or maxit
is reached.
The priors for the regression coefficients in the prevalence and incidence models can be controlled using beta.prior
, beta.prior.X
, sig.prior.X
, and tau.w
. Specifically:
-
beta.prior
determines the prior type for\beta_{xj}
(eithernormal
or Student-t
t
). -
beta.prior.X
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
\sigma
, withsig.prior.X
controlling the standard deviation.A zero-centered normal prior is assigned to
\beta_{wj}
, withtau.w
controlling its standard deviation (default: standard normal).
Sometimes model fitting can be improved by fixing the \sigma
parameter to a value, which is achieved through setting fix.sigma.X = TRUE
. Then, the value specified as sig.prior.X
is regarded as the correct value for \sigma
, akin to a point prior on this value. The functionality can also be used to obtain the exponential distribution, aking to a Markov model. For this choose dist='weibull'
, sig.prior.X = 1
, and fix.sigma.X=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 <- bayes.2S(...)
, it can be updated using mod_update <- bayes.2S(prev.run = mod)
. By default, ndraws
additional iterations are added unless otherwise specified via ndraws.update
. When updating, the number of discarded burn-in draws can be adjusted to half of the total draws (update.burnin = TRUE
) or remain at the initial number (update.burnin = FALSE
).
The Gibbs sampler requires starting values, which are obtained from an initial Bayesian interval-censored survival model using the specified dist
distribution. The jumping distribution variance and the number of MCMC draws for this initialization are controlled via ndraws.naive
and naive.run.prop.sd.X
. The default values typically suffice but may need adjustment if initialization fails (e.g., increasing ndraws.naive
or tuning naive.run.prop.sd.X
). If starting values are found but still lead to an infinite posterior at initialization, the error "Bad starting values" is returned. Then it usually sufficces to re-run bayes.2S
with an increased ndraws.naive
value.
Value
A list containing the following elements:
par.X.all |
An |
par.X.bi |
An |
X |
A matrix of posterior draws for the latent event times |
C |
A matrix of posterior draws for prevalence class membership |
ac.X |
A matrix with MCMC draws in rows and chains in columns, where each row indicates whether the Metropolis sampler accepted (1) or rejected (0) a sample. |
ac.X.cur |
Same as |
dat |
A data frame containing the last observed interval. |
priors |
A list of prior specifications for the model parameters, including |
runtime |
The total runtime of the MCMC sampler. |
Additionally, most input arguments are returned as part of the output for reference.
References
T. Klausch, B. I. Lissenberg-Witte, and V. M. Coupe (2024). "A Bayesian prevalence-incidence mixture model for screening outcomes with misclassification.", doi:10.48550/arXiv.2412.16065.
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, doi:10.2307/2669940.
Examples
library(BayesPIM)
# Generate data according to the Klausch et al. (2024) PIM
set.seed(2025)
dat <- gen.dat(kappa = 0.7, n = 1e3, theta = 0.2,
p = 1, p.discrete = 1,
beta.X = c(0.2, 0.2), beta.W = c(0.2, 0.2),
v.min = 20, v.max = 30, mean.rc = 80,
sigma.X = 0.2, mu.X = 5, dist.X = "weibull",
prob.r = 1)
# Initial model fit with fixed test sensitivity kappa (approx. 1-3 minutes runtime)
mod <- bayes.2S(Vobs = dat$Vobs,
Z.X = dat$Z,
Z.W = dat$Z,
r = dat$r,
kappa = 0.7,
update.kappa = FALSE,
ndraws = 1e4,
chains = 2,
prop.sd.X = 0.008,
parallel = TRUE,
dist.X = "weibull")
# Inspect results
mod$runtime # Runtime of the Gibbs sampler
plot(trim.mcmc(mod$par.X.all, thining = 10)) # MCMC chains including burn-in, also see ?trim.mcmc
plot(trim.mcmc(mod$par.X.bi, thining = 10)) # MCMC chains excluding burn-in
apply(mod$ac.X, 2, mean) # Acceptance rates per chain
gelman.diag(mod$par.X.bi) # Gelman convergence diagnostics
# Model updating
mod_update <- bayes.2S(prev.run = mod) # Adds ndraws additional MCMC draws
mod_update <- bayes.2S(prev.run = mod,
ndraws.update = 1e3) # Adds ndraws.update additional MCMC draws
# Example with kappa estimated/updated
mod2 <- bayes.2S(Vobs = dat$Vobs,
Z.X = dat$Z,
Z.W = dat$Z,
r = dat$r,
kappa = 0.7,
update.kappa = TRUE,
kappa.prior = c(0.7, 0.1), # Beta prior, mean = 0.7, s.d. = 0.1
ndraws = 1e4,
chains = 2,
prop.sd.X = 0.008,
parallel = TRUE,
dist.X = "weibull")
# Inspect results
mod2$runtime # runtime of Gibbs sampler
plot( trim.mcmc( mod2$par.X.all, thining = 10) ) # kappa returned as part of the mcmc.list