gibbs: Run a Gibbs sampler on one or multiple types of Bayesian...
In CNPBayes: Bayesian mixture models for copy number polymorphisms

Description Usage Arguments Details Value See Also Examples

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.

1
2
3

gibbs(model = c("SB", "MB", "SBP", "MBP", "TBM"), dat, mp, hp.list,
  batches, k_range = c(1, 4), max_burnin = 32000, min_GR = 1.2,
  top = 2, df = 100, min_effsize = 500, maplabel, mprob)

`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 `McmcParams`
`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 `dat` indicating the batch in which the sample was processed
`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.
`min_GR`	length-1 numeric vector specifying the maximum Gelman-Rubin (GR) statistic. If the GR statistic is above this value, the marginal likelihood will not be estimated.
`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,
        min_GR=1.2,
        top=2,
        df=100))