1

Has anyone managed to obtain a Monte Carlo error for a parameter when running bayesian model un R2OpenBugs?

It is provided in a standard output of OpenBugs, but when run under R2OpenBugs, the log file doesn't have MC error.Is there a way to ask R2OpenBugs to calculate MC error? Or maybe there is a way to calculate it manually? Please, let me know if you heard of any way to do that. Thank you!

Here is the standard log output of R2OpenBugs:

$stats
              mean      sd val2.5pc    median val97.5pc sample
beta0      1.04700 0.13250   0.8130   1.03800   1.30500   1500
beta1     -0.31440 0.18850  -0.6776  -0.31890   0.03473   1500
beta2     -0.05437 0.05369  -0.1648  -0.05408   0.04838   1500
deviance 588.70000 7.87600 575.3000 587.50000 606.90000   1500

$DIC
       Dbar  Dhat   DIC    pD
t     588.7 570.9 606.5 17.78
total 588.7 570.9 606.5 17.78
marina_esp
  • 13
  • 3

2 Answers2

0

A simple way to calculate Monte Carlo standard error (MCSE) is to divide the standard deviation of the chain by the square root of the effective number of samples. The standard deviation is provided in your output, but the effective sample size should be given as n.eff (the rightmost column) when you print the model output - or at least that is the impression I get from:

https://cran.r-project.org/web/packages/R2OpenBUGS/vignettes/R2OpenBUGS.pdf

I don't use OpenBugs any more so can't easily check for you, but there should be something there that indicates the effective sample size (this is NOT the same as the number of iterations you have sampled, as it also takes into account the loss of information due to correlation within the chains).

Otherwise you can obtain it yourself by extracting the raw MCMC chains and then either computing the effective sample size using the coda package (?coda::effectiveSize) or just use LaplacesDemon::MCSE to calculate the Monte Carlo standard error directly. For more information see:

https://rdrr.io/cran/LaplacesDemon/man/MCSE.html

Note that some people (including me!) would suggest focusing on the effective sample size directly rather than looking at the MCSE, as the old "rule of thumb" that MCSE should be less than 5% of the sample standard deviation is equivalent to saying that the effective sample size should be at least 400 (1/0.05^2). But opinions do vary :)

Matt Denwood
  • 2,537
  • 1
  • 12
  • 16
  • 1
    Hey Matt, many thanks to your answer - indeed, I found the effective number of parameters in the output of the R2OpenBugs runs. I saw also that coda can provide effective sample size and time-series error that I can use instead of MC error, so problem solved! Thanks! – marina_esp Jul 17 '22 at 11:24
0

The MCMC-error is named Time-series SE, and can be found in the statistics section of the summary of the coda object:

library(R2OpenBUGS)
library(coda)
my_result <- bugs(...., codaPg = TRUE)
my_coda <- read.bugs(my_result)
summary(my_coda$statistics)
MartineJ
  • 591
  • 5
  • 16