Add DIC

[itsdfish]

Hi all,

I worked with Rob to create a library of model comparison metrics at StatsModelComparisons. This package has a generic interface for using model comparison metrics. This PR provides a method for computing DIC from a chain object. Please let me know if there are any edits you would like me to make.

[itsdfish]

It looks like this failed due to a compatibility issue. MCMCChains tests on Julia 1.0.5, but StanDump requires Julia 1.1 or higher. What is the best path forward?

[devmotion]

IMO it is a bit problematic that the method requires a parameter :lp with the log-likelihood. The name is not standardized, i.e., it might not exist or it might denote something different such as the log joint probability or something else, depending on the sampling algorithm. Probably this was also the reason for the current implementation which requires a log-likelihood function.

Maybe it would be better to remove dic from MCMCChains and refer users to other packages such as StatsModelComparisons. Alternatively, one could fix the current implementation and add a simple example in the docstring, e.g., for the case where :lp is actually the log-likelihood.

[itsdfish]

Part of the rationale for creating StatsModelComparisons was to make model comparison and MCMC output more modular. One approach would be to add methods to each participating package to preprocess and interface with the primary method defined in StatsModelComparisons. An alternative approach would be to define some of those connecting methods in StatsModelComparisons, in which case, you can refer users to StatsModelComparisons for model comparison. The third option is to define your own methods for DIC, WAIC etc. in MCMCChains to do all the calculations. I think option 1 or 2 is preferable because it makes the ecosystem more standardized, modular and reduces redundancy, but of course you should choose the option that you think fits best with the goals of MCMCChains.

Thanks for pointing out the issue with lp. I also saw log_density in the chain output. It's not clear to me why there are two variables that give the same values in most cases. Would it be better to change lp to log_density? I am going to ping @goedman so he can participate in these discussions.

[devmotion]

I totally agree that it would be good to make the statistical methods more modular - this was also the main motivation for changing (most) implementations in MCMCChains to operate on arrays instead of Chain objects recently.

Thanks for pointing out the issue with lp. I also saw log_density in the chain output. It's not clear to me why there are two variables that give the same values in most cases. Would it be better to change lp to log_density? I am going to ping @goedman so he can participate in these discussions.

Using a different name does not resolve the general issue here. MCMCChains does not enforce any specific parameter name for log-likelihood values (if they exist at all) and lets you rename all parameters arbitrarily. So anything that is based on the evaluation of a specific parameter name is likely to be broken - it will error if it does not exist and, even worse, it will yield completely incorrect results if a parameter of this name exists but does not denote log-likelihood values.

[itsdfish]

Ok. I understand the issue now. I was suffering from a bit of myopia. I usually use MCMCChains with specific samplers in Turing. So the posterior log likelihoods have a specific name. In that case, the simplest solution is probably adding an additional argument for the variable name:

function dic(chain::Chains, post_log_like::Symbol)
    lps = get(chain, post_log_like).post_log_like[:] |> Array
    return dic(lps)
end

Do you think something like that would work?

[devmotion]

Yep, I think that sounds good. I would just suggest

function dic(chain::Chains, loglik::Symbol)
    lps = Array(chain[:, loglik, :])
    return dic(lps)
end

instead.

[itsdfish]

Thanks. I like that better. I modified your method so that lps is a vector rather than a matrix.

[devmotion]

MCMCChains tests on Julia 1.0.5, but StanDump requires Julia 1.1 or higher. What is the best path forward?

Do you know why StanDump requires Julia 1.1?

[itsdfish]

I'm not sure. Let me find out.

[tpapp]

As noted in tpapp/StanDump.jl#7, I don't quite remember, but if you really need 1.0, PRs are welcome.

That said, do you really need 1.0?

[devmotion]

I personally don't need 1.0 (or anything apart from the latest version) but I think it would be unfortunate to drop Julia 1.0 if it can be supported without much effort. Currently, it does not lead to any additional maintenance burden. AFAIK Julia 1.6 might become the new LTS version only after Julia 1.7 is released, so I'm not sure how long it will take until a new LTS version is available.

[tpapp]

it does not lead to any additional maintenance burden

Let's see — again, PRs are welcome.

Generally I think that so much happened since 1.0 in terms of new features that it does not make sense to support it.

[goedman]

The current version of StatsModelComparisons.jl uses StanSample.jl for testing purposes. Wouldn't that be an ongoing source of issues if StatsModelComparisons.jl is made part of testing MCMCChains.jl.

In the past there was a link between Turing and CmdStan for comparisons. I'm no longer sure that's still the case, but if so, that could also be problematic as it would mix and match CmdStan and StanSample.

If there are good reasons to go this route I would suggest to redo the StatsModelComparisons.jl tests using Turing.jl or DynamicHMC.jl.

In the PR, in src/MCMCChains.jl line 27: import ModelComparisons: dic uses the initial proposed package name.

[itsdfish]

"In the PR, in src/MCMCChains.jl line 27: import ModelComparisons: dic uses the initial proposed package name."

Good catch. I fixed the import statement.

"The current version of StatsModelComparisons.jl uses StanSample.jl for testing purposes."

Would this issue be resolved if you specified StanSample in extras instead of compat? It was not clear to me based on the documentation for Pkg.jl.

[devmotion]

The current version of StatsModelComparisons.jl uses StanSample.jl for testing purposes.

I.e., it is used only in the tests? Then it should not be listed among the dependencies but in an [extras] section (https://julialang.github.io/Pkg.jl/v1/creating-packages/#Test-specific-dependencies-in-Julia-1.0-and-1.1).

[itsdfish]

@devmotion, I noticed that too. Would the conflict be resolved if listed under extras?

[devmotion]

Yes, since the packages listed under extras are not pulled in and installed when the package is installed - they are only used for testing.

[goedman]

Just to make sure I understand this right. For the import statement in MCMCChains.jl, StatsModelsComparisons.jl needs to be in the dependencies of MCMCChains right?

In StatsModelsComparisons.jl I can turn all examples into tests and then make StanSample.jl a Test dependency.

Interesting!

[devmotion]

I opened a PR that also uses RecipeBase instead of StatsPlots: StatisticalRethinkingJulia/StatsModelComparisons.jl#4

[itsdfish]

Rob, I believe that is the case.

[goedman]

Very nice changes David! Thanks a lot. I'll remove the mc_path. At some point I ran in troubles with making it a constant but that was a few years ago. For this package it is not needed.

Alvaro definitely used the R and Python sources to do the translations after permission by the author (Aki Vehtari I think used Matlab). I will indeed update the LICENSE.md file.

[devmotion]

I added some suggestions.

[devmotion]

Looks good to me 👍 I am not sure what's the best way to test the implementation but maybe we could at least call it with some random data to make sure it is not completely broken?

[itsdfish]

Yeah. That is a good idea. I also added a test to ensure that it fails as expected too. dic was no longer visible after your changes. I added a using and export statement. I'm not sure how you want to handle that. importing dic should prevent other variables in StatsModelComparisons from being exported in the MCMCChains.

[devmotion]

I noticed that in the mean time the license of StatsModelComparisons was changed to GPL-3 - does this mean we would have to relicense MCMCChains.dic (or maybe even other parts?) as well? This seems quite unfortunate, in particular since StatsModelComparisons.dic is a completely straightforward Julia implementation and only StatsModelComparisons.psisloo was derived from a GPL-3 licensed MATLAB implementation.

[itsdfish]

Based on our previous conversation, I was under the impression that the GPL license would only apply to the code for psisloo methods. @goedman, did your plans change?

[goedman]

Yes, this recent change to only auto-merge packages with OSI-approved licenses was a surprise to me as well. A problem is that TuringLang is a product while StatisticalRethinkingJulia is just for learning purposes. I have derived the WAIC part from the SR implementation so legally it is probably also GPL-ish. But I'm definitely open to change the license file if that is possible and would work for TuringLang.

[devmotion]

@cpfiffer What's your opinion about the license issue?

[cpfiffer]

Dunno, I'm kind of the wrong person to ask -- I don't know much about licensing issues.

[devmotion]

Unfortunately, I also don't know much about these problems. I would like to fix dic and merge this PR, and I think it is great if there's a general interface for methods such as dic that we can implement instead of rolling our own implementations. However, my gut feeling is that the MCMCChains.dic would have to be licensed under GPLv3 as well since StatsModelComparisons uses this license (as far as I understand, we don't have to relicense the other parts of MCMCChains). Even though it's just a single function I think it would be unfortunate - the implementation of dic is quite trivial and, as far as I understand, in principle it could also use the MIT license in StatsModelComparisons.

I also just noticed that StatsBase already defines aic, bic, and aicc (it seems they are not overloaded or reexported in StatsModelComparisons even though the package depends on StatsBase but reimplemented - @itsdfish is this intentional?) but it seems dic is missing. Maybe one should add dic to StatsBase? Then we could skip the dependency on StatsModelComparisons completely for now and avoid all licensing issues.

[itsdfish]

@devmotion, I like your idea of submitting a PR to StatsBase for DIC. I did not know that they implemented other model comparison metrics.

I will look at the algorithm for PSISLOO to see if it would be feasible to implement. It looks much more involved, so it is probably a long shot for me.