gibbs {CNPBayes} | R Documentation |
Model types: SB (SingleBatchModel): hierarchical model with mixture component-specific means and variances; MB (MultiBatchModel): hierarchical model with mixture component- and batch-specific means and variances; SBP (SingleBatchPooled): similar to SB model but with a pooled MBP (MultiBatchPooled): similar to MB model but with a pooled estimate of the variance (across mixture components) for each batch.
gibbs(model = c("SB", "MB", "SBP", "MBP"), dat, mp, hp.list, batches, k_range = c(1, 4), max_burnin = 32000, top = 2, df = 100, min_effsize = 500)
model |
a character vector indicating which models to fit (any combination of 'SB', 'MB', 'SBP', and 'MBP') |
dat |
numeric vector of the summary copy number data for each sample at a single CNP (e.g., the median log R ratio for each sample) |
mp |
an object of class |
hp.list |
a list of hyperparameters for each of the different models. If missing, this list will be generated automatically with default hyperparameters that work well for copy number data |
batches |
an integer vector of the same length as |
k_range |
a length-two numeric vector providing the minimum and maximum number of components to model. For example, c(1, 3) will fit mixture models with 1, 2, and 3 components. |
max_burnin |
the maximum number of burnin iterations. See details. |
top |
a length-one numeric vector indicating how many of the top models to return. |
df |
length-1 numeric vector for t-distribution degrees of freedom |
min_effsize |
length-1 numeric vector specifying the minimum effective size of the MCMC simulations. If below this value, the marginal likelihood will not be estimated. |
For each model specified, a Gibbs sampler will be initiated
The number of mixture models fit depends on k_range
and
A list of models of length top
sorted by decreasing
gelman.diag
effectiveSize
marginalLikelihood
See ggChains
set.seed(100) nbatch <- 3 k <- 3 means <- matrix(c(-2.1, -2, -1.95, -0.41, -0.4, -0.395, -0.1, 0, 0.05), nbatch, k, byrow = FALSE) sds <- matrix(0.15, nbatch, k) sds[, 1] <- 0.3 N <- 1000 truth <- simulateBatchData(N = N, batch = rep(letters[1:3], length.out = N), p = c(1/10, 1/5, 1 - 0.1 - 0.2), theta = means, sds = sds) ## clearly not enough simulations mp <- McmcParams(iter=10, burnin=5) ## this should be much higher > 1,000 max_burnin <- 10 ## suppress the warnings from the short chains suppressWarnings(gibbs(model=c("SB", "SBP", "MB", "MBP"), dat=y(truth), mp=mp, batches=batch(truth), k_range=c(1, 4), max_burnin=max_burnin, top=2, df=100))