1. Overview
This vignette walks through Bayesian estimation of the
three-parameter Exponentiated Danish (ED) submodel with
bayes_betadanish(). The same call with
submodel = FALSE fits the full four-parameter Beta-Danish
distribution.
From version 1.0.0 the default sampler is the multi-chain adaptive
Metropolis algorithm used in the thesis on which the package is based.
It is written in base R, so no additional package is needed. Each chain
adapts the scale and then the shape of its proposal during the burn-in
and stops adapting before any draw is retained. The fit then reports the
potential scale reduction factor (R-hat) and the effective sample size
of every parameter, and flags the result as converged only when every
R-hat is at most 1.05 and every effective sample size is at least
400.
The single-chain sampler of earlier versions is still available as
method = "metrop1R", which requires the MCMCpack and coda
packages.
2. ED submodel on bladder cancer remission times
library(BetaDanish)
data("remission")
fit_bayes <- bayes_betadanish(
time = remission$time,
status = remission$status,
submodel = TRUE,
burnin = 2000, mcmc = 4000,
n_chains = 4,
seed = 1
)
print(fit_bayes)
#>
#> Beta-Danish Bayesian Fit (3-Parameter ED Submodel)
#> Adaptive Metropolis: 4 chain(s), 2000 burn-in + 4000 iterations each, thin 1
#>
#> parameter mean sd q2.5 q25 median q75 q97.5
#> b 5.772700 5.821500 2.2308000 3.414400 4.416100 6.002500 19.88100
#> c 1.529600 0.280560 1.0734000 1.328300 1.495700 1.701300 2.19610
#> k 0.069555 0.045676 0.0071347 0.037042 0.062239 0.089256 0.18831
#> HPD_lower HPD_upper
#> 1.569500 11.59100
#> 1.033400 2.07100
#> 0.001507 0.15827
#>
#> Convergence diagnostics (log scale):
#> parameter Rhat ESS acceptance_rate
#> b 1.0355 471.1 0.229
#> c 1.0128 621.0 0.229
#> k 1.0303 422.5 0.229
#>
#> Converged (R-hat <= 1.05, ESS >= 400): yes
The chains start at the maximum-likelihood estimate (the first chain)
and at random perturbations of it (the others), so agreement between
them is evidence that the sampler has forgotten where it started.
3. Posterior diagnostics
fit_bayes$diagnostics
#> parameter Rhat ESS acceptance_rate
#> 1 b 1.035524 471.0715 0.2287083
#> 2 c 1.012839 620.9782 0.2287083
#> 3 k 1.030278 422.5018 0.2287083
fit_bayes$converged
#> [1] TRUE
plot(fit_bayes)
Each chain is drawn in its own colour in the trace panels. Chains
that overlap and show no drift are what a converged sampler looks like;
the density panels pool all chains and mark the 95% highest posterior
density interval.
4. Posterior survival vs Kaplan-Meier
post_mean <- setNames(fit_bayes$summary$mean, fit_bayes$summary$parameter)
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")
5. Bayesian vs MLE comparison
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[names(coef(fit_mle))], 4)
)
#> MLE Bayes_post_mean
#> b 4.2304 5.7727
#> c 1.5256 1.5296
#> k 0.0648 0.0696
6. The four-parameter model
For the full model the posterior of the product a * c is
usually far better determined than those of a and
c separately, because near the lower tail the two enter the
density almost only through their product. The summary therefore carries
an extra row, a_times_c. The four-parameter posterior needs
longer chains than the ED one; check converged before using
it.
fit_bd <- bayes_betadanish(remission$time, remission$status,
submodel = FALSE, burnin = 20000, mcmc = 40000,
seed = 1)
fit_bd$summary
fit_bd$converged