Skip to content
3 changes: 0 additions & 3 deletions docs/src/api.md
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@ eti!
AbstractELPDResult
PSISLOOResult
elpd_estimates
information_criterion
loo
```

Expand All @@ -36,12 +35,10 @@ loo
```@docs
ModelComparisonResult
compare
model_weights
```

The following model weighting methods are available
```@docs
AbstractModelWeightsMethod
BootstrappedPseudoBMA
PseudoBMA
Stacking
Expand Down
5 changes: 2 additions & 3 deletions src/PosteriorStats.jl
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,10 @@ export PSIS, PSISResult, psis, psis!

# LOO-CV
export AbstractELPDResult, PSISLOOResult
export elpd_estimates, information_criterion, loo
export elpd_estimates, loo

# Model weighting and comparison
export AbstractModelWeightsMethod, BootstrappedPseudoBMA, PseudoBMA, Stacking, model_weights
export BootstrappedPseudoBMA, PseudoBMA, Stacking
export ModelComparisonResult, compare

# Summary statistics
Expand All @@ -49,7 +49,6 @@ export eti, eti!, hdi, hdi!
export loo_pit, r2_score

const DEFAULT_CI_PROB = 0.89f0
const INFORMATION_CRITERION_SCALES = (deviance=-2, log=1, negative_log=-1)

include("utils.jl")
include("show_prettytable.jl")
Expand Down
32 changes: 20 additions & 12 deletions src/compare.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,15 @@ The ELPD is estimated by Pareto smoothed importance sampling leave-one-out cross

# Keywords

- `weights_method::AbstractModelWeightsMethod=Stacking()`: the method to be used to weight
the models. See [`model_weights`](@ref) for details
- `method=Stacking()`: the method to be used to weight the models
(see [Yao2018](@citet) for details):
+ [`Stacking`](@ref): Stacking of predictive distributions. The default and recommended
approach, as it performs well even when the true data generating process is not
included among the candidate models.
+ [`BootstrappedPseudoBMA`](@ref): pseudo-Bayesian Model averaging using Akaike-type
weighting, where the weights are stabilized using the Bayesian bootstrap.
+ [`PseudoBMA`](@ref): pseudo-Bayesian Model averaging using Akaike-type weighting
(not recommended).
- `sort::Bool=true`: Whether to sort models by decreasing ELPD.

# Returns
Expand Down Expand Up @@ -54,7 +61,7 @@ Compare the same models from pre-computed PSIS-LOO results and computing
```jldoctest compare; setup = :(using Random; Random.seed!(23))
julia> elpd_results = mc.elpd_result;

julia> compare(elpd_results; weights_method=BootstrappedPseudoBMA())
julia> compare(elpd_results; method=BootstrappedPseudoBMA())
ModelComparisonResult with BootstrappedPseudoBMA weights
rank elpd se_elpd elpd_diff se_elpd_diff weight p se_p
non_centered 1 -31 1.5 0 0.0 0.51 0.9 0.32
Expand All @@ -64,17 +71,18 @@ ModelComparisonResult with BootstrappedPseudoBMA weights
# References

- [Spiegelhalter2002](@cite) Spiegelhalter et al. J. R. Stat. Soc. B 64 (2002)
- [Yao2018](@cite) Yao et al. Bayesian Analysis 13, 3 (2018)
"""
function compare(
inputs;
weights_method::AbstractModelWeightsMethod=Stacking(),
method::AbstractModelWeightsMethod=Stacking(),
model_names=_indices(inputs),
sort::Bool=true,
)
length(model_names) === length(inputs) ||
throw(ArgumentError("Length of `model_names` must match length of `inputs`"))
elpd_results = map(_maybe_loo, inputs)
weights = model_weights(weights_method, elpd_results)
weights = model_weights(method, elpd_results)
perm = _sortperm(elpd_results; by=x -> elpd_estimates(x).elpd, rev=true)
i_elpd_max = first(perm)
elpd_max_i = elpd_estimates(elpd_results[i_elpd_max]; pointwise=true).elpd
Expand All @@ -89,7 +97,7 @@ function compare(
se_elpd_diff = map(last, se_elpd_diff_and)
rank = _assimilar(elpd_results, (1:length(elpd_results))[perm])
result = ModelComparisonResult(
model_names, rank, elpd_diff, se_elpd_diff, weights, elpd_results, weights_method
model_names, rank, elpd_diff, se_elpd_diff, weights, elpd_results, method
)
sort || return result
return _permute(result, perm)
Expand Down Expand Up @@ -117,13 +125,13 @@ struct ModelComparisonResult{E,N,R,W,ER,M}
elpd_diff::E
"Standard error of the ELPD difference"
se_elpd_diff::E
"Model weights computed with `weights_method`"
"Model weights computed with `method`"
weight::W
"""`AbstactELPDResult`s for each model, which can be used to access useful stats like
ELPD estimates, pointwise estimates, and Pareto shape values for PSIS-LOO"""
elpd_result::ER
"Method used to compute model weights with [`model_weights`](@ref)"
weights_method::M
"Method used to compute model weights"
method::M
end

#### custom tabular show methods
Expand All @@ -140,7 +148,7 @@ function _show(io::IO, mime::MIME, r::ModelComparisonResult; kwargs...)
cols = Tables.columnnames(r)[2:end]
table = NamedTuple{cols}(Tables.columntable(r))

weights_method_name = _typename(r.weights_method)
method_name = _typename(r.method)
weights = table.weight
digits_weights = ceil(Int, -log10(maximum(weights))) + 1
weight_formatter = PrettyTables.ft_printf(
Expand All @@ -150,7 +158,7 @@ function _show(io::IO, mime::MIME, r::ModelComparisonResult; kwargs...)
io,
mime,
table;
title="ModelComparisonResult with $(weights_method_name) weights",
title="ModelComparisonResult with $(method_name) weights",
row_labels,
extra_formatters=(weight_formatter,),
kwargs...,
Expand All @@ -160,7 +168,7 @@ end
function _permute(r::ModelComparisonResult, perm)
return ModelComparisonResult(
(_permute(getfield(r, k), perm) for k in fieldnames(typeof(r))[1:(end - 1)])...,
r.weights_method,
r.method,
)
end

Expand Down
30 changes: 0 additions & 30 deletions src/elpdresult.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@ each, from which other relevant estimates can be computed.

Subtypes implement the following functions:
- [`elpd_estimates`](@ref)
- [`information_criterion`](@ref)
"""
abstract type AbstractELPDResult end

Expand All @@ -29,35 +28,6 @@ Return the (E)LPD estimates from the `result`.
"""
function elpd_estimates end

"""
$(FUNCTIONNAME)(elpd, scale::Symbol)

Compute the information criterion for the given `scale` from the `elpd` estimate.

`scale` must be one of `$(keys(INFORMATION_CRITERION_SCALES))`.

See also: [`loo`](@ref)
"""
function information_criterion(estimates, scale::Symbol)
scale_value = INFORMATION_CRITERION_SCALES[scale]
return scale_value * estimates.elpd
end

"""
$(FUNCTIONNAME)(result::AbstractELPDResult, scale::Symbol; pointwise=false)

Compute information criterion for the given `scale` from the existing ELPD `result`.

`scale` must be one of `$(keys(INFORMATION_CRITERION_SCALES))`.

If `pointwise=true`, then pointwise estimates are returned.
"""
function information_criterion(
result::AbstractELPDResult, scale::Symbol; pointwise::Bool=false
)
return information_criterion(elpd_estimates(result; pointwise), scale)
end

function _lpd_pointwise(log_likelihood, dims)
ndraws = prod(Base.Fix1(size, log_likelihood), dims)
lpd = LogExpFunctions.logsumexp(log_likelihood; dims)
Expand Down
11 changes: 10 additions & 1 deletion src/loo.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,13 @@ end
Compute the Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO).
[Vehtari2017, LOOFAQ](@cite)

!!! details Relationship to information criteria

PSIS-LOO estimates the expected log pointwise predictive density (ELPD) and effective
number of parameters (p). The ELPD can be converted to the "deviance" scale sometimes
used for information criteria by multiplying it by `-2`. See [LOOFAQ](@citet) for more
details.

`log_likelihood` must be an array of log-likelihood values with shape
`(chains, draws[, params...])`.

Expand All @@ -46,7 +53,9 @@ Compute the Pareto-smoothed importance sampling leave-one-out cross-validation (
[`MCMCDiagnosticTools.ess`](@extref).
- `kwargs`: Remaining keywords are forwarded to [`PSIS.psis`](@extref).

See also: [`PSISLOOResult`](@ref)
# Returns

- [`PSISLOOResult`](@ref): A container for the PSIS-LOO results.

# Examples

Expand Down
58 changes: 0 additions & 58 deletions src/model_weights.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,64 +9,6 @@ Subtypes implement [`model_weights`](@ref)`(method, elpd_results)`.
"""
abstract type AbstractModelWeightsMethod end

"""
model_weights(elpd_results; method=Stacking())
model_weights(method::AbstractModelWeightsMethod, elpd_results)

Compute weights for each model in `elpd_results` using `method`.

`elpd_results` is a `Tuple`, `NamedTuple`, or `AbstractVector` with
[`AbstractELPDResult`](@ref) entries. The weights are returned in the same type of
collection.

[`Stacking`](@ref) is the recommended approach, as it performs well even when the true data
generating process is not included among the candidate models. See [Yao2018](@citet) for
details.

See also: [`AbstractModelWeightsMethod`](@ref), [`compare`](@ref)

# Examples

Compute [`Stacking`](@ref) weights for two models:

```jldoctest model_weights; filter = [r"└.*", r"(\\d+\\.\\d{3})\\d*" => s"\\1"]
julia> using ArviZExampleData

julia> models = (
centered=load_example_data("centered_eight"),
non_centered=load_example_data("non_centered_eight"),
);

julia> elpd_results = map(models) do idata
log_like = PermutedDimsArray(idata.log_likelihood.obs, (2, 3, 1))
return loo(log_like)
end;
┌ Warning: 1 parameters had Pareto shape values 0.7 < k ≤ 1. Resulting importance sampling estimates are likely to be unstable.
└ @ PSIS ~/.julia/packages/PSIS/...

julia> model_weights(elpd_results; method=Stacking()) |> pairs
pairs(::NamedTuple) with 2 entries:
:centered => 3.50546e-31
:non_centered => 1.0
```

Now we compute [`BootstrappedPseudoBMA`](@ref) weights for the same models:

```jldoctest model_weights; setup = :(using Random; Random.seed!(94))
julia> model_weights(elpd_results; method=BootstrappedPseudoBMA()) |> pairs
pairs(::NamedTuple) with 2 entries:
:centered => 0.492513
:non_centered => 0.507487
```

# References

- [Yao2018](@cite) Yao et al. Bayesian Analysis 13, 3 (2018)
"""
function model_weights(elpd_results; method::AbstractModelWeightsMethod=Stacking())
return model_weights(method, elpd_results)
end

# Akaike-type weights are defined as exp(-AIC/2), normalized to 1, which on the log-score
# IC scale is equivalent to softmax
akaike_weights!(w, elpds) = LogExpFunctions.softmax!(w, elpds)
Expand Down
10 changes: 5 additions & 5 deletions test/compare.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ end
@test mc1.elpd_diff.non_centered == 0.0
@test mc1.elpd_diff.centered > 0
@test mc1.weight == NamedTuple{(:non_centered, :centered)}(
model_weights(eight_schools_loo_results)
PosteriorStats.model_weights(Stacking(), eight_schools_loo_results)
)
@test mc1.elpd_result ==
NamedTuple{(:non_centered, :centered)}(eight_schools_loo_results)
Expand All @@ -41,11 +41,11 @@ end
end

@testset "keywords are forwarded" begin
mc2 = compare(eight_schools_loo_results; weights_method=PseudoBMA())
mc2 = compare(eight_schools_loo_results; method=PseudoBMA())
@test !_isequal(mc2, compare(eight_schools_loo_results))
@test mc2.weights_method === PseudoBMA()
@test mc2.method === PseudoBMA()
mc3 = compare(eight_schools_loo_results; sort=false)
for k in filter(!=(:weights_method), propertynames(mc1))
for k in filter(!=(:method), propertynames(mc1))
if k === :name
@test getproperty(mc3, k) == reverse(getproperty(mc1, k))
else
Expand Down Expand Up @@ -98,7 +98,7 @@ end
end

@testset "show" begin
mc5 = compare(eight_schools_loo_results; weights_method=PseudoBMA())
mc5 = compare(eight_schools_loo_results; method=PseudoBMA())
@test sprint(show, "text/plain", mc1) == """
ModelComparisonResult with Stacking weights
rank elpd se_elpd elpd_diff se_elpd_diff weight p se_p
Expand Down
11 changes: 0 additions & 11 deletions test/loo.jl
Original file line number Diff line number Diff line change
Expand Up @@ -32,17 +32,6 @@ using Test
@test loo_result.psis_result.reff == pointwise.reff
@test loo_result.psis_result.pareto_shape == pointwise.pareto_shape
end
@testset "information criterion" begin
@test information_criterion(loo_result, :log) == estimates.elpd
@test information_criterion(loo_result, :negative_log) == -estimates.elpd
@test information_criterion(loo_result, :deviance) == -2 * estimates.elpd
@test information_criterion(loo_result, :log; pointwise=true) ==
pointwise.elpd
@test information_criterion(loo_result, :negative_log; pointwise=true) ==
-pointwise.elpd
@test information_criterion(loo_result, :deviance; pointwise=true) ==
-2 * pointwise.elpd
end
end
end
# @testset "keywords forwarded" begin
Expand Down
Loading
Loading