Gibbs sampler

Description

Gibbs Sampler Algorithm to fit a three component normal mixture to z-scores

Usage

bacon(teststatistics = NULL, effectsizes = NULL, standarderrors = NULL,
  niter = 5000L, nburnin = 2000L, nbins = 1000, trim = 0.999,
  level = 0.05, verbose = FALSE, priors = list(sigma = list(alpha = 1.28,
  beta = 0.36), mu = list(lambda = c(0, 3, -3), tau = c(1000, 100, 100)),
  epsilon = list(gamma = c(90, 5, 5))))

Arguments

Value

object of class-Bacon

Author(s)

mvaniterson

References

Implementation is based on a version from Zhihui Liu https://macsphere.mcmaster.ca/handle/11375/9368

Examples

##simulate some test-statistic from a normal mixture
##and run bacon
y <- rnormmix(2000, c(0.9, 0, 1, 0, 4, 1))
bc <- bacon(y)
##extract all estimated mixture parameters
estimates(bc)
##extract inflation
inflation(bc)
##extract bias
bias(bc)

##extract bias and inflation corrected test-statistics
head(tstat(bc))

##inspect the Gibbs Sampling output
traces(bc)
posteriors(bc)
fit(bc)

##simulate multiple sets of test-statistic from a normal mixture
##and run bacon
y <- matrix(rnormmix(10*2000, c(0.9, 0, 1, 0, 4, 1)), ncol=10)
bc <- bacon(y)
##extract all estimated mixture parameters
estimates(bc)
##extract only the inflation
inflation(bc)
##extract only the bias
bias(bc)
##extract bias and inflation corrected P-values
head(pval(bc))
##extract bias and inflation corrected test-statistics
head(tstat(bc))

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(bacon)
Loading required package: ggplot2
Loading required package: BiocParallel
Loading required package: ellipse
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/bacon/bacon.Rd_%03d_medium.png", width=480, height=480)
> ### Name: bacon
> ### Title: Gibbs sampler
> ### Aliases: bacon
> 
> ### ** Examples
> 
> ##simulate some test-statistic from a normal mixture
> ##and run bacon
> y <- rnormmix(2000, c(0.9, 0, 1, 0, 4, 1))
> bc <- bacon(y)
> ##extract all estimated mixture parameters
> estimates(bc)
           p.0        p.1        p.2        mu.0     mu.1      mu.2  sigma.0
[1,] 0.9149751 0.04540628 0.03961866 0.008803744 3.072176 -2.980066 0.973314
      sigma.1  sigma.2
[1,] 3.215098 2.422526
> ##extract inflation
> inflation(bc)
 sigma.0 
0.973314 
> ##extract bias
> bias(bc)
       mu.0 
0.008803744 
> 
> ##extract bias and inflation corrected test-statistics
> head(tstat(bc))
           [,1]
[1,]  3.3547193
[2,]  0.2891248
[3,] -1.2378830
[4,]  0.6932090
[5,] -0.5214156
[6,] -2.7739474
> 
> ##inspect the Gibbs Sampling output
> traces(bc)
> posteriors(bc)
> fit(bc)
> 
> ##simulate multiple sets of test-statistic from a normal mixture
> ##and run bacon
> y <- matrix(rnormmix(10*2000, c(0.9, 0, 1, 0, 4, 1)), ncol=10)
> bc <- bacon(y)
Detected 2 workers!
 Running in parallel!
> ##extract all estimated mixture parameters
> estimates(bc)
            p.0        p.1        p.2          mu.0     mu.1      mu.2
 [1,] 0.9087300 0.05157525 0.03969479 -0.0125123187 3.003531 -3.016948
 [2,] 0.9247407 0.02812204 0.04713722  0.0073809920 2.999192 -2.983687
 [3,] 0.9178518 0.04106253 0.04108569 -0.0103019001 2.979784 -3.011203
 [4,] 0.9154095 0.03733342 0.04725705 -0.0056951483 2.998013 -3.025140
 [5,] 0.9040199 0.04641846 0.04956166 -0.0040141435 2.998946 -3.050163
 [6,] 0.9265751 0.03584039 0.03758449 -0.0008153726 2.939128 -2.908333
 [7,] 0.8925215 0.05671603 0.05076246  0.0146175373 2.978800 -2.937664
 [8,] 0.9162281 0.05184305 0.03192882 -0.0064406061 3.094536 -2.984730
 [9,] 0.9063187 0.05008929 0.04359197 -0.0065676185 2.964585 -2.963472
[10,] 0.9081254 0.04529965 0.04657499 -0.0119304411 3.097835 -2.949008
        sigma.0  sigma.1  sigma.2
 [1,] 1.0297355 3.646010 2.227626
 [2,] 1.0272668 1.725541 2.815489
 [3,] 1.0234439 2.652566 2.248209
 [4,] 1.0197880 2.659098 2.463979
 [5,] 0.9412091 2.040122 2.335229
 [6,] 0.9944553 2.372087 2.726762
 [7,] 0.9628074 3.157276 2.767606
 [8,] 1.0335693 2.845803 2.272625
 [9,] 1.0266529 2.871041 2.521263
[10,] 0.9848894 2.589971 2.842602
> ##extract only the inflation
> inflation(bc)
 [1] 1.0297355 1.0272668 1.0234439 1.0197880 0.9412091 0.9944553 0.9628074
 [8] 1.0335693 1.0266529 0.9848894
> ##extract only the bias
> bias(bc)
 [1] -0.0125123187  0.0073809920 -0.0103019001 -0.0056951483 -0.0040141435
 [6] -0.0008153726  0.0146175373 -0.0064406061 -0.0065676185 -0.0119304411
> ##extract bias and inflation corrected P-values
> head(pval(bc))
           [,1]       [,2]         [,3]         [,4]         [,5]      [,6]
[1,] 0.05815205 0.35749362 2.395684e-17 7.103640e-01 0.0089762012 0.1171387
[2,] 0.50013233 0.82825475 2.335452e-01 6.136641e-05 0.9270022367 0.1037776
[3,] 0.10725455 0.44524292 1.755921e-01 9.118047e-01 0.5567070873 0.2786014
[4,] 0.37213817 0.45281238 2.677628e-01 9.542784e-01 0.7704204351 0.7515907
[5,] 0.90099644 0.07387008 3.606889e-01 6.578322e-01 0.0003132259 0.9407800
[6,] 0.21864743 0.15715430 7.421502e-01 4.958916e-01 0.3169277481 0.7387116
          [,7]       [,8]        [,9]     [,10]
[1,] 0.9247569 0.81388132 0.002816869 0.5855144
[2,] 0.3218910 0.38453697 0.444972775 0.5770296
[3,] 0.2835238 0.91621232 0.290975694 0.9595532
[4,] 0.1618951 0.24184238 0.930388611 0.2323650
[5,] 0.6164188 0.05245107 0.039363885 0.5992374
[6,] 0.4893227 0.70743328 0.992906454 0.4845341
> ##extract bias and inflation corrected test-statistics
> head(tstat(bc))
           [,1]       [,2]       [,3]        [,4]        [,5]       [,6]
[1,] -1.8945500 -0.9201514 -8.4727977  0.37136730  2.61295922  1.5668973
[2,] -0.6742815 -0.2169405 -1.1912761 -4.00749404 -0.09161714  1.6268092
[3,] -1.6106568 -0.7633697 -1.3544525  0.11076249  0.58773961 -1.0834666
[4,] -0.8924754  0.7507346 -1.1082293  0.05733488 -0.29182499  0.3165426
[5,]  0.1244027  1.7874172  0.9140532 -0.44290812 -3.60410873 -0.0742895
[6,] -1.2301327 -1.4147074  0.3290073 -0.68096821 -1.00079123  0.3335601
           [,7]       [,8]         [,9]       [,10]
[1,]  0.0944435 -0.2354219 -2.987046455  0.54534774
[2,]  0.9905795 -0.8695672 -0.763822921  0.55772876
[3,]  1.0724371  0.1052060 -1.055983629 -0.05071428
[4,] -1.3987262  1.1703942  0.087355914 -1.19428878
[5,]  0.5009322  1.9394130  2.060362650 -0.52549755
[6,] -0.6913865  0.3753055  0.008890559  0.69902862
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>