## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
library(BayesTSM)

# Generate data
set.seed(1)
dat = gendat(
             n = 1000,  # Sample size
             p = 2,     # Number of normally distributed covariates
             sigma.X = 0.3,         # True scale parameter of X
             mu.X    = 2,           # True intercept parameter of X
             beta.X  = c(0.5,0.5),  # True slope parameters of all covariates of X
             sigma.S = 0.5,         # True scale parameter of S
             mu.S    = 1,           # True intercept parameter of S
             beta.S  = c(0.5,0.5),  # True slope parameters of all covariates of S
             dist.X  = 'weibull',   # Distribution of X
             dist.S  = 'weibull',   # Distribution of S
             v.min   = 1,           # Minimum time between screening moments
             v.max   = 5,           # Maximum time between screening moments
             Tmax    = 2e2,         # Maximum number of screening times
             mean.rc = 10           # Mean time to right censoring, exponential distribution
)

## -----------------------------------------------------------------------------
head(dat)

## -----------------------------------------------------------------------------
table(dat$d)

## ----eval = FALSE-------------------------------------------------------------
# # Run bayestsm Gibbs sampler with data augmentation and slice sampling of the parameters
# 
# d = dat$d
# L = dat$L
# R = dat$R
# Z = dat[, c("Z.1", "Z.2")]
# 
# mod_slice = bayestsm(
#   d = d,
#   L = L,
#   R = R,
#   Z.X = Z,
#   Z.S = Z,
#   mc = 1e4,
#   warmup = 5e2,
#   thinning = 10,
#   chains = 4,
#   update_till_convergence = FALSE,
#   MH = FALSE,
#   dist.X = "weibull",
#   dist.S = "weibull",
#   seed_chains = 1:4
# )

## ----eval = FALSE-------------------------------------------------------------
# plot(mod_slice, plot.X = F)

## ----mcmc-traceplot, echo = FALSE, out.width = "100%", fig.cap = "Trace plots for the MCMC chains of the s-model (after thinning, 10000 raw draws)."----
knitr::include_graphics("figures/mcmc-traceplot.png")

## ----eval = FALSE-------------------------------------------------------------
# # Updating previous run
# mod_slice_update = bayestsm(
#   prev.run = mod_slice, # pass previous bayestsm object
#   mc_update = 2e4
# )

## ----eval = FALSE-------------------------------------------------------------
# # Automated updating till convergence
# mod_slice_weibull = bayestsm(
#   d = d,
#   L = L,
#   R = R,
#   Z.X = Z,
#   Z.S = Z,
#   mc = 1e4,
#   warmup = 5e2,
#   thinning = 10,
#   chains = 4,
#   update_till_convergence = TRUE,
#   mc_update = 1e4,
#   MH = FALSE,
#   dist.X = "weibull",
#   dist.S = "weibull"
# )

## ----eval = FALSE-------------------------------------------------------------
# mod_slice_weibull$runtime

## ----eval = FALSE-------------------------------------------------------------
# # Lognormal model
# mod_slice_lognormal = bayestsm(
#    d = d,
#    L = L,
#    R = R,
#    Z.X = Z,
#    Z.S = Z,
#    mc = 1e4,
#    warmup = 5e2,
#    thinning = 10,
#    chains = 4,
#    update_till_convergence = TRUE,
#    mc_update = 1e4,
#    MH = FALSE,
#    dist.X = "lognormal",
#    dist.S = "lognormal",
#    seed_chains = 5:8
#  )

## ----eval = FALSE-------------------------------------------------------------
# # Exponential model
# mod_slice_exponential = bayestsm(
#    d = d,
#    L = L,
#    R = R,
#    Z.X = Z,
#    Z.S = Z,
#    mc = 1e4,
#    warmup = 5e2,
#    thinning = 10,
#    chains = 4,
#    update_till_convergence = TRUE,
#    mc_update = 1e4,
#    MH = FALSE,
#    dist.X = "weibull",
#    dist.S = "weibull",
#    fix.sigma.X = T, # Fix sigma.X at sig.prior.X (default 1)
#    fix.sigma.S = T, # Fix sigma.S at sig.prior.S (default 1)
#    seed_chains = 9:12
#  )

## ----eval=F-------------------------------------------------------------------
# get_IC(mod_slice_weibull, warmup =500, cores = NULL)

## ----eval=F-------------------------------------------------------------------
# get_IC(mod_slice_lognormal, warmup =500, cores = NULL)

## ----eval=F-------------------------------------------------------------------
# get_IC(mod_slice_exponential, warmup =500, cores = NULL)

## ----eval=FALSE---------------------------------------------------------------
# plot(mod_slice_weibull)

## ----eval=FALSE---------------------------------------------------------------
# plot(mod_slice_weibull, warmup = 500)

## ----eval = FALSE-------------------------------------------------------------
# summary(mod_slice_weibull, warmup = 500)

## ----eval = FALSE-------------------------------------------------------------
# summary(mod_slice_weibull$par.X)

## ----eval = FALSE-------------------------------------------------------------
# summary( trim.mcmc( mod_slice_weibull$par.X, burnin = 500) )

## ----eval = FALSE-------------------------------------------------------------
# ppCIF( mod_slice_weibull,
#        type = 'quantiles',
#        warmup = 500,
#        quant = c(5, 10) )

## ----eval = FALSE-------------------------------------------------------------
# fix_Z.X = c(1, 1)
# fix_Z.S = c(1, 1)

## ----eval = FALSE-------------------------------------------------------------
# fix_Z.X = c(1, NA)
# fix_Z.S = c(1, NA)

## ----eval = FALSE-------------------------------------------------------------
# ppCIF( mod_slice_weibull,
#        type = 'quantiles',
#        warmup = 500,
#        quant = c(5, 10),
#        fix_Z.X = c(1, NA),
#        fix_Z.S = c(1, NA)
# )

## ----eval = FALSE-------------------------------------------------------------
# pp_grid <- ppCIF( mod_slice_weibull, type = 'quantiles',  warmup = 500 )
# plot(pp_grid, xlim=c(0,50))

## ----ppCIF, echo = FALSE, out.width = "100%", fig.cap = "Posterior predictive CIFs obtained via default quantile grid."----
knitr::include_graphics("figures/ppCIF_quantiles.png")

## ----eval = FALSE-------------------------------------------------------------
# pq_grid <- ppCIF( mod_slice_weibull, type = 'percentiles',  warmup = 500 )
# plot(pq_grid, xlim=c(0,50))

## ----pqCIF, echo = FALSE, out.width = "100%", fig.cap = "Posterior predictive CIFs obtained via default percentile grid."----
knitr::include_graphics("figures/ppCIF_percentiles.png")

## ----eval = FALSE-------------------------------------------------------------
# mod_MH = bayestsm(
#   d = d,
#   L = L,
#   R = R,
#   Z.X = Z,
#   Z.S = Z,
#   mc = 5e5,
#   warmup = 1e5,
#   thinning = 100,
#   chains = 4,
#   update_till_convergence = TRUE,
#   MH = TRUE,
#   dist.X = "weibull",
#   dist.S = "weibull",
#   seed_chains = 1:4
# )

## ----eval = F-----------------------------------------------------------------
# summary(mod_MH)

## ----eval = F-----------------------------------------------------------------
# mod_MH$runtime

## ----eval = FALSE-------------------------------------------------------------
# ess_per_second_slice <- mod_slice_weibull$convergence$eff_s /
#                         as.numeric(mod_slice_weibull$runtime, units = "secs")
# 
# ess_per_second_MH <- mod_MH$convergence$eff_s /
#                      as.numeric(mod_MH$runtime, units = "secs")
# 
# rbind(ess_per_second_slice, ess_per_second_MH)

## ----eval = F-----------------------------------------------------------------
# log_aft_prior <- function(eta, tau = 4, sig.prior = 1, beta.prior = "t") { ... }

## ----eval = FALSE-------------------------------------------------------------
#   medLR = median(c(L[is.infinite(R)], R[is.finite(R)]))
#   L   = L/medLR
#   R   = R/medLR