Runs a Metropolis Coupled MCMC (MC\(^3\)) sampler in order to estimate the joint posterior distribution of the model.
cure_rate_MC3(formula, data, nChains = 12, mcmc_cycles = 15000,
alpha = NULL,nCores = 1, sweep = 5, a_g = 1, b_g = 1,
a_l = 2.1, b_l = 1.1, mu_b = NULL,
Sigma = NULL, g_prop_sd = 0.045,
lambda_prop_scale = 0.03, b_prop_sd = NULL,
initialValues = NULL, plot = TRUE, adjust_scales = FALSE,
verbose = FALSE, tau_mala = 1.5e-05, mala = 0.15,
promotion_time = list(family = "weibull",
prior_parameters = matrix(rep(c(2.1, 1.1), 2), byrow = TRUE, 2, 2),
prop_scale = c(0.1, 0.2)), single_MH_in_f = 0.2, c_under = 1e-9)An object of class bayesCureModel, i.e. a list with the following entries
Object of class mcmc (see the coda package), containing the generated MCMC sample for the target chain. The column names correspond to
g_mcmc Sampled values for parameter \(\gamma\)
lambda_mcmc Sampled values for parameter \(\lambda\)
alpha1_mcmc... Sampled values for parameter \(\alpha_1\)... of the promotion time distribution \(F(\cdot;\alpha_1,\ldots,\alpha_d)\). The subsequent \(d-1\) columns contain the sampled values for all remaining parameters, \(\alpha_2,...,\ldots,\alpha_d\), where \(d\) depens on the family used in promotion_time.
b0_mcmc Sampled values for the constant term of the regression (present only in the case where the design matrix X contains a column of 1s).
b1_mcmc... Sampled values for the regression coefficient for the first explanatory variable, and similar for all subsequent columns.
A data frame with the simulated binary latent status for each censored item.
The complete log-likelihood for the target chain.
the acceptance rate of proposed swappings between adjacent MCMC chains.
The complete log-likelihood for all chains.
the input data, model specification and selected prior parameters values.
the logarithm of the (non-augmented) posterior distribution (after integrating the latent cured-status parameters out), up to a normalizing constant.
The Maximum A Posterior estimate of parameters.
Bayesian Information Criterion of the fitted model.
Akaike Information Criterion of the fitted model.
The Cox-Snell residuals of the fitted model.
The maximum log-likelihood value for the target chain.
The number of parameters of the model.
The initial values per chain.
an object of class formula: a symbolic description of the model to be fitted. The left-hand side should be a Surv object, a class inherited from the survival package. The binary censoring indicators are interpreted as a time-to-event (1) or as a censoring time (0).
a data frame containing all variable names included in formula.
Positive integer corresponding to the number of heated chains in the MC\(^3\) scheme.
Length of the generated MCMC sample. Default value: 15000. Note that each MCMC cycle consists of sweep (see below) usual MCMC iterations.
A decreasing sequence in \([1,0)\) of nChains temperatures (or heat values). The first value should always be equal to 1, which corresponds to the target posterior distribution (that is, the first chain).
The number of cores used for parallel processing. In case where nCores > 1, the nChains will be processed in parallel using nCores cores. This may speed up significantly the run-time of the algorithm in Linux or macOS, but it is not suggested in Windows.
The number of usual MCMC iterations per MCMC cycle. Default value: 10.
Parameter \(a_{\gamma}\) of the prior distribution of \(\gamma\).
Parameter \(b_{\gamma}\) of the prior distribution of \(\gamma\).
Shape parameter \(a_{\lambda}\) of the Inverse Gamma prior distribution of \(\lambda\).
Scale parameter \(b_{\lambda}\) of the Inverse Gamma prior distribution of \(\lambda\).
Mean (\(\boldsymbol\mu\)) of the multivariate normal prior distribution of regression coefficients. Should be a vector whose length is equal to \(k\), i.e. the number of columns of the design matrix X. Default value: rep(0, k).
Covariance matrix of the multivariate normal prior distribution of regression coefficients.
The scale of the proposal distribution for single-site updates of the \(\gamma\) parameter.
The scale of the proposal distribution for single-site updates of the \(\lambda\) parameter.
The scale of the proposal distribution for the update of the \(\beta\) parameter (regression coefficients).
A list of initial values for each parameter (optional).
Plot MCMC sample on the run. Default: TRUE.
Boolean. If TRUE the MCMC sampler runs an initial phase of a small number of iterations in order to tune the scale of the proposal distributions in the Metropolis-Hastings steps.
Print progress on the terminal if TRUE.
Scale of the Metropolis adjusted Langevin diffussion proposal distribution.
A number between \([0,1]\) indicating the proportion of times the sampler attempts a MALA proposal. Thus, the probability of attempting a typical Metropolis-Hastings move is equal to 1 - mala.
A list with details indicating the parametric family of distribution describing the promotion time and corresponding prior distributions. See `details`.
The probability for attempting a series of single site updates in the typical Metropolis-Hastings move. Otherwise, with probability 1 - single_MH_in_f a Metropolis-Hastings move will attempt to update all continuous parameters simultaneously. It only makes sense when mala < 1.
A small positive number (much smaller than 1) which is used as a threshold in the CDF of the promotion time for avoiding underflows in the computation of the log-likelihood function. Default value: 1e-9.
Panagiotis Papastamoulis
It is advised to scale all continuous explanatory variables in the design matrix, so their sample mean and standard deviations are equal to 0 and 1, respectively. No missing or infinite values are accepted.
The promotion_time argument should be a list containing the following entries
family Character string specifying the family of distributions \(\{F(\cdot;\boldsymbol\alpha);\boldsymbol\alpha\in\mathcal A\}\) describing the promotion time.
prior_parameters Values of hyper-parameters in the prior distribution of the parameters \(\boldsymbol\alpha\).
prop_scale The scale of the proposal distributions for each parameter in \(\boldsymbol\alpha\).
K Optional. The number of mixture components in case where a mixture model is fitted, that is, when setting distribution to either 'gamma_mixture' or 'user_mixture'.
dirichlet_concentration_parameter Optional. Relevant only in the case of the 'gamma_mixture' or 'user_mixture'. Positive scalar (typically, set to 1) determining the (common) concentration parameter of the Dirichlet prior distribution of mixing proportions.
The family specifies the distributional family of promotion time and corresponds to a character string with available choices: 'exponential', 'weibull', 'gamma', 'logLogistic', 'gompertz', 'lomax', 'dagum', 'gamma_mixture'. User defined promotion time distributions and finite mixtures of them can be also fitted using the options family = 'user' and family = 'user_mixture', respectively. See the 'Note' below for further details.
The joint prior distribution of \(\boldsymbol\alpha = (\alpha_1,\ldots,\alpha_d)\) factorizes into products of inverse Gamma distributions for all (positive) parameters of \(F\). Moreover, in the case of 'gamma_mixture' or 'user_mixture', the joint prior also consists of another term to the Dirichlet prior distribution on the mixing proportions.
The prop_scale argument should be a vector with length equal to the length of vector \(d\) (number of elements in \(\boldsymbol\alpha\)), containing (positive) numbers which correspond to the scale of the proposal distribution. Note that these scale parameters are used only as initial values in case where adjust_scales = TRUE.
Papastamoulis and Milienos (2024). Bayesian inference and cure rate modeling for event history data. TEST doi: 10.1007/s11749-024-00942-w.
Papastamoulis and Milienos (2024b). bayesCureRateModel: Bayesian Cure Rate Modeling for Time to Event Data in R. arXiv:2409.10221
Plummer M, Best N, Cowles K, Vines K (2006). "CODA: Convergence Diagnosis and Output Analysis for MCMC." R News, 6(1), 7-11.
Therneau T (2024). A Package for Survival Analysis in R. R package version 3.7-0, https://CRAN.R-project.org/package=survival.
cure_rate_mcmc
# simulate toy data just for cran-check purposes
set.seed(10)
n = 4
# censoring indicators
stat = rbinom(n, size = 1, prob = 0.5)
# covariates
x <- matrix(rnorm(2*n), n, 2)
# observed response variable
y <- rexp(n)
# define a data frame with the response and the covariates
my_data_frame <- data.frame(y, stat, x1 = x[,1], x2 = x[,2])
fit1 <- cure_rate_MC3(survival::Surv(y, stat) ~ x1 + x2,
data = my_data_frame,
promotion_time = list(family = 'weibull'),
nChains = 2, nCores = 1,
mcmc_cycles = 3, sweep = 2)
Run the code above in your browser using DataLab