<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/using-mcmc-output-to-efficiently-estimate-bayes-factors/"></div></div><div><a data-url="http://www.r-bloggers.com/using-mcmc-output-to-efficiently-estimate-bayes-factors/" data-counturl="http://www.r-bloggers.com/using-mcmc-output-to-efficiently-estimate-bayes-factors/" data-count="horizontal" data-text="Using MCMC output to efficiently estimate Bayes factors" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://xianblog.wordpress.com/2016/05/19/using-mcmc-output-to-efficiently-estimate-bayes-factors/"> R – Xi'an's Og</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><strong>A</strong>s I was checking for software to answer <a target="_blank" href="http://stats.stackexchange.com/a/210196/7224">a query on X validated</a> about generic Bayes factor derivation, I came across an R software called <a target="_blank" href="https://cran.r-project.org/package=BayesFactor">BayesFactor</a>, which only applies in regression settings and relies on the <a target="_blank" href="https://xianblog.wordpress.com/2009/11/16/on-resolving-the-savage-dickey-paradox/">Savage-Dickey representation</a> of the Bayes factor</p>
<p><img title="B_{01}=\dfrac{f(y|\theta^0)}{m(y)}=\dfrac{\pi(\theta^0|y)}{\pi(\theta^0)}" alt="B_{01}=\dfrac{f(y|\theta^0)}{m(y)}=\dfrac{\pi(\theta^0|y)}{\pi(\theta^0)}" class="latex" src="https://s0.wp.com/latex.php?latex=B_%7B01%7D%3D%5Cdfrac%7Bf%28y%7C%5Ctheta%5E0%29%7D%7Bm%28y%29%7D%3D%5Cdfrac%7B%5Cpi%28%5Ctheta%5E0%7Cy%29%7D%7B%5Cpi%28%5Ctheta%5E0%29%7D&bg=000000&%23038;fg=B0B0B0&%23038;s=0"></p>
<p>when the null hypothesis writes as θ=θ⁰ (and possibly additional nuisance parameters with [roughly speaking] an independent prior). As we discussed in <a target="_blank" href="https://xianblog.wordpress.com/2010/07/12/savage-dickey-published/">our paper with Jean-Michel Marin</a> [which got ignored by large!], this representation of the Bayes factor is based on picking a very specific version of the prior, or more exactly of three prior densities. Assuming such versions are selected, I wonder at the performances of this approximation, given that it involves approximating the marginal posterior at θ⁰….</p>
<blockquote>
<p><em><span>“To ensure that the Bayes factor we compute using the Savage–Dickey ratio is the the ratio of marginal densities that we intend, the condition (…) is easily met by models which specify priors in which the nuisance parameters are independent of the parameters of interest.” Morey et al. (2011)</span></em></p>
</blockquote>
<p>First, when reading <a target="_blank" href="http://www.sciencedirect.com/science/article/pii/S0022249611000666">Morey at al. (2011)</a>, I realised (a wee bit late!) that Chib’s method is nothing but a version of the Savage-Dickey representation when the marginal posterior can be estimated in a parametric (Rao-Blackwellised) way. However, outside hierarchical models based on conjugate priors such parametric approximations are intractable and non-parametric versions must be invoked instead, which necessarily degrades the quality of the method. A degradation that escalates with the dimension of the parameter θ. In addition, I am somewhat perplexed by the use of a Rao-Blackwell argument in the setting of the Dickey-Savage representation. Indeed this representation assumes that</p>
<p><img title="\pi_1(\psi|theta_0)=\pi_0(\ps) \quad\text{or}\quad \pi_1(\theta_0,\psi)=\pi_1(\theta_0)\pi_0(\psi)" alt="\pi_1(\psi|theta_0)=\pi_0(\ps) \quad\text{or}\quad \pi_1(\theta_0,\psi)=\pi_1(\theta_0)\pi_0(\psi)" class="latex" src="https://s0.wp.com/latex.php?latex=%5Cpi_1%28%5Cpsi%7Ctheta_0%29%3D%5Cpi_0%28%5Cps%29+%5Cquad%5Ctext%7Bor%7D%5Cquad+%5Cpi_1%28%5Ctheta_0%2C%5Cpsi%29%3D%5Cpi_1%28%5Ctheta_0%29%5Cpi_0%28%5Cpsi%29&bg=000000&%23038;fg=B0B0B0&%23038;s=0"></p>
<p>which means that [the specific version of] the conditional density of θ⁰ given ψ should not depend on the nuisance parameter. But relying on a Rao-Blackwellisation leads to estimate the marginal posterior via full conditionals. Of course, θ given ψ and y may depend on ψ, but still… <a target="_blank" href="http://www.sciencedirect.com/science/article/pii/S0022249611000666">Morey at al. (2011)</a> advocate the recourse to Chib’s formula as optimal but this obviously requires the full conditional to be available. They acknowledge this point as moot, since it is sufficient from their perspective to specify a conjugate prior. They consider this to be a slight modification of the model (p.377). However, I see the evaluation of an estimated density at a <em>single</em> (I repeat, <em>single</em>!) point as being the direst part of the method as it is clearly more sensitive to approximations that the evaluation of a whole integral, since the later incorporates an averaging effect by definition. Hence, even if this method was truly available for all models, I would be uncertain of its worth when compared with other methods, except the harmonic mean estimator of course!</p>
<p>On the side, <a target="_blank" href="http://www.sciencedirect.com/science/article/pii/S0022249611000666">Morey at al. (2011)</a> study a simple one-sample t test where they use an improper prior on the nuisance parameter σ, under both models. While the Savage-Dickey representation is correct in this special case, I fail to see why the identity would apply in every case under an improper prior. In particular, independence does not make sense with improper priors. The authors also indicate the possible use of this Bayes factor approximation for encompassing models. At first, I thought this could be most useful in our testing by mixture framework where we define an encompassing model as a mixture. However, I quickly realised that using a Beta Be(a,a) prior on the weight α with a&lt;1 leads to an infinite density value at both zero and one, hence cannot be compatible with a Savage-Dickey representation of the Bayes factor.</p>
<p>Filed under: <a target="_blank" href="https://xianblog.wordpress.com/category/books/">Books</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/statistics/r-statistics/">R</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/statistics/">Statistics</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/university-life/">University life</a> Tagged: <a target="_blank" href="https://xianblog.wordpress.com/tag/bayes-factors/">Bayes factors</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/bayesfactor/">BayesFactor</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/chibs-approximation/">Chib’s approximation</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/dickey-savage-ratio/">Dickey-Savage ratio</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/r/">R</a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/xianblog.wordpress.com/30587/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/xianblog.wordpress.com/30587/" border="0"></a> <img height="1" alt="" width="1" src="https://pixel.wp.com/b.gif?host=xianblog.wordpress.com&blog=5051449&%23038;post=30587&%23038;subd=xianblog&%23038;ref=&%23038;feed=1" border="0"></p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/using-mcmc-output-to-efficiently-estimate-bayes-factors/"></div></div><div><a data-url="http://www.r-bloggers.com/using-mcmc-output-to-efficiently-estimate-bayes-factors/" data-counturl="http://www.r-bloggers.com/using-mcmc-output-to-efficiently-estimate-bayes-factors/" data-count="vertical" data-text="Using MCMC output to efficiently estimate Bayes factors" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://xianblog.wordpress.com/2016/05/19/using-mcmc-output-to-efficiently-estimate-bayes-factors/"> R – Xi'an's Og</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/fieuDZ6VUbA"></a></div>