# Bayesian Model Selection

**URL:** <https://discourse.julialang.org/t/bayesian-model-selection/55725>\
**Category:** Statistics\
**Tags:** turing\
**Created:** [February 21, 2021, 4:34pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725 "2021-02-21T16:34:11Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![ANasc](https://avatars.discourse-cdn.com/v4/letter/a/41988e/32.png) [@ANasc](https://discourse.julialang.org/u/ANasc)\
**Post date:** [February 21, 2021, 4:34pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/1 "2021-02-21T16:34:11Z")

</div>

Hi All,

This question might reveal my limited knowledge of Bayesian estimation, but I am having difficulty figuring out how to perform basic Bayesian model selection in Turing. The documentation for MCMCChains.jl has a short section on DIC, which I’m not very familiar with. How does one calculate Bayes factors, for example? I can’t even figure out how to get model Evidence from the output of sample(). (I’m using Metropolis-Hastings algorithm). Any tips or directions of where to look would be greatly appreciated.

Thanks

---

<div class="post-metadata">

**Author:** ![rikh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rikh/32/204104_2.png) [@rikh](https://discourse.julialang.org/u/rikh)\
**Post date:** [February 22, 2021, 9:15am UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/2 "2021-02-22T09:15:25Z")

</div>

> [@ANasc](#):
>
> How does one calculate Bayes factors, for example?

I don’t know much about the metrics for model selection, but do know that various statisticians argue against Bayes factors (Tendeiro & Kiers, [2019](https://doi.org/10.1037/met0000221); Gelman et al., [2020](https://doi.org/https://doi.org/10.1201/b16018)).

As far as I understand, McElreath ([2020](https://doi.org/10.1201/9780429029608)) and Gelman et al. ([2020](https://doi.org/https://doi.org/10.1201/b16018)) advise to do model comparison based on predictive performance measures such as cross-validation and AIC. McElreath ([2020](https://doi.org/10.1201/9780429029608)) argues that these metrics are quite similar.

But, I think that it also depends on where you want to publish. My co-authors aren’t comfortable with cross-validation, so I will stick to a multiverse analysis (Steegen et al., [2016](https://doi.org/10.1177%2F1745691616658637)) with a multilevel model to reduce overfitting and clear graphs to show what the model does.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [February 22, 2021, 9:51am UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/3 "2021-02-22T09:51:08Z")

</div>

As far as I can tell, there are two options currently, First, there is a Julia interface to [ArviZ](https://arviz-devs.github.io/ArviZ.jl/v0.2/api/), which is written in Python. Second, there is a [package](https://github.com/StatisticalRethinkingJulia/ParetoSmoothedImportanceSampling.jl) associated with the book Statistical Rethinking that implements WAIC and LOO. The downside to the second option is that it does not currently have diagnostics for the model selection criteria.

---

<div class="post-metadata">

**Author:** ![ANasc](https://avatars.discourse-cdn.com/v4/letter/a/41988e/32.png) [@ANasc](https://discourse.julialang.org/u/ANasc)\
**Post date:** [February 23, 2021, 5:48am UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/4 "2021-02-23T05:48:37Z")

</div>

Thank you both for the links, these are very useful. I will brush up on which approaches are most appropriate for my case and see if I can implement something.

I’ll also look more closely at DIC in MCMCChains. I can’t seem to get it to work. It gives me an error that I haven’t used the correct inputs for logpdf(), and I’m not sure exactly what needs to be in there (for both the distribution and x). I’m assuming it’s asking for the likelihood function (mine is BinomialLogit). Here is the function:

````nohighlight
lpfun ​=​ ​function​ ​f​(chain​::​Chains​) ​#​ function to compute the logpdf values​
    niter, nparams, nchains ​=​ ​size​(chain)
    lp ​=​ ​zeros​(niter ​+​ nchains) ​#​ resulting logpdf values​
    ​for​ i ​=​ ​1​:​nparams
        lp ​+=​ ​map​(p ​->​ ​logpdf​( ​...​ , x), ​Array​(chain[:,i,:]))
    ​end​
    ​return​ lp
​end​
DIC, pD ​=​ ​dic​(chn, lpfun)```
````

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [February 24, 2021, 6:42pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/5 "2021-02-24T18:42:39Z")

</div>

> [@ANasc](#):
>
> `lp ​+=​ ​map​(p ​->​ ​logpdf​( ​...​ , x), ​Array​(chain[:,i,:]))`

I think the problem is that you are supposed to add your distribution and data to the code. For example,

```julia
lpfun ​=​ ​function​ ​f​(chain​::​Chains​) ​#​ function to compute the logpdf values​
    niter, nparams, nchains ​=​ ​size​(chain)
    lp ​=​ ​zeros​(niter ​+​ nchains) ​#​ resulting logpdf values​
    ​for​ i ​=​ ​1​:​nparams
        lp ​+=​ ​map​(p ​->​ ​logpdf​(BinomialLogit(p...), your_data), ​Array​(chain[:,i,:]))
    ​end​
    ​return​ lp
​end​
DIC, pD ​=​ ​dic​(chn, lpfun)

```

However, I think this API might be difficult for users. It might be better to do something like this:

```julia
DIC, pD ​=​ ​dic​(chn, distribution_object, the_data)

```

Under my proposal, the function dic would create the `lpfun` function and perform the calculations as normal. This way the user does not need to modify functions.

@cpfiffer, would this be a better way to compute DIC or did I misunderstand the documentation?

---

<div class="post-metadata">

**Author:** ![cpfiffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cpfiffer/32/208747_2.png) [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)\
**Post date:** [February 24, 2021, 6:53pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/6 "2021-02-24T18:53:23Z")

</div>

Yeah, I think that’s probably a better way to do it. It seems like a pretty trivial thing to add that extension into MCMCChains – maybe a good first issue for somebody.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [February 24, 2021, 7:54pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/7 "2021-02-24T19:54:18Z")

</div>

I tried to change the API, but I encountered several issues with the `lpfun` function. It almost seems like it was never fully tested and implemented. I opened an [issue](https://github.com/TuringLang/MCMCChains.jl/issues/267).

---

<div class="post-metadata">

**Author:** ![ANasc](https://avatars.discourse-cdn.com/v4/letter/a/41988e/32.png) [@ANasc](https://discourse.julialang.org/u/ANasc)\
**Post date:** [February 25, 2021, 3:37pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/8 "2021-02-25T15:37:07Z")

</div>

@Christopher_Fisher Thanks for working on this, I really appreciate it.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [February 25, 2021, 7:03pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/9 "2021-02-25T19:03:02Z")

</div>

`No problem. The new plan is to add DIC to [PSS](https://github.com/StatisticalRethinkingJulia/ParetoSmoothedImportanceSampling.jl) tomorrow or Saturday.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [February 27, 2021, 6:55pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/10 "2021-02-27T18:55:16Z")

</div>

@ANasc, just a quick status update: [DIC](https://github.com/StatisticalRethinkingJulia/ModelComparisons.jl/pulls) was added to a new repo called `ModelComparisons`, along with BIC and AIC. This repo also has WAIC and PSILOO. It will take a few days for registration to be approved. I will also submit a change to [MCMCChains](https://github.com/itsdfish/MCMCChains.jl) so that you can calculate DIC directly from a Chain object. If you need to use `DIC` immediately, you can switch to my forks for `MCMCChains` and `ModelComparisons` until new releases are official. I plan on adding functionality to hook `MCMCChains` into WAIC and `PSILOO` in the next week or so.

---

<div class="post-metadata">

**Author:** ![ANasc](https://avatars.discourse-cdn.com/v4/letter/a/41988e/32.png) [@ANasc](https://discourse.julialang.org/u/ANasc)\
**Post date:** [February 28, 2021, 5:16am UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/11 "2021-02-28T05:16:29Z")

</div>

Amazing! Thanks so much. I look forward to trying it out. I think this will be a very useful addition for a lot of people.

---

<div class="post-metadata">

**Author:** ![Storopoli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/storopoli/32/209278_2.png) [@Storopoli](https://discourse.julialang.org/u/Storopoli)\
**Post date:** [March 3, 2021, 7:43am UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/12 "2021-03-03T07:43:54Z")

</div>

Great Job @Christopher_Fisher! Excited to see `PSIS-LOO` coming to life in `MCMCChains.jl`!

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [April 16, 2021, 9:33pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/13 "2021-04-16T21:33:55Z")

</div>

> [@Christopher\_Fisher](#):
>
> . I plan on adding functionality to hook `MCMCChains` into WAIC and `PSILOO` in the next week or so.

Wow. Very excited about this! Thanks slot.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [April 22, 2021, 1:26pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/14 "2021-04-22T13:26:51Z")

</div>

We encountered a few issues during the process. One is a licensing issue for PSILOO. My plan is to try to implement a new version as a work around. I am also trying to figure out the best place to have this functionality. I will follow up with an update once I know more.

---

<div class="post-metadata">

**Author:** ![Storopoli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/storopoli/32/209278_2.png) [@Storopoli](https://discourse.julialang.org/u/Storopoli)\
**Post date:** [April 22, 2021, 2:03pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/15 "2021-04-22T14:03:52Z")

</div>

Isn’t it under Creative Commons or something? Can’t we give deserved and mandatory attribution (citation) and use it?

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [April 22, 2021, 2:47pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/16 "2021-04-22T14:47:16Z")

</div>

The PSILOO code is under a GNU license. Other developers are reluctant to use it because under some interpretations, the GNU licence would apply to dependencies and supercede the MIT license.

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [May 1, 2021, 11:01am UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/17 "2021-05-01T11:01:10Z")

</div>

I think Turin.jl can be used gor hypothesis testing, but I’m completely lost by the documentation.

---

<div class="post-metadata">

**Author:** ![Storopoli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/storopoli/32/209278_2.png) [@Storopoli](https://discourse.julialang.org/u/Storopoli)\
**Post date:** [May 2, 2021, 8:12am UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/18 "2021-05-02T08:12:50Z")

</div>

You can try with the Savage-Dickey method: [Bayesian hypothesis testing for psychologists: A tutorial on the Savage–Dickey method - ScienceDirect](https://www.sciencedirect.com/science/article/abs/pii/S0010028509000826)

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [May 2, 2021, 8:58pm UTC](https://discourse.julialang.org/t/bayesian-model-selection/55725/19 "2021-05-02T20:58:15Z")

</div>

Looks, good, is there a package version, or do I need to put it together. I take it I can use MCMC, and it looks like this uses Mann Whitney U Test?