# Bayesian model comparison / selection
**URL:** https://discourse.julialang.org/t/bayesian-model-comparison-selection/80420
**Category:** Statistics
**Tags:** question, turing, bayesian-inference
**Created:** [May 3, 2022, 11:54am UTC](https://discourse.julialang.org/t/bayesian-model-comparison-selection/80420 "2022-05-03T11:54:53Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@uferc](https://discourse.julialang.org/u/uferc)
#### Post date: [May 3, 2022, 11:54am UTC](https://discourse.julialang.org/t/bayesian-model-comparison-selection/80420/1 "2022-05-03T11:54:53Z")
Hello everyone,
my question revolves around the topic of Bayesian model comparison and henceforce selection.
I’ve build 3 alternative models (each NUTS(0.65), 4 chains, 8000 draws) and fit my data with the Turing.jl lib.
However I could not find a way in Turing to test loo characteristics and came across some discussions in the internet - from where I implemented the following for each model using the ArviZ interface:
```julia
ℓ = Turing.pointwise_loglikelihoods(param_mod, chains)
ℓ_mat = reduce(hcat, values(ℓ));
ℓ_arr = reshape(ℓ_mat, 1, size(ℓ_mat)...); # (chain_idx, sample_idx, parameter_idx)
data = ArviZ.from_mcmcchains(
chains,
library = "Turing",
log_likelihood = Dict("y" => ℓ_arr)
)
criterion = ArviZ.loo(data)
```
Finally I build a Dict and fed it into the compare function:
```julia
compare_dict = Dict(
"model 1" => criterion1,
"model 2" => criterion2,
"model 3" => criterion3,
)
compare("compare_dict", ic="loo")
```
I also tried to feed in the ArviZ interface data – which doesn’t work neither. I always get an none informative error message:
```julia
LoadError: PyError ($(Expr(:escape, :(ccall(#= C:\Users\xx\.julia\packages\PyCall\7a7w0\src\pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw)))))
AttributeError('Encountered error in ic computation of compare.')
```
Here in the documentation [Reference · ArviZ.jl](https://julia.arviz.org/v0.2/reference/#ArviZ.compare-Tuple) there is only a link "See documentation for [`arviz.compare`](https://arviz-devs.github.io/arviz/generated/arviz.compare.html)" provided where I find just as little information about the function as on the previous page.
From previous discussions here in the forum I know about the StatisticalRethinkingJulia / ParetoSmoothedImportanceSampling.jl package but I couldn’t find a compare function.
From the discussion here [Bayesian Model Selection - #8 by ANasc](https://discourse.julialang.org/t/bayesian-model-selection/55725/8), I got that the ModelComparisons package was in the end not added to the MCMCChains repo…?
Any recommendation how to approach this problem? How can I find out how the compare-function works, which argument it takes in which format? I there any other package that includes model comparison?
I am very happy about any recommendation!
---
### Author:  [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)
#### Post date: [May 3, 2022, 12:10pm UTC](https://discourse.julialang.org/t/bayesian-model-comparison-selection/80420/2 "2022-05-03T12:10:33Z")
I’m afraid that I don’t have an answer to your immediate question. However, another approach is to use ParetoSmooth.jl. This might be relevant: [Home · ParetoSmooth.jl](https://beta.turing.ml/ParetoSmooth.jl/stable/#ParetoSmooth.ModelComparison)
---
### Author:  [@sethaxen](https://discourse.julialang.org/u/sethaxen)
#### Post date: [May 3, 2022, 1:08pm UTC](https://discourse.julialang.org/t/bayesian-model-comparison-selection/80420/3 "2022-05-03T13:08:23Z")
> [@uferc](#):
>
> Finally I build a Dict and fed it into the compare function:
>
> ```julia
> compare_dict = Dict(
> "model 1" => criterion1,
> "model 2" => criterion2,
> "model 3" => criterion3,
> )
> compare("compare_dict", ic="loo")
>
> ```
>
> I also tried to feed in the ArviZ interface data – which doesn’t work neither. I always get an none informative error message:
>
> ```julia
> LoadError: PyError ($(Expr(:escape, :(ccall(#= C:\Users\xx\.julia\packages\PyCall\7a7w0\src\pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw)))))
> AttributeError('Encountered error in ic computation of compare.')
>
> ```
Ah, thanks, this is a bug. I have [opened an issue](https://github.com/arviz-devs/ArviZ.jl/issues/183) to track this. ArviZ converts the `arviz.ELPDData` returned by `arviz.loo` to a `DataFrame` on the Julia side, but the Python `compare` function expects an `ELPDData` input, and currently ArviZ doesn’t convert it back.
There are 2 quick workarounds.
First, we’ll load an example `InferenceData`:
```julia
julia> using ArviZ
julia> idata = load_arviz_data("radon");
```
Option 1: pass a `Dict` of `InferenceData` to `compare`:
```julia
julia> compare_dict1 = Dict(
"model 1" => idata,
"model 2" => idata,
"model 3" => idata,
);
julia> compare(compare_dict1; ic="loo")
3×10 DataFrame
Row │ name rank loo p_loo d_loo weight se dse warning loo_scale
│ String Int64 Float64 Float64 Float64 Float64 Float64 Float64 Bool String
─────┼────────────────────────────────────────────────────────────────────────────────────────────
1 │ model 1 0 -1027.12 26.7643 0.0 0.333333 28.8476 0.0 false log
2 │ model 3 1 -1027.12 26.7643 0.0 0.333333 28.8476 0.0 false log
3 │ model 2 2 -1027.12 26.7643 0.0 0.333333 28.8476 0.0 false log
```
Option 2: bypass ArviZ.jl’s conversion to `DataFrame` by calling the Python function directly:
```julia
julia> loo_result = ArviZ.arviz.loo(idata)
PyObject Computed from 2000 by 919 log-likelihood matrix
Estimate SE
elpd_loo -1027.12 28.85
p_loo 26.76 -
------
Pareto k diagnostic values:
Count Pct.
(-Inf, 0.5] (good) 919 100.0%
(0.5, 0.7] (ok) 0 0.0%
(0.7, 1] (bad) 0 0.0%
(1, Inf) (very bad) 0 0.0%
julia> compare_dict2 = Dict(
"model 1" => loo_result,
"model 2" => loo_result,
"model 3" => loo_result,
);
julia> compare(compare_dict2; ic="loo")
3×10 DataFrame
Row │ name rank loo p_loo d_loo weight se dse warning loo_scale
│ String Int64 Float64 Float64 Float64 Float64 Float64 Float64 Bool String
─────┼────────────────────────────────────────────────────────────────────────────────────────────
1 │ model 1 0 -1027.12 26.7643 0.0 0.333333 28.8476 0.0 false log
2 │ model 3 1 -1027.12 26.7643 0.0 0.333333 28.8476 0.0 false log
3 │ model 2 2 -1027.12 26.7643 0.0 0.333333 28.8476 0.0 false log
```
> [@](#):
>
> Here in the documentation [Reference · ArviZ.jl](https://julia.arviz.org/v0.2/reference/#ArviZ.compare-Tuple) there is only a link "See documentation for [`arviz.compare`](https://arviz-devs.github.io/arviz/generated/arviz.compare.html)" provided where I find just as little information about the function as on the previous page.
Whenever ArviZ.jl wraps a Python function with no modification, rather than manually convert all docs to Julia and try to keep them up-to-date, we refer the users to the Python docs. In the future these will all be pure Julia functions with their own documentation.
---
### Author:  [@uferc](https://discourse.julialang.org/u/uferc)
#### Post date: [May 3, 2022, 2:40pm UTC](https://discourse.julialang.org/t/bayesian-model-comparison-selection/80420/4 "2022-05-03T14:40:20Z")
Thank you so much @sethaxen for the fast reply! I tried the second workaround and got my result. Thanks 🤗