get.ppd.2S {BayesPIM} | R Documentation |
Posterior Predictive Cumulative Incidence Function
Description
Computes the posterior predictive cumulative incidence function (CIF) from a bayes.2S prevalence-incidence mixture model. The function can return quantiles corresponding to user-specified percentiles (i.e., time points at which the cumulative probability reaches certain thresholds), or vice versa (percentiles at user-specified quantiles). Additionally, it allows for marginal or conditional CIFs of either the mixture population (including prevalence as a point-mass at time zero) or the non-prevalent (healthy) subpopulation.
Usage
get.ppd.2S(
mod,
fix_Z.X = NULL,
fix_Z.W = NULL,
pst.samples = 1000,
perc = seq(0, 1, 0.01),
type = "x",
ppd.type = "percentiles",
quant = NULL
)
Arguments
mod |
A fitted prevalence-incidence mixture model of class |
fix_Z.X |
Either |
fix_Z.W |
Same as |
pst.samples |
Integer; number of posterior samples to draw when computing the
posterior predictive CIF. Must not exceed the total available posterior samples
in |
perc |
A numeric vector of cumulative probabilities (i.e., percentiles in (0,1))
for which time points are returned when |
type |
Character; |
ppd.type |
Character; |
quant |
A numeric vector of time points for which the function returns
cumulative probabilities when |
Details
For a prevalence-incidence mixture model, some fraction of the population may already have experienced the event (prevalent cases) at baseline, while the remaining (healthy) fraction has not. This function estimates the CIF in two ways:
-
type = "xstar"
(mixture CIF): Includes a point-mass at time zero representing baseline prevalence, with incidence beginning thereafter. -
type = "x"
(non-prevalent CIF): Excludes prevalent cases, so it only shows incidence among the initially healthy subpopulation.
You may request a marginal CIF by setting both fix_Z.X = NULL
and
fix_Z.W = NULL
, thus integrating over all covariates. Alternatively, a
conditional CIF can be obtained by partially or fully specifying fixed
covariate values in fix_Z.X
(and optionally fix_Z.W
) while integrating
out the unspecified covariates (NA
entries).
The function operates in two main modes:
-
ppd.type = "quantiles"
: Given a set ofperc
(cumulative probabilities), returns corresponding quantiles (time points). -
ppd.type = "percentiles"
: Given a set ofquant
(time points), returns corresponding percentiles (cumulative probabilities).
Value
A list
with some or all of the following elements:
med.cdf
-
If
ppd.type = "quantiles"
,med.cdf
contains the median quantiles (time points) across posterior samples for each percentile inperc
.If
ppd.type = "percentiles"
,med.cdf
contains the median cumulative probabilities across posterior samples for each time point inquant
.
med.cdf.ci
-
A 2-row matrix with the 2.5% and 97.5% posterior quantiles for the estimated
med.cdf
, reflecting uncertainty. quant
-
If
ppd.type = "percentiles"
, this is a copy of the inputquant
(time points). perc
-
If
ppd.type = "quantiles"
, this is a copy of the inputperc
(cumulative probabilities).
Examples
# 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)
# Fit a bayes.2S model (example with moderate ndraws = 2e4)
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")
###################
# (1) Provide percentiles, get back quantiles (times)
###################
cif_nonprev <- get.ppd.2S(mod, pst.samples = 1e3, type = "x",
ppd.type = "quantiles", perc = seq(0, 1, 0.01))
cif_mix <- get.ppd.2S(mod, pst.samples = 1e3, type = "xstar",
ppd.type = "quantiles", perc = seq(0, 1, 0.01))
# Plot: Non-prevalent stratum CIF vs. mixture CIF (marginal)
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1,2))
plot(cif_nonprev$med.cdf, cif_nonprev$perc, type = "l", xlim = c(0,300), ylim = c(0,1),
xlab = "Time", ylab = "Cumulative Incidence")
lines(cif_nonprev$med.cdf.ci[1,], cif_nonprev$perc, lty = 2)
lines(cif_nonprev$med.cdf.ci[2,], cif_nonprev$perc, lty = 2)
plot(cif_mix$med.cdf, cif_mix$perc, type = "l", xlim = c(0,300), ylim = c(0,1),
xlab = "Time", ylab = "Cumulative Incidence")
lines(cif_mix$med.cdf.ci[1,], cif_mix$perc, lty = 2)
lines(cif_mix$med.cdf.ci[2,], cif_mix$perc, lty = 2)
###################
# (2) Provide quantiles (times), get back percentiles (cumulative probabilities)
###################
cif2_nonprev <- get.ppd.2S(mod, pst.samples = 1e3, type = "x",
ppd.type = "percentiles", quant = 1:300)
cif2_mix <- get.ppd.2S(mod, pst.samples = 1e3, type = "xstar",
ppd.type = "percentiles", quant = 1:300)
# Plot: Non-prevalent vs. mixture CIF using times on the x-axis
plot(cif2_nonprev$quant, cif2_nonprev$med.cdf, type = "l", xlim = c(0,300), ylim = c(0,1),
xlab = "Time", ylab = "Cumulative Incidence")
lines(cif2_nonprev$quant, cif2_nonprev$med.cdf.ci[1,], lty = 2)
lines(cif2_nonprev$quant, cif2_nonprev$med.cdf.ci[2,], lty = 2)
plot(cif2_mix$quant, cif2_mix$med.cdf, type = "l", xlim = c(0,300), ylim = c(0,1),
xlab = "Time", ylab = "Cumulative Incidence")
lines(cif2_mix$quant, cif2_mix$med.cdf.ci[1,], lty = 2)
lines(cif2_mix$quant, cif2_mix$med.cdf.ci[2,], lty = 2)
###################
# (3) Conditional CIFs by fixing some covariates
###################
cif_mix_m1 <- get.ppd.2S(mod, fix_Z.X = c(-1, NA), pst.samples = 1e3,
type = "xstar", ppd.type = "quantiles", perc = seq(0,1,0.01))
cif_mix_0 <- get.ppd.2S(mod, fix_Z.X = c(0, NA), pst.samples = 1e3,
type = "xstar", ppd.type = "quantiles", perc = seq(0,1,0.01))
cif_mix_p1 <- get.ppd.2S(mod, fix_Z.X = c(1, NA), pst.samples = 1e3,
type = "xstar", ppd.type = "quantiles", perc = seq(0,1,0.01))
# Plot: mixture CIF for three different values of the first covariate
par(mfrow = c(1,1))
plot(cif_mix_m1$med.cdf, cif_mix_m1$perc, type = "l", xlim = c(0,300), ylim = c(0,1),
xlab = "Time", ylab = "Cumulative Incidence", col=1)
lines(cif_mix_0$med.cdf, cif_mix_m1$perc, col=2)
lines(cif_mix_p1$med.cdf, cif_mix_m1$perc, col=3)
par(oldpar)