# A Deep Learning Method for Comparing Bayesian Hierarchical Models

Lasse Elsemüller
Institute of Psychology
Heidelberg University
<EMAIL_ADDRESS>&Martin Schnuerch
Department of Psychology
University of Mannheim
<EMAIL_ADDRESS>&Paul-Christian Bürkner
Department of Statistics
TU Dortmund University
<EMAIL_ADDRESS>&Stefan T. Radev
Cluster of Excellence STRUCTURES
Heidelberg University
<EMAIL_ADDRESS>
###### Abstract

Bayesian model comparison (BMC) offers a principled approach for assessing the
relative merits of competing computational models and propagating uncertainty
into model selection decisions. However, BMC is often intractable for the
popular class of hierarchical models due to their high-dimensional nested
parameter structure. To address this intractability, we propose a deep
learning method for performing BMC on any set of hierarchical models which can
be instantiated as probabilistic programs. Since our method enables amortized
inference, it allows efficient re-estimation of posterior model probabilities
and fast performance validation prior to any real-data application. In a
series of extensive validation studies, we benchmark the performance of our
method against the state-of-the-art bridge sampling method and demonstrate
excellent amortized inference across all BMC settings. We then showcase our
method by comparing four hierarchical evidence accumulation models that have
previously been deemed intractable for BMC due to partly implicit likelihoods.
Additionally, we demonstrate how transfer learning can be leveraged to enhance
training efficiency. We provide reproducible code for all analyses and an
open-source implementation of our method.

##  1 Introduction

Hierarchical or multilevel models (HMs) play an increasingly important
methodological role in the social and cognitive sciences [134, 192]. HMs
embody probabilistic and structural information about nested data occurring
frequently in various settings, such as educational research [209],
experimental psychology [211], epidemiology [161] or astrophysics [159], to
name just a few. Crucially, HMs can often extract more information from rich
data structures than their non-hierarchical counterparts (e.g., aggregate
analyses), while retaining a relatively high intrinsic interpretability of
their structural components (i.e., parameters). Moreover, viewed as formal
instantiations of scientific hypotheses, HMs can be employed to systematically
assign preferences to these hypotheses by means of formal model comparison.
For example, [153] proposed a powerful framework based on Bayesian HMs for
formulating and testing competing theoretical positions on quantitative vs.
qualitative individual differences.

We consider Bayesian model comparison (BMC) as a principled framework for
comparing and ranking competing HMs via Occam’s razor [162, 169, 168].
However, standard BMC is analytically intractable for non-trivial HMs, as it
requires marginalization over high-dimensional parameter spaces. Moreover, BMC
for complex HMs without explicit likelihoods (i.e., HMs available only as
randomized simulators) becomes increasingly hopeless and precludes many
interesting applications in the rapidly expanding field of simulation-based
inference [127].

In this work, we propose to tackle the problem of BMC for arbitrarily complex
HMs from a simulation-based perspective using deep learning. In particular, we
build on the BayesFlow framework [186, 184] for simulation-based Bayesian
inference and propose a novel hierarchical neural network architecture for
approximating Bayes factors (BFs) and posterior model probabilities (PMPs) for
any collection of HMs.

Our neural approach circumvents the steps of explicitly fitting all models and
marginalizing over the parameter space of each model. Thus, it is applicable
to both HMs with explicit likelihood functions and HMs accessible only through
Monte Carlo simulations (i.e., with implicit likelihood functions). Moreover,
our neural networks come with an efficient way to compute their calibration
error [151], which provides an important diagnostic for self-consistency.
Lastly, trained networks can be adapted to related tasks, substantially
reducing the computational burden when dealing with demanding simulators.

The remainder of this paper is organized as follows. In Section 2, we
introduce the theoretical background and related work on (hierarchical) BMC.
We then present the rationale and details of our deep learning method in
Section 3. In Sections 4.1 and 4.2, we present two validation studies of the
proposed method: One that includes toy models for illustrative purposes and
one that includes two popular classes of models from the field of cognitive
psychology. In Section 4.3, we then apply our method to compare hierarchical
diffusion decision models with partly intractable likelihoods on a real data
set. Finally, Section 5 summarizes our contributions and discusses future
perspectives.

##  2 Theoretical Background

###  2.1 Bayesian Hierarchical Modeling

In order to streamline statistical analyses, researchers rely on assumptions
about the probabilistic structure or symmetry of the assumed data-generating
process. For instance, the canonical IID assumption in psychological modeling
states that (multivariate) observations are independent of each other and
sampled from the same latent probability distribution [198, 180].

However, more complex dependencies may arise in a variety of contexts. For
instance, if there are repeated measurements per participant or participants
belong to different natural groups (e.g., school classes, working groups), the
respective observations exhibit higher correlations within those clusters than
across them. Ignoring this nested structure in statistical analyses may result
in biased conclusions [198]. Bayesian HMs formalize this structural knowledge
by assuming that observations are sampled from a multilevel generative process
[136].

For instance, the generative recipe for a two-level Bayesian HM can be written
as:

$\displaystyle\boldsymbol{\eta}$ $\displaystyle\sim p(\boldsymbol{\eta})$ (1)
$\displaystyle\boldsymbol{\theta}_{m}$ $\displaystyle\sim
p(\boldsymbol{\theta}\,\,|\,\boldsymbol{\eta})\text{ for }m=1,\dots,M$ (2)
$\displaystyle\boldsymbol{x}_{mn}$ $\displaystyle\sim
p(\boldsymbol{x}\,\,|\,\,\boldsymbol{\theta}_{m})\text{ for }n=1,\dots,N_{m},$
(3)

where $\boldsymbol{\eta}$ denotes the group-level parameters,
$\boldsymbol{\theta}_{m}$ denotes the individual parameters in group $m$ and
$\boldsymbol{x}_{mn}$ represents the $n$-th observation in group $m$. Such a
model suggests the following (non-unique) factorization of the joint
distribution:

$p(\boldsymbol{\eta},\\{\boldsymbol{\theta}_{m}\\},\\{\boldsymbol{x}_{mn}\\})=\\\
p(\boldsymbol{\eta})\prod_{m=1}^{M}p\left(\boldsymbol{\theta}_{m}\,\vline\,\boldsymbol{\eta}\right)\prod_{n=1}^{N_{m}}p\left(\boldsymbol{x}_{mn}\,|\,\boldsymbol{\theta}\right).$
(4)

The set notation $\\{\boldsymbol{\theta}_{m}\\}$ and
$\\{\boldsymbol{x}_{mn}\\}$ implies that the number of groups and observations
in each group can vary across simulations, data sets and experiments and that
these quantities are exchangeable.

HMs can be considered as a compromise between a separate analysis of each
group (no-pooling) that neglects the information contained in the rest of the
data and an aggregate analysis of the data (complete pooling) that loses the
distinction between intra-group and inter-group variability [160]. The partial
pooling of information induced by HMs leads to more stable and accurate
individual estimates through the shrinkage properties of multilevel priors,
whereby single estimates inform each other [136, 120].

(a) Marginal Likelihoods

(b) Posterior Model Probabilities

Figure 1: Hypothetical BMC setting with a simple model $\mathcal{M}_{1}$ and a
more complex model $\mathcal{M}_{2}$. (a) The complex model which accounts for
a broader range of observations needs to spread its marginal likelihood to
cover its larger generative scope. It does so at the cost of diminished
sharpness. Thus, even though observation $x_{1}$ is well within its generative
scope, the simpler model $\mathcal{M}_{1}$ yields a higher marginal likelihood
and is therefore preferred. In contrast, observation $x_{0}$ has a higher
marginal likelihood under model $\mathcal{M}_{2}$, as it is very unlikely to
be generated by the simpler model $\mathcal{M}_{1}$. (b) The corresponding
posterior model probabilities (PMPs) given a uniform model prior.

Despite having desirable properties, hierarchical modeling comes at the cost
of increased complexity and computational demands. These increased demands
make it hard or even impossible to compare competing HMs within the
probabilistic framework of BMC. Before we highlight these challenges, we first
describe the basics of BMC for non-hierarchical models.

###  2.2 Bayesian Model Comparison

The starting point of BMC is a collection of $J$ competing generative models
$\mathcal{M}=\\{\mathcal{M}_{1},\mathcal{M}_{2},\dots,\mathcal{M}_{J}\\}$.
Each $\mathcal{M}_{j}$ is associated with a prior
$p\left(\boldsymbol{\theta}_{j}\,|\,\mathcal{M}_{j}\right)$ on the parameters
$\boldsymbol{\theta}_{j}$ and a generative mechanism, which is either defined
analytically through a (tractable) likelihood density function
$p\left(\boldsymbol{x}\,|\,\boldsymbol{\theta}_{j},\mathcal{M}_{j}\right)$ or
realized as a Monte Carlo simulation program
$g_{j}(\boldsymbol{\theta},\boldsymbol{z})$ with random states
$\boldsymbol{z}$. Together, the prior and the likelihood define the Bayesian
joint model

$\displaystyle
p\left(\boldsymbol{\theta}_{j},\boldsymbol{x}\,|\,\mathcal{M}_{j}\right)=p\left(\boldsymbol{\theta}_{j}\,|\,\mathcal{M}_{j}\right)p\left(\boldsymbol{x}\,|\,\boldsymbol{\theta}_{j},\mathcal{M}_{j}\right),$
(5)

which is also tacitly defined for simulator-based models by marginalizing the
joint distribution
$p\left(\boldsymbol{x},\boldsymbol{z}\,|\,\boldsymbol{\theta}_{j},\mathcal{M}_{j}\right)$
over all possible execution paths (i.e., random states) of the simulation
program to obtain the implicit likelihood

$\displaystyle
p\left(\boldsymbol{x}\,|\,\boldsymbol{\theta}_{j},\mathcal{M}_{j}\right)=\int
p\left(\boldsymbol{x},\boldsymbol{z}\,|\,\boldsymbol{\theta}_{j},\mathcal{M}_{j}\right)\text{d}\boldsymbol{z}.$
(6)

This integral is typically intractable for complex simulators [127], which
makes it impossible to evaluate the likelihood and use standard Bayesian
methods for parameter inference or model comparison.

The likelihood function, be it explicit or implicit, is a key object in
Bayesian inference. When the parameters $\boldsymbol{\theta}$ are
systematically varied and the data $\boldsymbol{x}$ held constant, the
likelihood quantifies the relative fit of each model instantiation (defined by
a fixed configuration $\boldsymbol{\theta}$) to the observed data.

When we marginalize the Bayesian joint model (Eq. 5) over its parameter space,
we obtain the marginal likelihood or Bayesian evidence [169, see]:

$\displaystyle p\left(\boldsymbol{x}\,|\,\mathcal{M}_{j}\right)=\int
p\left(\boldsymbol{x}\,|\,\boldsymbol{\theta}_{j},\mathcal{M}_{j}\right)p\left(\boldsymbol{\theta}_{j}\,|\,\mathcal{M}_{j}\right)\text{d}\boldsymbol{\theta}_{j}.$
(7)

The marginal likelihood can be interpreted as the probability that we would
generate data $\boldsymbol{x}$ from model $\mathcal{M}_{j}$ when we randomly
sample from the model’s parameter prior
$p\left(\boldsymbol{\theta}_{j}\,|\,\mathcal{M}_{j}\right)$. Moreover, the
marginal likelihood is a central quantity for prior predictive hypothesis
testing or model selection [162, 182, 191]. It is well-known that the marginal
likelihood encodes a notion of Occam’s razor arising from the basic principles
of probability [162]. Thus, the marginal likelihood provides a foundation for
the widespread use of Bayes factors [155, BFs;] or posterior model
probabilities [126, PMPs;] for BMC.

The relative evidence for a pair of models can be computed through the ratio
of marginal likelihoods for the two competing models $\mathcal{M}_{j}$ and
$\mathcal{M}_{k}$,

$\displaystyle\text{BF}_{jk}=\frac{p\left(\boldsymbol{x}\,|\,\mathcal{M}_{j}\right)}{p\left(\boldsymbol{x}\,|\,\mathcal{M}_{k}\right)}.$
(8)

This ratio is called Bayes factor (BF) and is widely used for quantifying
pairwise model preference in Bayesian settings [162, 155]. Accordingly, a
$\text{BF}_{jk}$ $>1$ indicates preference for model $j$ over model $k$ given
available data $\boldsymbol{x}$. Alternatively, one can directly focus on the
(marginal) posterior probability of a model $\mathcal{M}_{j}$,

$p\left(\mathcal{M}_{j}\,|\,\boldsymbol{x}\right)=\frac{p(\boldsymbol{x}\,|\,\mathcal{M}_{j})\,p(\mathcal{M}_{j})}{\sum_{j=1}^{J}p\left(\boldsymbol{x}\,|\,\mathcal{M}_{j}\right)\,p(\mathcal{M}_{j})},$
(9)

where $p(\mathcal{M}_{j})$ is a categorical (typically uniform) prior
distribution encoding a researcher’s prior beliefs regarding the plausibility
of each considered model. This prior distribution is then updated with the
information contained in the marginal likelihood
$p(\boldsymbol{x}\,|\,\mathcal{M}_{j})$ to obtain the corresponding posterior
model probability (PMP), $p(\mathcal{M}_{j}\,|\,\boldsymbol{x})$. Occasionally
in the text, we will refer to the vector of PMPs for all $J$ models as
$\boldsymbol{\pi}$ and to the individual PMPs as $\pi_{j}$. The ratio of two
PMPs, known as posterior odds, is in turn connected to the Bayes factor via
the corresponding model priors:

$\frac{p(\mathcal{M}_{j}\,|\,\boldsymbol{x})}{p(\mathcal{M}_{k}\,|\,\boldsymbol{x})}=\frac{p\left(\boldsymbol{x}\,|\,\mathcal{M}_{j}\right)}{p\left(\boldsymbol{x}\,|\,\mathcal{M}_{k}\right)}\times\frac{p(\mathcal{M}_{j})}{p(\mathcal{M}_{k})}.$
(10)

Despite its intuitive appeal, the marginal likelihood (and thus BFs and PMPs)
represents a well-known and widely appreciated source of intractability in
Bayesian workflows, since it typically involves a multi-dimensional integral
(Eq. 7) over potentially unbounded parameter spaces [168, 146]. Furthermore,
the marginal likelihood becomes doubly intractable when the likelihood
function is itself not available (e.g., in simulation-based settings), thereby
making the comparison of such models a challenging and sometimes, up to this
point, hopeless endeavor.

Unsurprisingly, estimating the marginal likelihood (Eq. 7) in the context of
hierarchical models becomes even more challenging, since the number of
parameters over which we need to perform marginalization grows dramatically
(i.e., parameters at all hierarchical levels enter the computation). These
computational demands render the probabilistic comparison of HMs based on BFs
or PMPs analytically intractable even for relatively simple models with
explicit (analytical) likelihoods. Therefore, researchers need to resort to
costly, approximate methods which typically only work for models with explicit
likelihoods [137, 173, 146].

###  2.3 Approximate Bayesian Model Comparison

####  2.3.1 Explicit Likelihoods

The most efficient approximate methods to date require all candidate models to
possess explicitly available likelihood functions. For the most simple
scenario in which two HMs are nested (e.g., through an equality constraint on
a parameter), the Savage-Dickey density ratio [131] provides a convenient
approximation of the BF [217]. Typically, however, the candidate models are
not nested but exhibit notable structural differences. Thus, a general-purpose
method is needed to encompass the entire plethora of model comparison
scenarios arising in practical applications.

A more general method, and the current state-of-the-art for comparing HMs in
psychological and cognitive modeling [150, 145, 195], is given by bridge
sampling [115, 174]. Bridge sampling has enabled comparisons within families
of complex process models, such as multinomial processing trees [150, MPTs;]
or evidence accumulation models [145, EAMs;], and serves as a simple add-on
for Markov chain Monte Carlo (MCMC) based Bayesian workflows.

Crucially, bridge sampling relies on the posterior draws generated by an MCMC
sampler [122, e.g., Stan;] to efficiently approximate the marginal likelihood
of each respective model [146]. Note, however, that bridge sampling requires
considerably more random draws for stable results than standard parameter
estimation [147, usually about an order of magnitude more;]. Moreover, the
approximation quality of bridge sampling is dependent on the convergence of
the MCMC chains [145]. Finally, there are no strong theoretical guarantees
that the approximations are unbiased and accurately reflect the true marginal
likelihoods [195].

####  2.3.2 Implicit Likelihoods

With the rise of complex, high-resolution models, intractable likelihood
functions (i.e., functions that do not admit a closed form or are too costly
to evaluate) become more and more common in statistical modeling. Such models
are not limited to psychology and cognitive science [210, 180], but are also
common in fields such as neuroscience [141], epidemiology [185], population
genetics [183] or astrophysics [157]. Despite the common term likelihood-free,
simulator-based models still possess an implicitly defined likelihood (see
Section 2.2) from which we can obtain random draws through Monte Carlo
simulations. This enables model comparison through simulation-based methods,
usually by means of approximate Bayesian computation [183, 170, 175, ABC;].

Traditional (rejection-based) ABC methods for BMC repeatedly simulate data
sets from the specified generative models, retaining only those simulations
that are sufficiently similar to the empirical data. To enable the calculation
of this (dis-)similarity even in high-dimensional cases, the information
contained in the simulated data sets is reduced by computing hand-crafted
summary statistics, such as the mean and variance [202, 129]. The resulting
acceptance rates of the candidate models represent the approximations of their
PMPs [170, 175].

Even for non-hierarchical models, ABC methods are known to be notoriously
inefficient and highly dependent on the concrete choice of summary statistics
[170, 127]. This choice is even more challenging for HMs, as modelers now have
to retain an optimal amount of information on multiple levels. Moreover, the
rapidly growing number of summary statistics reduces the probability that a
simulated data set is similar enough to the empirical data, which vastly
increases the number of required simulations [113, 170].

Regardless of the number of summary statistics, their manual computation
carries the danger of insufficiently summarizing the simulations and thereby
producing biased approximations [170, a phenomenon known as curse of
insufficiency;]. While many improvements of rejection-based ABC have been
proposed, most notably ABC-MCMC [171, 207], ABC-SMC [200], as well as Gibbs
ABC [208] for Bayesian hierarchical modeling in particular [123, 135, see
also], these advancements are still limited by their dependence on hand-
crafted summary statistics or kernel density estimation methods.

Figure 2: Our proposed hierarchical neural network architecture for encoding
permutation invariance in the transformation of nested, two-level data into
posterior model probabilities. A first invariant module $\Sigma_{I}^{(1)}$
reduces all $N_{m}$ observations within each of the $M$ groups to a single
intermediary embedding vector $\widetilde{\boldsymbol{x}}_{m}$. For
readability of the figure, we display $N_{m}=N$ as constant across each group.
A second invariant module $\Sigma_{I}^{(2)}$ reduces all intermediary
embedding vectors to a hierarchical embedding vector $\boldsymbol{z}$, which
gets passed through an inference network to arrive at the final vector
$\widehat{\boldsymbol{\pi}}$ of approximated posterior model probabilities.

Recent developments, such as ABC-RF [183], combine ABC with machine learning
methods to build more expressive approximators for BMC problems. Accordingly,
model comparison is treated as a supervised learning problem – the simulated
data encompasses a training set for a machine learning algorithm that learns
to recognize the true generative model from which the data set was simulated.
The machine learning approach reduces the inefficiency problem that haunts
rejection-based ABC methods, but does not alleviate the curse of insufficiency
[170].

###  2.4 Bayesian Model Comparison with Neural Networks

Recently, [184] explored a method for simulation-based BMC using specialized
neural networks. The authors proposed to jointly train two specialized neural
networks using Monte Carlo simulations from each candidate model in
$\mathcal{M}$: a summary network and an evidential network. The goal of the
summary network is to extract maximally informative (in the optimal case,
sufficient) summary statistics from complex data sets. The goal of the
evidential network is to approximate PMPs as accurately as possible and,
optionally, to quantify their epistemic uncertainty.

Importantly, simulation-based training of neural networks enables amortized
inference for both implicit and explicit likelihood models. Amortization is a
property that ensures rapid inference for an arbitrary amount of data sets
after a potentially high computational investment for simulation and training
[186, 184, 176]. As a consequence, the calibration [151, 203] or the
inferential adequacy [194, 195] of an amortized Bayesian method are
embarrassingly easy to validate in practice.

In contrast, non-amortized methods, such as ABC-MCMC [207] or ABC-SMC [200]
need to repeat all computations from scratch for each observed data set.
Thereby, it is often infeasible to assess their calibration or inferential
adequacy in the pre-data phase of a Bayesian workflow [138].

Unfortunately, the evidential method proposed by [184] is not applicable to
HMs due to their nested probabilistic structure which cannot be tackled via
previous summary networks. This severely limits the applicability of the
method in quantitative research, where hierarchical models have been advocated
as a default choice [172, 192, 165]. In the following, we describe how to
extend the original method to enable amortized BMC for HMs.

##  3 Method

At its core, our method involves a multilevel permutation invariant neural
network which is aligned to the probabilistic symmetry of the underlying HMs
(see Figure 2 for a visualization). We hold that any method which does not
rely on ad hoc summary statistics should take this probabilistic symmetry
(e.g., exchangeability) into account in order to ensure the structural
faithfulness of its approximations. Moreover, respecting the probabilistic
symmetry implied by a generative model cannot only make simulation-based
training easier but also suggests a particular architecture for building
neural Bayesian approximators.

###  3.1 Permutation Invariance

Permutation invariance is the functional equivalent of the probabilistic
notion of exchangeability [136, 118], which roughly states that the order of
random variables should not influence their joint probability.

To illustrate this point, consider the model in Eq. 4, which has two
exchangeable levels by design, indexed by $m\in\\{1\dots,M\\}$ and
$n\in\\{1,\dots,N_{m}\\}$. In a setting familiar to social scientists, we
might have $M$ individuals, each of whom provides $N_{m}$ (multivariate)
responses on some scale or in repeated trials of an experiment. Now, suppose
that we want to compare a set of HMs
$\mathcal{M}=\\{\mathcal{M}_{1},\dots,\mathcal{M}_{J}\\}$ of the form given by
Eq. 4 that might differ in various ways (e.g., different prior/hyperprior
assumptions or disparate likelihoods). Due to the structure of the models, the
PMPs $p(\mathcal{M}\,|\,\\{\boldsymbol{x}_{mn}\\})$ depend on neither the
ordering of the individuals nor the ordering of their responses (which also
holds true for the corresponding BFs).

More precisely, if $\mathbb{S}(\cdot)$ is an arbitrary permutation of an index
set, then

$p(\mathcal{M}\,|\,\\{\boldsymbol{x}_{mn}\\})=p(\mathcal{M}\,|\,\mathbb{S}(\\{\boldsymbol{x}_{mn}\\}))$
(11)

for any $\mathbb{S}(\cdot)$ acting on
$\\{1\dots,M\\}\times\\{1,\dots,N_{1}\\}\times\dots\times\\{1,\dots,N_{M}\\}$
where $\times$ denotes the Cartesian product of two (index) sets. Note that
this notation implies that only permuting each $m$ and permuting each $n$
within, but not across each group $m$ is allowed. The property of permutation
invariance is immediately obvious from the right-hand side of Eq. 4 that
involves two nested products (products being permutation invariant
transformations when seen as functions operating on sets). Naturally, learning
permutation invariance directly from data or simulations is hardly feasible
with standard neural networks, even for non-nested data. Indeed, for non-
hierarchical generative models, [184] propose to use composite permutation
invariant networks as employed by [220]. In the following section, we
generalize this architectural concept to the hierarchical setting.

###  3.2 Hierarchical Invariant Neural Network Architecture

Permutation invariant networks differ from standard feedforward networks in
that they can process inputs of different sizes and encode the probabilistic
symmetry of the data directly (i.e., remove the need to learn the symmetry
implicitly during training by supervised learning alone).

For the purpose of BMC with HMs, we realize a hierarchical permutation
invariant function via a stack of invariant modules $\Sigma_{I}^{(l)}$ for
each hierarchical level $l=1,\dots,L$ of the Bayesian model (see Figure 2).
Each invariant module performs an equivariant non-linear transformation
$h^{(l)}_{1}$ acting on the individual data points, followed by a pooling
operator (e.g., sum or max) and a further non-linear transformation
$h^{(l)}_{2}$ acting on the pooled data.

In order to preserve hierarchical symmetry, we apply each $\Sigma_{I}^{(l)}$
independently to each nested sequence of data points. To make this point
concrete, consider the two-level model given by Eq. 4 and let data point
$\boldsymbol{x}_{mn}$ denote the multivariate response of person $m$ in trial
$n$ of some data collection experiment. Accordingly, the first invariant
module $\Sigma_{I}^{(1)}$ operates by reducing the trial data
$\\{\boldsymbol{x}_{n}\\}_{m}$ of each person $m$ to a single person-vector
$\widetilde{\boldsymbol{x}}_{m}$ of fixed size:

$\widetilde{\boldsymbol{x}}_{m}=\Sigma_{I}^{(1)}\left(\\{\boldsymbol{x}_{n}\\}_{m}\right)=h^{(1)}_{2}\left(\sum_{n=1}^{N_{m}}h^{(1)}_{1}(\boldsymbol{x}_{mn})\right),$
(12)

where $h_{1}$ and $h_{2}$ are implemented as simple feedforward neural
networks with trainable parameters suppressed for clarity. The second
invariant module $\Sigma_{I}^{(2)}$ then compresses all person vectors to a
final vector $\boldsymbol{z}$ of fixed size:

$\boldsymbol{z}=\Sigma_{I}^{(2)}\left(\\{\boldsymbol{x}_{m}\\}\right)=h^{(2)}_{2}\left(\sum_{m=1}^{M}h^{(2)}_{1}(\widetilde{\boldsymbol{x}}_{m})\right).$
(13)

In this way, the architecture becomes completely independent of the number of
persons $M$ or number of trials per person $N_{m}$, which could vary
arbitrarily across persons. The vector $\boldsymbol{z}$, whose dimensionality
represents a tunable hyperparameter, can be interpreted as encoding learned
summary statistics for the BMC task at hand (to be discussed shortly).
Moreover, it is easy to see that $\boldsymbol{z}$ is independent of the
ordering of persons or the ordering of trials within persons, as necessitated
by the model formulation in Eq. 4. Thus, the composition
$\Sigma_{I}^{(2)}\circ\Sigma_{I}^{(1)}(\\{\boldsymbol{x}_{mn}\\})$ reduces a
hierarchical data set with two levels to a single vector $\boldsymbol{z}$
which respects the probabilistic symmetry implied by the particular
hierarchical model formulation.

###  3.3 Increasing the Capacity of Invariant Networks

Encoding an entire hierarchical data set $\\{\boldsymbol{x}_{mn}\\}$ into a
single vector $\boldsymbol{z}$ forces the composite neural network to perform
massive data compression, creating a potential information bottleneck. For
complex generative models, this task can become rather challenging and will
depend highly on the representational capacity of the neural network (i.e.,
its ability to extract informative data set embeddings). Fortunately, we can
enhance the simple architecture described in the preceding paragraph by using
insights from [220] and [118].

In order to increase the capacity of the previously introduced invariant
transformation, we can stack together multiple equivariant modules
$\Sigma_{E}^{(l)}$. Each equivariant module implements a combination of
equivariant and invariant transformations. For instance, focusing on our two-
level model example (Eq. 4), the transformations at level $1$ for each person
$m$ are now given by:

$\displaystyle\widetilde{\boldsymbol{x}}_{m}$
$\displaystyle=h^{(1)}_{2}\left(\sum_{n=1}^{N_{m}}h^{(1)}_{1}(\boldsymbol{x}_{mn})\right)$
(14) $\displaystyle\widetilde{\boldsymbol{x}}_{mn}$
$\displaystyle=h_{3}^{(1)}(\left[\boldsymbol{x}_{mn},\widetilde{\boldsymbol{x}}_{m}\right])\text{\quad
for \quad}n=1,\dots,N_{m},$ (15)

where $h_{3}$ is also implemented as a simple feedforward neural network. In
this way, each intermediary output $\widetilde{\boldsymbol{x}}_{mn}$ of the
equivariant module now contains information from all data points, so the
network can learn considerably more flexible transformations. Moreover, we can
stack $K$ equivariant modules followed by an invariant module, in order to
obtain a deep invariant module, which for the first hierarchical level ($l=1$)
takes the following form:

$\widetilde{\boldsymbol{x}}_{m}=(\Sigma_{I}^{(1)}\circ\Sigma_{E}^{(K,1)}\circ\cdot\cdot\cdot\circ\Sigma_{E}^{(1,1)})(\\{\boldsymbol{x}_{n}\\}_{m}).$
(16)

Compared to the simple invariant module from Eq. 12, the deep invariant module
involves a larger number of computations but allows the network to learn more
expressive representations. Accordingly, the transformation for the second
hierarchical level ($l=2$), which yields the final summary representation
$\boldsymbol{z}$, is given by:

$\boldsymbol{z}=(\Sigma_{I}^{(2)}\circ\Sigma_{E}^{(K^{\prime},2)}\circ\cdot\cdot\cdot\circ\Sigma_{E}^{(1,2)})(\\{\widetilde{\boldsymbol{x}}_{m}\\}),$
(17)

where the number of equivariant modules $K^{\prime}$ for level $2$ can differ
from the number of equivariant modules $K$ for level $1$. In our experiments,
reported in Section 4, we observe a clear advantage of using deep invariant
networks over their simple counterparts. Furthermore, for two-level models, we
find that the performance of the networks is largely insensitive to the choice
of $K$ or $K^{\prime}$.

###  3.4 Learning the Model Comparison Problem

In order to get from the learned summary representation $\boldsymbol{z}$ to an
approximation of the analytic PMPs $\widehat{\boldsymbol{\pi}}$, we apply a
final neural classifier (i.e., the inference network)
$\mathcal{I}(\boldsymbol{z})=\widehat{\boldsymbol{\pi}}$, as visualized in
Figure 2. We deviate from the Dirichlet-based setting in [184], since we found
that implementing the inference network as a standard softmax classifier [142]
provides slightly better calibration and leads to more stable training in the
specific context of HMs.

Denoting the entire hierarchical neural network as
$f_{\boldsymbol{\phi}}(\\{\boldsymbol{x}\\})=\widehat{\boldsymbol{\pi}}$ and
an arbitrary hierarchical data set as $\\{\boldsymbol{x}\\}$, we aim to
minimize the expected logarithmic loss

$\min_{\boldsymbol{\phi}}\mathbb{E}_{p(\mathcal{M},\\{\boldsymbol{x}\\})}\left[-\sum_{j=1}^{J}\mathbb{I}_{\mathcal{M}_{j}}\cdot\log
f_{\boldsymbol{\phi}}(\\{\boldsymbol{x}\\})_{j}\right],$ (18)

where $\boldsymbol{\phi}$ represents the vector of trainable neural network
parameters (e.g., weights and biases), $\mathbb{I}_{\mathcal{M}_{j}}$ is the
indicator function for the “true” model. The expectation runs over the joint
generative (mixture) distribution of all models
$p(\mathcal{M},\\{\boldsymbol{x}\\})$, which we access through Monte Carlo
simulations. Since the logarithmic loss is a strictly proper loss [140], it
drives the outputs of $f_{\boldsymbol{\phi}}(\\{\boldsymbol{x}\\})$ to
estimate the actual PMPs $p(\mathcal{M}\,|\,\\{\boldsymbol{x}\\})$ as best as
possible. Thus, perfect convergence in theory guarantees that the network
outputs the analytically correct PMPs which asymptotically select the “true”
model in the closed world or the model that minimizes the Kullback-Leibler
divergence to the “true” data generating process in the open world [112].

In practice, we approximate Eq. 18 over a training set of $B$ simulations from
the competing HMs. Each entry $b$ for $b=1,...,B$ in this training set
represents a hierarchical data set $\\{\boldsymbol{x}^{(b)}\\}$ itself along
with a corresponding one-hot encoded vector for the “true” model index
$\mathcal{M}^{(b)}_{j}$. The latter denotes the model from which the data set
was generated and serves as the “ground truth” for supervised learning.

Similarly to [184], our neural method encodes an implicit preference for
simpler HMs (i.e., Occam’s razor) inherent in all marginal likelihood-based
methods [169, see]. Since our simulation-based training approximates an
expectation over the marginal likelihoods of all HMs
$p(\mathcal{M})\,p(\boldsymbol{x}\,|\,\mathcal{M})$, data sets generated by a
simpler HM will tend to be more similar compared to those generated by a more
complex one (cf. Figure 1). Thus, data sets that are plausible under both HMs
will be generated more often by the simpler model than by the more complex
model. A sufficiently expressive neural network will capture this behavior by
assigning a higher PMP for the simpler model111Assuming equal prior model
probabilities., thereby capturing complexity differences arising directly from
the generative behavior of the HMs.

Finally, to increase training efficiency when working under a limited
simulation budget, we also explore a novel pre-training method inspired by
transfer learning [205, 114]. First, we train the networks on data sets with a
reduced number of exchangeable units (e.g., reducing the number of
observations at level $l=1$). This procedure accelerates training since it
uses fewer simulator calls and the forward pass through the networks becomes
cheaper. In a second step, we generate data with a realistic number of
exchangeable units. Crucially, since we can use the pre-trained network from
step one as a better-than-random initialization, we need considerably fewer
simulations than if we trained the network from scratch. Indeed, our real-data
application in Experiment 4.3 demonstrates the utility of this training
method.

##  4 Experiments

In this section, we first conduct two simulation studies in which we
extensively test the approximation performance of our hierarchical neural
method. We start with a comparison of two nested toy HMs in Section 4.1,
followed by a comparison of two complex non-nested HMs of cognition in Section
4.2. For both validation studies, we test our method internally by examining
the calibration of the approximated PMPs. Additionally, we validate our method
externally by benchmarking its performance against the current state-of-the-
art for comparing HMs, namely, bridge sampling [137, 147]. To enable this
challenging benchmark, we limit our validation studies to the comparison of
models with explicit likelihoods to which bridge sampling is applicable.

Finally, in a real-data application, we use our deep learning method to
compare four hierarchical EAMs of response time data in Section 4.3. Two of
these models have no analytic likelihood, which makes the entire BMC setup
intractable with current state-of-the-art methods (e.g., bridge sampling).
Moreover, with this example, we also address the utility of a novel EAM, the
Lévy flight model [215], that has previously been impossible to investigate
directly using Bayesian HMs.

For all experiments, we assume uniform model priors $p(\mathcal{M}_{j})=1/J$.
All computations are performed on a single-GPU machine with an NVIDIA RTX 3070
graphics card and an AMD Ryzen 5 5600X processor. The reported computation
times are measured as wall-clock times. Details on the implementation of our
neural networks and the employed training procedures are provided in Appendix
A. Code for reproducing all results from this paper is freely available at
https://github.com/bayesflow-org/Hierarchical-Model-Comparison. Additionally,
our proposed method is implemented in the BayesFlow Python library for
amortized Bayesian workflows [187].

###  4.1 Validation Study 1: Hierarchical Normal Models

In this first experiment, we examine a simple and controllable model
comparison setup to examine the behavior of our method under various
conditions, before moving on to more complex scenarios. Inspired by [144], we
compare two hierarchical normal models $\mathcal{M}_{1}$ and
${\mathcal{M}_{2}}$ that share the same hierarchical structure

$\displaystyle\tau^{2}$ $\displaystyle\sim\text{Normal}_{+}(0,1)$ (19)
$\displaystyle\sigma^{2}$ $\displaystyle\sim\text{Normal}_{+}(0,1)$ (20)
$\displaystyle\theta_{m}$
$\displaystyle\sim\text{Normal}(\mu,\sqrt{\tau^{2}})\text{ for }m=1,\dots,M$
(21) $\displaystyle x_{mn}$
$\displaystyle\sim\text{Normal}(\theta_{m},\sqrt{\sigma^{2}})\text{ for
}n=1,\dots,N_{m},$ (22)

with Normal${}_{+}(\cdot)$ denoting a zero-truncated normal distribution. The
models differ with respect to the parameter $\mu$ that describes the location
of the individual-level parameters ${\theta_{m}}$: Whereas $\mathcal{M}_{1}$
assumes the location of ${\theta_{m}}$ to be fixed at $0$, the more flexible
$\mathcal{M}_{2}$ allows for $\mu$ to vary

$\displaystyle\mathcal{M}_{1}\text{: }\mu$ $\displaystyle=0$ (23)
$\displaystyle\mathcal{M}_{2}\text{: }\mu$
$\displaystyle\sim\text{Normal}(0,1).$ (24)

####  4.1.1 Calibration

The most important properties of an approximate inference method are the
trustworthiness of its results and, more pragmatically, whether we can
diagnose the lack of trustworthiness in a given application. A useful proxy
for trustworthiness is the calibration of a probabilistic classifier, which
measures how closely the predicted probabilities of outcomes match their true
underlying probabilities [151, 195].

However, computing the calibration of a BMC procedure is hardly feasible in a
non-amortized setting, since it involves applying the method to a large number
of simulated data sets. For bridge sampling, for example, that would imply re-
fitting the models via MCMC and running bridge sampling on at least hundreds,
if not thousands of simulated data sets. The calibration of our networks, on
the other hand, can be determined almost immediately after training due to
their amortization property [184].

In the following experiments, we assess the calibration of our networks
visually (via calibration curves) and numerically (via a measure of
calibration error). For generating a calibration curve [130, 181], we first
sort the predicted PMPs $\widehat{\pi}^{(s)}_{j}$ on $S$ simulated data sets
$s=1,\ldots,S$, which we then partition into $I$ equally spaced probability
bins $i=1,\ldots,I$ (we use $I=15$ bins for all validation experiments). For
each model $j$ and each bin $i$ containing a set $\mathcal{B}_{ij}$ of
predicted model indices, we compute the mean prediction for the model
(predicted probability, PP) and the actual fraction of this model being true
(true probability, TP) as follows:

$\displaystyle\text{PP}(\mathcal{B}_{ij})$
$\displaystyle:=\frac{1}{\lvert\mathcal{B}_{ij}\rvert}\sum_{b\in\mathcal{B}_{ij}}\widehat{\pi}^{(b)}_{j},$
(25) $\displaystyle\text{TP}(\mathcal{B}_{ij})$
$\displaystyle:=\frac{1}{\lvert\mathcal{B}_{ij}\rvert}\sum_{b\in\mathcal{B}_{ij}}\mathbb{I}_{\mathcal{M}^{(b)}_{j}},$
(26)

where $\mathbb{I}$ again denotes the indicator function for the “true model”.
These two quantities varying over the bins form the $X$\- and $Y$-axis of a
calibration curve. A well-calibrated model comparison method with an agreement
in each bin (as indicated by a diagonal line) thus yields approximations that
reflect the true probabilities of the compared models [151]. We further
summarize this information via the Expected Calibration Error [178, ECE;] as a
single number bounded between $0$ and $1$, which we estimate by averaging the
individual deviations between predicted and true probability in each bin:

$\widehat{\text{ECE}}_{j}:=\sum_{i=1}^{I}\frac{\lvert\mathcal{B}_{ij}\rvert}{S}\bigg{\lvert}\text{PP}(\mathcal{B}_{ij})-\text{TP}(\mathcal{B}_{ij})\bigg{\rvert}.$
(27)

If follows from Eq. 27 that a perfect ECE can be achieved by always predicting
indifferent probabilities (e.g., $\widehat{\pi}_{1}=\widehat{\pi}_{2}=.5$ when
comparing two models). We therefore complement our calibration assessment by
measuring the accuracy of recovery, for which we dichotomize the predicted
PMPs $\widehat{\pi}^{(s)}_{j}$ on $S$ simulated data sets into one-vs-rest
model predictions $\widehat{\mathcal{M}}^{(s)}_{j}$:

$\text{Acc}_{j}:=\frac{1}{S}\,\mathbb{I}_{\widehat{\mathcal{M}}^{(s)}_{j}=\mathcal{M}^{(s)}_{j}}.$
(28)

Thus, in our BMC context, accuracy roughly is to ECE what sharpness is to
posterior calibration in Bayesian parameter estimation [121, 124].222We focus
on the accuracy since we use a uniform model prior $p(\mathcal{M})$, but other
metrics of predictive performance, such as the logarithmic scoring rule, would
have been expedient as well.

(a) Median calibration curve and confidence intervals (CIs) for data sets of
$M=50$ groups with $N_{m}=50$ observations within each group.

(b) Median expected calibration errors (ECEs) and confidence intervals for
data sets of $M=50$ groups with differing numbers of observations $N_{m}$
within each group.

Figure 3: Validation study 1: Calibration results for (a) the neural network
trained on fixed data set sizes and (b) the neural network trained on data
sets with varying numbers of observations. Medians and confidence intervals
(CIs) are computed over $25$ repetitions.

##### Fixed data set sizes

In the first calibration experiment, we examine the performance of our method
for the most simple application case of learning a model comparison problem on
a specific (fixed) data set size. Here, all data sets simulated for training
and validating the network consist of $M=50$ groups and $N_{m}=50$
observations for each group $m=1,\dots,M$.

We train the network for $10,000$ backpropagation steps, taking $6$ minutes.
Subsequently, we calculate its calibration on $5,000$ held-out validation data
sets and repeat this process $25$ times to obtain stable results with
uncertainty quantification. 3(a) depicts the resulting median calibration
curve. Its close alignment to the dashed diagonal line (representing perfect
calibration) indicates that the PMP approximations are well-calibrated (median
ECE over all repetitions of $\widehat{\text{ECE}}=0.014$). The curve’s
coverage of the full range of predicted probabilities and the median accuracy
of $\widehat{\text{Acc}}=.89$ confirm that the excellent calibration does not
stem from indifferent predictions. The subsequent comparison of our method to
bridge sampling suggests that this accuracy is indeed close to the upper bound
imposed by the aleatoric uncertainty in the model-implied data.

##### Data sets with varying numbers of observations

We now train our hierarchical network to approximate BMC over a range of
hierarchical data sets with varying numbers of observations within groups
$N_{m}$. This amortization over observation sizes would provide a substantial
efficiency gain if a researcher desires to compare HMs on multiple data sets
with differing $N_{m}$, as only a single network would have to be trained for
all data sets.333Note that we refer to variability between data sets. We
describe an approach for handling within data set variability of nested trials
in Section 4.3. In our validation setup, each simulated data set still
consists of $M=50$ groups, but now the number of observations within those
groups varies in $N_{m}=1,\ldots,100$.

We train the network for $20,000$ training steps, taking $13$ minutes. At each
training step, we draw the number of observations for the current batch of
simulations from a discrete uniform distribution
$N_{m}\sim\text{Uniform}_{D}(1,100)$. For each $N_{m}$ used during training,
we evaluate the calibration $25$ times on $5,000$ held-out simulated
validation data sets. This repetition procedure allows us to quantify the
uncertainty of our ECE estimates.

3(b) plots the median ECE values for each observation size. The neural network
achieves high calibration with a median ECE over all observation sizes (and
repetitions) of $\widehat{\text{ECE}}=0.012$. Moreover, the unsystematic
pattern of the median curve and the homoscedastic variation between the
observation sizes indicate that the network has learned the model comparison
task equally well for all settings (with the ECE only rising slightly for the
poorly identifiable $N_{m}=1$ setting). Together, the low calibration error
and the accurate model predictions (median accuracy
$\widehat{\text{Acc}}=.88$) indicate that our method incurs no trade-off
between calibration and accuracy. We additionally observe no bias towards a
model in all but the smallest observation sizes (see Figure 9 for accuracy and
bias examinations in all settings).

(a) Expected Calibration Error (ECE).

(b) Accuracy of recovery.

Figure 4: Validation study 1: Results for the neural network trained and
tested over variable data set sizes.

##### Data sets with varying numbers of groups and observations

In the third calibration experiment, we test the ability of the network to
learn a model comparison problem over a range of data sets with varying
numbers of groups $M$ and varying observations per group $N_{m}$. This
training scheme allows for amortized model comparison on multiple data sets
with different sizes, which can be especially useful for a priori sample size
determination on simulated data. Additionally, the trained network can be
stored and reused on future data sets with yet-unknown sample sizes. For this
experiment, training and validation data sets are simulated with
$M=1,\ldots,100$ groups and $N_{m}=1,\ldots,100$ observations, resulting in a
vast variability of data set sizes between $1$ up to $10,000$ data points.

Given the complexity of the learning task, we now train the network for
$40,000$ training steps, taking $36$ minutes. At each training step, we draw
the number of groups and observations from discrete uniform distributions
$M\sim\text{Uniform}_{D}(1,100)$ and $N_{m}\sim\text{Uniform}_{D}(1,100)$. We
estimate calibration on $5,000$ held-out simulations for each combination of
$M$ and $N_{m}$. As this implies simulating $50,000,000$ data sets, we forego
the repetition procedure employed in the previous experiments.

Figure 4 depicts the calibration and accuracy results for all combinations of
$M$ and $N_{m}$. We observe low ECEs for the vast majority of settings in 4(a)
(median ECE over all settings of $\widehat{\text{ECE}}=0.013$). In other
words, the trained network is capable of generating highly calibrated PMPs
over a broad range of data set sizes. Moreover, the BMC results are sensitive
to the number of nested observations $N_{m}$, but not to the number of groups
$M$, in our experimental setups. The only systematic drop in calibration
occurs for data sets containing just a few nested observations ($N_{m}\leq
5$). Considering that we observed better calibration for this low number of
observations in a network trained on data sets with varying $N_{m}$ (see
3(b)), we surmise that the drop in the edge areas in 4(a) arises from the
challenging learning task over vastly different data set sizes [128, a
phenomenon known as amortization gap;]. The overall low (i.e., good) ECEs for
all cases but the poorly identifiable $N_{m}=1$ setting suggest that the
networks’ approximations are generally trustworthy. This is further confirmed
by 4(b), where the observable accuracy pattern assures that this high
calibration does not arise from a trade-off with predictive performance.
Despite the demanding amortization setting, the network achieves an excellent
median accuracy of $\widehat{\text{Acc}}=.88$, similar to the earlier
experiments. We also find no indication of bias in any of the test settings
except the $N_{m}=1$ setting (see Figure 10). Marginal diagnostic plots for
all metrics are provided in Figure 11.

####  4.1.2 Bridge Sampling Comparison

After validating the general trustworthiness of our method, we now benchmark
it against the current gold standard for comparing HMs, namely, bridge
sampling, as implemented by [147]. As the non-amortized nature of bridge
sampling restricts the feasible number of test sets, we conduct the
benchmarking on $100$ test sets which are simulated equally from
$\mathcal{M}_{1}$ and $\mathcal{M}_{2}$. All simulated data sets consist of
$M=50$ groups and $N_{m}=50$ observations per group. The fixed sample sizes of
the test sets allow us to compare the two most distinct networks from Section
4.1.1 to bridge sampling: The fixed network that is trained for this specific
sample size and the more complex variable network that is trained for
amortized model comparison over variable sample sizes between $M=1,\ldots,100$
groups and $N_{m}=1,\ldots,100$ observations per group.

For bridge sampling, we first run four parallel MCMC chains with a warm-up
period of $1,000$ draws and $49,000$ post-warm-up posterior draws per chain in
Stan [122, 201]. We assess convergence through a visual inspection of the MCMC
chains and an assessment of the $\widehat{R}$, bulk ESS and tail ESS metrics
[213]. Afterwards, we use the posterior draws to approximate PMPs and BFs with
the bridgesampling R package [147]. We confirm the sufficiency of the total of
$196,000$ posterior draws by assessing the variability between multiple runs
as in [195], which yields highly similar results. Further insights via our
calibration diagnostics are precluded by bridge sampling being a non-amortized
method.

##### Approximation performance

Table 1: Validation study 1: Performance metrics for the comparison between
hierarchical normal models.

| Accuracy | MAE | RMSE | Log-Score | SBC
---|---|---|---|---|---
Bridge sampling | 0.86 (0.03) | 0.19 (0.03) | 0.32 (0.03) | 0.32 (0.06) | -0.02 (0.04)
Fixed network | 0.84 (0.04) | 0.19 (0.03) | 0.32 (0.03) | 0.32 (0.06) | -0.01 (0.04)
Variable network | 0.86 (0.03) | 0.19 (0.03) | 0.32 (0.03) | 0.31 (0.06) | -0.01 (0.04)

33footnotetext: Note. Bootstrapped mean values and standard errors (in
parentheses) are presented. We use $1000$ bootstrap versions of the test data
sets and estimate the standard errors from the bootstrap standard deviations
of the metrics.

As we compare approximate PMPs, we can use a number of complementary metrics
commonly employed to evaluate the quality of probabilistic predictions. First,
we quantify the fraction of times the correct model $\mathcal{M}^{(s)}_{j}$
underlying a simulated data set $s$ was detected, that is, the accuracy of
recovery (see Equation 28). Second, we assess the Mean Absolute Error (MAE) to
investigate the average deviation of the approximated model probabilities
$\widehat{\pi}^{(s)}_{j}$ from a perfect classification:

$\text{MAE}_{j}:={\frac{1}{S}\sum_{s=1}^{S}\left|\widehat{\pi}^{(s)}_{j}-\mathbb{I}_{\mathcal{M}_{j}}^{(s)}\right|}.$
(29)

Third, we measure the Root Mean Squared Error (RMSE), which places particular
emphasis on large prediction errors, to detect whether one method produces
highly incorrect approximations more frequently than the other:

$\text{RMSE}_{j}:=\sqrt{\frac{1}{S}\sum_{s=1}^{S}\left(\widehat{\pi}^{(s)}_{j}-\mathbb{I}_{\mathcal{M}_{j}}^{(s)}\right)^{2}}.$
(30)

Fourth, we calculate the Log-Score following the logarithmic scoring rule:

$\text{LogScore}_{j}:=-\frac{1}{S}\sum_{s=1}^{S}\left[\mathbb{I}_{\mathcal{M}_{j}}^{(s)}\cdot\text{log}\widehat{\pi}^{(s)}_{j}\right].$
(31)

Its property as a strictly proper scoring rule implies that it is
asymptotically minimized if and only if the approximate probabilities equal
the true probabilities [140]. Lastly, we measure simulation-based calibration
[203, SBC;] as adapted by [195] for model inference by the difference between
the prior probability for a model and its average posterior probability in the
test sets:

$\text{SBC}_{j}:=p(\mathcal{M}_{j})-\frac{1}{S}\sum_{s=1}^{S}\widehat{\pi}^{(s)}_{j}.$
(32)

We evaluate all metrics for $\mathcal{M}_{2}$, so that a bias towards
$\mathcal{M}_{1}$ is indicated by positive SBC values and a bias towards
$\mathcal{M}_{2}$ by negative SBC values.

Table 1 depicts the comparison results for our experimental setting. All
metrics show equal performances for bridge sampling and the two neural network
variants, with any differences being well within the range of the standard
errors.

##### Approximation convergence

In the following, we analyze the degree of convergence between the two methods
at the level of individual data sets. We explore this visually by contrasting
the PMP and (natural logarithmic) BF approximations of bridge sampling with
the two neural network variants in Figure 5. We observe that the two methods’
PMP approximations agree for the easy cases where the true underlying model is
clearly classifiable. Thus, discrepancies between the two methods arise mainly
for data sets with predicted PMPs close to $\widehat{\pi}=0.5$. Even for the
data sets with the largest discrepancies, the two methods do not map to
qualitatively different decisions: $\widehat{\pi}^{(\text{bridge})}_{2}=.67$
and $\widehat{\pi}^{(\text{neural})}_{2}=.79$ for the fixed network,
$\widehat{\pi}^{(\text{bridge})}_{2}=.32$ and
$\widehat{\pi}^{(\text{neural})}_{2}=.25$ for the variable network. Most
importantly, we detect no systematic pattern in these deviations.

Figure 5: Validation study 1: Comparison of approximation results obtained via
bridge sampling vs. the neural network trained on fixed data set sizes (left)
and the neural network trained on variable data set sizes (right). For
visibility purposes, the Bayes factor plots include only those 73 data sets
for which bridge sampling approximated a $\mathrm{BF}_{21}<1,000,000$ (plots
with all data sets are provided in Appendix B).

As BFs represent the ratio of marginal likelihoods, they allow for a closer
inspection of the degree of agreement between the methods in those edge cases
with PMPs close to $0$ or $1$. We observe a close convergence for data sets
classified as stemming from $\mathcal{M}_{1}$. Considering the predictions
favoring $\mathcal{M}_{2}$, there are discrepancies for data sets with log BFs
$>9.49$. Since this corresponds to BFs $>13,000$ and PMPs $>.9999$, it is not
visible in the PMP approximation plots. We obtain such extreme results only
for $\mathcal{M}_{2}$, as this model allows for deviations of the group level
parameters’ location from $0$ and enables the occurrence of extreme evidence
in its favor. The divergence in this area of extreme evidence emerges most
likely from the loss function employed for training the neural networks: The
logarithmic loss obtained from a minuscule deviation of the PMP from $1$ is
near $0$, which results in a negligible incentive for further optimization of
the network’s weights. We could reject a competing explanation based on
limited floating-point precision, since training with an increased floating-
point precision from $32$-bit to $64$-bit resulted in identical patterns. For
visibility purposes, we exclude the 27 data sets for which bridge sampling
approximated a BF $>1,000,000$ for the BF plots in Figure 5, all continuing
the observed plateau pattern. Plots with all 100 data sets are provided in
Appendix B.

The divergence we encountered provides insights into the technical nature of
our method but only arises in cases of extreme evidence. Thus, it is far from
altering the substantive conclusions derived from the simulated BMC setting.
Considering the convergence between the two methods in the realm of practical
relevance, we can conclude that our method produces highly similar
approximations to bridge sampling in this scenario.

##### Approximation time

Both bridge sampling and our deep learning method can be divided into two
computational phases. For bridge sampling, the first phase consists of drawing
from the posterior parameter distributions (taking $52$ seconds per data set
on average). Bridge sampling itself takes place in the second phase (taking
$38$ seconds on average). Notably, in contrast to amortized inference with
neural networks, both phases need to be repeated for each (simulated or
observed) data set. Taking the initial compilation time of $42$ seconds into
account, bridge sampling consequently took $152$ minutes for BMC on our $100$
test data sets.

For the neural networks, the first phase (training) is resource-intensive
(taking $6$ minutes for the fixed network and $36$ minutes for the variable
network). The second phase (inference) is then performed in near real-time
(inference on all $100$ test data sets took $0.0004$ seconds for the fixed
network and $0.007$ seconds for the variable network) and thus amortizes the
training cost over multiple applications. For the simple HMs compared here,
the amortization gains of our networks over bridge sampling come into effect
after performing BMC on $4$ (fixed network) or $24$ (variable network) data
sets.

We acknowledge our likely suboptimal choices of computational steps for the
bridge sampling workflow or the neural networks and hence wish to stress the
general patterns of non-amortized vs. amortized methods demonstrated here. In
general, we expect an advantage of bridge sampling in terms of efficiency in
situations where only one or a few data sets are available and obtaining a
large number of posterior draws is feasible. The demonstrated amortization
property of our method might not be so relevant for inference on a single
hierarchical data set, but it becomes crucial for performing calibration or
recovery studies, which necessitate multiple re-fits of the same model [195].

(a)

(b)

(c)

Figure 6: Validation study 2: Results for the comparison between hierarchical
SDT and MPT models. Table 2: Validation study 2: Performance metrics for the
comparison between hierarchical SDT and MPT models.

| Accuracy | MAE | RMSE | Log Score | SBC
---|---|---|---|---|---
Bridge sampling | 0.95 (0.02) | 0.1 (0.02) | 0.22 (0.03) | 0.16 (0.04) | -0.01 (0.04)
Neural network | 0.95 (0.02) | 0.1 (0.02) | 0.21 (0.03) | 0.16 (0.04) | 0.00 (0.04)

33footnotetext: Note. Bootstrapped mean values and standard errors (in
parentheses) are presented. We use $1000$ bootstrap versions of the test data
sets and estimate the standard errors from the bootstrap standard deviations
of the metrics.

###  4.2 Validation Study 2: Hierarchical SDT vs. MPT Models

We now extend our validation experiments from the simple setup with nested HMs
to the comparison of non-nested HMs of cognition. In this simulation study, we
examine the ability of our method to distinguish between data sets generated
either from an HM based on signal detection theory [143, SDT model;] or a
hierarchical multinomial processing tree model [189, MPT model;]. For
illustrative purposes, we embed our simulation study within an old-new
recognition scenario, where participants indicate whether or not a stimulus
was previously presented to them.

We ensure a challenging model comparison setting via three design aspects:
First, we specify both models to possess a similar generative behavior, that
is, hardly distinguishable prior predictive distributions of hit rates and
false alarm rates (prior predictive plots are provided in Appendix C). Second,
data sets of old-new recognition typically contain low information as they
only consist of binary variables indicating the stimulus type and response,
respectively. Third, we further amplify the information sparsity of the data
sets by choosing a particularly small size for all data sets of $M=25$
simulated participants and $N_{m}=50$ observations per participant.

A major difference between the compared cognitive model classes lies in the
assumption of a continuous latent process by the SDT model and discrete
processes (or states) by the MPT model. Our specification of the SDT model
follows the hierarchical formulation of the standard equal-variance model by
[190]. As the competing MPT model, we specify a hierarchical latent-trait two-
high-threshold model [164], which, in contrast to the SDT model, explicitly
models correlations between its parameters. We follow the convention of
restricting the parameters that describe the probability of recognizing a
previously presented stimulus as old and a distractor stimulus as new to be
equal, $D_{O}=D_{N}$, to render the MPT model identifiable [199, 133]. Our
prior choices for the parameters of both models are described in Appendix C.

We train the neural network for $50,000$ training steps. As in Section 4.1, we
first leverage the amortization property of our method to inspect its
calibration for the current model comparison task.
LABEL:fig:calibration_SDTMPT shows that the trained neural network generates
well-calibrated PMP approximations (median ECE over $25$ repetitions of
$\widehat{\text{ECE}}=0.009$).

Next, we assess whether the observed calibration of the network translates
into a competitive performance relative to bridge sampling. The benchmarking
setup ($50$ simulated data sets from each model) and the implementation of the
bridge sampling workflow follow the procedure described in Section 4.1.2.

The classification metrics depicted in Table 2 reveal the excellent
performance of both methods, despite the challenging BMC scenario. We further
observe a high degree of convergence between approximate PMPs derived by the
two methods (cf. LABEL:fig:BS_vs_NN_SDTMPT_PMPs). Again, we find discrepancies
between bridge sampling and our method in areas of extreme evidence (see
Figure 14 for log BFs). As depicted in LABEL:fig:BS_vs_NN_SDTMPT_comp_times,
obtaining PMP approximations for the $100$ test data sets took more than $6$
hours for bridge sampling and $55$ minutes for the neural network. For this
comparison of more complex cognitive models than in Section 4.1.2, the
amortization advantage of our method emerges when analyzing $15$ or more data
sets. Note that this advantage would quickly show up in validation studies
involving multiple model re-fits (e.g., bootstrap, sensitivity analysis or
cross-validation).

All validation experiments so far have been set up in an $\mathcal{M}$-closed
setting, with validation data simulated from the set $\mathcal{M}$ of models
under consideration [116]. Therefore, as a final validation, we test whether
our method also behaves sensibly in an $\mathcal{M}$-open setting, where none
of the models generated the test data. For this, we simulate $100$ noise data
sets with the same hierarchical structure as before but generate the binary
values for stimulus types and responses from a Bernoulli distribution with
$p=0.5$. Our neural method agrees with bridge sampling by assigning very high
PMPs to the SDT model for all noise data sets
($\bar{\pi}^{(\text{bridge})}_{SDT}=.999965$;
$\bar{\pi}^{(\text{neural})}_{SDT}=.999958$). Correspondingly, the deviations
between both methods are minimal. We thus observe a close alignment between
bridge sampling and our neural method in both a well-specified and a
misspecified scenario. This tentative result suggests that our amortized
estimates are faithful approximations not only in an $\mathcal{M}$-closed but
also an $\mathcal{M}$-open setting, at least for this BMC scenario.

The converging results from the two validation studies demonstrate that our
neural method generates well-calibrated and accurate PMP approximations.
Despite our method only accessing the likelihood function indirectly via
simulations, it can successfully compete with bridge sampling, which has
direct access to the likelihood function.

###  4.3 Application: Hierarchical Evidence Accumulation Models

In the following, we showcase the utility of our method by comparing complex
hierarchical EAMs in a real-data situation where likelihood-based methods such
as bridge sampling would not be applicable. More precisely, we seek to test
the explanatory power of different stochastic diffusion model formulations
proposed by [215] for experimental response time data.

The so-called Lévy flight model increases the flexibility of the standard
Wiener diffusion model [188] but renders its likelihood function intractable
with standard numerical approximations [216]. The complete incorporation of
all information through hierarchical modeling and the realization of BMC has
consequently been infeasible so far. Thus, in a recent study, [219] had to
resort to a separate computation of the Bayesian Information Criterion (BIC)
for each participant with subsequent aggregation. We aim to extend the study
of [219] by comparing fully hierarchical EAMs through PMPs and BFs. Moreover,
we intend to answer the question formulated by [219] as to whether the
superior performance of the more complex models in their study stems from an
insufficient punishment of model flexibility by the BIC. In addition to
addressing a substantive research question in this application, we also
demonstrate multiple advantages of our deep learning method on empirical data:

  * •

Compare HMs with intractable likelihoods: As our method is simulation-based,
including models with intractable likelihood functions in the comparison set
does not alter its feasibility.

  * •

Adequately model nested data: Our method alleviates computational challenges
that prevent modelers from adequately capturing the information contained in
nested data structures through HMs.

  * •

Re-use trained networks via fine-tuning: We accelerate the training of our
neural network by pre-training it on less complex simulated data and
subsequently fine-tuning it on simulated data resembling the actual
experimental setting.

  * •

Handle missing data: We train a neural network that can handle varying amounts
of missing data by randomly masking simulated data during the training
process.

  * •

Validate a trained network on simulated data: The amortized nature of our
method allows for extensive validation of a trained network prior to its
application to empirical data.

(a)

(b)

Figure 7: Real-data application: Validation results for the evidence
accumulation models on $2,000$ simulated data sets per model.

####  4.3.1 Model Specification

For this application, we consider a Lévy flight model with non-Gaussian noise
[215]. The Lévy flight process is driven by the following stochastic ordinary
differential equation:

$\displaystyle dx$ $\displaystyle=v\,dt+\sigma d\xi$ (33) $\displaystyle\xi$
$\displaystyle\sim\text{AlphaStable}(\alpha,\mu=0,\sigma=\frac{1}{\sqrt{2}},\beta=0),$
(34)

which represents a Lévy walk characterized by a fat-tailed stable noise
distribution. 444An earlier version of this work used the original formulation
by [215], which sets $\sigma=1$. For the special case of $\alpha=2.0$, which
is equivalent to the Wiener diffusion model, $\sigma=1$ leads to an unusual
diffusion constant (standard deviation of Gaussian noise) of $\sqrt{2}$,
whereas $\sigma=\frac{1}{\sqrt{2}}$ ensures the conventional diffusion
constant of $1$. Notably, model comparison results are highly sensitive to the
choice of $\sigma$. In the above equation, $x$ denotes the accumulated
(perceptual) evidence, $v$ denotes the rate of accumulation and $\alpha$
controls the tail exponent of the noise variate $\xi$. [215] and [219] argue
that the more abrupt changes in the information accumulation process that this
model allows for could provide a better description of human decision-making
than a Gaussian noise. The addition of Lévy noise renders the standard
numerical approximation of the diffusion model likelihood intractable [216].
Consequently, neither standard MCMC nor bridge sampling are applicable for
Bayesian parameter estimation and BMC, respectively.

There is an ongoing debate about the inclusion of additional parameters that
account for inter-trial variability in the diffusion model parameters: While
they can provide a better model fit, the estimation of inter-trial variability
parameters is often difficult and can result in unstable results [119, 166].
Thus, [219] also compared basic (without inter-trial variability parameters)
to full (with inter-trial variability parameters) versions of the drift-
diffusion and Lévy flight model.

Consequently, the set of candidate models considered here consists of four
EAMs with increasing flexibility (i.e., the scope of possible data patterns
that they can generate):

  * •

$\mathcal{M}_{1}$, the most parsimonious basic diffusion model with the
parameter $v$ describing the mean rate of information uptake, the parameter
$a$ describing the threshold at which a decision is made, the parameter
$z_{r}$ describing a bias of the starting point towards one decision
alternative and the parameter $t_{0}$ describing the non-decision time, that
is, the time spent encoding the stimulus and executing the decision.

  * •

$\mathcal{M}_{2}$, the basic Lévy flight model, in which the assumption of a
Wiener diffusion process with Gaussian noise is replaced by the above
introduced Lévy flight process. The additional free parameter $\alpha$ governs
the tail behavior of the noise distribution. The setting $\alpha=2$ is
equivalent to a Gaussian distribution, whereas $\alpha=1$ reduces to a Cauchy
distribution.

  * •

$\mathcal{M}_{3}$, the full diffusion model, which extends $\mathcal{M}_{1}$
with the parameters $s_{v_{m}}$, $s_{z_{m}}$ and $s_{t_{m}}$ that denote the
variability (i.e., standard deviations) of drift rate, starting point bias and
non-decision time, respectively, between trials.

  * •

$\mathcal{M}_{4}$, the full Lévy flight model, that possesses the largest
flexibility by including inter-trial variability parameters as well as the
flexible Lévy noise distribution controlled by $\alpha$.

Parameter priors and prior predictive checks are provided in Appendix D.1.

####  4.3.2 Data

The reanalyzed data set by [219] contains $40$ participants who completed a
total of $900$ trials of binary decision tasks (color discrimination and
lexical decision) each. On average, $3.17\%$ of trials per participant were
excluded due to extremely short or long reaction times.

Figure 8: Real-data application: Model posteriors on the empirical data set
with uncertainty under different data perturbations. We use $100$ bootstrap
samples for the bootstrapped results.

####  4.3.3 Simulation-Based Training

Since simulating data from EAMs can be challenging, especially when they
include non-Gaussian noise, we leverage the advantage that neural networks are
capable of transfer learning as described in Section 3.4. Transfer learning
describes the utilization of representations that had been previously learned
by a neural network in a particular task for a new, related task [179, e.g.,].
In this way, neural networks can be applied in small data settings (e.g., a
limited simulation budget) by re-using the training knowledge encoded from
structurally similar (possibly big data) settings.

For the purpose of model comparison, we first pre-train the network for $20$
epochs (passes over the whole training data) on $10,000$ simulated data sets
per model. These data sets resemble the empirical data in that they consist of
$40$ simulated participants, but differ in that the number of trials is
reduced by a factor of 9 ($100$ instead of $900$ trials per participant).
Afterwards, we fine-tune the network for additional $30$ epochs on $2,000$
simulated data sets per model that match the empirical data set with $40$
simulated participants and $900$ trials per participant. Thereby, we
considerably reduce the computational demand of the training process. We
further speed up the training phase by simulating all data prior to the
training of the network in the high-performance programming language Julia
[117]. Pre-training took 10 minutes for the simulations and 11 minutes for
training the networks. Fine-tuning took 18 minutes for the simulations and 16
minutes for training the networks, resulting in a total of 55 minutes for the
training phase.

To fully adapt the network to the characteristics of the empirical data, we
also simulate missing data during fine-tuning. In each training epoch, we
generate a random binary mask $\boldsymbol{f}$ coding the simulated missing
values. We sample the number of masked trials from a (discretized) normal
distribution truncated between $1$ and the number of trials, $900$. The
distributions’ mean and standard deviation match the amount and variability of
missing trials in the empirical data. We then perform an element-wise
multiplication $\tilde{\boldsymbol{x}}=\boldsymbol{x}\otimes\boldsymbol{f}$
and feed the “contaminated” data $\tilde{\boldsymbol{x}}$ to the network. This
procedure results in a robust network that can process various proportions of
missing data. We find rank stability of our results in the presence of up to
$25\%$ missing data in Appendix D.2.

####  4.3.4 Results

Before applying our trained network to the empirical data, we validate it on
$2,000$ simulated data sets per model. First, the individual calibration
curves in LABEL:fig:levy_validation_confusion show excellent calibration for
all models with $\widehat{\text{ECEs}}$ close to $0$. The calibration curves
now consist of $10$ instead of $15$ intervals to obtain stable results despite
the smaller amount of validation data sets per model. Second, we evaluate the
accuracy of recovery and patterns of misclassification through the confusion
matrix depicted in LABEL:fig:levy_validation_calibration. The confusion matrix
confirms that the excellent calibration of the network does not stem from
chance performance. It also reveals that the selection of the “true” model
becomes more difficult with increasing model complexity, which is a direct
consequence of the Occam’s razor property inherent in BMC (cf. Figure 1).

Table 3: Real-data application: Bayes factors (BFs) and posterior model probabilities (PMPs) estimated from data by [219]. The preferred model is indicated by an asterisk. | $\mathcal{M}_{1}$ | $\mathcal{M}_{2}$ | $\mathcal{M}_{3}$ | $\mathcal{M}_{4}$
---|---|---|---|---
$\text{BF}_{j3}$ | 9.63e-05 | 0.01 | * | 0.27
$\text{BF}_{3j}$ | 1.04e+04 | 78 | * | 3.71
PMP | 7.51e-05 | 1.00e-02 | 0.78 | 0.21

Table 3 presents the model comparison results on the empirical data set.
Additionally, Figure 8 displays the model posteriors under different data
perturbations. Consistent with the results of the non-hierarchical BIC
approach by [219], we find little evidence for both the basic diffusion model
$\mathcal{M}_{1}$ and the basic Lévy flight model $\mathcal{M}_{2}$. This
implies that the additional complexity of allowing parameters to vary between
trials in $\mathcal{M}_{3}$ and $\mathcal{M}_{4}$ is, even under the strict
penalization of prior-predictive flexibility in BMC, outweighed by better
model fit. Also in agreement with [219], we observe evidence for both
$\mathcal{M}_{3}$ and $\mathcal{M}_{4}$, but, in contrast to [219], our
results slightly favor the full diffusion model $\mathcal{M}_{3}$ over the
full Lévy flight model $\mathcal{M}_{4}$. Figure 8 confirms both the slight
advantage of $\mathcal{M}_{3}$ over $\mathcal{M}_{4}$ and the substantial
uncertainty associated with these results.

##  5 Discussion

Nested data are ubiquitous in the quantitative sciences, including
psychological and cognitive research [134]. Yet, to avoid dealing with the
complex dependencies resulting from these data, researchers often resort to
simpler analyses, ignoring potentially important structural information.
Hierarchical models (HMs) provide a flexible way to represent the multilevel
structure of nested data, but this flexibility can make Bayesian model
comparison a daunting undertaking.

In this work, we proposed a powerful remedy to this problem: Building on the
BayesFlow framework [186], we developed a neural network architecture that
enables approximate BMC for arbitrarily complex HMs. In two simulation
studies, we showed that our deep learning method is well-calibrated and
performs as accurately as bridge sampling, which is the current state-of-the-
art for comparing HMs with simple likelihoods. Moreover, in a subsequent real-
data application, we compared the relatively new Lévy flight model with
existing evidence accumulation models. Thus, we argue that our method is well-
suited to enhance the applicability of (complex) HMs in psychological
research. Below, we summarize the key properties and limitations of our method
while also outlining future research directions.

###  5.1 Amortized Inference

Our method offloads the computational demands for comparing HMs onto the
training phase of a custom neural network, allowing for near real-time model
comparison using the trained network. The resulting amortization offers
several advantages over non-amortized methods.

First, it enables thorough validation of a trained network on thousands of
simulated data sets, allowing large-scale simulation-based diagnostics to
become an integral part of the BMC workflow [195, 138]. Second, the trained
and validated network can be used not only for point estimates of BFs or PMPs
on empirical data but also for exploring the robustness of the results against
multiple data perturbations, as showcased in our real-data application.

Third, we demonstrated the feasibility of amortizing over variable data set
sizes in our first validation study. This is particularly advantageous in the
context of HMs since nested data sets often contain multiple exchangeable
levels with variable sizes (e.g., different numbers of clusters, participants
and observations). Analyzing multiple hierarchical data sets with variable
sizes only requires a single network that has seen different data set sizes
during training. The same network could also be used for various simulation
studies, such as the challenging task of designing maximally informative
experiments in a hierarchical BMC setting [177, 156].

Lastly, we showed that researchers do not even need to consider all possible
shapes of future data sets when training such a network, as they can use
transfer learning to efficiently adapt a trained network to a related setting.
Beyond allowing more flexibility in reusing networks across experiments,
researchers or even fields, transfer learning can also considerably reduce the
computational demands associated with comparing complex HMs. As demonstrated
in our real-data application, a network can be pre-trained on simulated data
sets with reduced size and fine-tuned afterwards on sizes matching the
empirical data.

###  5.2 Independence From Explicit Likelihoods

Unlike other popular methods for performing BMC on HMs, such as the Savage-
Dickey density ratio or bridge sampling, our method is not constrained by the
availability of an explicit likelihood function for all competing models. As
long as the models in question can be implemented as simulators, the neural
network can be trained to perform BMC on these models. The value of such a
method is evident, as it decouples the substantive task of model specification
from concerns about the feasibility of estimation methods.

Statistical models are instantiations of substantive knowledge or hypotheses.
As such, we argue that model specification should not be unduly restricted by
considerations of computational tractability – a sentiment that is closely
related to what [152] call the “specification-first-principle”. Our proposed
deep learning method satisfies this principle, as model specification may be
guided exclusively by substantive arguments with few concerns about
tractability. Thus, we believe that our method makes a contribution to the
recent upsurge of innovative psychological models [139, 154, 125] by allowing
for an efficient assessment of their incremental value in a hierarchical
setting.

###  5.3 Limitations and Outlook

One of the main challenges of approximate methods and, more broadly,
statistical inference is ensuring the faithfulness of the obtained results.
The outlined possibilities for validating the network and examining the
robustness of the results are important contributions of our method but come
with open questions. Concerning the validation of the network, framing model
comparison as a supervised learning problem allows us to draw from the rich
literature on classification performance metrics. Nevertheless, determining a
“good-enough” score for an approximate BMC method remains challenging, as the
optimally possible performance is application-specific and usually unknown.

Concerning the application of the network to empirical data, we showed in
Validation Study 2 that our method produces, at least in this scenario,
reasonable results when confronted with data not stemming from the models
under consideration. Moreover, our robustness checks are a practical proxy for
measuring the reliability of BMC results in a closed-world setting. However,
these checks cannot possibly capture the (lack of) absolute evidence for an
HM: As a relative method, BMC may indicate that one model fits the data
_better_ than a set of competing models, but it does not provide any measure
of how well (or poorly) the model itself approximates the underlying data-
generating process. A promising direction to address this limitation could be
the combination of our method with the recently proposed meta-uncertainty
framework for BMC [196], which can be greatly accelerated with amortized deep
learning methods. This combination could provide a principled delineation of
different uncertainty sources, enabling the detection of _model
misspecification_ cases where none of the competing HMs can explain the
observed data. Still, further research is needed to determine whether meta-
uncertainty can provide reliable evidence for the open vs. closed world
assumption in the context of HMs and prevent the dangers that simulation gaps
(i.e., as induced by model misspecification) pose for simulation-based
inference [158].

Since BMC is a marginal likelihood (i.e., prior predictive) approach, the
priors should be informed by scientific theory and will thus have a decisive
influence on the results [212]. We do not intend to re-iterate the ongoing
discussion about this property of BMC [152, 214, 148, 149], but want to
highlight a specific difficulty that arises for HMs: Parameter priors of an HM
are connected via multilevel dependencies, increasing the risk that poor prior
choices lead to non-intended model behavior [193, for a recent discussion of
this problem in cognitive modeling, see]. Therefore, prior predictive checks
and prior sensitivity analyses become especially important when conducting BMC
on competing HMs. While transfer learning reduces the computational demands of
retraining a neural network for sensitivity analyses, another avenue for
future research would be the amortization over different prior choices,
enabling immediate prior sensitivity assessment.

Finally, it should be noted that the version of our method explored here can
only compare HMs assuming exchangeable data at each hierarchical level.
Although the majority of HMs in social science research follow this
probabilistic symmetry, some researchers may want to compare non-exchangeable
HMs, for example, to study within-person dynamics [197, 167, 132].
Fortunately, the modularity of our method allows easy adaptation of the neural
network architecture to handle non-exchangeable HMs. To compare hierarchical
time series models with temporal dependencies at the lowest level, for
instance, the first invariant module could be exchanged for a recurrent
network, as proposed in [184]. Thus, future research could extend and validate
our method in BMC settings involving non-exchangeable HMs.

## Acknowledgments

Lasse Elsemüller was previously affiliated with the Department of Psychology,
University of Mannheim, and Paul-Christian Bürkner with the Cluster of
Excellence SimTech, University of Stuttgart. The authors thank Lukas
Schumacher for helpful comments on this manuscript.

LE and MS were supported by a grant from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation; GRK 2277) to the research training group
Statistical Modeling in Psychology (SMiP). LE was additionally supported by
the Google Cloud Research Credits program with the award GCP19980904. PCB was
supported by the Deutsche Forschungsgemeinschaft under Germany’s Excellence
Strategy – EXC-2075 - 390740016 (the Stuttgart Cluster of Excellence SimTech).
STR was supported by the Deutsche Forschungsgemeinschaft under Germany’s
Excellence Strategy – EXC-2181 - 390900948 (the Heidelberg Cluster of
Excellence STRUCTURES).

## References

  * [1] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S Corrado, Andy Davis, Jeffrey Dean and Matthieu Devin  “TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems”  In _arXiv preprint arXiv:1603.04467_ , 2015
  * [2] Andrew Barron, Mark J Schervish and Larry Wasserman  “The consistency of posterior distributions in nonparametric problems”  In _The Annals of Statistics_ 27.2 Institute of Mathematical Statistics, 1999, pp. 536–561
  * [3] Mark A Beaumont  “Approximate Bayesian computation in evolution and ecology”  In _Annual review of ecology, evolution, and systematics_ 41 Annual Reviews, 2010, pp. 379–406
  * [4] Yoshua Bengio, Jérôme Louradour, Ronan Collobert and Jason Weston  “Curriculum learning”  In _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 41–48
  * [5] Charles H Bennett  “Efficient estimation of free energy differences from Monte Carlo data”  In _Journal of Computational Physics_ 22.2 Elsevier, 1976, pp. 245–268
  * [6] José M Bernardo and Adrian FM Smith  “Bayesian Theory”  John Wiley & Sons, 1994
  * [7] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B Shah  “Julia: A fresh approach to numerical computing”  In _SIAM review_ 59.1 SIAM, 2017, pp. 65–98
  * [8] Benjamin Bloem-Reddy and Yee Whye Teh  “Probabilistic Symmetries and Invariant Neural Networks.”  In _J. Mach. Learn. Res._ 21, 2020, pp. 90–1
  * [9] Udo Boehm, Jeffrey Annis, Michael J. Frank, Guy E. Hawkins, Andrew Heathcote, David Kellen, Angelos-Miltiadis Krypotos, Veronika Lerche, Gordon D. Logan, Thomas J. Palmeri, Don van Ravenzwaaij, Mathieu Servant, Henrik Singmann, Jeffrey J. Starns, Andreas Voss, Thomas V. Wiecki, Dora Matzke and Eric-Jan Wagenmakers  “Estimating across-trial variability parameters of the Diffusion Decision Model: Expert advice and recommendations”  In _Journal of Mathematical Psychology_ 87, 2018, pp. 46–75
  * [10] Paul-Christian Bürkner  “brms: An R package for Bayesian multilevel models using Stan”  In _Journal of statistical software_ 80, 2017, pp. 1–28
  * [11] Paul-Christian Bürkner, Maximilian Scholz and Stefan T Radev  “Some models are useful, but how do we know which ones? Towards a unified Bayesian model taxonomy”  In _arXiv preprint arXiv:2209.02439_ , 2022
  * [12] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li and Allen Riddell  “Stan: A probabilistic programming language”  In _Journal of statistical software_ 76.1 Columbia Univ., New York, NY (United States); Harvard Univ., Cambridge, MA …, 2017
  * [13] Grégoire Clarté, Christian P Robert, Robin J Ryder and Julien Stoehr  “Componentwise approximate Bayesian computation via Gibbs-like steps”  In _Biometrika_ 108.3 Oxford University Press, 2021, pp. 591–607
  * [14] Lucas Clarté, Bruno Loureiro, Florent Krzakala and Lenka Zdeborová  “A study of uncertainty quantification in overparametrized high-dimensional models”  In _arXiv preprint arXiv:2210.12760_ , 2022
  * [15] Anne GE Collins and Amitai Shenhav  “Advances in modeling learning and decision-making in neuroscience”  In _Neuropsychopharmacology_ 47.1 Nature Publishing Group, 2022, pp. 104–118
  * [16] Peter Congdon  “Bayesian model choice based on Monte Carlo estimates of posterior model probabilities”  In _Computational statistics & data analysis_ 50.2 Elsevier, 2006, pp. 346–357
  * [17] Kyle Cranmer, Johann Brehmer and Gilles Louppe  “The frontier of simulation-based inference”  In _Proceedings of the National Academy of Sciences_ 117.48 National Acad Sciences, 2020, pp. 30055–30062
  * [18] Chris Cremer, Xuechen Li and David Duvenaud  “Inference suboptimality in variational autoencoders”  In _International Conference on Machine Learning_ , 2018, pp. 1078–1086  PMLR
  * [19] Katalin Csilléry, Michael GB Blum, Oscar E Gaggiotti and Olivier François  “Approximate Bayesian computation (ABC) in practice”  In _Trends in Ecology & Evolution_ 25.7, 2010, pp. 410–418
  * [20] Morris H DeGroot and Stephen E Fienberg  “The comparison and evaluation of forecasters”  In _Journal of the Royal Statistical Society: Series D (The Statistician)_ 32.1-2 Wiley Online Library, 1983, pp. 12–22
  * [21] James M Dickey and BP Lientz  “The weighted likelihood ratio, sharp hypotheses about chances, the order of a Markov chain”  In _The Annals of Mathematical Statistics_ JSTOR, 1970, pp. 214–226
  * [22] Charles C Driver and Manuel C Voelkle  “Hierarchical Bayesian continuous time dynamic modeling.”  In _Psychological Methods_ 23.4 American Psychological Association, 2018, pp. 774
  * [23] Edgar Erdfelder, Tina-Sarah Auer, Benjamin E Hilbig, André Aßfalg, Morten Moshagen and Lena Nadarevic  “Multinomial processing tree models: A review of the literature.”  In _Zeitschrift für Psychologie/Journal of Psychology_ 217.3 Hogrefe & Huber Publishers, 2009, pp. 108
  * [24] Simon Farrell and Stephan Lewandowsky  “Computational modeling of cognition and behavior”  Cambridge University Press, 2018
  * [25] Alexander Fengler, Lakshmi N Govindarajan, Tony Chen and Michael J Frank  “Likelihood approximation networks (LANs) for fast inference of simulation models in cognitive neuroscience”  In _Elife_ 10 eLife Sciences Publications Limited, 2021, pp. e65074
  * [26] Andrew Gelman  “Multilevel (hierarchical) modeling: what it can and cannot do”  In _Technometrics_ 48.3 Taylor & Francis, 2006, pp. 432–435
  * [27] Andrew Gelman and Xiao-Li Meng  “Simulating normalizing constants: From importance sampling to bridge sampling to path sampling”  In _Statistical science_ JSTOR, 1998, pp. 163–185
  * [28] Andrew Gelman, Aki Vehtari, Daniel Simpson, Charles C Margossian, Bob Carpenter, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner and Martin Modrák  “Bayesian workflow”  In _arXiv preprint arXiv:2011.01808_ , 2020
  * [29] Amin Ghaderi-Kangavari, Jamal Amani Rad and Michael D Nunez  “A general integrative neurocognitive modeling framework to jointly describe EEG and decision-making on single trials”  In _Computational Brain & Behavior_ Springer, 2023, pp. 1–60
  * [30] Tilmann Gneiting and Adrian E Raftery  “Strictly proper scoring rules, prediction, and estimation”  In _Journal of the American statistical Association_ 102.477 Taylor & Francis, 2007, pp. 359–378
  * [31] Pedro J Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William F Podlaski, Sara A Haddad and Tim P Vogels  “Training deep neural density estimators to identify mechanistic models of neural dynamics”  In _Elife_ 9 eLife Sciences Publications Limited, 2020, pp. e56261
  * [32] Will Grathwohl, Kuan-Chieh Wang, Jörn-Henrik Jacobsen, David Duvenaud, Mohammad Norouzi and Kevin Swersky  “Your classifier is secretly an energy based model and you should treat it like one”  In _arXiv preprint arXiv:1912.03263_ , 2019
  * [33] David Marvin Green and John A Swets  “Signal detection theory and psychophysics”  Wiley New York, 1966
  * [34] Quentin F. Gronau  “Hierarchical Normal Example (Stan)”  In _CRAN_ , 2021  URL: https://cran.csiro.au/web/packages/bridgesampling/vignettes/bridgesampling_example_stan.html
  * [35] Quentin F Gronau, Alexandra Sarafoglou, Dora Matzke, Alexander Ly, Udo Boehm, Maarten Marsman, David S Leslie, Jonathan J Forster, Eric-Jan Wagenmakers and Helen Steingroever  “A tutorial on bridge sampling”  In _Journal of mathematical psychology_ 81 Elsevier, 2017, pp. 80–97
  * [36] Quentin F Gronau, Henrik Singmann and Eric-Jan Wagenmakers  “bridgesampling: An R package for estimating normalizing constants”  In _arXiv preprint arXiv:1710.08162_ , 2017
  * [37] Quentin F Gronau and Eric-Jan Wagenmakers  “Limitations of Bayesian leave-one-out cross-validation for model selection”  In _Computational brain & behavior_ 2.1 Springer, 2019, pp. 1–11
  * [38] Quentin F Gronau and Eric-Jan Wagenmakers  “Rejoinder: More limitations of Bayesian leave-one-out cross-validation”  In _Computational Brain & Behavior_ 2.1 Springer, 2019, pp. 35–47
  * [39] Quentin F Gronau, Eric-Jan Wagenmakers, Daniel W Heck and Dora Matzke  “A simple method for comparing complex models: Bayesian model comparison for hierarchical multinomial processing tree models using Warp-III bridge sampling”  In _Psychometrika_ 84.1 Springer, 2019, pp. 261–284
  * [40] Quentin F Gronau, Andrew Heathcote and Dora Matzke  “Computing Bayes factors for evidence-accumulation models using Warp-III bridge sampling”  In _Behavior research methods_ 52.2 Springer, 2020, pp. 918–937
  * [41] Chuan Guo, Geoff Pleiss, Yu Sun and Kilian Q Weinberger  “On calibration of modern neural networks”  In _International Conference on Machine Learning_ , 2017, pp. 1321–1330  PMLR
  * [42] Julia M Haaf and Jeffrey N Rouder  “Developing constraint in Bayesian mixed models.”  In _Psychological methods_ 22.4 American Psychological Association, 2017, pp. 779
  * [43] Julia M Haaf, Fayette Klaassen and Jeffrey N Rouder  “Bayes factor vs. Posterior-Predictive Model Assessment: Insights from Ordinal Constraints”  In _PsyArXiv preprint_ PsyArXiv, 2021
  * [44] Andrew Heathcote and Dora Matzke  “Winner takes all! What are race models, and why and how should psychologists use them?”  In _Current Directions in Psychological Science_ 31.5 SAGE Publications Sage CA: Los Angeles, CA, 2022, pp. 383–394
  * [45] Daniel W Heck and Edgar Erdfelder  “Maximizing the expected information gain of cognitive modeling via design optimization”  In _Computational Brain & Behavior_ 2.3 Springer, 2019, pp. 202–209
  * [46] Daniel W Heck, Udo Boehm, Florian Böing-Messing, Paul-Christian Bürkner, Koen Derks, Zoltan Dienes, Qianrao Fu, Xin Gu, Diana Karimova and Henk AL Kiers  “A review of applications of the bayes factor in psychological research.”  In _Psychological Methods_ American Psychological Association, 2022
  * [47] Joeri Hermans, Nilanjan Banik, Christoph Weniger, Gianfranco Bertone and Gilles Louppe  “Towards constraining warm dark matter with stellar streams through neural simulation-based inference”  In _Monthly Notices of the Royal Astronomical Society_ 507.2 Oxford University Press, 2021, pp. 1999–2011
  * [48] Joeri Hermans, Arnaud Delaunoy, François Rozet, Antoine Wehenkel and Gilles Louppe  “Averting a crisis in simulation-based inference”  In _stat_ 1050, 2021, pp. 14
  * [49] Samuel R Hinton, TM Davis, Alex G Kim, D Brout, Christopher Brian D’Andrea, Richard Kessler, J Lasker, Chris Lidman, E Macaulay and A Möller  “Steve: A hierarchical Bayesian model for supernova cosmology”  In _The Astrophysical Journal_ 876.1 IOP Publishing, 2019, pp. 15
  * [50] Joop J Hox, Mirjam Moerbeek and Rens Van de Schoot  “Multilevel analysis: Techniques and applications”  Routledge, 2017
  * [51] Abdollah Jalilian and Jorge Mateu  “A hierarchical spatio-temporal model to analyze relative risk variations of COVID-19: a focus on Spain, Italy and Germany”  In _Stochastic Environmental Research and Risk Assessment_ 35.4 Springer, 2021, pp. 797–812
  * [52] Robert E Kass and Adrian E Raftery  “Bayes factors”  In _Journal of the american statistical association_ 90.430 Taylor & Francis, 1995, pp. 773–795
  * [53] Diederik P. Kingma and Jimmy Lei Ba  “Adam: A method for stochastic optimization”  In _3rd International Conference on Learning Representations_ , 2015, pp. 1–15
  * [54] Karl Christoph Klauer  “Hierarchical multinomial processing tree models: A latent-trait approach”  In _Psychometrika_ 75.1 Springer, 2010, pp. 70–98
  * [55] Michael D Lee  “How cognitive modeling can benefit from hierarchical Bayesian models”  In _Journal of Mathematical Psychology_ 55.1 Elsevier, 2011, pp. 1–7
  * [56] Veronika Lerche and Andreas Voss  “Model complexity in diffusion modeling: Benefits of making the model more parsimonious”  In _Frontiers in Psychology_ 7.1324, 2016, pp. 1–14
  * [57] Tom Lodewyckx, Francis Tuerlinckx, Peter Kuppens, Nicholas B Allen and Lisa Sheeber  “A hierarchical state space approach to affective dynamics”  In _Journal of mathematical psychology_ 55.1 Elsevier, 2011, pp. 68–83
  * [58] Sanae Lotfi, Pavel Izmailov, Gregory Benton, Micah Goldblum and Andrew Gordon Wilson  “Bayesian Model Selection, the Marginal Likelihood, and Generalization”  In _arXiv preprint arXiv:2202.11678_ , 2022
  * [59] David MacKay  “Information theory, inference and learning algorithms”  Cambridge university press, 2003
  * [60] Jean-Michel Marin, Pierre Pudlo, Arnaud Estoup and Christian Robert  “Likelihood-free model choice”  ChapmanHall/CRC Press Boca Raton, FL, 2018
  * [61] Paul Marjoram, John Molitor, Vincent Plagnol and Simon Tavaré  “Markov chain Monte Carlo without likelihoods”  In _Proceedings of the National Academy of Sciences_ 100.26 National Acad Sciences, 2003, pp. 15324–15328
  * [62] Richard McElreath  “Statistical rethinking: A Bayesian course with examples in R and Stan”  ChapmanHall/CRC, 2020
  * [63] Xiao-Li Meng and Stephen Schilling  “Warp bridge sampling”  In _Journal of Computational and Graphical Statistics_ 11.3 Taylor & Francis, 2002, pp. 552–586
  * [64] Xiao-Li Meng and Wing Hung Wong  “Simulating ratios of normalizing constants via a simple identity: a theoretical exploration”  In _Statistica Sinica_ JSTOR, 1996, pp. 831–860
  * [65] Ulf Kai Mertens, Andreas Voss and Stefan Radev  “ABrox—A user-friendly Python module for approximate Bayesian computation with a focus on model comparison”  In _PloS one_ 13.3 Public Library of Science San Francisco, CA USA, 2018, pp. e0193981
  * [66] Merijn Mestdagh, Stijn Verdonck, Kristof Meers, Tim Loossens and Francis Tuerlinckx  “Prepaid parameter estimation without likelihoods”  In _PLoS computational biology_ 15.9 Public Library of Science San Francisco, CA USA, 2019, pp. e1007181
  * [67] Jay I Myung and Mark A Pitt  “Optimal experimental design for model discrimination.”  In _Psychological review_ 116.3 American Psychological Association, 2009, pp. 499
  * [68] Mahdi Pakdaman Naeini, Gregory Cooper and Milos Hauskrecht  “Obtaining Well Calibrated Probabilities Using Bayesian Binning”  In _Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence_ AAAI Press, 2015, pp. 2901–2907
  * [69] Hong-Wei Ng, Viet Dung Nguyen, Vassilios Vonikakis and Stefan Winkler  “Deep learning for emotion recognition on small datasets using transfer learning”  In _Proceedings of the 2015 ACM on international conference on multimodal interaction_ , 2015, pp. 443–449
  * [70] Bruno Nicenboim, Daniel J Schad and Shravan Vasishth  “An introduction to Bayesian data analysis for cognitive science”, 2022  URL: https://vasishth.github.io/bayescogsci/book/
  * [71] Alexandru Niculescu-Mizil and Rich Caruana  “Predicting good probabilities with supervised learning”  In _Proceedings of the 22nd international conference on Machine learning_ , 2005, pp. 625–632
  * [72] Anthony O’Hagan  “Fractional Bayes factors for model comparison”  In _Journal of the Royal Statistical Society: Series B (Methodological)_ 57.1 Wiley Online Library, 1995, pp. 99–118
  * [73] Pierre Pudlo, Jean-Michel Marin, Arnaud Estoup, Jean-Marie Cornuet, Mathieu Gautier and Christian P Robert  “Reliable ABC model choice via random forests”  In _Bioinformatics_ 32.6 Oxford University Press, 2016, pp. 859–866
  * [74] Stefan T Radev, Marco D’Alessandro, Ulf K Mertens, Andreas Voss, Ullrich Köthe and Paul-Christian Bürkner  “Amortized bayesian model comparison with evidential deep learning”  In _IEEE Transactions on Neural Networks and Learning Systems_ IEEE, 2021
  * [75] Stefan T Radev, Frederik Graw, Simiao Chen, Nico T Mutters, Vanessa M Eichel, Till Bärnighausen and Ullrich Köthe  “OutbreakFlow: Model-based Bayesian inference of disease outbreak dynamics with invertible neural networks and its application to the COVID-19 pandemics in Germany”  In _PLoS computational biology_ 17.10 Public Library of Science San Francisco, CA USA, 2021, pp. e1009472
  * [76] Stefan T Radev, Ulf K Mertens, Andreas Voss, Lynton Ardizzone and Ullrich Köthe  “BayesFlow: Learning complex stochastic models with invertible neural networks”  In _IEEE transactions on neural networks and learning systems_ IEEE, 2020
  * [77] Stefan T Radev, Marvin Schmitt, Lukas Schumacher, Lasse Elsemüller, Valentin Pratz, Yannik Schälte, Ullrich Köthe and Paul-Christian Bürkner  “BayesFlow: Amortized Bayesian Workflows With Neural Networks”  In _arXiv preprint arXiv:2306.16015_ , 2023
  * [78] Roger Ratcliff, Philip L Smith, Scott D Brown and Gail McKoon  “Diffusion decision model: Current issues and history”  In _Trends in cognitive sciences_ 20.4 Elsevier, 2016, pp. 260–281
  * [79] David M Riefer and William H Batchelder  “Multinomial modeling and the measurement of cognitive processes.”  In _Psychological Review_ 95.3 American Psychological Association, 1988, pp. 318
  * [80] Jeffrey N Rouder and Jun Lu  “An introduction to Bayesian hierarchical models with an application in the theory of signal detection”  In _Psychonomic bulletin & review_ 12.4 Springer, 2005, pp. 573–604
  * [81] Jeffrey N Rouder and Richard D Morey  “Default Bayes factors for model selection in regression”  In _Multivariate Behavioral Research_ 47.6 Taylor & Francis, 2012, pp. 877–903
  * [82] Jeffrey N Rouder, Richard D Morey and Michael S Pratte  “Bayesian hierarchical models of cognition.”  In _New handbook of mathematical psychology: Foundations and methodology_ Cambridge University Press, 2017, pp. 504–551
  * [83] Alexandra Sarafoglou, Beatrice G Kuhlmann, Frederik Aust and Julia M Haaf  “Theory-Informed Refinement of Bayesian Hierarchical MPT Modeling”  In _PsyArXiv preprint_ PsyArXiv, 2022
  * [84] Daniel J Schad, Michael Betancourt and Shravan Vasishth  “Toward a principled Bayesian workflow in cognitive science.”  In _Psychological methods_ 26.1 American Psychological Association, 2021, pp. 103
  * [85] Daniel J Schad, Bruno Nicenboim, Paul-Christian Bürkner, Michael Betancourt and Shravan Vasishth  “Workflow techniques for the robust use of bayes factors.”  In _Psychological Methods_ American Psychological Association, 2022
  * [86] Marvin Schmitt, Stefan T Radev and Paul-Christian Bürkner  “Meta-Uncertainty in Bayesian Model Comparison”  In _arXiv preprint arXiv:2210.07278_ , 2022
  * [87] Lukas Schumacher, Paul-Christian Bürkner, Andreas Voss, Ullrich Köthe and Stefan T Radev  “Neural Superstatistics: A Bayesian Method for Estimating Dynamic Models of Cognition”  In _arXiv preprint arXiv:2211.13165_ , 2022
  * [88] H. Singmann and D. Kellen  “An introduction to mixed models for experimental psychology”  In _New methods in Cognitive Psychology_ New York: Routledge, 2019, pp. 4–31
  * [89] Henrik Singmann and David Kellen  “MPTinR: Analysis of multinomial processing tree models in R”  In _Behavior Research Methods_ 45.2 Springer, 2013, pp. 560–575
  * [90] Scott A Sisson, Yanan Fan and Mark M Tanaka  “Sequential monte carlo without likelihoods”  In _Proceedings of the National Academy of Sciences_ 104.6 National Acad Sciences, 2007, pp. 1760–1765
  * [91] Stan Development Team  “Stan Modeling Language Users Guide and Reference Manual” Version 2.21.0, 2019  URL: https://mc-stan.org
  * [92] Mikael Sunnåker, Alberto Giovanni Busetto, Elina Numminen, Jukka Corander, Matthieu Foll and Christophe Dessimoz  “Approximate bayesian computation”  In _PLoS computational biology_ 9.1 Public Library of Science, 2013, pp. e1002803
  * [93] Sean Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari and Andrew Gelman  “Validating Bayesian inference algorithms with simulation-based calibration”  In _arXiv preprint arXiv:1804.06788_ , 2018
  * [94] Tijmen Tieleman and Geoffrey Hinton  “Lecture 6.5-rmsprop, coursera: Neural networks for machine learning”  In _University of Toronto, Technical Report_ 6, 2012
  * [95] Lisa Torrey and Jude Shavlik  “Transfer learning”  In _Handbook of research on machine learning applications and trends: algorithms, methods, and techniques_ IGI global, 2010, pp. 242–264
  * [96] N.-Han Tran, Leendert Maanen, Andrew Heathcote and Dora Matzke  “Systematic Parameter Reviews in Cognitive Modeling: Towards a Robust and Cumulative Characterization of Psychological Processes in the Diffusion Decision Model”  In _Frontiers in Psychology_ 11, 2021
  * [97] Brandon M Turner and Per B Sederberg  “A generalized, likelihood-free method for posterior estimation”  In _Psychonomic bulletin & review_ 21.2 Springer, 2014, pp. 227–250
  * [98] Brandon M Turner and Trisha Van Zandt  “Hierarchical approximate Bayesian computation”  In _Psychometrika_ 79.2 Springer, 2014, pp. 185–209
  * [99] Esther Ulitzsch, Matthias Davier and Steffi Pohl  “A multiprocess item response model for not-reached items due to time limits and quitting”  In _Educational and Psychological Measurement_ 80.3 Sage Publications Sage CA: Los Angeles, CA, 2020, pp. 522–547
  * [100] Iris Van Rooij, Mark Blokpoel, Johan Kwisthout and Todd Wareham  “Cognition and intractability: A guide to classical and parameterized complexity analysis”  Cambridge University Press, 2019
  * [101] Joachim Vandekerckhove, Francis Tuerlinckx and Michael D Lee  “Hierarchical diffusion models for two-choice response times.”  In _Psychological methods_ 16.1 American Psychological Association, 2011, pp. 44
  * [102] Wolf Vanpaemel  “Prior sensitivity in theory testing: An apologia for the Bayes factor”  In _Journal of Mathematical Psychology_ 54.6 Elsevier, 2010, pp. 491–498
  * [103] Aki Vehtari, Daniel P Simpson, Yuling Yao and Andrew Gelman  “Limitations of “Limitations of Bayesian leave-one-out cross-validation for model selection””  In _Computational Brain & Behavior_ 2.1 Springer, 2019, pp. 22–27
  * [104] Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter and Paul-Christian Bürkner  “Rank-normalization, folding, and localization: an improved R for assessing convergence of MCMC (with discussion)”  In _Bayesian analysis_ 16.2 International Society for Bayesian Analysis, 2021, pp. 667–718
  * [105] Andreas Voss and Jochen Voss  “Fast-dm: A free program for efficient diffusion model analysis”  In _Behavior research methods_ 39.4 Springer, 2007, pp. 767–775
  * [106] Andreas Voss, Veronika Lerche, Ulf Mertens and Jochen Voss  “Sequential sampling models with variable boundaries and non-normal noise: A comparison of six models”  In _Psychonomic bulletin & review_ 26.3 Springer, 2019, pp. 813–832
  * [107] Eric-Jan Wagenmakers, Tom Lodewyckx, Himanshu Kuriyal and Raoul Grasman  “Bayesian hypothesis testing for psychologists: A tutorial on the Savage–Dickey method”  In _Cognitive psychology_ 60.3 Elsevier, 2010, pp. 158–189
  * [108] Thomas V Wiecki, Imri Sofer and Michael J Frank  “HDDM: Hierarchical Bayesian estimation of the Drift-Diffusion Model in Python”  In _Frontiers in Neuroinformatics_ 7 Frontiers, 2013
  * [109] Eva Marie Wieschen, Andreas Voss and Stefan Radev  “Jumping to conclusion? a lévy flight model of decision making”  In _The Quantitative Methods for Psychology_ 16.2, 2020, pp. 120–132
  * [110] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov and Alexander J Smola  “Deep sets”  In _Advances in neural information processing systems_ 30, 2017

## References

  * [111] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S Corrado, Andy Davis, Jeffrey Dean and Matthieu Devin  “TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems”  In _arXiv preprint arXiv:1603.04467_ , 2015
  * [112] Andrew Barron, Mark J Schervish and Larry Wasserman  “The consistency of posterior distributions in nonparametric problems”  In _The Annals of Statistics_ 27.2 Institute of Mathematical Statistics, 1999, pp. 536–561
  * [113] Mark A Beaumont  “Approximate Bayesian computation in evolution and ecology”  In _Annual review of ecology, evolution, and systematics_ 41 Annual Reviews, 2010, pp. 379–406
  * [114] Yoshua Bengio, Jérôme Louradour, Ronan Collobert and Jason Weston  “Curriculum learning”  In _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 41–48
  * [115] Charles H Bennett  “Efficient estimation of free energy differences from Monte Carlo data”  In _Journal of Computational Physics_ 22.2 Elsevier, 1976, pp. 245–268
  * [116] José M Bernardo and Adrian FM Smith  “Bayesian Theory”  John Wiley & Sons, 1994
  * [117] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B Shah  “Julia: A fresh approach to numerical computing”  In _SIAM review_ 59.1 SIAM, 2017, pp. 65–98
  * [118] Benjamin Bloem-Reddy and Yee Whye Teh  “Probabilistic Symmetries and Invariant Neural Networks.”  In _J. Mach. Learn. Res._ 21, 2020, pp. 90–1
  * [119] Udo Boehm, Jeffrey Annis, Michael J. Frank, Guy E. Hawkins, Andrew Heathcote, David Kellen, Angelos-Miltiadis Krypotos, Veronika Lerche, Gordon D. Logan, Thomas J. Palmeri, Don van Ravenzwaaij, Mathieu Servant, Henrik Singmann, Jeffrey J. Starns, Andreas Voss, Thomas V. Wiecki, Dora Matzke and Eric-Jan Wagenmakers  “Estimating across-trial variability parameters of the Diffusion Decision Model: Expert advice and recommendations”  In _Journal of Mathematical Psychology_ 87, 2018, pp. 46–75
  * [120] Paul-Christian Bürkner  “brms: An R package for Bayesian multilevel models using Stan”  In _Journal of statistical software_ 80, 2017, pp. 1–28
  * [121] Paul-Christian Bürkner, Maximilian Scholz and Stefan T Radev  “Some models are useful, but how do we know which ones? Towards a unified Bayesian model taxonomy”  In _arXiv preprint arXiv:2209.02439_ , 2022
  * [122] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li and Allen Riddell  “Stan: A probabilistic programming language”  In _Journal of statistical software_ 76.1 Columbia Univ., New York, NY (United States); Harvard Univ., Cambridge, MA …, 2017
  * [123] Grégoire Clarté, Christian P Robert, Robin J Ryder and Julien Stoehr  “Componentwise approximate Bayesian computation via Gibbs-like steps”  In _Biometrika_ 108.3 Oxford University Press, 2021, pp. 591–607
  * [124] Lucas Clarté, Bruno Loureiro, Florent Krzakala and Lenka Zdeborová  “A study of uncertainty quantification in overparametrized high-dimensional models”  In _arXiv preprint arXiv:2210.12760_ , 2022
  * [125] Anne GE Collins and Amitai Shenhav  “Advances in modeling learning and decision-making in neuroscience”  In _Neuropsychopharmacology_ 47.1 Nature Publishing Group, 2022, pp. 104–118
  * [126] Peter Congdon  “Bayesian model choice based on Monte Carlo estimates of posterior model probabilities”  In _Computational statistics & data analysis_ 50.2 Elsevier, 2006, pp. 346–357
  * [127] Kyle Cranmer, Johann Brehmer and Gilles Louppe  “The frontier of simulation-based inference”  In _Proceedings of the National Academy of Sciences_ 117.48 National Acad Sciences, 2020, pp. 30055–30062
  * [128] Chris Cremer, Xuechen Li and David Duvenaud  “Inference suboptimality in variational autoencoders”  In _International Conference on Machine Learning_ , 2018, pp. 1078–1086  PMLR
  * [129] Katalin Csilléry, Michael GB Blum, Oscar E Gaggiotti and Olivier François  “Approximate Bayesian computation (ABC) in practice”  In _Trends in Ecology & Evolution_ 25.7, 2010, pp. 410–418
  * [130] Morris H DeGroot and Stephen E Fienberg  “The comparison and evaluation of forecasters”  In _Journal of the Royal Statistical Society: Series D (The Statistician)_ 32.1-2 Wiley Online Library, 1983, pp. 12–22
  * [131] James M Dickey and BP Lientz  “The weighted likelihood ratio, sharp hypotheses about chances, the order of a Markov chain”  In _The Annals of Mathematical Statistics_ JSTOR, 1970, pp. 214–226
  * [132] Charles C Driver and Manuel C Voelkle  “Hierarchical Bayesian continuous time dynamic modeling.”  In _Psychological Methods_ 23.4 American Psychological Association, 2018, pp. 774
  * [133] Edgar Erdfelder, Tina-Sarah Auer, Benjamin E Hilbig, André Aßfalg, Morten Moshagen and Lena Nadarevic  “Multinomial processing tree models: A review of the literature.”  In _Zeitschrift für Psychologie/Journal of Psychology_ 217.3 Hogrefe & Huber Publishers, 2009, pp. 108
  * [134] Simon Farrell and Stephan Lewandowsky  “Computational modeling of cognition and behavior”  Cambridge University Press, 2018
  * [135] Alexander Fengler, Lakshmi N Govindarajan, Tony Chen and Michael J Frank  “Likelihood approximation networks (LANs) for fast inference of simulation models in cognitive neuroscience”  In _Elife_ 10 eLife Sciences Publications Limited, 2021, pp. e65074
  * [136] Andrew Gelman  “Multilevel (hierarchical) modeling: what it can and cannot do”  In _Technometrics_ 48.3 Taylor & Francis, 2006, pp. 432–435
  * [137] Andrew Gelman and Xiao-Li Meng  “Simulating normalizing constants: From importance sampling to bridge sampling to path sampling”  In _Statistical science_ JSTOR, 1998, pp. 163–185
  * [138] Andrew Gelman, Aki Vehtari, Daniel Simpson, Charles C Margossian, Bob Carpenter, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner and Martin Modrák  “Bayesian workflow”  In _arXiv preprint arXiv:2011.01808_ , 2020
  * [139] Amin Ghaderi-Kangavari, Jamal Amani Rad and Michael D Nunez  “A general integrative neurocognitive modeling framework to jointly describe EEG and decision-making on single trials”  In _Computational Brain & Behavior_ Springer, 2023, pp. 1–60
  * [140] Tilmann Gneiting and Adrian E Raftery  “Strictly proper scoring rules, prediction, and estimation”  In _Journal of the American statistical Association_ 102.477 Taylor & Francis, 2007, pp. 359–378
  * [141] Pedro J Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William F Podlaski, Sara A Haddad and Tim P Vogels  “Training deep neural density estimators to identify mechanistic models of neural dynamics”  In _Elife_ 9 eLife Sciences Publications Limited, 2020, pp. e56261
  * [142] Will Grathwohl, Kuan-Chieh Wang, Jörn-Henrik Jacobsen, David Duvenaud, Mohammad Norouzi and Kevin Swersky  “Your classifier is secretly an energy based model and you should treat it like one”  In _arXiv preprint arXiv:1912.03263_ , 2019
  * [143] David Marvin Green and John A Swets  “Signal detection theory and psychophysics”  Wiley New York, 1966
  * [144] Quentin F. Gronau  “Hierarchical Normal Example (Stan)”  In _CRAN_ , 2021  URL: https://cran.csiro.au/web/packages/bridgesampling/vignettes/bridgesampling_example_stan.html
  * [145] Quentin F Gronau, Andrew Heathcote and Dora Matzke  “Computing Bayes factors for evidence-accumulation models using Warp-III bridge sampling”  In _Behavior research methods_ 52.2 Springer, 2020, pp. 918–937
  * [146] Quentin F Gronau, Alexandra Sarafoglou, Dora Matzke, Alexander Ly, Udo Boehm, Maarten Marsman, David S Leslie, Jonathan J Forster, Eric-Jan Wagenmakers and Helen Steingroever  “A tutorial on bridge sampling”  In _Journal of mathematical psychology_ 81 Elsevier, 2017, pp. 80–97
  * [147] Quentin F Gronau, Henrik Singmann and Eric-Jan Wagenmakers  “bridgesampling: An R package for estimating normalizing constants”  In _arXiv preprint arXiv:1710.08162_ , 2017
  * [148] Quentin F Gronau and Eric-Jan Wagenmakers  “Limitations of Bayesian leave-one-out cross-validation for model selection”  In _Computational brain & behavior_ 2.1 Springer, 2019, pp. 1–11
  * [149] Quentin F Gronau and Eric-Jan Wagenmakers  “Rejoinder: More limitations of Bayesian leave-one-out cross-validation”  In _Computational Brain & Behavior_ 2.1 Springer, 2019, pp. 35–47
  * [150] Quentin F Gronau, Eric-Jan Wagenmakers, Daniel W Heck and Dora Matzke  “A simple method for comparing complex models: Bayesian model comparison for hierarchical multinomial processing tree models using Warp-III bridge sampling”  In _Psychometrika_ 84.1 Springer, 2019, pp. 261–284
  * [151] Chuan Guo, Geoff Pleiss, Yu Sun and Kilian Q Weinberger  “On calibration of modern neural networks”  In _International Conference on Machine Learning_ , 2017, pp. 1321–1330  PMLR
  * [152] Julia M Haaf, Fayette Klaassen and Jeffrey N Rouder  “Bayes factor vs. Posterior-Predictive Model Assessment: Insights from Ordinal Constraints”  In _PsyArXiv preprint_ PsyArXiv, 2021
  * [153] Julia M Haaf and Jeffrey N Rouder  “Developing constraint in Bayesian mixed models.”  In _Psychological methods_ 22.4 American Psychological Association, 2017, pp. 779
  * [154] Andrew Heathcote and Dora Matzke  “Winner takes all! What are race models, and why and how should psychologists use them?”  In _Current Directions in Psychological Science_ 31.5 SAGE Publications Sage CA: Los Angeles, CA, 2022, pp. 383–394
  * [155] Daniel W Heck, Udo Boehm, Florian Böing-Messing, Paul-Christian Bürkner, Koen Derks, Zoltan Dienes, Qianrao Fu, Xin Gu, Diana Karimova and Henk AL Kiers  “A review of applications of the bayes factor in psychological research.”  In _Psychological Methods_ American Psychological Association, 2022
  * [156] Daniel W Heck and Edgar Erdfelder  “Maximizing the expected information gain of cognitive modeling via design optimization”  In _Computational Brain & Behavior_ 2.3 Springer, 2019, pp. 202–209
  * [157] Joeri Hermans, Nilanjan Banik, Christoph Weniger, Gianfranco Bertone and Gilles Louppe  “Towards constraining warm dark matter with stellar streams through neural simulation-based inference”  In _Monthly Notices of the Royal Astronomical Society_ 507.2 Oxford University Press, 2021, pp. 1999–2011
  * [158] Joeri Hermans, Arnaud Delaunoy, François Rozet, Antoine Wehenkel and Gilles Louppe  “Averting a crisis in simulation-based inference”  In _stat_ 1050, 2021, pp. 14
  * [159] Samuel R Hinton, TM Davis, Alex G Kim, D Brout, Christopher Brian D’Andrea, Richard Kessler, J Lasker, Chris Lidman, E Macaulay and A Möller  “Steve: A hierarchical Bayesian model for supernova cosmology”  In _The Astrophysical Journal_ 876.1 IOP Publishing, 2019, pp. 15
  * [160] Joop J Hox, Mirjam Moerbeek and Rens Van de Schoot  “Multilevel analysis: Techniques and applications”  Routledge, 2017
  * [161] Abdollah Jalilian and Jorge Mateu  “A hierarchical spatio-temporal model to analyze relative risk variations of COVID-19: a focus on Spain, Italy and Germany”  In _Stochastic Environmental Research and Risk Assessment_ 35.4 Springer, 2021, pp. 797–812
  * [162] Robert E Kass and Adrian E Raftery  “Bayes factors”  In _Journal of the american statistical association_ 90.430 Taylor & Francis, 1995, pp. 773–795
  * [163] Diederik P. Kingma and Jimmy Lei Ba  “Adam: A method for stochastic optimization”  In _3rd International Conference on Learning Representations_ , 2015, pp. 1–15
  * [164] Karl Christoph Klauer  “Hierarchical multinomial processing tree models: A latent-trait approach”  In _Psychometrika_ 75.1 Springer, 2010, pp. 70–98
  * [165] Michael D Lee  “How cognitive modeling can benefit from hierarchical Bayesian models”  In _Journal of Mathematical Psychology_ 55.1 Elsevier, 2011, pp. 1–7
  * [166] Veronika Lerche and Andreas Voss  “Model complexity in diffusion modeling: Benefits of making the model more parsimonious”  In _Frontiers in Psychology_ 7.1324, 2016, pp. 1–14
  * [167] Tom Lodewyckx, Francis Tuerlinckx, Peter Kuppens, Nicholas B Allen and Lisa Sheeber  “A hierarchical state space approach to affective dynamics”  In _Journal of mathematical psychology_ 55.1 Elsevier, 2011, pp. 68–83
  * [168] Sanae Lotfi, Pavel Izmailov, Gregory Benton, Micah Goldblum and Andrew Gordon Wilson  “Bayesian Model Selection, the Marginal Likelihood, and Generalization”  In _arXiv preprint arXiv:2202.11678_ , 2022
  * [169] David MacKay  “Information theory, inference and learning algorithms”  Cambridge university press, 2003
  * [170] Jean-Michel Marin, Pierre Pudlo, Arnaud Estoup and Christian Robert  “Likelihood-free model choice”  ChapmanHall/CRC Press Boca Raton, FL, 2018
  * [171] Paul Marjoram, John Molitor, Vincent Plagnol and Simon Tavaré  “Markov chain Monte Carlo without likelihoods”  In _Proceedings of the National Academy of Sciences_ 100.26 National Acad Sciences, 2003, pp. 15324–15328
  * [172] Richard McElreath  “Statistical rethinking: A Bayesian course with examples in R and Stan”  ChapmanHall/CRC, 2020
  * [173] Xiao-Li Meng and Stephen Schilling  “Warp bridge sampling”  In _Journal of Computational and Graphical Statistics_ 11.3 Taylor & Francis, 2002, pp. 552–586
  * [174] Xiao-Li Meng and Wing Hung Wong  “Simulating ratios of normalizing constants via a simple identity: a theoretical exploration”  In _Statistica Sinica_ JSTOR, 1996, pp. 831–860
  * [175] Ulf Kai Mertens, Andreas Voss and Stefan Radev  “ABrox—A user-friendly Python module for approximate Bayesian computation with a focus on model comparison”  In _PloS one_ 13.3 Public Library of Science San Francisco, CA USA, 2018, pp. e0193981
  * [176] Merijn Mestdagh, Stijn Verdonck, Kristof Meers, Tim Loossens and Francis Tuerlinckx  “Prepaid parameter estimation without likelihoods”  In _PLoS computational biology_ 15.9 Public Library of Science San Francisco, CA USA, 2019, pp. e1007181
  * [177] Jay I Myung and Mark A Pitt  “Optimal experimental design for model discrimination.”  In _Psychological review_ 116.3 American Psychological Association, 2009, pp. 499
  * [178] Mahdi Pakdaman Naeini, Gregory Cooper and Milos Hauskrecht  “Obtaining Well Calibrated Probabilities Using Bayesian Binning”  In _Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence_ AAAI Press, 2015, pp. 2901–2907
  * [179] Hong-Wei Ng, Viet Dung Nguyen, Vassilios Vonikakis and Stefan Winkler  “Deep learning for emotion recognition on small datasets using transfer learning”  In _Proceedings of the 2015 ACM on international conference on multimodal interaction_ , 2015, pp. 443–449
  * [180] Bruno Nicenboim, Daniel J Schad and Shravan Vasishth  “An introduction to Bayesian data analysis for cognitive science”, 2022  URL: https://vasishth.github.io/bayescogsci/book/
  * [181] Alexandru Niculescu-Mizil and Rich Caruana  “Predicting good probabilities with supervised learning”  In _Proceedings of the 22nd international conference on Machine learning_ , 2005, pp. 625–632
  * [182] Anthony O’Hagan  “Fractional Bayes factors for model comparison”  In _Journal of the Royal Statistical Society: Series B (Methodological)_ 57.1 Wiley Online Library, 1995, pp. 99–118
  * [183] Pierre Pudlo, Jean-Michel Marin, Arnaud Estoup, Jean-Marie Cornuet, Mathieu Gautier and Christian P Robert  “Reliable ABC model choice via random forests”  In _Bioinformatics_ 32.6 Oxford University Press, 2016, pp. 859–866
  * [184] Stefan T Radev, Marco D’Alessandro, Ulf K Mertens, Andreas Voss, Ullrich Köthe and Paul-Christian Bürkner  “Amortized bayesian model comparison with evidential deep learning”  In _IEEE Transactions on Neural Networks and Learning Systems_ IEEE, 2021
  * [185] Stefan T Radev, Frederik Graw, Simiao Chen, Nico T Mutters, Vanessa M Eichel, Till Bärnighausen and Ullrich Köthe  “OutbreakFlow: Model-based Bayesian inference of disease outbreak dynamics with invertible neural networks and its application to the COVID-19 pandemics in Germany”  In _PLoS computational biology_ 17.10 Public Library of Science San Francisco, CA USA, 2021, pp. e1009472
  * [186] Stefan T Radev, Ulf K Mertens, Andreas Voss, Lynton Ardizzone and Ullrich Köthe  “BayesFlow: Learning complex stochastic models with invertible neural networks”  In _IEEE transactions on neural networks and learning systems_ IEEE, 2020
  * [187] Stefan T Radev, Marvin Schmitt, Lukas Schumacher, Lasse Elsemüller, Valentin Pratz, Yannik Schälte, Ullrich Köthe and Paul-Christian Bürkner  “BayesFlow: Amortized Bayesian Workflows With Neural Networks”  In _arXiv preprint arXiv:2306.16015_ , 2023
  * [188] Roger Ratcliff, Philip L Smith, Scott D Brown and Gail McKoon  “Diffusion decision model: Current issues and history”  In _Trends in cognitive sciences_ 20.4 Elsevier, 2016, pp. 260–281
  * [189] David M Riefer and William H Batchelder  “Multinomial modeling and the measurement of cognitive processes.”  In _Psychological Review_ 95.3 American Psychological Association, 1988, pp. 318
  * [190] Jeffrey N Rouder and Jun Lu  “An introduction to Bayesian hierarchical models with an application in the theory of signal detection”  In _Psychonomic bulletin & review_ 12.4 Springer, 2005, pp. 573–604
  * [191] Jeffrey N Rouder and Richard D Morey  “Default Bayes factors for model selection in regression”  In _Multivariate Behavioral Research_ 47.6 Taylor & Francis, 2012, pp. 877–903
  * [192] Jeffrey N Rouder, Richard D Morey and Michael S Pratte  “Bayesian hierarchical models of cognition.”  In _New handbook of mathematical psychology: Foundations and methodology_ Cambridge University Press, 2017, pp. 504–551
  * [193] Alexandra Sarafoglou, Beatrice G Kuhlmann, Frederik Aust and Julia M Haaf  “Theory-Informed Refinement of Bayesian Hierarchical MPT Modeling”  In _PsyArXiv preprint_ PsyArXiv, 2022
  * [194] Daniel J Schad, Michael Betancourt and Shravan Vasishth  “Toward a principled Bayesian workflow in cognitive science.”  In _Psychological methods_ 26.1 American Psychological Association, 2021, pp. 103
  * [195] Daniel J Schad, Bruno Nicenboim, Paul-Christian Bürkner, Michael Betancourt and Shravan Vasishth  “Workflow techniques for the robust use of bayes factors.”  In _Psychological Methods_ American Psychological Association, 2022
  * [196] Marvin Schmitt, Stefan T Radev and Paul-Christian Bürkner  “Meta-Uncertainty in Bayesian Model Comparison”  In _arXiv preprint arXiv:2210.07278_ , 2022
  * [197] Lukas Schumacher, Paul-Christian Bürkner, Andreas Voss, Ullrich Köthe and Stefan T Radev  “Neural Superstatistics: A Bayesian Method for Estimating Dynamic Models of Cognition”  In _arXiv preprint arXiv:2211.13165_ , 2022
  * [198] H. Singmann and D. Kellen  “An introduction to mixed models for experimental psychology”  In _New methods in Cognitive Psychology_ New York: Routledge, 2019, pp. 4–31
  * [199] Henrik Singmann and David Kellen  “MPTinR: Analysis of multinomial processing tree models in R”  In _Behavior Research Methods_ 45.2 Springer, 2013, pp. 560–575