## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----load_data, echo=FALSE----------------------------------------------------
load("vignette_data.rdata") # if it's in the same directory as your .Rmd
## -----------------------------------------------------------------------------
library(BayesPIM)
## ----eval = FALSE-------------------------------------------------------------
# # 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
# )
## -----------------------------------------------------------------------------
head(dat$Vobs)
## ----eval = FALSE-------------------------------------------------------------
# # An initial model fit with fixed test sensitivity kappa (approx. 1-3 minutes, depending on machine)
# 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 = 4,
# prop.sd.X = 0.008,
# parallel = TRUE,
# dist.X = 'weibull'
# )
#
# # Searching starting values by one naive run
# # Starting Gibbs sampler with 3 chains and 5000 iterations.
# # Now doing main run.
# # Starting Gibbs sampler with 4 chains and 10000 iterations.
## ----eval=FALSE---------------------------------------------------------------
# # Inspect results
# mod$runtime # runtime of Gibbs sampler
# plot( trim.mcmc( mod$par.X.all, thining = 10) ) # MCMC chains including burn-in
# plot( trim.mcmc( mod$par.X.bi, thining = 10) ) # MCMC chains excluding burn-in
# summary( mod$par.X.bi)
# apply(mod$ac.X, 2, mean) # Acceptance rates per chain
# gelman.diag(mod$par.X.bi) # Gelman convergence diagnostics
## ----eval = FALSE-------------------------------------------------------------
# # An initial model fit with a moderate number of ndraws (here 1e3)
# mod.ini = bayes.2S( Vobs = dat$Vobs,
# Z.X = dat$Z,
# Z.W = dat$Z,
# r= dat$r,
# kappa = 0.7,
# update.kappa = F,
# ndraws= 1e3,
# chains = 4,
# prop.sd.X = 0.005,
# parallel = T,
# dist.X = 'weibull'
# )
#
# # Searching starting values by one naive run
# # Starting Gibbs sampler with 3 chains and 5000 iterations.
# # Now doing main run.
# # Starting Gibbs sampler with 4 chains and 1000 iterations.
#
# # Running the automated search
# search.sd <- search.prop.sd(m = mod.ini)
#
# # Iteration 1
# # Acceptance rate was: 0.564
# # prop.sd.X is set to 0.011
# # Iteration 2
# # Acceptance rate was: 0.417
# # prop.sd.X is set to 0.019
# # Iteration 3
# # Acceptance rate was: 0.113
# # prop.sd.X is set to 0.011
# # Iteration 4
# # Acceptance rate was: 0.216
# # Success. Doubling number of MCMC draws: 2000
# # Iteration 5
# # Acceptance rate was: 0.223
# # Success. Doubling number of MCMC draws: 4000
# # Iteration 6
# # Acceptance rate was: 0.222
# # Finished calibrating proposal variance.
## -----------------------------------------------------------------------------
print(search.sd$prop.sd.X)
## ----eval=FALSE---------------------------------------------------------------
# # Model updating
# mod_update = bayes.2S( prev.run = mod ) # ndraws additional MCMC draws
#
# # Updating previous MCMC run.
# # Starting Gibbs sampler with 4 chains and 10000 iterations.
#
# mod_update = bayes.2S( prev.run = mod, ndraws.update = 1e3 ) # ndraws.update additional MCMC draws
#
# # Updating previous MCMC run.
# # Starting Gibbs sampler with 4 chains and 1000 iterations.
#
## ----eval = FALSE-------------------------------------------------------------
# # Get posterior predictive marginal CIF, we provide percentiles and obtain quantiles
# 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))
## ----fig.width=8, fig.height=4, eval = FALSE----------------------------------
# # Comparison 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, ty = 'l', ylim=c(0,1),
# xlim=c(0,300), xlab = 'Time', ylab ='Cum. prob.', main = 'Non-prevalence CIF')
# lines(cif_nonprev$med.cdf.ci[1,], cif_nonprev$perc, lty=2, col=2)
# lines(cif_nonprev$med.cdf.ci[2,], cif_nonprev$perc, lty=2, col=2)
# plot(cif_mix$med.cdf, cif_nonprev$perc, ty = 'l', ylim=c(0,1),
# xlim=c(0,300), xlab = 'Time', ylab ='Cum. prob.', main = 'Mixture CIF')
# lines(cif_mix$med.cdf.ci[1,], cif_nonprev$perc, lty=2, col=2)
# lines(cif_mix$med.cdf.ci[2,], cif_nonprev$perc, lty=2, col=2)
## ----echo=FALSE, out.width='90%'----------------------------------------------
# Actually show the saved image:
knitr::include_graphics("cif_marg.png")
## ----eval=FALSE---------------------------------------------------------------
# # Alternatively, we can provide quantiles and obtain percentiles
# 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)
## ----eval=FALSE---------------------------------------------------------------
# get.ppd.2S(mod, pst.samples = 1e3, type = 'xstar', ppd.type = "percentiles", quant = c(0,100))
## -----------------------------------------------------------------------------
print(quants)
## ----eval=FALSE---------------------------------------------------------------
# # Conditional CIFs, example conditional mixture CIF integrating over one covariate
# 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))
## ----fig.width=4, fig.height=4, eval=FALSE------------------------------------
# # Plot of CIFs for three levels of the first covariate
# par(mfrow = c(1,1))
# plot(cif_mix_m1$med.cdf, cif_mix_m1$perc, ty = 'l', ylim=c(0,1),
# xlim=c(0,300), xlab = 'Time', ylab ='Cum. prob.', main = 'Conditional mixture CIF')
# 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)
## ----echo=FALSE, out.width='50%'----------------------------------------------
# Actually show the saved image:
knitr::include_graphics("cif_cond.png")
## ----eval=FALSE---------------------------------------------------------------
# # An exponential model
# mod_exp = bayes.2S( Vobs = dat$Vobs,
# Z.X = dat$Z,
# Z.W = dat$Z,
# r= dat$r,
# kappa = 0.7,
# update.kappa = F,
# ndraws= 1e4,
# chains = 4,
# prop.sd.X = 0.008,
# parallel = T,
# fix.sigma.X = T,
# sig.prior.X = 1,
# dist.X = 'weibull'
# )
#
# # Searching starting values by one naive run
# # Starting Gibbs sampler with 3 chains and 5000 iterations.
# # Now doing main run.
# # Starting Gibbs sampler with 4 chains and 10000 iterations.
#
# # Get information criteria
# get.IC_2S(mod, samples = 1e3)
#
# # > IC_weib
# # WAIC1 WAIC2 DIC
# # [1,] 2083.599 2083.694 2083.359
#
# get.IC_2S(mod_exp, samples = 1e3)
#
# # > IC_exp
# # WAIC1 WAIC2 DIC
# # [1,] 2392.869 2392.945 2393.931
#