1. Multinomial Processing Trees

Data for MPTs consist of categorical responsesFootnote 4 from several participants to a set of items. MPTs are based on the assumption that these responses follow a multinomial distribution. MPTs reparametrize the category probabilities of the multinomial distribution in terms of the model parameters that represent the probabilities of latent cognitive processes (Riefer & Batchelder, Reference Riefer and Batchelder1988).

Consider the pair-clustering MPT depicted in Fig. 1. The model was developed for the measurement of the storage and retrieval processes that determine the recall of semantically related word pairs (Batchelder & Riefer, Reference Batchelder and Riefer1980). A typical pair-clustering study involves a free recall memory experiment, where participants are presented with a list of study words in a word-by-word fashion. The study list consists of two types of items: semantically related word pairs such as knife–fork and words without a category partner (i.e., singletons), such as dog. After the study phase, participants are required to recall as many of the study words as they can. Typically, semantically related word pairs are recalled consecutively as a “pair-cluster.”

The model represents the interplay between the hypothesized latent cognitive processes in a rooted tree structure. The pair-clustering MPT features $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K = 2$$\end{document} independent category systems. Each category system corresponds to a separate multinomial distribution: one for word pairs ( $k=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = 1$$\end{document} ) and one for singletons ( $k=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = 2$$\end{document} ). The category probabilities in each system are modeled using a separate subtree with a finite number of branches.

Figure 1. Pair-clustering MPT. Available at https://tinyurl.com/yb7bma4e under CC license https://creativecommons.org/licenses/by/2.0/.

Each branch of a subtree corresponds to a specific sequence of processing stages and terminates in one of $L_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_k$$\end{document} possible response categories denoted as $C_{\mathrm{kl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{kl}$$\end{document} , where $l=1,\dots ,L_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l = 1, \ldots , L_k$$\end{document} indexes the lth of $L_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_k$$\end{document} possible responses in subtree k. In the pair-clustering MPT, the recall of word pairs is scored into $L_1=4$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_1 = 4$$\end{document} categories: (1) Both words of the pair are recalled consecutively ( $C_{11}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{11}$$\end{document} ); (2) both words are recalled but not consecutively ( $C_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{12}$$\end{document} ); (3) only one word is recalled ( $C_{13}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{13}$$\end{document} ); (4) no word is recalled ( $C_{14}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{14}$$\end{document} ). The recall of singletons is scored into $L_2=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2 = 2$$\end{document} response categories: (1) The word is recalled ( $C_{21}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{21}$$\end{document} ); (2) the word is not recalled ( $C_{22}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{22}$$\end{document} ).

The response category probabilities are expressed as a function of the MPT parameters, ${\theta }_p\in (0,1)\forall p\in \{1,2,\dots ,P\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _p \in (0, 1) \; \; \forall p \in \{1, 2, \ldots , P\}$$\end{document} , which can be collected in a vector $\theta =({\theta }_1,{\theta }_2,\dots ,{\theta }_P)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta } = (\theta _1, \theta _2, \ldots , \theta _P)$$\end{document} . The pair-clustering MPT features four parameters: $\theta =(c,r,u,a)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta } = (c, r, u, a)$$\end{document} . The cluster-storage parameter c corresponds to the probability that a word pair is stored as a cluster in memory. The cluster-retrieval parameter r corresponds to the conditional probability that a clustered word pair is retrieved from memory during the test phase. The model assumes that stored and retrieved word clusters are always recalled consecutively. The storage-retrieval parameter u corresponds to the conditional probability that a member of a word pair is stored and retrieved, given that the word pair was not clustered. The model makes the simplifying assumption that words from unclustered pairs are never recalled consecutively. The singleton storage-retrieval parameter a corresponds to the probability that a singleton is stored and retrieved. In many applications, researchers impose the constraint that $a=u$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = u$$\end{document} .

The response category probabilities are obtained as follows. First, we obtain the probability of each branch that terminates in a given response category. Let $B_{\mathrm{klm}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_{klm}$$\end{document} denote the mth of $M_{\mathrm{kl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{kl}$$\end{document} branches that terminate in response category $C_{\mathrm{kl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{kl}$$\end{document} . The probability of branch $B_{\mathrm{klm}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_{klm}$$\end{document} is obtained by traversing the tree from root to leaf and multiplying the encountered parameters:

(4) $\begin{array}{c}\text{Pr}(B_{\mathrm{klm}}\mid \theta )=\prod _{p=1}^P{\theta }_p^{v_{\mathrm{klmp}}}{(1-{\theta }_p)}^{w_{\mathrm{klmp}}},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {Pr}(B_{klm} \mid \varvec{\theta }) = \prod ^P_{p=1}\theta ^{v_{klmp}}_{p} (1-\theta _p)^{w_{klmp}}, \end{aligned}$$\end{document}

where $v_{\mathrm{klmp}}\ge 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${v_{klmp}} \ge 0$$\end{document} and $w_{\mathrm{klmp}}\ge 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${w_{klmp}}\ge 0$$\end{document} are the number of nodes on branch $B_{\mathrm{klm}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_{klm}$$\end{document} that are related to parameter ${\theta }_p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _p$$\end{document} , $p=1,\dots ,P$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=1, \ldots , P$$\end{document} , and $1-{\theta }_p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1-\theta _p$$\end{document} , respectively. Second, we sum the probabilities of the $M_{\mathrm{kl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{kl}$$\end{document} branches that terminate in $C_{\mathrm{kl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{kl}$$\end{document} :

(5) $\begin{array}{c}\text{Pr}(C_{\mathrm{kl}}\mid \theta )=\sum _{m=1}^{M_{\mathrm{kl}}}\text{Pr}(B_{\mathrm{klm}}\mid \theta ).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {Pr}(C_{kl} \mid \varvec{\theta }) = \sum _{m = 1}^{M_{kl}} \text {Pr}(B_{klm} \mid \varvec{\theta }). \end{aligned}$$\end{document}

For instance, the probability of response category $C_{14}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{14}$$\end{document} is given by $\text{Pr}(C_{14}\mid \theta )=c(1-r)+(1-c){(1-u)}^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Pr}(C_{14} \mid \varvec{\theta }) = c \; (1 - r) + (1 -c) \; (1 - u)^2$$\end{document} .

The probability of the observed response frequencies across category systems denoted by $n=(n_{11},\dots ,n_{1L_1},\dots ,n_{K1},\dots ,n_{K{L}_K})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{n} = (n_{11},\ldots , n_{1L_1}, \ldots , n_{K1}, \ldots , n_{KL_K})$$\end{document} , where $n_{\mathrm{kl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{kl}$$\end{document} is the observed response frequency for category $l=1,\dots ,L_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l = 1,\ldots , L_k$$\end{document} in category system (subtree) $k=1,\dots ,K$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = 1,\ldots , K$$\end{document} , is given by a product-multinomial distribution:

(6) $\begin{array}{c}\text{Pr}(N=n\mid \theta )=\prod _{k=1}^K\left(\frac{J_k!}{n_{k1}!\times n_{k2}!\times ...\times n_{k{L}_k}!},\prod _{l=1}^{L_k},{\left(\text{Pr},(,C_{\mathrm{kl}},\mid ,\theta ,)\right)}^{n_{\mathrm{kl}}}\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {Pr}(\varvec{N} = \varvec{n} \mid \varvec{\theta }) = \prod ^K_{k=1} \left\{ {\frac{J_k!}{n_{k1}! \times n_{k2}! \times ... \times n_{kL_k}!}} \prod ^{L_k}_{l=1}\left[ \text {Pr}(C_{kl} \mid \varvec{\theta })\right] ^{n_{kl}}\right\} , \end{aligned}$$\end{document}

where $J_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_k$$\end{document} denotes the number of items in category system k (see also Klauer, Reference Klauer2010; Matzke et al., Reference Matzke, Dolan, Batchelder and Wagenmakers2015).

1.1. Bayesian Hierarchical MPTs: The Latent-Trait Approach

Bayesian hierarchical approaches explicitly model heterogeneity in participants by introducing a group-level distribution from which the participant-level parameters are drawn (e.g., Gelman & Hill, Reference Gelman and Hill2007; Gill, Reference Gill2002; Lee, Reference Lee2011; Lee & Wagenmakers, Reference Lee and Wagenmakers2013; Rouder & Lu, Reference Rouder and Lu2005).Footnote 5 Here we focus on Klauer’s (Reference Klauer2010) latent-trait approach that relies on a multivariate normal group-level distribution to describe the between-subject variability and the correlations between the participant-level parameters.

To model participant heterogeneity, observed responses are aggregated over items, but not over participants, resulting in a vector of category frequencies for each participant i: $n_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{n}_i$$\end{document} , $i=1,2,\dots ,I$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i = 1, 2, \ldots , I$$\end{document} , where I is the total number of participants. Each participant obtains a participant-specific parameter vector ${\theta }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_i$$\end{document} of length P.

The latent-trait approach assumes that the probit-transformed participant-level parameter vectors ${\theta }_i^{{}^{\prime }}={\Phi }^{-1}\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_{i}^{'} = \Phi ^{-1}(\varvec{\theta }_{i})$$\end{document} follow a P-dimensional multivariate normal distribution with mean vector $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} and covariance matrix $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} : ${\theta }_i^{{}^{\prime }}\sim N_P(\mu ,\Sigma )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \varvec{\theta }_{i}^{'} \sim \mathcal {N}_P(\varvec{\mu }, \varvec{\Sigma })$$\end{document} . The probit transformation ${\Phi }^{-1}\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi ^{-1}(\varvec{\theta }_i)$$\end{document} is defined component-wise, where ${\Phi }^{-1}(·)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi ^{-1}(\cdot )$$\end{document} corresponds to the inverse of the cumulative distribution function of the normal distribution. Priors are assigned to $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} and $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} . We follow earlier implementations of the latent-trait approach and assign independent standard normal distributions to the P components of $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} (Heck, Arnold, & Arnold, Reference Heck, Arnold and Arnold2018a; Matzke et al., Reference Matzke, Dolan, Batchelder and Wagenmakers2015). This choice corresponds to uniform priors on the probability scale for the grand means. For the covariance matrix $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} , a convenient prior choice would be an inverse Wishart prior with degrees of freedom $\nu =P+1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu = P + 1$$\end{document} and identity scale matrix. This setting leads to uniform priors on the correlation parameters; however, this choice is constraining on the standard deviation parameters. Although changing the degrees of freedom $\nu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document} affords more flexibility for modeling the standard deviations, it comes at the cost of constraining the prior on the correlation parameters (Gelman & Hill, Reference Gelman and Hill2007).

This dilemma can be circumvented by using a scaled inverse Wishart prior as introduced by Gelman and Hill (Reference Gelman and Hill2007) and proposed in the context of hierarchical MPT modeling by Klauer (Reference Klauer2010). Compared to a regular inverse Wishart prior, the scaled version has the advantage that it allows one to model the standard deviations more flexibly while retaining the desirable uniform prior on the correlation parameters. The scaled inverse Wishart prior is based on the following decomposition of the covariance matrix $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} :

(7) $\begin{array}{c}\Sigma =\text{Diag}\left(\xi \right)Q\text{Diag}\left(\xi \right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\Sigma } = \text {Diag}(\varvec{\xi }) \, \varvec{Q} \, \text {Diag}(\varvec{\xi }), \end{aligned}$$\end{document}

where $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} is a vector of P scaling parameters and $Q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Q}$$\end{document} corresponds to the $P\times P$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P \times P$$\end{document} unscaled covariance matrix. The scaled inverse Wishart prior is obtained by placing a regular inverse Wishart prior on the unscaled covariance matrix $Q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Q}$$\end{document} and a suitable prior on the vector of scaling parameters $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} .

We follow Klauer (Reference Klauer2010) and assign $Q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Q}$$\end{document} an inverse Wishart prior with degrees of freedom $\nu =P+1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu = P + 1$$\end{document} and scale matrix $I_P$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{I}_P$$\end{document} (i.e., $P\times P$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P \times P$$\end{document} identity matrix). For the P components of $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} , we follow Heck et al. (Reference Heck, Arnold and Arnold2018a) and use independent uniform priors that range from zero to ten. These choices correspond to relatively diffuse priors for the standard deviations of the random effects on the probit scale and uniform priors for the correlations between the random effects.

Note that these prior distributions have been proposed in a context of parameter estimation, where the exact choice of the prior is irrelevant as long as sufficiently informative data are available. In contrast, in the context of model comparison, the priors have an important and lasting effect: As shown in Eq. (3), the marginal likelihood is obtained by taking a weighted average of the probability of the data across all possible parameter settings where the weights correspond to the parameters’ prior density. We argue that the standard normal and uniform priors for the grand means and the correlations, respectively, provide a reasonable default setting also from the perspective of model comparison. The choice of the prior for $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} is less straightforward. We report the results corresponding to the default setting of the recently developed MPT software package TreeBUGS (Heck et al., Reference Heck, Arnold and Arnold2018a), but we probed the robustness of our conclusions with a sensitivity analysis using ${\xi }_p\sim \text{Uniform}(0,{\xi }_{\text{max}})\forall p\in \{1,2,\dots ,P\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _p \sim \text {Uniform}(0, \xi _{\text {max}}) \, \forall p \in \{1, 2, \ldots , P\}$$\end{document} , with ${\xi }_{\text{max}}=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 2$$\end{document} instead of ${\xi }_{\text{max}}=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 10$$\end{document} , a prior that was chosen based on the implied group-level distributions on the probability scale. As the conclusions were unaffected by the choice of the upper bound, the results of the sensitivity analysis are mentioned only briefly and are presented in more detail in Supplemental Materials available at https://osf.io/rycg6/.

Under these prior settings, the probit-transformed participant-level MPT parameter vectors can be written as:

(8) $\begin{array}{c}{\theta }_i^{{}^{\prime }}=\mu +\xi \odot {\omega }_i,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\theta }_{i}^{'} = \varvec{\mu } + \varvec{\xi } \odot \varvec{\omega }_i, \end{aligned}$$\end{document}

where ${\omega }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\omega }_i$$\end{document} is the P-dimensional vector with the unscaled random effects for participant i and $\odot$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\odot $$\end{document} denotes the Hadamard product (i.e., entry-wise multiplication, e.g., Liu & Trenkler, Reference Liu and Trenkler2008). The unscaled random effects are drawn from a P-dimensional zero-centered multivariate normal distribution with covariance matrix $Q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Q}$$\end{document} : ${\omega }_i\sim N_P(0,Q)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\omega }_i \sim \mathcal {N}_P(\varvec{0},\varvec{Q})$$\end{document} .

Note that the model is overparameterized: $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} and $Q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Q}$$\end{document} cannot be interpreted separately. Similarly, the unscaled random effects ${\omega }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\omega }_i$$\end{document} cannot be interpreted on their own but need to be combined with the scaling parameter vector $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} to form the random effects of interest. The scaling parameters $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} , the unscaled covariance matrix $Q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Q}$$\end{document} , and the unscaled random effects ${\omega }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\omega }_i$$\end{document} are not of interest in themselves and are simply an artifact of using a flexible scaled inverse Wishart prior on $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} : The parameters of interest are ${\theta }_i^{{}^{\prime }}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_{i}^{'}$$\end{document} , $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} , and $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} . Therefore, the scaled inverse Wishart prior can be regarded as a form of parameter expansion (e.g., Gelman & Hill, Reference Gelman and Hill2007) which has been reported to speed up convergence when fitting the model using Markov chain Monte Carlo sampling (MCMC; e.g., Gamerman & Lopes, Reference Gamerman and Lopes2006).

The reader is referred to Klauer (Reference Klauer2010) and Matzke et al. (Reference Matzke, Dolan, Batchelder and Wagenmakers2015) for a more detailed description of the latent-trait approach. Parameter estimation may proceed using MCMC sampling implemented in standard Bayesian statistical software such as JAGS (Plummer, Reference Plummer, Hornik, Leisch and Zeileis2003) or Stan (Stan Development Team, 2016).

1.2. Computing the Marginal Likelihood

The marginal likelihood for latent-trait MPTs is given by:Footnote 6

(9) $\begin{array}{cc}\text{Pr}(N=n)=& \int ...\int \prod _{i=1}^I\left(\stackrel{\text{individual-level}}{\overbrace{\text{Pr}(N_i=n_i\mid \mu ,\xi ,{\omega }_i)}},\stackrel{\text{group-level}}{\overbrace{p({\omega }_i\mid Q)}}\right)\stackrel{\text{priors}}{\overbrace{p\left(Q\right)p\left(\mu \right)p\left(\xi \right)}}\text{d}Q\text{d}\mu \text{d}\xi \text{d}{\omega }_1...\text{d}{\omega }_I\\ =& \int ...\int \prod _{i=1}^I[\underset{\text{Pr}(N_i=n_i\mid \mu ,\xi ,{\omega }_i)}{\underbrace{{\prod }_{k=1}^K\{\frac{J_k!}{n_{ik1}!\times n_{ik2}!\times ...\times n_{ik{L}_k}!}{\prod }_{l=1}^{L_k}{\left(\text{Pr},(,C_{\mathrm{kl}},\mid ,\mu ,,,\xi ,,,{\omega }_i,)\right)}^{n_{\mathrm{ikl}}}\}}}\\ & \times \underset{p({\omega }_i\mid Q)}{\underbrace{{\left(2,\pi \right)}^{-\frac{P}{2}}{\left(Q\right)}^{-\frac{1}{2}}exp\{-\frac{1}{2}{\omega }_i^{\top }{Q}^{-1}{\omega }_i\}}}]\\ & \times \underset{p\left(Q\right)}{\underbrace{\frac{1}{2^{\frac{\nu P}{2}}{\Gamma }_P\left(\frac{\nu }{2}\right)}{\left(Q\right)}^{-\frac{\nu +P+1}{2}}exp\left(-,\frac{1}{2},\text{tr},\left(Q^{-1}\right)\right)}}\\ & \times \underset{p\left(\mu \right)}{\underbrace{{\left(2,\pi \right)}^{-\frac{P}{2}}exp\{-\frac{1}{2}{\mu }^{\top }\mu \}}}\underset{p\left(\xi \right)}{\underbrace{{\left({\xi }_{\text{max}}\right)}^{-P}}}\text{d}Q\text{d}\mu \text{d}\xi \text{d}{\omega }_1...\text{d}{\omega }_I,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {Pr}(\varvec{N} = \varvec{n}) =&\int ... \int \prod _{i = 1}^{I} \left[ \overbrace{\text {Pr}(\varvec{N}_i = \varvec{n}_i \mid \varvec{\mu }, \varvec{\xi }, \varvec{\omega }_i)}^{\text {individual-level}} \overbrace{p(\varvec{\omega }_i \mid \varvec{Q})}^{\text {group-level}} \right] \overbrace{p(\varvec{Q})p(\varvec{\mu }) p(\varvec{\xi })}^{\text {priors}} \text {d}\varvec{Q} \text {d}\varvec{\mu } \text {d}\varvec{\xi } \text {d}\varvec{\omega }_1 ... \text {d}\varvec{\omega }_I \nonumber \\ =&{\int } ... {\int } \prod _{i = 1}^{I} \Bigg [\underbrace{\prod ^K_{k=1} \Bigg \{{\frac{J_k!}{n_{ik1}! \times n_{ik2}! \times ... \times n_{ikL_k}!}} \prod ^{L_k}_{l=1}\left[ \text {Pr}(C_{kl} \mid \varvec{\mu }, \varvec{\xi }, \varvec{\omega }_i)\right] ^{n_{ikl}}\Bigg \}}_{\text {Pr}(\varvec{N}_i = \varvec{n}_i \mid \varvec{\mu }, \varvec{\xi }, \varvec{\omega }_i)} \nonumber \\&\quad \times \underbrace{\left( 2\pi \right) ^{-\frac{P}{2}} \left| \varvec{Q}\right| ^{-\frac{1}{2}} \exp \bigg \{-\frac{1}{2} \varvec{\omega }_i^\top \varvec{Q}^{-1}\varvec{\omega }_i\bigg \}}_{p(\varvec{\omega }_i \mid \varvec{Q})}\Bigg ] \nonumber \\&\quad \times \underbrace{\frac{1}{2^{\frac{\nu P}{2}}\Gamma _P(\frac{\nu }{2})} \left| \varvec{Q}\right| ^{-\frac{\nu +P+1}{2}}\exp \left\{ -\frac{1}{2}{\text {tr}}\left( \varvec{Q}^{-1}\right) \right\} }_{p(\varvec{Q})} \nonumber \\&\quad \times \underbrace{\left( 2\pi \right) ^{-\frac{P}{2}} \exp \bigg \{-\frac{1}{2} \varvec{\mu }^\top \varvec{\mu }\bigg \}}_{p(\varvec{\mu })} \underbrace{\left( \xi _{\text {max}}\right) ^{-P}}_{p(\varvec{\xi })} \text {d}\varvec{Q} \text {d}\varvec{\mu } \text {d}\varvec{\xi } \text {d}\varvec{\omega }_1 ... \text {d}\varvec{\omega }_I, \end{aligned}$$\end{document}

where ${\Gamma }_P\left(a\right)={\pi }^{P(P-1)/4}{\prod }_{j=1}^P\Gamma \left(a,+,\frac{1-j}{2}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{P}(a)=\pi ^{{P(P-1)/4}}\prod _{{j=1}}^{P}\Gamma \left( a+\frac{1-j}{2}\right) $$\end{document} and $\Gamma \left(z\right)={\int }_0^{\infty }{x}^{z-1}{e}^{-x}\text{d}x$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma (z)=\int _{0}^{\infty }x^{z-1}e^{-x}\,\text {d}x$$\end{document} are the multivariate and regular gamma function, respectively. In this parametrization, we do not need to explicitly integrate out the participant-level parameter vectors ${\theta }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_i$$\end{document} since they are functions of $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} , $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} , and ${\omega }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\omega }_i$$\end{document} (see Eq. (8)).

We exploit the fact that the covariance matrix $Q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Q}$$\end{document} in Eq. (9) can be integrated out in closed form (see also, Overstall & Forster, Reference Overstall and Forster2010); a detailed derivation is provided in Supplemental Materials. The marginal likelihood is then given by:

(10) $\begin{array}{cc}\text{Pr}(N=n)=& \int ...\int \prod _{i=1}^I\left(\prod _{k=1}^K,\{,\frac{J_k!}{n_{ik1}!\times n_{ik2}!\times ...\times n_{ik{L}_k}!},\prod _{l=1}^{L_k},{\left(\text{Pr},(,C_{\mathrm{kl}},\mid ,\mu ,,,\xi ,,,{\omega }_i,)\right)}^{n_{\mathrm{ikl}}},\}\right)\\ & \times \frac{{\Gamma }_P\left(\frac{\nu +I}{2}\right)}{{\Gamma }_P\left(\frac{\nu }{2}\right)}\frac{{\pi }^{-\frac{\mathrm{IP}}{2}}}{{\left({\Omega }^{\top }\Omega +I_P\right)}^{\frac{\nu +I}{2}}}\times {\left(2,\pi \right)}^{-\frac{P}{2}}exp\{-\frac{1}{2}{\mu }^{\top }\mu \}\\ & \times {\left({\xi }_{\text{max}}\right)}^{-P}\text{d}\mu \text{d}\xi \text{d}{\omega }_1...\text{d}{\omega }_I,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {Pr}(\varvec{N} = \varvec{n})= & {} {\int } ... {\int } \prod _{i = 1}^{I} \left[ \prod ^K_{k=1} \Bigg \{{\frac{J_k!}{n_{ik1}! \times n_{ik2}! \times ... \times n_{ikL_k}!}} \prod ^{L_k}_{l=1}\left[ \text {Pr}(C_{kl} \mid \varvec{\mu }, \varvec{\xi }, \varvec{\omega }_i)\right] ^{n_{ikl}}\Bigg \}\right] \nonumber \\&\quad \times \frac{\Gamma _P(\frac{\nu + I}{2})}{\Gamma _P(\frac{\nu }{2})} \frac{\pi ^{-\frac{I P}{2}}}{\left| { \varvec{\Omega }^\top \varvec{\Omega } + {\varvec{I}_P}}\right| ^{\frac{\nu + I}{2}}} \times \left( 2\pi \right) ^{-\frac{P}{2}} \exp \bigg \{-\frac{1}{2} \varvec{\mu }^\top \varvec{\mu }\bigg \}\nonumber \\&\quad \times \left( \xi _{\text {max}}\right) ^{-P} \text {d}\varvec{\mu } \text {d}\varvec{\xi } \text {d}\varvec{\omega }_1 ... \text {d}\varvec{\omega }_I, \end{aligned}$$\end{document}

where $\Omega$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Omega }$$\end{document} is an $I\times P$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I \times P$$\end{document} matrix of the P-dimensional random effects vectors ${\omega }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\omega }_i$$\end{document} of the I participants. Even after integrating out $Q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Q}$$\end{document} , the expression for the marginal likelihood is still a high-dimensional integral (i.e., $P(I+2)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(I + 2)$$\end{document} dimensions); the challenge is to find a method which yields accurate estimates of this integral.

2. Warp-III Bridge Sampling for MPTs

We propose to use Warp-III bridge sampling (Meng & Schilling, Reference Meng and Schilling2002; Meng & Wong, Reference Meng and Wong1996; Overstall, Reference Overstall2010), an advanced version of bridge sampling, to evaluate the high-dimensional integral in Eq. (10). Bridge sampling is a general method for estimating normalizing constantsFootnote 7, a problem that is not only encountered in Bayesian inference, but also in likelihood-based approaches (Gelman & Meng, Reference Gelman and Meng1998). We first outline the basic principles of bridge sampling and then present the details of the advanced Warp-III method. The reader is referred to the recent tutorial by Gronau et al. (Reference Gronau, Sarafoglou, Matzke, Ly, Boehm and Marsman2017) for a detailed explanation of the general bridge sampling approach.

Let $\zeta =(\mu ,\xi ,{\omega }_1,\dots ,{\omega }_I)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\zeta } = (\varvec{\mu }, \varvec{\xi }, \varvec{\omega }_1, \ldots , \varvec{\omega }_I)$$\end{document} be the vector of quantities that must be integrated out to obtain the marginal likelihood, so that

(11) $\begin{array}{c}\begin{array}{cc}\text{Pr}(N=n)=& \int \text{Pr}(N=n\mid \zeta )p\left(\zeta \right)\text{d}\zeta .\end{array}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \text {Pr}(\varvec{N} = \varvec{n}) =&\int \text {Pr}(\varvec{N} = \varvec{n} \mid \varvec{\zeta }) \, p(\varvec{\zeta }) \text {d}\varvec{\zeta }. \end{aligned} \end{aligned}$$\end{document}

General bridge sampling is based on the following identity:

(12) $\begin{array}{c}1=\frac{\int \stackrel{\text{bridge function}}{\overbrace{h\left(\zeta \right)}}p(\zeta \mid N=n)\stackrel{\text{proposal distribution}}{\overbrace{g\left(\zeta \right)}}\text{d}\zeta }{\int h\left(\zeta \right)\underset{\text{posterior distribution}}{\underbrace{p(\zeta \mid N=n)}}g\left(\zeta \right)\text{d}\zeta },\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} 1 = \frac{ {\int } \overbrace{h(\varvec{\zeta })}^{\text {bridge function}} \, p(\varvec{\zeta } \mid \varvec{N} = \varvec{n}) \, \overbrace{g(\varvec{\zeta })}^{\text {proposal distribution}} \text {d}\varvec{\zeta }}{{\int } h(\varvec{\zeta }) \underbrace{p(\varvec{\zeta } \mid \varvec{N} = \varvec{n})}_{\text {posterior distribution}} g(\varvec{\zeta }) \text {d}\varvec{\zeta }}, \end{aligned}$$\end{document}

where $p(\zeta \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\varvec{\zeta } \mid \varvec{N} = \varvec{n})$$\end{document} is the posterior distribution of $\zeta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\zeta }$$\end{document} , $g\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\zeta })$$\end{document} is the probability density function of a proposal distribution, and $h\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h(\varvec{\zeta })$$\end{document} is a function such that $0<\left(\int ,h,\left(,\zeta ,\right),,p,(,\zeta ,\mid ,N,=,n,),,g,\left(,\zeta ,\right),\text{d},\zeta \right)<\infty$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0< \left| \int h(\varvec{\zeta }) \, p(\varvec{\zeta } \mid \varvec{N} = \varvec{n}) \, g(\varvec{\zeta }) \text {d}\varvec{\zeta }\right| < \infty $$\end{document} . It follows from Eq. (12) that

(13) $\begin{array}{c}\text{Pr}(N=n)=\frac{\int h\left(\zeta \right)\text{Pr}(N=n\mid \zeta )p\left(\zeta \right)g\left(\zeta \right)\text{d}\zeta }{\int h\left(\zeta \right)g\left(\zeta \right)p(\zeta \mid N=n)\text{d}\zeta }=\frac{E_{g\left(\zeta \right)}\left(h,\left(,\zeta ,\right),,\text{Pr},(,N,=,n,\mid ,\zeta ,),,p,\left(,\zeta ,\right)\right)}{E_{p(\zeta \mid N=n)}\left(h,\left(,\zeta ,\right),,g,\left(,\zeta ,\right)\right)}.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {Pr}(\varvec{N} = \varvec{n}) = \frac{{\int } h(\varvec{\zeta }) \, \text {Pr}(\varvec{N} = \varvec{n} \mid \varvec{\zeta }) \, p(\varvec{\zeta }) \, g(\varvec{\zeta }) \text {d}\varvec{\zeta }}{{\int } h(\varvec{\zeta }) \, g(\varvec{\zeta }) \, p(\varvec{\zeta } \mid \varvec{N} = \varvec{n}) \text {d}\varvec{\zeta }} = \frac{\mathbb {E}_{g(\varvec{\zeta })} \left[ h(\varvec{\zeta }) \, \text {Pr}(\varvec{N} = \varvec{n} \mid \varvec{\zeta }) \, p(\varvec{\zeta })\right] }{\mathbb {E}_{p(\varvec{\zeta } \mid \varvec{N} = \varvec{n})} \left[ h(\varvec{\zeta }) \, g(\varvec{\zeta })\right] }. \end{aligned}$$\end{document}

The bridge sampling estimate of the marginal likelihood is then obtained by sampling from $g\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\zeta })$$\end{document} and $p(\zeta \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\varvec{\zeta } \mid \varvec{N} = \varvec{n})$$\end{document} and then using Monte Carlo approximations to estimate the expected values.

The optimal choice of $h\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h(\varvec{\zeta })$$\end{document} , one that minimizes the relative mean-squared error of the estimator, is given by:

(14) $\begin{array}{c}{h}_o\left(\zeta \right)\propto {\left(s_1,,\text{Pr},(N=n\mid \zeta ),,p,\left(\zeta \right),+,s_2,,\text{Pr},(N=n),,g,\left(\zeta \right)\right)}^{-1},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} h_o(\varvec{\zeta }) \propto \left[ s_1 \, \text {Pr}(\varvec{N} = \varvec{n} \mid \varvec{\zeta }) \, p(\varvec{\zeta }) + s_2 \, \text {Pr}(\varvec{N} = \varvec{n}) \, g(\varvec{\zeta }) \right] ^{-1}, \end{aligned}$$\end{document}

where $s_i=\frac{D_i}{D_1+D_2},i\in \{1,2\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i = \frac{D_i}{D_1 + D_2}, \; i \in \{1,2\}$$\end{document} , where $D_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_1$$\end{document} and $D_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_2$$\end{document} denote the number of draws from $p(\zeta \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\varvec{\zeta } \mid \varvec{N} = \varvec{n})$$\end{document} and $g\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\zeta })$$\end{document} , respectively, used to approximate the expected values (Meng & Wong, Reference Meng and Wong1996). We set $D_1=D_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_1 = D_2$$\end{document} . Note that $h_o$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_o$$\end{document} is only optimal if the draws from the posterior distribution are independent which is not the case with MCMC procedures. To account for this fact, we replace $D_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_1$$\end{document} in defining the weights $s_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_1$$\end{document} and $s_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_2$$\end{document} by the effective sample size obtained using the codaR package (Plummer, Best, Cowles, & Vines, Reference Plummer, Best, Cowles and Vines2006).Footnote 8 As $h_o\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_o(\varvec{\zeta })$$\end{document} depends on $\text{Pr}(N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Pr}(\varvec{N} = \varvec{n})$$\end{document} , the very quantity we want to estimate, we follow Meng and Wong (Reference Meng and Wong1996) and use an iterative scheme to update an initial guess of the marginal likelihood until convergence:Footnote 9

(15) $\begin{array}{c}\widehat{\text{Pr}}{(N=n)}^{(t+1)}=\frac{\frac{1}{D_2}\sum _{r=1}^{D_2}\frac{l_{2,r}}{s_1{l}_{2,r}+s_2\widehat{\text{Pr}}{(N=n)}^{\left(t\right)}}}{\frac{1}{D_1}\sum _{j=1}^{D_1}\frac{1}{s_1{l}_{1,j}+s_2\widehat{\text{Pr}}{(N=n)}^{\left(t\right)}}},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \hat{\text {Pr}}(\varvec{N} = \varvec{n})^{(t+1)} = \frac{\frac{1}{D_2}\sum \limits _{r = 1}^{D_2}\frac{l_{2, r}}{s_1 l_{2, r} + s_2 \hat{\text {Pr}}(\varvec{N} = \varvec{n})^{(t)}}}{\frac{1}{D_1}\sum \limits _{j = 1}^{D_1}\frac{1}{s_1 l_{1, j} + s_2 \hat{\text {Pr}}(\varvec{N} = \varvec{n})^{(t)}}}, \end{aligned}$$\end{document}

where $l_{1,j}=\frac{\text{Pr}(N=n\mid {{\zeta }^{\ast }}_j)p\left({{\zeta }^{\ast }}_j\right)}{g\left({{\zeta }^{\ast }}_j\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l_{1,j} = \frac{\text {Pr}(\varvec{N} = \varvec{n} \mid \varvec{\zeta ^*}_j) \, p(\varvec{\zeta ^*}_j)}{g(\varvec{\zeta ^*}_j)}$$\end{document} , $l_{2,r}=\frac{\text{Pr}(N=n\mid {\tilde{\zeta }}_r)p\left({\tilde{\zeta }}_r\right)}{g\left({\tilde{\zeta }}_r\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l_{2,r} = \frac{\text {Pr}(\varvec{N} = \varvec{n} \mid \varvec{\tilde{\zeta }}_r) \, p(\varvec{\tilde{\zeta }}_r)}{g(\varvec{\tilde{\zeta }}_r)}$$\end{document} , $\{{{\zeta }^{\ast }}_1,\dots ,{{\zeta }^{\ast }}_{D_1}\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\varvec{\zeta ^*}_1, \ldots , \varvec{\zeta ^*}_{D_1}\}$$\end{document} are $D_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_1$$\end{document} draws from $p(\zeta \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\varvec{\zeta } \mid \varvec{N} = \varvec{n})$$\end{document} , and $\{{\tilde{\zeta }}_1,\dots ,{\tilde{\zeta }}_{D_2}\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\varvec{\tilde{\zeta }}_1, \ldots , \varvec{\tilde{\zeta }}_{D_2}\}$$\end{document} are $D_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_2$$\end{document} draws from $g\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\zeta })$$\end{document} .

A remaining question is how to choose $g\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\zeta })$$\end{document} . The precision of the bridge sampling estimator is governed by the number of samples from $g\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\zeta })$$\end{document} and the overlap between $g\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\zeta })$$\end{document} and $p(\zeta \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\varvec{\zeta } \mid \varvec{N} = \varvec{n})$$\end{document} (Meng & Wong, Reference Meng and Wong1996). Therefore, $g\left(\zeta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\zeta })$$\end{document} should closely resemble the posterior distribution. For instance, we may choose a multivariate normal distribution for g with mean vector and covariance matrix that match the corresponding quantities of the posterior samples. Although the multivariate normal approach works well in many applications (e.g., Gronau et al., Reference Gronau, Sarafoglou, Matzke, Ly, Boehm and Marsman2017; Overstall & Forster, Reference Overstall and Forster2010), it can be inefficient when the posterior distribution is skewed.

Warp-III improves upon the multivariate normal bridge sampling approach by matching, not only the first two, but also the third moment (i.e., skewness) of g and the posterior distribution. Consequently, in case there is no skewness, Warp-III results in estimates with the same precision as the ones from the simpler multivariate normal approach. However, crucially, in the presence of skewness, Warp-III is able to match g and the posterior distribution more closely which results in a higher precision of the marginal likelihood estimates compared to the simpler approach. How much of an improvement Warp-III is over the simpler multivariate normal approach may depend on the particular example at hand.

Figure 2. Matching the proposal and posterior distribution with warping. Histograms show the posterior distribution; density lines show the standard normal proposal distribution. Available at https://tinyurl.com/y7owvsz3 under CC license https://creativecommons.org/licenses/by/2.0/.

In Warp-III, g is fixed to a multivariate standard normal distribution. The posterior distribution is then manipulated—“warped”—so that its mean vector, covariance matrix, and skew match g. Crucially, the warped posterior distribution retains the normalizing constant of the posterior distribution. Figure 2 illustrates the rationale of the Warp-III transformation for the univariate case. The histogram in the upper-left panel shows hypothetical “unbounded” posterior samples that can range across the entire real line; the solid line shows the standard normal proposal distribution g. The overlap between the two distributions is clearly suboptimal. Bridge sampling applied to these two distributions can be thought of as “Warp-0” because the posterior distribution is not modified. The upper-right panel illustrates “Warp-I”: Subtracting the mean of the posterior samples from all posterior samples matches the first moment of the distributions. The lower-right panel illustrates “Warp-II”: Dividing the zero-centered posterior samples by their standard deviation matches the first two moments of the distributions. This approach is practically equivalent to the multivariate normal bridge sampling approach described above. Lastly, the lower-left panel illustrates Warp-III: Randomly assigning a minus sign to the standardized posterior samples matches also the third moment of the distributions.

Warp-III assumes that all components of the parameter vector can range across the entire real line. In the context of latent-trait MPTs, this assumption is not fulfilled since ${\xi }_p\in (0,{\xi }_{\text{max}})\forall p\in \{1,\dots ,P\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _p \in (0, \xi _{\text {max}}) \; \forall p \in \{1, \ldots , P\}$$\end{document} . We therefore transform $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} so that ${\xi }_{\text{trans}}={\Phi }^{-1}\left(\frac{\xi }{{\xi }_{\text{max}}}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }_{\text {trans}} = \Phi ^{-1}\left( \frac{\varvec{\xi }}{\xi _{\text {max}}}\right) $$\end{document} with Jacobian ${\left({\xi }_{\text{max}}\right)}^P{N}_P({\xi }_{\text{trans}};0,I_P)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( \xi _{\text {max}}\right) ^{P} \mathcal {N}_P(\varvec{\xi }_{\text {trans}}; \varvec{0}, \varvec{I}_P)$$\end{document} , where $N_P(x;y,Z)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {N}_P(\varvec{x}; \varvec{y}, \varvec{Z})$$\end{document} denotes the probability density function of a P-dimensional normal distribution with mean vector $y$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{y}$$\end{document} and covariance matrix $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Z}$$\end{document} which is evaluated for the vector $x$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}$$\end{document} .Footnote 10 Let $\psi =(\mu ,{\xi }_{\text{trans}},{\omega }_1,\dots ,{\omega }_I)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\psi } = (\varvec{\mu }, \varvec{\xi }_{\text {trans}} , \varvec{\omega }_1, \ldots , \varvec{\omega }_I)$$\end{document} denote the resulting parameter vector where all components are on the real line.

Warp-III is then based on applying the following stochastic transformation to $\psi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\psi }$$\end{document} :

(16) $\begin{array}{c}\eta =\underset{\text{symmetry}}{\underbrace{b}}\times \underset{\text{covariance}I}{\underbrace{R^{-1}}}\times \underset{\text{mean}0}{\underbrace{(\psi -v)}},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\eta } = \underbrace{b}_{\textstyle \text {symmetry}} \times \underbrace{\varvec{R}^{-1}}_{\textstyle \text {covariance } \varvec{I}} \times \;\; \underbrace{(\varvec{\psi } - \varvec{v})}_{\textstyle \text {mean } \varvec{0}}, \end{aligned}$$\end{document}

where $b\sim \text{Bernoulli}\left(0.5\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b \sim \text {Bernoulli}(0.5)$$\end{document} on $\{-1,1\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{-1, 1\}$$\end{document} and $v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{v}$$\end{document} corresponds to the expected value of $\psi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\psi }$$\end{document} (i.e., the mean vector). The matrix $R$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{R}$$\end{document} is obtained via the Cholesky decomposition of the covariance matrix of $\psi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\psi }$$\end{document} , denoted as $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{S}$$\end{document} , thus, $S=R{R}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{S} = \varvec{R} \varvec{R}^\top $$\end{document} . In practice, $v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{v}$$\end{document} and $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{S}$$\end{document} are unknown and must be approximated using the posterior samples. Note that Eq. (16) simply generalizes the intuition illustrated in Fig. 2 for the univariate case to the general case with multiple parameters.

Due to the Bernoulli random variable b, the warped posterior density has the form of a mixture density (see also Overstall, Reference Overstall2010, p. 70):

(17) $\begin{array}{c}\begin{array}{cc}{p}_{\eta }(\eta \mid N=n)& =\frac{\left(R\right)}{2}\left(\frac{{\tilde{p}}_{\psi }(v-R\eta \mid N=n)}{\text{Pr}(N=n)},+,\frac{{\tilde{p}}_{\psi }(v+R\eta \mid N=n)}{\text{Pr}(N=n)}\right)\\ & =\frac{{\tilde{p}}_{\eta }(\eta \mid N=n)}{\text{Pr}(N=n)},\end{array}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} p_{\varvec{\eta }}(\varvec{\eta } \mid \varvec{N} = \varvec{n})&= \frac{\left| \varvec{R}\right| }{2} \left[ \frac{\tilde{p}_{\varvec{\psi }}(\varvec{v} - \varvec{R}\varvec{\eta } \mid \varvec{N} = \varvec{n})}{\text {Pr}(\varvec{N} = \varvec{n})} + \frac{\tilde{p}_{\varvec{\psi }}(\varvec{v} + \varvec{R}\varvec{\eta } \mid \varvec{N} = \varvec{n})}{\text {Pr}(\varvec{N} = \varvec{n})}\right] \\&= \frac{\tilde{p}_{\varvec{\eta }}(\varvec{\eta } \mid \varvec{N} = \varvec{n})}{\text {Pr}(\varvec{N} = \varvec{n})}, \end{aligned} \end{aligned}$$\end{document}

where ${\tilde{p}}_{\eta }(\eta \mid N=n)=\frac{\left(R\right)}{2}\left({\tilde{p}}_{\psi },(v-R\eta \mid N=n),+,{\tilde{p}}_{\psi },(v+R\eta \mid N=n)\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{p}_{\varvec{\eta }}(\varvec{\eta } \mid \varvec{N} = \varvec{n}) = \frac{\left| \varvec{R}\right| }{2} \left[ \tilde{p}_{\varvec{\psi }}(\varvec{v} - \varvec{R}\varvec{\eta } \mid \varvec{N} = \varvec{n}) + \tilde{p}_{\varvec{\psi }}(\varvec{v} + \varvec{R}\varvec{\eta } \mid \varvec{N} = \varvec{n})\right] $$\end{document} denotes the un-normalized warped posterior distribution and ${\tilde{p}}_{\psi }(·\mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{p}_{\varvec{\psi }}(\cdot \mid \varvec{N} = \varvec{n})$$\end{document} denotes the un-normalized posterior distribution that has been transformed to the real line (but not warped). This proves that the warped posterior distribution retains the normalizing constant of the original posterior distribution.

The Warp-III estimator of the marginal likelihood is then derived by using the warped posterior distribution $p_{\eta }(\eta \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\varvec{\eta }}(\varvec{\eta } \mid \varvec{N} = \varvec{n})$$\end{document} instead of $p(\zeta \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\varvec{\zeta } \mid \varvec{N} = \varvec{n})$$\end{document} in Eq. (12). Equation (13) shows that this results in a ratio of two expected values, where the numerator is an expected value with respect to the multivariate standard normal proposal distribution $g\left(\eta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\eta })$$\end{document} and the denominator is an expected value with respect to the warped posterior distribution $p_{\eta }(\eta \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\varvec{\eta }}(\varvec{\eta } \mid \varvec{N} = \varvec{n})$$\end{document} . Hence, we could obtain an estimate of the marginal likelihood by first warping the posterior samples using Eq. (16), then sampling from the proposal distribution, and applying the iterative updating scheme in Eq. (15).

However, in line with the literature (e.g., Sinharay & Stern, Reference Sinharay and Stern2005), we rewrite the expected value in the denominator of Eq. (13) in terms of the unbounded posterior samples that are transformed to the real line but are not warped; a derivation is provided in Supplemental Materials. The estimate of the marginal likelihood is then obtained by applying the iterative scheme in Eq. (15) using:

(18) $\begin{array}{c}{l}_{1,j}=\frac{\frac{\left(R\right)}{2}\left({\tilde{p}}_{\psi },(2v-{{\psi }^{\ast }}_j\mid N=n),+,{\tilde{p}}_{\psi },({{\psi }^{\ast }}_j\mid N=n)\right)}{g\left(R^{-1},\left({{\psi }^{\ast }}_j,-,v\right)\right)},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} l_{1,j} = \frac{\frac{\left| \varvec{R}\right| }{2} \left[ \tilde{p}_{\varvec{\psi }}(2\varvec{v} - \varvec{\psi ^*}_j \mid \varvec{N} = \varvec{n}) + \tilde{p}_{\varvec{\psi }}(\varvec{\psi ^*}_j \mid \varvec{N} = \varvec{n})\right] }{g\left( \varvec{R}^{-1}\left( \varvec{\psi ^*}_j - \varvec{v}\right) \right) }, \end{aligned}$$\end{document}

and

(19) $\begin{array}{c}{l}_{2,r}=\frac{\frac{\left(R\right)}{2}\left({\tilde{p}}_{\psi },(v-R\widetilde{{\eta }_r}\mid N=n),+,{\tilde{p}}_{\psi },(v+R\widetilde{{\eta }_r}\mid N=n)\right)}{g\left({\tilde{\eta }}_r\right)},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} l_{2,r} =\frac{\frac{\left| \varvec{R}\right| }{2} \left[ \tilde{p}_{\varvec{\psi }}(\varvec{v} - \varvec{R}\tilde{\varvec{\eta }_r} \mid \varvec{N} = \varvec{n}) + \tilde{p}_{\varvec{\psi }}(\varvec{v} + \varvec{R}\tilde{\varvec{\eta }_r} \mid \varvec{N} = \varvec{n})\right] }{g(\tilde{\varvec{\eta }}_r)}, \end{aligned}$$\end{document}

where $\{{{\psi }^{\ast }}_1,\dots ,{{\psi }^{\ast }}_{D_1}\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\varvec{\psi ^*}_1, \ldots , \varvec{\psi ^*}_{D_1}\}$$\end{document} are $D_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_1$$\end{document} draws from $p_{\psi }(\psi \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\varvec{\psi }}(\varvec{\psi } \mid \varvec{N} = \varvec{n})$$\end{document} and $\{{\tilde{\eta }}_1,\dots ,{\tilde{\eta }}_{D_2}\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\varvec{\tilde{\eta }}_1, \ldots , \varvec{\tilde{\eta }}_{D_2}\}$$\end{document} are $D_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_2$$\end{document} draws from the proposal distribution $g\left(\eta \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\varvec{\eta })$$\end{document} . Furthermore, ${\tilde{p}}_{\psi }(\psi \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{p}_{\varvec{\psi }}(\varvec{\psi }\mid \varvec{N} = \varvec{n})$$\end{document} denotes the un-normalized posterior density of the unbounded posterior samples; it is therefore written in terms of ${\xi }_{\text{trans}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }_{\text {trans}}$$\end{document} and is adjusted by the Jacobian term:Footnote 11

(20) $\begin{array}{cc}{\tilde{p}}_{\psi }(\psi \mid N=n)=& \prod _{i=1}^I\left(\prod _{k=1}^K,\{,\frac{J_k!}{n_{ik1}!\times n_{ik2}!\times ...\times n_{ik{L}_k}!},\prod _{l=1}^{L_k},{\left(\text{Pr},(,C_{\mathrm{kl}},\mid ,\mu ,,,{\xi }_{\text{trans}},,,{\omega }_i,)\right)}^{n_{\mathrm{ikl}}},\}\right)\\ & \times \frac{{\Gamma }_P\left(\frac{\nu +I}{2}\right)}{{\Gamma }_P\left(\frac{\nu }{2}\right)}\frac{{\pi }^{-\frac{\mathrm{IP}}{2}}}{{\left({\Omega }^{\top }\Omega +I_P\right)}^{\frac{\nu +I}{2}}}\times {\left(2,\pi \right)}^{-\frac{P}{2}}exp\{-\frac{1}{2}{\mu }^{\top }\mu \}\\ & \times {\left(2,\pi \right)}^{-\frac{P}{2}}exp\{-\frac{1}{2}{\xi }_{\text{trans}}^{\top }{\xi }_{\text{trans}}\}.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \tilde{p}_{\varvec{\psi }}(\varvec{\psi }\mid \varvec{N} = \varvec{n})= & {} \prod _{i = 1}^{I} \left[ \prod ^K_{k=1} \Bigg \{{\frac{J_k!}{n_{ik1}! \times n_{ik2}! \times ... \times n_{ikL_k}!}} \prod ^{L_k}_{l=1}\left[ \text {Pr}(C_{kl} \mid \varvec{\mu }, \varvec{\xi }_{\text {trans}}, \varvec{\omega }_i)\right] ^{n_{ikl}}\Bigg \}\right] \nonumber \\&\quad \times \frac{\Gamma _P(\frac{\nu + I}{2})}{\Gamma _P(\frac{\nu }{2})} \frac{\pi ^{-\frac{I P}{2}}}{\left| { \varvec{\Omega }^\top \varvec{\Omega }+ {\varvec{I}_P}}\right| ^{\frac{\nu + I}{2}}} \times \left( 2\pi \right) ^{-\frac{P}{2}} \exp \bigg \{-\frac{1}{2} \varvec{\mu }^\top \varvec{\mu }\bigg \}\nonumber \\&\quad \times \left( 2\pi \right) ^{-\frac{P}{2}} \exp \bigg \{-\frac{1}{2} \varvec{\xi }_{\text {trans}}^\top \varvec{\xi }_{\text {trans}}\bigg \}. \end{aligned}$$\end{document}

Note that rewriting the expected value in terms of ${\tilde{p}}_{\psi }(\psi \mid N=n)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{p}_{\varvec{\psi }}(\varvec{\psi }\mid \varvec{N} = \varvec{n})$$\end{document} is only a technical nicety. This approach is identical to applying the Warp-III transformation to the posterior samples and then using the iterative scheme with the warped posterior density and a multivariate standard normal proposal distribution.

3. Empirical Examples

3.1 Example 1: Nested Model Comparison

We re-analyzed the pair-clustering data set reported in Riefer, Knapp, Batchelder, Bamber, and Manifold (Reference Riefer, Knapp, Batchelder, Bamber and Manifold2002) using the hierarchical latent-trait approach.Footnote 12 Experiment 4 examined the memory of patients with brain damage due to prolonged alcoholism in comparison with a control group of alcoholic patients without indications of brain damage. The participants attempted to memorize the same list of 20 categorically related word pairs in a series of six study-test trials.Footnote 13 For demonstration purposes, we focused on the free recall performance of the 21 control participants. Specifically, we investigated whether the model parameters change from the first to the second trial indicating a change in the storage and retrieval processes as a function of practice using posterior model probabilities and Bayesian model averaging.

3.1.1. Model Specification

To model differences in parameters, we augmented Eq. (8) with a parameter vector that captures the difference in parameters between the two trials: $\delta =({\delta }_c,{\delta }_r,{\delta }_u)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\delta } = (\delta _c, \delta _r, \delta _u)$$\end{document} . The probit-transformed parameter vectors of participant i for the first trial ( ${\theta }_{1,i}^{{}^{\prime }}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_{1,i}^{'}$$\end{document} ) and the second trial ( ${\theta }_{2,i}^{{}^{\prime }}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_{2,i}^{'}$$\end{document} ) are then obtained as follows:

(21) $\begin{array}{c}\begin{array}{cc}{\theta }_{1,i}^{{}^{\prime }}& =\stackrel{\begin{array}{c}\text{group mean}\\ \text{for first trial}\end{array}}{\overbrace{\mu -\frac{\delta }{2}}}+\xi \odot {\omega }_i,\\ {\theta }_{2,i}^{{}^{\prime }}& =\underset{\begin{array}{c}\text{group mean}\\ \text{for second trial}\end{array}}{\underbrace{\mu +\frac{\delta }{2}}}+\xi \odot {\omega }_i.\end{array}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \varvec{\theta }_{1,i}^{'}&= \; \overbrace{\varvec{\mu } - \frac{\varvec{\delta }}{2}}^{\begin{array}{c} \text {group mean} \\ \text {for first trial} \end{array}} \;\,\,+ \;\;\; \varvec{\xi } \odot \varvec{\omega }_i,\\ \varvec{\theta }_{2,i}^{'}&= \underbrace{\varvec{\mu } + \frac{\varvec{\delta }}{2}}_{\begin{array}{c} \text {group mean} \\ \text {for second trial} \end{array}} + \;\;\; \varvec{\xi } \odot \varvec{\omega }_i. \end{aligned} \end{aligned}$$\end{document}

For an alternative approach to modeling within-subject differences in model parameters, the reader is referred to Rouder, Lu, Morey, Sun, and Speckman (Reference Rouder, Lu, Morey, Sun and Speckman2008).

Table 1 shows the $2^3=8$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2^3 = 8$$\end{document} nested models that implement the eight sets of possible parameter constraints. $M_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_1$$\end{document} allows all three parameters to vary between trials so that $\delta =({\delta }_c,{\delta }_r,{\delta }_u)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\delta } = (\delta _c, \delta _r, \delta _u)$$\end{document} . In contrast, $M_8$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_8$$\end{document} posits that none of the parameters vary between trials so that $\delta =(0,0,0)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\delta } = (0, 0, 0)$$\end{document} . Models $M_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_2$$\end{document} to $M_7$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_7$$\end{document} are between these extremes and allow either one or two parameters to vary between trials.

Table 1. Overview of the eight nested models for the analysis of the first two trials of the pair-clustering data set reported in Riefer et al. (Reference Riefer, Knapp, Batchelder, Bamber and Manifold2002).

Note. $M_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_1$$\end{document} allows all three parameters to vary between trials, and $M_8$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_8$$\end{document} posits that none of the parameters vary between trials. Models $M_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_2$$\end{document} to $M_7$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_7$$\end{document} are between these extremes.

We used independent zero-centered normal priors for the components of $\delta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\delta }$$\end{document} . We explored a narrow ( ${\sigma }_{\delta }^{\text{narrow}}\approx 0.52$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma ^\text {narrow}_{\delta } \approx 0.52$$\end{document} ), medium ( ${\sigma }_{\delta }^{\text{medium}}\approx 0.84$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma ^\text {medium}_{\delta } \approx 0.84$$\end{document} ), and a wide ( ${\sigma }_{\delta }^{\text{wide}}\approx 1.28$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma ^\text {wide}_{\delta } \approx 1.28$$\end{document} ) zero-centered normal prior to assess the sensitivity of the results to the width of the test-relevant prior distribution. As shown in Supplemental Materials, the standard deviations ${\sigma }_{\delta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\delta }$$\end{document} were chosen to correspond to small, medium, and large effects on the probability scale centered around 0.5. Priors for the remaining parameters followed the specification described earlier.

We estimated the posterior distribution of the model parameters using JAGS by adapting the script provided by Matzke et al. (Reference Matzke, Dolan, Batchelder and Wagenmakers2015). The JAGS code is available in Supplemental Materials. We ran three MCMC chains with overdispersed start values, discarded the first 4000 posterior samples as burn in, and retained only every 20th sample to reduce autocorrelation. Results reported below are based on a total of 90,000 posterior samples. Convergence of the MCMC chains was assessed by visual inspection and the $\hat{R}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{R}$$\end{document} statistic ( $\hat{R}<1.05$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{R}<1.05$$\end{document} for all parameters; Gelman & Rubin, Reference Gelman and Rubin1992).

Figure 3. Posterior distributions of the probit group-level means (plotted on the probability scale) from the full model $M_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_1$$\end{document} for the analysis of the first two trials of the pair-clustering data reported in Riefer et al. (Reference Riefer, Knapp, Batchelder, Bamber and Manifold2002). The solid lines correspond to the posteriors for the first trial and the dotted lines to the posteriors for the second trial. Available at https://tinyurl.com/y9a33l4t under CC license https://creativecommons.org/licenses/by/2.0/.

Figure 3 shows the resulting posterior distributions of the probit group-level means from the full model $M_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_1$$\end{document} ; the parameters were transformed back to the probability scale. The posteriors were computed using the medium prior setting ( ${\sigma }_{\delta }^{\text{medium}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma ^\text {medium}_{\delta }$$\end{document} )—results obtained with the narrow and wide prior were highly similar and are not displayed. The plot of the posterior distributions based on the alternative prior choice for the elements of $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} (i.e., uniform priors with upper bound ${\xi }_{\text{max}}=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 2$$\end{document} instead of ${\xi }_{\text{max}}=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 10$$\end{document} ) was visually almost indistinguishable from the one presented here and has hence been relegated to Supplemental Materials. The cluster-storage c parameter did not change substantially, whereas the storage-retrieval u, and especially the cluster-retrieval r parameter, seemed to increase from the first trial to the second.

3.1.2. Computing Marginal Likelihoods with Warp-III

Equation (20) was adjusted to include the relevant prior distributions for the elements of $\delta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\delta }$$\end{document} . For each model, we split the 90,000 posterior samples in two equal parts (first and second half of the iterations per chain) and used the first part for estimating $R$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{R}$$\end{document} and $v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{v}$$\end{document} and the second part for the iterative updating scheme in Eq. (15) (Overstall & Forster, Reference Overstall and Forster2010). Hence, $D_1=D_2=45,000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_1 = D_2 = 45,000$$\end{document} . To assess the accuracy of the resulting estimates, we repeated this procedure 50 times.Footnote 14 We implemented the procedure in R (R Core Team, 2016). For efficiency, we parallelized the computations, and coded the computationally intensive elements in efficient C++ code which was called from within R using Rcpp (Eddelbuettel et al., Reference Eddelbuettel, François, Allaire, Chambers, Bates and Ushey2011). Using a standard personal computer and four CPU cores, computing the marginal likelihood for each repetition took less than one minute per model. The code is available in Supplemental Materials.

Figure 4. Posterior model probabilities (left panel) and posterior inclusion probabilities (right panel) for the analysis of the first two trials of the pair-clustering data reported in Riefer et al. (Reference Riefer, Knapp, Batchelder, Bamber and Manifold2002) obtained with Warp-III bridge sampling. In the left panel, the x-axis indicates which parameters were allowed to vary from the first to the second trial (e.g., $c-u$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c-u$$\end{document} corresponds to $M_3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_3$$\end{document} where r was fixed between trials). Gray symbols show the results of the 50 repetitions, and black symbols display the posterior model probabilities and posterior inclusion probabilities that are based on the median of the 50 estimated log marginal likelihoods. Circles show results obtained with the narrow prior, diamonds with the medium prior, and triangles with the wide prior. The dotted lines show the prior model probabilities and prior inclusion probabilities. Available at https://tinyurl.com/yaxbj9o6 under CC license https://creativecommons.org/licenses/by/2.0/.

3.1.3. Posterior Model Probabilities

To formally quantify evidence for the differences in parameters, we computed the posterior model probabilities of the eight models using the marginal likelihoods obtained with Warp-III. We assumed that all models were equally likely a priori. The left panel of Fig. 4 shows the posterior model probabilities for the narrow, medium, and wide prior settings. The plot of the posterior model probabilities based on the alternative prior choice for the elements of $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} (i.e., uniform priors with upper bound ${\xi }_{\text{max}}=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 2$$\end{document} instead of ${\xi }_{\text{max}}=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 10$$\end{document} ) was visually almost indistinguishable from the one presented here and has hence been relegated to Supplemental Materials. Formal model comparison confirmed the results of the visual inspection of the posterior distributions shown in Fig. 3: $M_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_2$$\end{document} , the model that allows for a difference in r and u, received the most support from the data. As expected, the width of the test-relevant prior $\delta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\delta }$$\end{document} influenced the value of the marginal likelihood, but it did not change the conclusions qualitatively. Warp-III provided accurate estimates of the posterior model probabilities as indicated by the small variability across the 50 repetitions (i.e., gray symbols). For this nested example, the posterior model probabilities can be also obtained using the Savage–Dickey density ratio representation of the Bayes factor (Dickey & Lientz, Reference Dickey and Lientz1970; Wagenmakers, Lodewyckx, Kuriyal, & Grasman, Reference Wagenmakers, Lodewyckx, Kuriyal and Grasman2010). As shown in Supplemental Materials, the Savage–Dickey procedure resulted in posterior model probabilities that were highly similar to the ones obtained with Warp-III.

3.1.4. Bayesian Model Averaging

Bayesian model averaging does not require researchers to commit to a single “best” model; it allows researchers to acknowledge uncertainty about the choice of the correct model (e.g., Hoeting et al., Reference Hoeting, Madigan, Raftery and Volinsky1999; Rouder et al., Reference Rouder, Morey, Verhagen, Swagman and Wagenmakers2017). This is achieved by considering the posterior inclusion probabilities of the parameters. Posterior inclusion probabilities quantify the model-averaged evidence for a change in a given parameter; they can be obtained by summing the posterior model probabilities of the models that allow the parameter to differ between the trials. For instance, the posterior inclusion probability of the c parameter is obtained by summing the posterior model probabilities of $M_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_1$$\end{document} , $M_3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_3$$\end{document} , $M_4$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_4$$\end{document} , and $M_6$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_6$$\end{document} . Posterior inclusion probabilities are then compared to the prior inclusion probabilities, in this case 0.5, which are obtained in an analogous manner but based on the prior model probabilities.Footnote 15 The right panel of Fig. 4 shows the posterior inclusion probabilities for the three prior settings. The plot of the posterior inclusion probabilities based on the alternative prior choice for the elements of $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} (i.e., uniform priors with upper bound ${\xi }_{\text{max}}=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 2$$\end{document} instead of ${\xi }_{\text{max}}=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 10$$\end{document} ) was visually almost indistinguishable from the one presented here and has hence been relegated to Supplemental Materials. The posterior inclusion probabilities of the r and u parameter are higher than the prior inclusion probabilities, indicating evidence for a difference in these parameters between trials. In contrast, the posterior inclusion probability of c is lower than the corresponding prior inclusion probability, indicating evidence for invariance between the trials. As before, the width of the $\delta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\delta }$$\end{document} prior does not change the conclusions qualitatively.

3.1.5. Substantive Contribution

The data from Riefer et al. (Reference Riefer, Knapp, Batchelder, Bamber and Manifold2002) have been analyzed in a number of articles. The original article analyzed the aggregated data (an approach known to suffer from limitations in case there is heterogeneity across participants, e.g., Klauer, Reference Klauer2006) and considered the p values of $G^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G^2$$\end{document} statistics to investigate whether parameters differ across trials. Smith and Batchelder (Reference Smith and Batchelder2010) re-analyzed a subset of the data using the hierarchical beta-MPT model (which specifies group-level beta distributions and thus differs from the latent-trait approach that we used).Footnote 16 To investigate whether parameters differ across trials, Smith and Batchelder (a) considered the posterior distribution of the difference between trials for the group-level mean parameters and (b) ran a classical paired sample t test on the individual-level parameter estimates. These approaches, however, do not allow one to quantify evidence for an invariance (i.e., a simpler model where some parameters do not differ across trials) on a continuous scale in a systematic way and, crucially, they do not allow one to disentangle “absence of evidence” (i.e., the data are uninformative) and “evidence of absence” (i.e., the data support a simpler model).Footnote 17 These shortcomings can be addressed by computing Bayes factors and posterior model and posterior inclusion probabilities. “Absence of evidence” can be inferred from Bayes factors close to one and posterior model and posterior inclusion probabilities close to the corresponding prior probabilities. In contrast, “evidence of absence” can be inferred from large Bayes factors in favor of the simpler model, and in situations when the posterior model probability of the simpler model is the highest or when the posterior inclusion probability is smaller than the prior inclusion probability.

Our Bayesian re-analysis suggests that there is strong evidence that the probability of retrieving word pairs that have been stored as a cluster (i.e., r) changed from the first to the second trial. Furthermore, there is evidence that the probability of storing and retrieving words that have not been stored as a cluster (i.e., u) differed between the two trials. Crucially, our approach also allowed us to conclude that there is some evidence that the probability of storing a word pair as a cluster (i.e., c) did not change from the first to the second trial (although this evidence is not that pronounced since the posterior inclusion probability for a difference in c is—depending on the prior choice—relatively close to the prior inclusion probability of .5). Another key improvement in our analysis over the above-mentioned analyses is the use of Bayesian model averaging. In this example, $M_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_2$$\end{document} received the highest posterior probability; however, $M_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_1$$\end{document} also received substantive posterior probability. Therefore, selecting a single best model (i.e., $M_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}_2$$\end{document} ) and basing final inference solely on this model might be suboptimal at best and misleading at worst. In contrast, when using the model-averaged posterior inclusion probabilities for drawing conclusions about which parameters differ between trials, one takes into account all models under consideration according to their plausibilities in light of the observed data.

Finally, note that one might argue that this data set is relatively small and is thus uninformative. However, one strength of the Bayesian approach is that it allows one to quantify whether the data are informative or not. For this example, the Bayesian results suggest that the data are in fact informative which is indicated by posterior model/inclusion probabilities that are quite different from the corresponding prior probabilities.

3.2. Example 2: Non-Nested Model Comparison

We re-analyzed data from Experiment 2 reported by Fazio et al. (Reference Fazio, Brashier, Payne and Marsh2015) who investigated the influence of knowledge on the illusory truth effect. The illusory truth effect refers to the phenomenon that, in the absence of knowledge about the truth status of a statement, repeated statements are easier to process and are judged more truthful than new statements. Fazio et al., however, provided evidence that participants tend to rely on the ease of processing (i.e., fluency) even when they have knowledge about the statement.

We re-analyzed data from 39 participants who indicated the truthfulness (i.e., “true”/“false”) of 176 statements, half of which were true and half of which were false. Half of the statements were likely to be known according to general knowledge norms (“known” statements) and half of them were likely to be unknown (“unknown” statements). An example of a true known statement is “The Pacific Ocean is the largest ocean on Earth.” An example of a false unknown statement is “Billy the Kid’s last name is Garrett.” To manipulate fluency, half of the statements were presented twice, once in the exposure phase and once in the truth-rating phase, whereas the other half was only presented in the truth-rating phase. Hence, the experiment had a 2 (truth status: true vs. false) $\times$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 2 (assumed knowledge: known vs. unknown) $\times$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 2 (repetition: repeated vs. not repeated) balanced within-subject design, and each cell of the design featured 22 statements.

Figure 5. Knowledge-conditional (top panel) and fluency-conditional (bottom panel) MPTs. Available at https://tinyurl.com/ya8sovfr under CC license https://creativecommons.org/licenses/by/2.0/.

3.2.1. Model Specification

Fazio et al. (Reference Fazio, Brashier, Payne and Marsh2015) constructed two MPTs to study the illusory truth effect. The knowledge-conditional model depicted in the top panel of Fig. 5 assumes that participants rely on knowledge when assessing truthfulness and only rely on fluency when they are unable to retrieve knowledge about the statement. Parameter k represents the probability of retrieving knowledge about the statement from memory. If knowledge is retrieved, participants are assumed to give the correct response (i.e., “true” for true statements and “false” for false statements). If no knowledge is retrieved with probability $1-k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - k$$\end{document} , participants rely on fluency with probability f and respond “true.” If participants do not rely on fluency with probability $1-f$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - f$$\end{document} , they guess “true” with probability g and “false” with probability $1-g$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - g$$\end{document} . Responses to true statements are scored into the categories $C_{11}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{11}$$\end{document} (correct “true” response) and $C_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{12}$$\end{document} (incorrect “false” response). Responses to false statements are scored into the categories $C_{21}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{21}$$\end{document} (incorrect “true” response) and $C_{22}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{22}$$\end{document} (correct “false” response). In contrast, the fluency-conditional model depicted in the bottom panel reflects the notion that participants mainly rely on fluency and only use knowledge in the absence of fluency. The models feature the same set of parameters, but they assume a different conditional probability structure.

For each model, we replicated the two subtrees four times (i.e., a total of eight subtrees per model) to accommodate the design of the experiment: The first replicate corresponded to known true and false statements that were not repeated, the second to known true and false statements that were repeated, the third to unknown true and false statements that were not repeated, and the fourth to unknown true and false statements that were repeated. Following Fazio et al. (Reference Fazio, Brashier, Payne and Marsh2015), we used separate knowledge parameters for known ( $k_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_k$$\end{document} ) and unknown ( $k_u$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_u$$\end{document} ) statements, and separate fluency parameters for repeated statements ( $f_r$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_r$$\end{document} ) and statements shown only once ( $f_n$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_n$$\end{document} ). The guessing parameter g was constrained to be equal across the four replicates. We implemented the models within the hierarchical latent-trait approach, using the prior specifications described earlier.

We estimated the posterior distribution of the model parameters using JAGS, ran three MCMC chains with overdispersed start values, discarded the first 4000 posterior samples as burn in, and retained only every 50th sample. Results reported below are based on a total of 180, 000 posterior samples. The posterior distributions of the group-level mean parameters are displayed in Supplemental Materials.

3.2.2. Computing Bayes Factors with Warp-III

For each model, we split the 180, 000 posterior samples in two equal parts (first and second half of the iterations per chain) and used the first part for estimating $R$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{R}$$\end{document} and $v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{v}$$\end{document} and the second part for the iterative updating scheme in Eq. (15) ( $D_1=D_2=90,000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_1 = D_2 = 90,000$$\end{document} ). Using a standard personal computer and four CPU cores, computing the marginal likelihood took approximately three minutes per model.

The resulting marginal likelihoods were used to compute the Bayes factor in favor of the fluency-conditional model over the knowledge-conditional model. To assess the accuracy of the resulting Bayes factor, we repeated this procedure 50 times. Estimates of the Bayes factor ranged from $1.3\times {10}^{42}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.3 \times 10^{42}$$\end{document} to $3.6\times {10}^{43}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3.6 \times 10^{43}$$\end{document} in favor of the fluency-conditional model. Estimates of the Bayes factor based on the alternative prior choice for the elements of $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} (i.e., uniform priors with upper bound ${\xi }_{\text{max}}=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 2$$\end{document} instead of ${\xi }_{\text{max}}=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{\text {max}} = 10$$\end{document} ) ranged from $1.7\times {10}^{41}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.7 \times 10^{41}$$\end{document} to $1.7\times {10}^{43}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.7 \times 10^{43}$$\end{document} in favor of the fluency-conditional model. In line with the conclusion drawn by Fazio et al. (Reference Fazio, Brashier, Payne and Marsh2015) based on the $G^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G^2$$\end{document} statistic, this result provides overwhelming evidence in favor of the fluency-conditional model.Footnote 18

Figure 6 displays the Warp-III Bayes factor estimates (on the log scale) in white as a function of the number of posterior samples used in the bridge sampling procedure.Footnote 19 As a comparison, the estimates based on the simpler multivariate normal bridge sampling approach are displayed in gray. As the number of posterior samples increases, the Bayes factor estimates become more precise. For this particular example, it is apparent that the Warp-III estimates are less variable than the estimates based on the simpler multivariate normal approach.

Figure 6. Log Bayes factor estimates in favor of the fluency-conditional (FC) model over the knowledge-conditional (KC) model as a function of the number of posterior samples. The Warp-III estimates are displayed in white, and the estimates based on the simpler multivariate normal approach are displayed in gray. Available at https://tinyurl.com/ydbfev7w under CC license https://creativecommons.org/licenses/by/2.0/.

3.2.3. Substantive Contribution

The authors of the original article analyzed the aggregated data (again, an approach known to be suboptimal in case there is heterogeneity across participants) and considered the $G^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G^2$$\end{document} statistics with corresponding p values. Based on the fact that the knowledge-conditional model had a larger, significant $G^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G^2$$\end{document} statistic compared to the fluency-conditional model that had a lower, nonsignificant $G^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G^2$$\end{document} statistic, the authors concluded that the knowledge-conditional model fit the data poorly and the fluency-conditional model fit the data well. Therefore, the authors favored the fluency-conditional model based on two binary accept–reject decisions. This makes it difficult to gauge the degree of support that the data provide in favor of the fluency-conditional model. The Bayes factor may be 10, or 100, or 1000—these are very different levels of evidence. In fact, our analysis shows that the Bayes factor is about $1.3\times {10}^{42}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.3 \times 10^{42}$$\end{document} – $3.6\times {10}^{43}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3.6 \times 10^{43}$$\end{document} in favor of the fluency-conditional model, which represents an overwhelming amount of evidence.

It could be argued that, since the compared models have the same number of parameters, comparing $G^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G^2$$\end{document} statistics may result in choosing the same model as based on considering AIC or BIC. AIC is asymptotically equivalent to cross-validation (Stone, Reference Stone1977) which is known to be inconsistent in the sense that, when the number of observations goes to infinity, the data-generating model will not be chosen with certainty (Shao, Reference Shao1993). In contrast, when using Bayes factors, model selection consistency is generally fulfilled (Bayarri, Berger, Forte, & García-Donato, Reference Bayarri, Berger, Forte and García-Donato2012). Although the BIC is a rough approximation of the Bayes factor, we believe that it is better to compute proper Bayes factors which are transparent with respect to the prior assumptions.

Finally, one might argue again that this data set is relatively small and is thus uninformative. However, the resulting Bayes factor is very different from 1, indicating that the data are in fact highly informative with respect to adjudicating between the fluency-conditional and the knowledge-conditional models.

5. Data Availability Statement

The datasets analyzed during the current study are available on the Open Science Framework: https://osf.io/rycg6/.