## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  warning = FALSE,
  message = FALSE
)
have_mcmc <- requireNamespace("MCMCpack", quietly = TRUE) &&
             requireNamespace("coda",     quietly = TRUE)
knitr::opts_chunk$set(eval = have_mcmc)

## -----------------------------------------------------------------------------
# library(BetaDanish)
# data("remission")
# fit_bayes <- bayes_betadanish(
#   time     = remission$time,
#   status   = remission$status,
#   submodel = TRUE,
#   burnin   = 2000, mcmc = 5000,
#   tune     = 0.5,
#   seed     = 1
# )
# print(fit_bayes)

## -----------------------------------------------------------------------------
# draws <- fit_bayes$draws
# op <- par(mfrow = c(3, 1), mar = c(3, 4, 2, 1))
# coda::traceplot(draws)
# par(op)

## -----------------------------------------------------------------------------
# post_mean <- summary(draws)$statistics[, "Mean"]
# b <- post_mean["b"]; c <- post_mean["c"]; k <- post_mean["k"]
# km <- survival::survfit(survival::Surv(time, status) ~ 1, data = remission)
# plot(km, conf.int = FALSE, xlab = "Time (months)",
#      ylab = "Survival probability",
#      main = "Posterior mean ED fit on remission data")
# t_grid <- seq(0.1, max(remission$time), length.out = 200)
# S_post <- pbetadanish(t_grid, a = 1, b = b, c = c, k = k,
#                       lower.tail = FALSE)
# lines(t_grid, S_post, col = "red", lwd = 2)
# legend("topright",
#        legend = c("Kaplan-Meier", "Posterior-mean ED"),
#        col = c("black", "red"), lty = 1, lwd = c(1, 2), bty = "n")

## -----------------------------------------------------------------------------
# fit_mle <- fit_betadanish(survival::Surv(time, status) ~ 1,
#                           data = remission, submodel = TRUE,
#                           n_starts = 3)
# cbind(
#   MLE  = round(coef(fit_mle), 4),
#   Bayes_post_mean = round(post_mean, 4)
# )

