Fits a Bayesian mixture model using Markov Chain Monte Carlo (MCMC) methods

Description

This method estimates the posterior distribution of a Bayesian mixture model using Markov Chain Monte Carlo (MCMC) methods and calculates measures of this distribution. The mixture model may consist of normal components (with a fixed expectation of 0), exponential components and gamma components, which may be mirrored in order to model negative values.

Usage

bayesMixModel(z, normNull=c(), expNeg=c(), expPos=c(), gamNeg=c(), gamPos=c(), sdNormNullInit=c(), rateExpNegInit=c(), rateExpPosInit=c(), shapeGamNegInit=c(), scaleGamNegInit=c(), shapeGamPosInit=c(), scaleGamPosInit=c(), piInit, classificationsInit, dirichletParInit=1, shapeDir=1, scaleDir=1, weightsPrior="FDD", sdAlpha, shapeNorm0=c(), scaleNorm0=c(), shapeExpNeg0=c(), scaleExpNeg0=c(), shapeExpPos0=c(), scaleExpPos0=c(), shapeGamNegAlpha0=c(), shapeGamNegBeta0=c(), scaleGamNegAlpha0=c(), scaleGamNegBeta0=c(), shapeGamPosAlpha0=c(), shapeGamPosBeta0=c(), scaleGamPosAlpha0=c(), scaleGamPosBeta0=c(), itb, nmc, thin, average="mean",sdShape)

Arguments

Details

The convergence of Markov chains must be assessed prior to an interpretation of results. Inspection of trace plots via plotChains is therefore urgently recommended. Iterations during which one of the chains has not yet reached stationarity should not be taken into account for analysis and can be excluded by setting an appropriate burn-in value itb. Autocorrelation between subsequent chain values can be reduced by thinning the chain, setting an appropriate value for thin. To ensure a sufficient number of iterations for the chains after the burn-in, nmc should be increased when the thinning is increased. The standard deviations of the proposal distribution in a Metropolis-Hastings step should be tuned to achieve a medium-level acceptance rate (e.g., 0.3-0.7): A very low acceptance rate would cause a long running time of the algorithm, while a very high acceptance rate typically leads to autocorrelation between the values of the chain. Acceptance is documented for each iteration in the chains slot of objects of class MixModelBayes-class.

Value

An object of class MixModelBayes-class storing results, data, priors, initial values and information about convergence.

Author(s)

Martin Schaefer (martin.schaefer@udo.edu)

See Also

plotChains, MixModelBayes-class

Examples

set.seed(1000)
z <- c(rnorm(1000, 0, 0.5), rnorm(1000, 0, 1))
mm <- bayesMixModel(z, normNull=1:2, sdNormNullInit=c(0.1, 0.2),
  piInit=c(1/2, 1/2), shapeNorm0=c(1, 1), scaleNorm0=c(1, 1),
  shapeExpNeg0=c(), scaleExpNeg0=c(), 
  shapeExpPos0=c(), scaleExpPos0=c(), sdAlpha=1, itb=100, nmc=1000, thin=10)
mm
plotComponents(mm)
plotChains(mm, chain="pi")

z <- c(rnorm(200, 0, 1), rnorm(200, 0, 5), rexp(200, 0.1), -rexp(200, 0.1))
mm <- bayesMixModel(z, normNull=1:2, gamNeg=3, gamPos=4,
  sdNormNullInit=c(1, 1),
  shapeGamNegInit=1, scaleGamNegInit=1, shapeGamPosInit=1, scaleGamPosInit=1,
  shapeNorm0=c(1,3), scaleNorm0=c(1,3), sdAlpha=1,
  shapeGamNegAlpha0=1, shapeGamNegBeta0=1,
  scaleGamNegAlpha0=1, scaleGamNegBeta0=1,
  shapeGamPosAlpha0=1, shapeGamPosBeta0=1,
  scaleGamPosAlpha0=1, scaleGamPosBeta0=1, sdShape=0.025, 
  itb=100, nmc=1000, thin=10)
mm
plotComponents(mm)
plotChains(mm, chain="pi")

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(epigenomix)
Loading required package: Biobase
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: S4Vectors
Loading required package: stats4

Attaching package: 'S4Vectors'

The following objects are masked from 'package:base':

    colMeans, colSums, expand.grid, rowMeans, rowSums

Loading required package: IRanges
Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: SummarizedExperiment
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/epigenomix/bayesMixModel.Rd_%03d_medium.png", width=480, height=480)
> ### Name: bayesMixModel
> ### Title: Fits a Bayesian mixture model using Markov Chain Monte Carlo
> ###   (MCMC) methods
> ### Aliases: bayesMixModel bayesMixModel,numeric-method
> 
> ### ** Examples
> 
> set.seed(1000)
> z <- c(rnorm(1000, 0, 0.5), rnorm(1000, 0, 1))
> mm <- bayesMixModel(z, normNull=1:2, sdNormNullInit=c(0.1, 0.2),
+   piInit=c(1/2, 1/2), shapeNorm0=c(1, 1), scaleNorm0=c(1, 1),
+   shapeExpNeg0=c(), scaleExpNeg0=c(), 
+   shapeExpPos0=c(), scaleExpPos0=c(), sdAlpha=1, itb=100, nmc=1000, thin=10)
[1] "Iteration: 50"
[1] "Iteration: 100"
[1] "Iteration: 150"
[1] "Iteration: 200"
[1] "Iteration: 250"
[1] "Iteration: 300"
[1] "Iteration: 350"
[1] "Iteration: 400"
[1] "Iteration: 450"
[1] "Iteration: 500"
[1] "Iteration: 550"
[1] "Iteration: 600"
[1] "Iteration: 650"
[1] "Iteration: 700"
[1] "Iteration: 750"
[1] "Iteration: 800"
[1] "Iteration: 850"
[1] "Iteration: 900"
[1] "Iteration: 950"
[1] "Iteration: 1000"
[1] "Iteration: 1050"
[1] "Iteration: 1100"
> mm

MixModel object
    Number of data points:  2000
    Number of components:   2
        1: NormNull
             mean = 0
             sd = 0.5501342
           weight pi = 0.6329708
           classified data points: 1666
        2: NormNull
             mean = 0
             sd = 1.081251
           weight pi = 0.3670292
           classified data points: 334

> plotComponents(mm)
> plotChains(mm, chain="pi")
> 
> z <- c(rnorm(200, 0, 1), rnorm(200, 0, 5), rexp(200, 0.1), -rexp(200, 0.1))
> mm <- bayesMixModel(z, normNull=1:2, gamNeg=3, gamPos=4,
+   sdNormNullInit=c(1, 1),
+   shapeGamNegInit=1, scaleGamNegInit=1, shapeGamPosInit=1, scaleGamPosInit=1,
+   shapeNorm0=c(1,3), scaleNorm0=c(1,3), sdAlpha=1,
+   shapeGamNegAlpha0=1, shapeGamNegBeta0=1,
+   scaleGamNegAlpha0=1, scaleGamNegBeta0=1,
+   shapeGamPosAlpha0=1, shapeGamPosBeta0=1,
+   scaleGamPosAlpha0=1, scaleGamPosBeta0=1, sdShape=0.025, 
+   itb=100, nmc=1000, thin=10)
[1] "Iteration: 50"
[1] "Iteration: 100"
[1] "Iteration: 150"
[1] "Iteration: 200"
[1] "Iteration: 250"
[1] "Iteration: 300"
[1] "Iteration: 350"
[1] "Iteration: 400"
[1] "Iteration: 450"
[1] "Iteration: 500"
[1] "Iteration: 550"
[1] "Iteration: 600"
[1] "Iteration: 650"
[1] "Iteration: 700"
[1] "Iteration: 750"
[1] "Iteration: 800"
[1] "Iteration: 850"
[1] "Iteration: 900"
[1] "Iteration: 950"
[1] "Iteration: 1000"
[1] "Iteration: 1050"
[1] "Iteration: 1100"
> mm

MixModel object
    Number of data points:  800
    Number of components:   4
        1: GamNeg
             shape = 0.8763034
             scale = 10.08081
           weight pi = 0.3723671
           classified data points: 266
        2: NormNull
             mean = 0
             sd = 1.885232
           weight pi = 0.001049356
           classified data points: 4
        3: NormNull
             mean = 0
             sd = 0.7871551
           weight pi = 0.2194311
           classified data points: 224
        4: GamPos
             shape = 1.142207
             scale = 6.474989
           weight pi = 0.4071524
           classified data points: 306

> plotComponents(mm)
> plotChains(mm, chain="pi")
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>