View source: R/curepwe_logml_map.R
| curepwe.logml.map | R Documentation |
Uses Markov chain Monte Carlo (MCMC) and bridge sampling to estimate the logarithm of the marginal likelihood of a CurePWE model under the meta-analytic predictive (MAP) prior. The MAP prior is equivalent to the prior induced by the Bayesian hierarchical model (BHM).
The arguments related to MCMC sampling are utilized to draw samples from the MAP prior. These samples are then used to compute the logarithm of the normalizing constant of the MAP prior using only historical data set.
curepwe.logml.map(
post.samples,
bridge.args = NULL,
iter_warmup = 1000,
iter_sampling = 1000,
chains = 4,
...
)
post.samples |
output from |
bridge.args |
a |
iter_warmup |
number of warmup iterations to run per chain. Defaults to 1000. See the argument |
iter_sampling |
number of post-warmup iterations to run per chain. Defaults to 1000. See the argument |
chains |
number of Markov chains to run. Defaults to 4. See the argument |
... |
arguments passed to |
The function returns a list with the following objects
"curepwe_bhm"
the estimated logarithm of the marginal likelihood
an object of class bridge or bridge_list containing the output from using bridgesampling::bridge_sampler()
to compute the logarithm of the normalizing constant of the Bayesian hierarchical model (BHM) using all data sets
an object of class bridge or bridge_list containing the output from using
bridgesampling::bridge_sampler() to compute the logarithm of the normalizing constant of the BHM using historical
data set
the minimum estimated bulk effective sample size of the MCMC sampling
the maximum Rhat
Gronau, Q. F., Singmann, H., and Wagenmakers, E.-J. (2020). bridgesampling: An r package for estimating normalizing constants. Journal of Statistical Software, 92(10).
if (instantiate::stan_cmdstan_exists()) {
if(requireNamespace("survival")){
library(survival)
data(E1684)
data(E1690)
## take subset for speed purposes
E1684 = E1684[1:100, ]
E1690 = E1690[1:50, ]
## replace 0 failure times with 0.50 days
E1684$failtime[E1684$failtime == 0] = 0.50/365.25
E1690$failtime[E1690$failtime == 0] = 0.50/365.25
E1684$cage = as.numeric(scale(E1684$age))
E1690$cage = as.numeric(scale(E1690$age))
data_list = list(currdata = E1690, histdata = E1684)
nbreaks = 3
probs = 1:nbreaks / nbreaks
breaks = as.numeric(
quantile(E1690[E1690$failcens==1, ]$failtime, probs = probs)
)
breaks = c(0, breaks)
breaks[length(breaks)] = max(10000, 1000 * breaks[length(breaks)])
d.bhm = curepwe.bhm(
formula = survival::Surv(failtime, failcens) ~ treatment + sex + cage + node_bin,
data.list = data_list,
breaks = breaks,
chains = 1, iter_warmup = 500, iter_sampling = 1000
)
curepwe.logml.map(
post.samples = d.bhm,
bridge.args = list(silent = TRUE),
chains = 1, iter_warmup = 1000, iter_sampling = 2000
)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.