Quantifying Uncertainty in Transdimensional Markov Chain Monte Carlo Using Discrete Markov Models

TL;DR

This paper introduces a novel method to quantify the uncertainty in model probabilities estimated by transdimensional MCMC, using a Markov model approach to improve accuracy and assess effective sample size.

Contribution

It proposes a new approach to estimate the precision of model probabilities in transdimensional MCMC using a transition matrix-based Markov model.

Findings

Accurately assesses uncertainty in posterior model probabilities

Provides estimates for effective sample size of MCMC output

Demonstrates effectiveness in two model-selection examples

Abstract

Bayesian analysis often concerns an evaluation of models with different dimensionality as is necessary in, for example, model selection or mixture models. To facilitate this evaluation, transdimensional Markov chain Monte Carlo (MCMC) relies on sampling a discrete indexing variable to estimate the posterior model probabilities. However, little attention has been paid to the precision of these estimates. If only few switches occur between the models in the transdimensional MCMC output, precision may be low and assessment based on the assumption of independent samples misleading. Here, we propose a new method to estimate the precision based on the observed transition matrix of the model-indexing variable. Assuming a first order Markov model, the method samples from the posterior of the stationary distribution. This allows assessment of the uncertainty in the estimated posterior model probabilities, model ranks, and Bayes factors. Moreover, the method provides an estimate for the effective sample size of the MCMC output. In two model-selection examples, we show that the proposed approach provides a good assessment of the uncertainty associated with the estimated posterior model probabilities.

Tables

Table 1. Table 1: Logistic regression data set by Healy ( 1988 ) .

Condition (A)	Antitoxin (B)	Death	Survival
More Severe	Yes	15	6
	No	22	4
Less Severe	Yes	5	15
	No	7	5

Table 2. Table 2: Estimated posterior model probabilities in percent.

	Kuo and Mallick (1998)				Carlin and Chib (1995)
Model	$\text{Mean}\left(\widehat{𝝅}\right)$	$\text{SD}\left(\widehat{𝝅}\right)$	${\overline{\text{SD}}}_{\text{iid}}$	${\overline{\text{SD}}}_{\text{Markov}}$	$\text{Mean}\left(\widehat{𝝅}\right)$	$\text{SD}\left(\widehat{𝝅}\right)$	${\overline{\text{SD}}}_{\text{iid}}$	${\overline{\text{SD}}}_{\text{Markov}}$
1	0.51	0.24	0.07	0.16	0.57	0.35	0.06	0.39
A	49.28	1.38	0.50	1.22	48.55	7.14	0.49	6.92
B	1.14	0.44	0.10	0.26	1.26	0.63	0.10	0.73
A+B	43.85	1.25	0.50	1.10	43.61	7.41	0.49	7.19
AB	5.22	0.37	0.22	0.34	6.00	3.38	0.21	3.82

Table 3. Table 3: Models with the highest posterior probability for the 2 6 superscript 2 6 2^{6} contingency table.

		Posterior model probabilities $𝝅$				Rank $\tau$
#	Model	$\text{Mean}\left(\widehat{𝝅}\right)$	$\text{SD}\left(\widehat{𝝅}\right)$	${\overline{\text{SD}}}_{\text{iid}}$	${\overline{\text{SD}}}_{\text{Markov}}$	$\text{Mean}\left(\tau \right)$	$\text{SD}\left(\tau \right)$	$\overline{\text{SD}\left(\tau \right)}$	$\tau =\#$	$\overline{P\left(\tau =\#\right)}$	$\tau \le 10$	$\overline{P\left(\tau \le 10\right)}$
1	CE	18.78	1.34	0.13	1.02	1.00	0.00	0.03	1.00	1.00	1.00	1.00
2	BE	11.92	0.94	0.11	0.84	2.00	0.00	0.04	1.00	1.00	1.00	1.00
3	BE + CE	7.12	1.11	0.09	0.43	3.34	0.61	0.37	.72	.78	1.00	1.00
4	BF + CE	6.57	1.20	0.08	0.52	3.94	0.84	0.42	.71	.75	1.00	1.00
5	BE + BF	4.20	0.85	0.07	0.41	5.42	1.59	0.21	.92	.93	.96	.99
6	CE + EF	2.77	0.50	0.06	0.33	6.80	1.71	0.58	.62	.65	.94	1.00
7	BE + BF + CE	2.53	0.60	0.05	0.24	8.24	5.64	0.54	.58	.66	.92	1.00
8	CE + ADE	1.88	0.30	0.05	0.25	8.72	1.35	0.80	.47	.56	.95	.95
9	BE + EF	1.76	0.38	0.04	0.26	9.43	3.21	0.88	.45	.54	.92	.93
10	BE + ADE	1.19	0.22	0.04	0.19	12.05	3.11	1.40	.32	.39	.39	.56

Equations

$${\mathbf{\pi }}^{\top }\mathbf{P}=1\cdot {\mathbf{\pi }}^{\top },$$

$${\mathbf{p}}_i\equiv (p_{i1},\dots ,p_{iI})\sim \mathcal{D}(ϵ,\dots ,ϵ).$$

$${\mathbf{p}}_i^{(r)}\sim \mathcal{D}(n_{i1}+ϵ,\dots ,n_{iI}+ϵ).$$

$$\tilde{\mathbf{\pi }}\mid \tilde{\mathbf{n}}\sim \mathcal{D}({\tilde{n}}_1,\dots ,{\tilde{n}}_I),$$

$${\hat{T}}_{\text{eff}}=\sum _{i=1}^I{\hat{\alpha }}_i-I^2ϵ.$$

$$y_{jl}$$

$$\log \left(\frac{p_{jl}}{1-p_{jl}}\right)$$

$$\log {\mu }_j=\varphi +{\mathbf{x}}_j^{\top }\mathbf{\beta },$$

$$L\left(\mathbf{\alpha }\right)=\prod _{r=1}^R\left[\frac{\Gamma \left({\sum }_i{\alpha }_i\right)}{{\prod }_i\Gamma ({\alpha }_i)}\prod _i{\left({\pi }_i^{(r)}\right)}^{{\alpha }_i-1}\right].$$

Full text

Keywords: reversible jump MCMC, product space MCMC, Bayesian model selection, posterior model probabilities, Bayes factor.

1 Introduction

Transdimensional Markov chain Monte Carlo (MCMC) methods provide an indispensable tool for the Bayesian analysis of models with varying dimensionality (Sisson, 2005). An important application is Bayesian model selection, where the aim is to estimate posterior model probabilities $p(\mathcal{M}_{i}\mid\bm{y})$ for a set of models $\mathcal{M}_{i}$, $i=1,\dots,I$ given the data $\bm{y}$ (Kass and Raftery, 1995). In order to ensure that the Markov chain converges to the correct stationary distribution, transdimensional MCMC methods such as reversible jump MCMC (Green, 1995) or the product space approach (Carlin and Chib, 1995) match the dimensionality of parameter spaces across different models (e.g., by adding parameters and link functions). Transdimensional MCMC methods have proven to be very useful for the analysis of many statistical models including capture-recapture models (Arnold et al., 2010), generalized linear models (Forster et al., 2012), factor models (Lopes and West, 2004), and mixture models (Frühwirth-Schnatter, 2001), and are widely used in substantive applications such as selection of phylogenetic trees (Opgen-Rhein et al., 2005), gravitational wave detection in physics (Karnesis, 2014), or cognitive models in psychology (Lodewyckx et al., 2011; Heck et al., 2017).

Crucially, transdimensional MCMC methods always include a discrete parameter $z$ with values in $1,\dots,I$ indexing the competing models. At iteration $t=1,\dots,T$, posterior samples are obtained for the indexing variable $z^{(t)}$ and the model parameters, which are usually continuous and differ in dimensionality (for a review, see Sisson, 2005). For instance, a Gibbs sampling scheme can be adopted (Barker and Link, 2013), in which the indexing variable $z$ and the continuous model parameters are updated in alternating order. Such a sampler switches between models depending on the current values of the continuous parameters, and then updates these parameters in light of the current model $\mathcal{M}_{i}$ conditionally on the value of $z^{(t)}=i$ (Barker and Link, 2013). Given convergence of the MCMC chain, the sequence $z^{(t)}$ follows a discrete stationary distribution with probabilities $\bm{\pi}=(\pi_{1},\dots,\pi_{I})^{\top}$. Due to the design of the sampler, these probabilities are identical to the posterior model probabilities of interest, $\pi_{i}=p(\mathcal{M}_{i}\mid\bm{y})$ and, given uniform model priors $p(\mathcal{M}_{i})=1/I$, also proportional to the marginal likelihoods $p(\bm{y}\mid\mathcal{M}_{i})$. Hence, transdimensional MCMC samplers can be used to directly estimate these posterior probabilities as the relative frequencies of samples $z^{(t)}$ falling into the $I$ categories, $\hat{\pi}_{i}=1/T\sum_{t}\mathbbm{I}(z^{(t)}=i)$, where $\mathbbm{I}$ is the indicator function. Due to the ergodicity of the Markov chain, this estimator is ensured to be asymptotically unbiased (Green, 1995; Carlin and Chib, 1995).

Usually, dependencies due to MCMC sampling are taken into account for continuous parameters (Jones et al., 2006; Flegal and Gong, 2015; Doss et al., 2014). In contrast, however, the estimate $\bm{\hat{\pi}}=(\hat{\pi}_{1},\dots,\hat{\pi}_{I})^{\top}$ based on the sequence of discrete samples $z^{(t)}$ is usually reported without quantifying estimation uncertainty due to MCMC sampling. Often, the samples $z^{(t)}$ are correlated to a substantial, but unknown, degree because of infrequent switching between models. This is illustrated in Figure 1, which shows a sequence of independent and correlated samples $z^{(t)}$ in Panels A and B, respectively. Inference about the stationary distribution $\bm{\pi}$ is more reliable in the first case compared to the second case, in which the autocorrelation reduces the amount of information available about $\bm{\pi}$ (cf. Section 3). The standard error $\text{SE}(\hat{\pi}_{i})=\sqrt{\hat{\pi}_{i}(1-\hat{\pi}_{i})/T}$ that assumes independent sampling will obviously underestimate the true variability of the estimate $\bm{\hat{\pi}}$ if samples are correlated (Green, 1995; Sisson, 2005). To obtain a measure of precision, Green (1995) proposed running several independent MCMC chains $c=1,\dots,C$ and computing the standard deviation of the estimates $\bm{\hat{\pi}}^{(c)}$ across these independent replications. However, for complex models, this method might require a substantial amount of additional computing time for burn-in and adaption and thus can be infeasible in practice.

Assessing the precision of the estimate $\bm{\hat{\pi}}$, which depends on the autocorrelation of the sequence of discrete MCMC samples $z^{(t)}$, is of major importance. In case of model selection, it must be ensured that the estimated posterior probabilities $p(\mathcal{M}_{i}\mid\bm{y})$ are sufficiently precise for drawing substantive conclusions. This issue is especially important when estimating a ratio of marginal probabilities, that is, the Bayes factor $B_{ij}=p(\bm{y}\mid\mathcal{M}_{i})/p(\bm{y}\mid\mathcal{M}_{j})$ (Jeffreys, 1961). Moreover, it is often of interest to compute the effective sample size defined as the number of independent samples that would provide the same amount of information as the given MCMC output for estimating $\bm{\pi}$ with $\bm{\hat{\pi}}$. Besides providing an intuitive measure of precision, a minimum effective sample size can serve as a principled and theoretically justified stopping rule for MCMC sampling (Gong and Flegal, 2016). However, standard methods of estimating the effective sample size (e.g., computing the spectral density at zero; Plummer et al., 2006; Heidelberger and Welch, 1981) are tailored to continuous parameters. When applied to the model-indexing variable $z^{(t)}$ of a transdimensional MCMC method, these methods neglect the discreteness of $z^{(t)}$. Depending on the specific numerical labels used for the different models (e.g., $(1,2,3,4)$ vs. $(1,4,2,3)$), spectral decomposition can lead to widely varying and arbitrary estimates for the effective sample size (see Section 4).

In summary, transdimensional MCMC is a very important and popular method for Bayesian inference (Sisson, 2005). However, little attention has been paid to the analysis of the resulting MCMC output, which requires that one takes into account the autocorrelation as well as the discrete nature of the model-indexing variable. As a solution, we propose to fit a discrete, first order Markov model to the MCMC output $z^{(t)}$ to assess the precision of the estimated stationary distribution $\bm{\hat{\pi}}$. Whereas several diagnostics have previously been proposed to assess the convergence of transdimensional MCMC samplers (e.g., Brooks and Giudici, 2000; Castelloe and Zimmerman, 2002; Brooks et al., 2003a; Sisson and Fan, 2007), we are unaware of any methods that quantify the precision of the point estimate $\bm{\hat{\pi}}$.

2 Method

2.1 A Discrete Markov Model for Transdimensional MCMC Output

The proposed method approximates the output of a transdimensional MCMC method (i.e., the sampled iterations $z^{(t)}$) by a discrete Markov model $\mathcal{M}^{\text{Markov}}$ with transition matrix $\bm{P}$. This model explicitly accounts for autocorrelation, which in turn allows quantifying estimation uncertainty for the discrete stationary distribution $\bm{\pi}$. The entries of $\bm{P}$ are defined as the transition probabilities $p_{ij}=P(z^{(t+1)}=j\mid z^{(t)}=i)$ for all $i,j=1,\dots,I$, with rows summing to one, $\sum_{j=1}^{I}p_{ij}=1$. According to the discrete Markov model, the probability distribution of the indexing variable $z^{(t)}$ at iteration $t$ is given by multiplying the transposed initial distribution $\bm{\pi}_{0}^{\top}$ by the transition matrix $t$ times, $P(z^{(t)}=i)=[\bm{\pi}_{0}^{\top}\bm{P}^{t}]_{i}$. The proposed method estimates the transition matrix $\bm{P}$ as a free parameter based on the sufficient statistic $\bm{N}$, the matrix of frequencies $n_{ij}$ counting the observed transitions from $z^{(t)}=i$ to $z^{(t+1)}=j$ (Anderson and Goodman, 1957).

Due to the construction of the transdimensional MCMC sampler, the discrete indexing variable $z^{(t)}$ follows a stationary distribution with a constant probability vector $\bm{\pi}$ (i.e., the posterior model probabilities of interest). Hence, when modeling the sequence $z^{(t)}$ with the discrete Markov model $\mathcal{M}^{\text{Markov}}$, this implies that the transition matrix $\bm{P}$ must satisfy the condition for stationarity

[TABLE]

meaning that the probability vector $\bm{\pi}$ is the left eigenvector of the matrix $\bm{P}$ with eigenvalue one (with $\bm{\pi}$ normalized to sum to one; Anderson and Goodman, 1957). Based on the model $\mathcal{M}^{\text{Markov}}$, an estimator for $\bm{\pi}$ is thus obtained by computing the eigenvector of $\bm{P}$ with eigenvalue one (Barker and Link, 2013).

However, we are less interested in a new estimator $\bm{\hat{\pi}}$ of the stationary distribution but rather in the precision of this estimate. To quantify estimation uncertainty, we thus fit the model $\mathcal{M}^{\text{Markov}}$ with $\bm{P}$ as a free parameter in a Bayesian framework by drawing posterior samples $\bm{P}^{(r)}$ ($r=1,\dots,R$ ). Similar to a parametric bootstrap, this Bayesian sampling approach has the advantage that we can easily quantify estimation uncertainty (i.e., the dispersion of the posterior distribution of $\bm{P}$) by computing descriptive statistics of the samples $\bm{P}^{(r)}$ (e.g., the standard deviation or credibility intervals). Moreover, we can directly quantify the estimation uncertainty of derived quantities such as the posterior model probabilities, model ranks, or Bayes factors (see Section 2.2). In the following, it is important to distinguish between the posterior distribution of $\bm{P}$ given the sufficient statistic $\bm{N}$, which quantifies the uncertainty of $\bm{P}$ due to estimation error of the transdimensional MCMC method, and the posterior distribution of the models given the empirical data, which is represented by the constant vector of probabilities $\bm{\pi}$ for a specific data set.

Next, we define a prior distribution for the parameter $\bm{P}$ of the model $\mathcal{M}^{\text{Markov}}$. Given that the transition matrix $\bm{P}$ includes one probability vector $\bm{p}_{i}$ for each row $i$, we assume independent Dirichlet distributions with parameter $\epsilon\geq 0$ for each row,

[TABLE]

Conditional on the MCMC output $\bm{N}$, the estimation uncertainty of $\bm{P}$ is approximated by drawing $R$ posterior samples $\bm{P}^{(r)}$. Since the Dirichlet prior is conjugate to the multinomial distribution, independent samples $\bm{P}^{(r)}$ can efficiently be drawn from the Dirichlet distribution with parameters

[TABLE]

Based on these samples, the estimation uncertainty of the stationary probabilities $\bm{\pi}$ is assessed by computing the (normalized) eigenvector with eigenvalue one for each sample $\bm{P}^{(r)}$ (Eq. 1). Algorithm 1 provides an overview of the computational steps of the proposed method as pseudo-code.

With regard to the prior parameter $\epsilon$, small values should be chosen to reduce its influence on the estimation of $\bm{P}$. In principle, the improper prior $\epsilon=0$ can be used, which minimizes the impact of the prior on the estimated stationary distribution. This improper prior also ensures that the results do not hinge on the set of models that could possibly be sampled, but were never actually observed in the sequence $z^{(t)}$. For such unsampled models, the corresponding rows and columns of the observed transition matrix $\bm{N}$ are filled with zeros. With $\epsilon=0$, the relevant eigenvector of the transition matrix $\bm{P}\mid\bm{N}$ is thus identical to that of a reduced matrix $\bm{P^{*}}\mid\bm{N^{*}}$ that includes only the transitions for the subset of models sampled in $z^{(t)}$. However, in our simulations, this improper Dirichlet prior proved to be numerically unstable and resulted in more variable point estimates than the standard i.i.d. estimate or the proper prior discussed next.

Here, we use the weakly informative prior $\epsilon=1/I$ as a default, which has an impact equivalent to one observation for each row of the observed transition matrix $\bm{N}$. By putting a small weight on all values of the transition matrix $\bm{P}$, this prior serves as a regularization of the posterior (Alvares et al., 2018). However, in scenarios where the number of models exceeds the number of iterations of the transdimensional MCMC method (i.e., $I\gg T$), such a regularization assigns substantial probability weight to models that are never observed in $z^{(t)}$. To limit the effect of the prior, we thus set $\epsilon=1/I^{*}$ only for those $I^{*}$ models that were observed in $z^{(t)}$ and $\epsilon=0$ for the remaining models. Besides reducing the impact of the prior, this choice has the computational advantage that one can draw posterior samples and compute eigenvectors for the reduced matrix $\bm{P^{*}}\mid\bm{N^{*}}$ that includes only the sampled models. In the two examples in Sections 4 and 5, this prior has proved to be numerically robust and resulted in point estimates close to the standard i.i.d. estimates.

As a third alternative, the prior can be adapted to the structure of specific transdimensional MCMC implementations, which only implement switches to a small subset of the competing models. For instance, in variable selection, regression parameters are often added or removed one at a time, resulting in a birth-death process (Stephens, 2000). For these kinds of samplers, the Dirichlet parameters $\epsilon_{ij}$ can be set to zero selectively. However, such adjustments will be dependent on the chosen MCMC sampling scheme. The default choice of $\epsilon=1/I^{*}$ for sampled models and $\epsilon=0$ for nonsampled models provides a good compromise of being very general and numerically robust, while having a small effect on the posterior. However, in general, the choice of $\epsilon$ becomes less influential as the number of MCMC samples increases (especially if the row sums of $\bm{N}$ are large).

2.2 Estimation Uncertainty

Based on the posterior samples $\bm{P}^{(r)}$ of the transition matrix and the derived model probabilities $\bm{\pi}^{(r)}$, it is straightforward to estimate the stationary distribution by the posterior mean $\bm{\hat{\pi}}$ (alternatively, the median or mode may be used). More importantly, however, estimation uncertainty due to the transdimensional MCMC method can directly be assessed by plotting the estimated posterior densities for each $\pi_{i}$. To quantify the precision of the estimate $\bm{\hat{\pi}}$, one can report posterior standard deviations or credibility intervals for the components $\hat{\pi}_{i}$. These component-wise summary statistics are most useful if the number of models $I$ is relatively small.

An important advantage of drawing posterior samples $\bm{\pi}^{(r)}$ in a Bayesian framework (instead of using asymptotic approximations for the standard error of $\hat{\bm{\pi}}$) is that one can directly quantify estimation uncertainty for other quantities of interest. For very large numbers of sampled models, the assessment of estimation uncertainty can be focused on the subset of $k$ models with the highest posterior model probabilities. Within the sampling approach, estimation uncertainty for the $k$ best-performing models can easily be assessed by computing ranks for each of the posterior samples $\bm{\pi}^{(r)}$. Then, the variability of these model ranks across the $R$ samples can be summarized, for instance, by the percentage of identical rank orders for the $k$ best-performing models, or the percentages of how often each model is included within the subset of the $k$ best-performing models (i.e., has a rank smaller or equal to $k$).

In case of model selection, dispersion statistics such as the posterior standard deviation are also of interest with respect to the Bayes factor $B_{ij}$ (Kass and Raftery, 1995). To judge the estimation uncertainty for the Bayes factor, one can evaluate the corresponding posterior distribution by computing the derived quantities $B_{ij}^{(r)}=\pi_{i}^{(r)}/\pi_{j}^{(r)}$ (given uniform prior model probabilities). Precision can also be assessed for model-averaging contexts when comparing subsets of models against each other (e.g., regression models including a specific effect vs. those not including it). Given such disjoint sets of model indices $M_{s}\subset\{1,\dots,I\}$, the posterior probability for each subset of models is directly obtained by summing the posterior samples $\pi_{i}^{(r)}$ for all $i\in M_{s}$. Note that it is invalid to aggregate across model subsets or to drop sampled models before applying the proposed Markov approach because functions of discrete Markov chains (e.g., collapsing the $I$ original states into a subset of $S$ states) are not Markovian in general (Burke and Rosenblatt, 1958).

2.3 Effective Sample Size

Besides quantifying estimation uncertainty, the posterior samples $\bm{\pi}^{(r)}$ can be used to estimate the effective sample size for the transdimensional MCMC output. For this purpose, we consider the benchmark model $\mathcal{M}^{\text{iid}}$ under the ideal scenario of drawing independent samples $\tilde{z}^{(t)}$ from the categorical distribution with probabilities $\bm{\tilde{\pi}}$. For this model, we assume an improper Dirichlet prior on the stationary distribution, $\bm{\tilde{\pi}}\sim\mathcal{D}(0,\dots,0)$ (whereas the Markov model assumes a Dirichlet prior on the transition probabilities). Since this prior is conjugate to the multinomial distribution, the posterior for the stationary distribution $\bm{\tilde{\pi}}$ is given by

[TABLE]

conditional on the observed frequencies ${\tilde{n}_{i}}=\sum_{t=1}^{T}\mathbb{I}(\tilde{z}^{(t)}=i)$. Note that the transition frequencies are rendered irrelevant in this i.i.d. model, since there are no dependencies in the sampled iterations $\tilde{z}^{(t)}$.

Given the dependent samples $z^{(t)}$ of a transdimensional MCMC chain, we can now compare the empirical posterior distribution of $\bm{\pi}$ estimated with the model $\mathcal{M}^{\text{Markov}}$ against the theoretically expected posterior distribution of $\bm{\tilde{\pi}}$ under the hypothetical model $\mathcal{M}^{\text{iid}}$. Essentially, we match the latter distribution to the former to estimate the effective sample size as the total number of independent samples $T_{\text{iid}}=\sum_{i}\tilde{n}_{i}$ that would result in a similar dispersion as that estimated by the Markov model. To estimate the $\tilde{n}_{i}$, the i.i.d. posterior distribution in Eq. 4 is fitted to the posterior distribution of the Markov model by estimating the shape parameters $\alpha_{1},\dots,\alpha_{I}$ of a Dirichlet distribution given the sampled $\bm{\pi}^{(r)}$ (which can be achieved by an efficient maximum-likelihood algorithm by Minka, 2000, see Appendix). Next, a comparison of the estimated Dirichlet parameters $\widehat{\alpha}_{i}$ with the conjugate posterior in Eq. 4 yields $\widehat{\tilde{n}}=\widehat{\alpha}_{i}$, which implies that the dispersion of the posterior model probabilities $\bm{\pi}^{(r)}$ is equivalent to having observed $\widehat{T}_{\text{iid}}=\sum_{i}\hat{\alpha}_{i}$ independent samples. However, the samples $\bm{\pi}^{(r)}$ are not only informed by the samples $z^{(t)}$ of the transdimensional MCMC sampler, but also by the prior distribution of the Markov model, which is irrelevant for estimating the effective sample size. Hence, to estimate the effective sample size for the transdimensional MCMC sampler, it is necessary to subtract the prior sample size $I^{2}\epsilon$ of the Markov model (cf. Eq. 2), which reflects the relative weight of the prior, since the Dirichlet shape parameter $\epsilon$ occurs $I$ times in each row of the $I\times I$ transition matrix $\bm{P}$ (Alvares et al., 2018). Overall, it follows that the effective sample size under the assumption of independent sampling from a multinomial distribution is estimated as

[TABLE]

Note that it is necessary to replace $I$ by $I^{*}$ in Eq. 5 if the Markov model uses only those $I^{*}$ models that were actually sampled in $z^{(t)}$. Importantly, the estimate $\widehat{T}_{\text{eff}}$ takes the discreteness of the indexing variable $z$ into account and does not depend on specific (but arbitrary) numerical values of the model indices.

2.4 Remarks

The proposed method quantifies estimation uncertainty by fitting a discrete Markov model to transdimensional MCMC output. For this purpose, a simplifying assumption is made that is not guaranteed to hold. Whereas samples of the full model space $(z^{(t)},\bm{\theta}^{(t)})$ necessarily follow a Markov process by construction, this does not imply that the samples $z^{(t)}$ follow a Markov chain marginally (Brooks et al., 2003b; Lodewyckx et al., 2011). In practice, the iterations of the model-indexing variable $z^{(t)}$ might have higher-order dependencies since transition probabilities depend on the exact state of the MCMC sampler in each of the models’ parameter spaces. However, in Sections 4 and 5 we show in two empirical examples that the proposed simplification (i.e., fitting a Markov chain of order one) is sufficient to account for autocorrelations in the samples $z^{(t)}$ in practice. Moreover, the approximation by a first-order Markov chain provides a trade-off between ignoring dependencies completely (i.e., assuming i.i.d. samples) and accounting for any higher-order dependencies (which will likely increase the computational burden especially for large numbers of models). Note that it is a common practice to rely on simplifying assumptions for the analysis of simulation output; for instance, a standard approach of estimating the effective sample size for continuous parameters assumes that the output sequence can be modeled as a covariance stationary process with a smooth log spectrum (Heidelberger and Welch, 1981).

The proposed method of fitting a discrete Markov model is very general and can be applied irrespective of specific transdimensional MCMC implementations. Moreover, it requires only the sampled sequence $z^{\text{(t)}}$ of the discrete parameter or the matrix $\bm{N}$ with the observed frequency of transitions. If output from multiple independent chains $c=1,\dots,C$ is available, the transition frequency matrices $\bm{N}^{(1)},\dots,\bm{N}^{(C)}$ can simply be summed before applying the method. This follows directly from Bayesian updating of the stationary distribution $\bm{\pi}$. Essentially, each chain provides independent evidence for the posterior of the transition matrix $\bm{P}$, which is reflected by using the sums $\sum_{c}n_{ij}^{(c)}$ for the conjugate Dirichlet prior in Eq. 3. Note that this feature can be used to compare the efficiency of many short versus few long MCMC chains.

In the R package MCMCprecision (Heck et al., 2018), we provide an implementation that relies on the efficient computation of eigenvectors in the C++ library Armadillo (Sanderson and Curtin, 2016), accessible in R via the package RcppArmadillo (Eddelbuettel and Sanderson, 2014). On a notebook with an Intel® i7-7700HQ processing unit, drawing $R=5,000$ samples from the posterior distribution for 10 (100) sampled models requires approximately 150 milliseconds (28 seconds). Similar to any MCMC or bootstrap approach, the choice of the number of samples $R$ depends on the summary statistic used to quantify uncertainty. Whereas more samples are required to approximate the density distribution (e.g., $R\geq 5,000$), less samples (e.g., $R\approx 1,000$) are sufficient to approximate the SD of the estimated posterior model probabilities. Since the samples $\bm{\pi}^{(r)}$ are independently drawn and SDs are usually sufficient to quantify uncertainty, the choice $R=1,000$ is often sufficient in practice (however, for the simulations below, we use $R=5,000$).

3 Illustration: Effect of Autocorrelation

Before applying the proposed method to actual output of transdimensional MCMC samplers, we first illustrate its use in an idealized setting, where the interest is in approximating the posterior model probabilities $\bm{\pi}=(.85,.13,.02)^{\top}$ by drawing random samples $z^{(t)}$. To investigate the effect of independent versus dependent sampling, we generated sequences $z^{(t)}$ from the Markov model $\mathcal{M}^{\text{Markov}}$ with the stationary distribution $\bm{\pi}$. To induce autocorrelation, we defined a mixture process for each iteration $t$. With probability $\beta$, the discrete indexing variable was identical to the current model, $z_{t+1}=z_{t}$. In contrast, with probability $1-\beta$, the value $z_{t+1}$ was sampled from the given stationary distribution $\bm{\pi}$. Thereby, increasing values of $\beta$ resulted in a larger autocorrelation of the sequence $z^{(t)}$ as illustrated for $\beta=0$ and $\beta=0.8$ in Figure 1A and 1B, respectively.

For varying levels of $\beta=0,0.1,\dots,0.8$, we sampled 500 replications with $T=1,000$ iterations each, applied the proposed method (with $R=5,000$) and computed the precision of the estimate $\bm{\hat{\pi}}$. The main interest is in the posterior SD and in the coverage probability, defined as the probability that the data-generating values $\bm{\pi}$ are in the 90% credibility interval defined by the 5% and 95% quantiles. As a benchmark, we also computed these summary statistics under the (false) assumption that the samples $z^{(t)}$ were independently drawn by fitting the model $\mathcal{M}^{\text{iid}}$ with the Dirichlet posterior distribution in Eq. 4. Note that the latter uncertainty estimate is equivalent to the standard Monte Carlo error that assumes independent sampling.

Figure 2 shows the results of this simulation. In Figure 2A, the three panels correspond to the estimation uncertainty (i.e., the posterior SD) of the three posterior model probabilities $\bm{\pi}=(\pi_{1},\pi_{2},\pi_{3})^{\top}$. The estimated posterior SD of the Markov model indicated increasing uncertainty for larger values of $\beta$, thus taking the increasing autocorrelation into account. In contrast, the model $\mathcal{M}^{\text{iid}}$ assumes independence a priori, and thus, the posterior uncertainty was independent of $\beta$. As a result of this, the corresponding 90% credibility interval was less likely to include the data-generating value $\bm{\pi}$ for increasing values of $\beta$ (see Figure 2B), whereas the Markov model provided an accurate description of the estimation uncertainty for any degree of dependence.

4 Variable Selection in Logistic Regression

In the following, we apply the proposed method to the problem of selecting variables in a logistic regression, an example introduced by Dellaportas et al. (2000) to highlight the implementation of transdimensional MCMC in BUGS (see also Dellaportas et al., 2002; Ntzoufras, 2002). Table 1 shows the frequencies of deaths and survivals conditional on severity and whether patients received treatment (i.e., antitoxin medication; Healy, 1988). To emphasize the importance of considering estimation uncertainty for the posterior model probabilities, we compare the efficiency of two transdimensional MCMC approaches, which can both be implemented in JAGS (Plummer, 2003).

The full logistic regression model assumes a binomial distribution $\mathcal{B}$ of the survival frequencies $y_{jl}$ and a linear model on the logit-transformed survival probabilities $p_{jl}$,

[TABLE]

where $n_{jl}$ are the total number of patients in condition $jl$ and $\bm{\beta}$ the regression coefficient for the effect-coded variables $a_{j}$, $b_{l}$, and $(ab)_{jl}$. Variable selection is required to choose between $I=5$ models: the intercept-only model I, the three main effect models A, B, and A+B, and the model AB that includes the interaction. For comparability, we use the same priors as Dellaportas et al. (2000) and assume a centered Gaussian prior with variance $\sigma^{2}=8$ for each regression parameter, $\beta_{k}\sim\mathcal{N}(0,8)$. Moreover, the model probabilities were set to be uniform, $p(\mathcal{M}_{i})=1/5$. Note that nonuniform prior probabilities might be used to protect against multiple testing issues (i.e., Bayes multiplicity; Scott and Berger, 2010).

One of the two implemented transdimensional MCMC approaches uses unconditional priors (Kuo and Mallick, 1998, KM98) and includes indicator variables $\gamma_{ik}\in\{0,1\}$ for each regression coefficient $\beta_{k}$ in model $\mathcal{M}_{i}$. The parameter $\bm{\gamma}_{i}$ determines which regression coefficients are included by removing some of the additive terms of the linear model in Equation 7. Details about the full and conditional posterior distributions are provided by Dellaportas et al. (2000, p. 7).

As a second transdimensional MCMC approach, we implemented the method of Carlin and Chib (1995; CC95), which stacks up all model parameters into a new parameter $\bm{\theta}=(z,\bm{\beta}_{1},\dots,\bm{\beta}_{I})$, where $\bm{\beta}_{i}$ is the vector of regression parameters of model $\mathcal{M}_{i}$. Thereby, this approach samples a total of 12 regression parameters along with the indexing variable $z$. Note that the method of Carlin and Chib (1995) uses pseudo-priors $p(\bm{\beta}_{i}\mid\mathcal{M}_{j})$, $i\neq j$, that do not influence the statistical inference about $p(\bm{y}\mid\mathcal{M}_{i})$ and $p(\bm{\beta}_{i}\mid\bm{y},\mathcal{M}_{i})$. However, these pseudo-priors determine the conditional proposal probabilities $p(z\mid\bm{y},\bm{\beta}_{1},\dots,\bm{\beta}_{I})$ of switching between the models and thereby affect the efficiency of the MCMC chain. In substantive applications, these pseudo-priors should be chosen to match the posterior $p(\bm{\beta}_{i}\mid\mathcal{M}_{i})$ in order to ensure high probabilities of switching between the models (cf. Carlin and Chib, 1995; Barker and Link, 2013). Here, however, we did not optimize the sampling scheme and used $\bm{\beta}_{ik}\mid\mathcal{M}_{j}\sim\mathcal{N}(0,8)$ for the pseudo-priors to illustrate that our method can correctly detect the lower precision resulting from this suboptimal choice.

Figure 3 shows the estimated posterior distribution ($R=5,000$) of the posterior model probabilities using one Markov chain with 11,000 iterations (including 1,000 burn-in samples). The vertical black lines show the reference values for $\bm{\pi}$, approximated with very high accuracy by the KM98 approach using eight independent chains and one million samples each. As expected, the (incorrect) assumption that $z^{(t)}$ are sampled independently resulted in overconfidence in the point estimates of the CC95 approach. For all models, the corresponding posterior distributions missed the correct value and did not identify this estimation uncertainty. This shows the importance of assessing the dependency in the samples $z^{(t)}$ in order to judge the estimation uncertainty for the estimated posterior model probabilities. As a remedy, the proposed Markov approach resulted in a posterior distribution that covered the target values with high probability. Moreover, the novel estimation method revealed that the KM98 implementation had a higher precision compared to the CC95 approach, which was likely due to the (intentionally not optimized) choice of the pseudo-priors in the latter method. Hence, the Markov model allows comparison of the estimation uncertainty of different transdimensional MCMC methods for the model probabilities $\bm{\pi}$.

To test the validity of the proposed method more rigorously, we replicated the previous analysis 500 times. Thereby, the estimated precision can be compared against the actual sampling variability of the estimated model probabilities. For both transdimensional MCMC methods, Table 2 shows the mean estimated model probabilities in percent. Across replications, the point estimates (posterior means) from the Markov and the i.i.d. approach were very similar with a median absolute difference of $0.03\%$ and $0.31\%$ for the KM98 and CC95 implementations, respectively. To judge whether the estimated precision (i.e., the mean posterior standard deviations $\overline{\text{SD}}_{\text{iid}}$ and $\overline{\text{SD}}_{\text{Markov}}$) is valid, Table 2 shows the empirical SD of the estimates $\bm{\hat{\pi}}$ across replications. The results show that the assumption of independent samples $z^{(t)}$ leads to an overconfident assessment of the precision for the estimated model probabilities, $\overline{\text{SD}}_{\text{iid}}\ll\text{SD}(\bm{\hat{\pi}})$, which is especially severe for the less efficient CC95 implementation. In contrast, the Markov approach provided good estimates of the actual estimation uncertainty, $\overline{\text{SD}}_{\text{Markov}}\approx\text{SD}(\bm{\hat{\pi}})$. Moreover, for the MCMC method by Carlin and Chib (1995), the larger SDs indicate a smaller efficiency compared to the unconditional prior approach by Kuo and Mallick (1998). This theoretically expected result is due to the suboptimal choice of pseudo-priors. However, note that this difference in efficiency may be overlooked when merely computing relative proportions based on the sampled indexing variable $z^{(t)}$ (i.e., when implicitly assuming independent samples).

The higher efficiency of the KM98 approach becomes even clearer when assessing the median of the estimated effective sample size, which was $2,043$ for the KM98 approach compared to only $65$ for the CC95 method. As discussed above, commonly used estimators of effective sample size for continuous parameters (e.g., Plummer et al., 2006) should not be applied to the discrete model-indexing variable $z$ because they depend on the arbitrary numerical labels used for the models. If such methods are applied nevertheless, the resulting estimate for the effective sample size cannot be interpreted because it is not invariant under permutations of the arbitrary model indices used for the discrete parameter $z$. To illustrate this, Figure 4 shows the distribution of the estimated effective sample size when applying the spectral decomposition available in the R package coda (Plummer et al., 2006) to all 120 permutations of the model indices $(1,\dots,5)$ for a fixed sequence $z^{(t)}$. Since this method incorrectly assumes that the discrete variable $z$ is continuous, the estimated effective sample size was not invariant, but varied considerably depending on the specific labeling of the models (gray histogram). In contrast, the proposed Markov approach resulted in a well-defined, invariant estimate $\widehat{T}_{\text{eff}}=1,921$ (vertical black line) by explicitly accounting for the discreteness of $z$.

Finally, we show that the posterior samples $\bm{\pi}^{(t)}$ of the model $\mathcal{M}^{\text{Markov}}$ can directly be used to assess the uncertainty of Bayes factor estimates. For instance, substantive applications could be interested in testing whether to include the interaction term of condition (A) and treatment (B) in a logistic regression model. Given the output of a single MCMC run with 10,000 samples, Figure 5 shows the resulting posterior distribution of the Bayes factor $B_{\text{A+B,AB}}$ in favor for the absence of an interaction. Similar to the posterior model probabilities, the i.i.d. approach resulted in overconfidence regarding the estimate and most of the probability mass excluded the correct value 8.51 (approximated with a precision of $\text{SD}=0.020$). In contrast, the Markov approach corrected for dependencies in the samples $z^{(t)}$ and included the correct value. The same pattern emerged across the 500 replications, that is, the mean estimated SD of the Bayes factor approximated the corresponding empirical SD of the Bayes factor estimates (KM98: 0.56 vs. 0.60; CC95: 74.7 vs. 114.3). When using transdimensional MCMC, Bayes factors cannot be expected to be reliably estimated if models are never or very infrequently sampled (e.g., Model 1 in Table 2). For instance, the Bayes factor $B_{\text{A},\text{B}}\approx 43.8$ was estimated very imprecisely even in the KM98 approach (mean SD = 13.0; empirical SD = 24.3). To obtain more precise Bayes factor estimates in the presence of infrequently sampled models, it is recommended to rerun the transdimensional MCMC chain including only the two relevant models of interest (Lodewyckx et al., 2011).

5 Log-Linear Models for a $2^{6}$ Contingency Table

The application of the proposed method is also feasible in realistic scenarios with hundreds of sampled models. To illustrate this, we reanalyzed the $2^{6}$ complete contingency table by Edwards and Havránek (1985), which includes six risk factors for coronary heart disease (i.e., smoking, strenuous mental work, strenuous physical work, systolic blood pressure, ratio of $\alpha$ and $\beta$ lipoproteins, and family anamnesis of coronary heart disease). We are interested in finding the most parsimonious log-linear model, which accounts for the cell frequencies $y_{j}$ of cell $j$ ($j=1,\dots,2^{6}$) by assuming a Poisson distribution with mean $\mu_{j}$ and

[TABLE]

where $\phi$ is the intercept, $\bm{\beta}$ the vector of regression parameters, and $\bm{x}_{j}^{\top}$ the (transposed) design vector, which selects the elements of $\bm{\beta}$ included for modeling cell $j$. Here, we consider the class of hierarchical log-linear models that only allow the inclusion of an interaction if the corresponding marginal effects and lower interaction terms are included in the model as well (e.g., Overstall and King, 2014b).

To select between all 7.8 million possible hierarchical log-linear models (Dellaportas and Forster, 1999), we used the reversible jump algorithm proposed by Forster et al. (2012), which is implemented in the R package conting (Overstall and King, 2014a). Assuming a unit information prior (Ntzoufras et al., 2003), we sampled 100,000 iterations, discarded 10,000 as burn-in, and applied the proposed Markov chain method by drawing $R=5,000$ samples for the posterior model probabilities of the $I^{*}$ sampled models. To assess whether the estimated uncertainty accurately quantifies sampling variability, we ran 200 replications initialized with randomly chosen models.

Across replications, 5,805 models were sampled (on average, 562.7 per replication). Table 3 shows the results for the 10 models with the highest posterior probabilities. All of these 10 models included the six main effects (A: smoking, B: strenuous mental work, C: strenuous physical work, D: systolic blood pressure, E: ratio of $\alpha$ and $\beta$ lipoproteins, F: family anamnesis of coronary heart disease) and the first-order interactions AC, AD, AE, BC, and DE, but differed with respect to including the remaining interactions. Despite the large number of iterations, the estimation uncertainty (i.e., the posterior SD) of the posterior model probabilities was relatively large, indicating that the samples $z^{(t)}$ were autocorrelated to a substantial degree. This is also reflected by the effective sample size, which was estimated to be $\widehat{T}_{\text{eff}}=4,259$ on average ($\text{SD}=181$), approximately $5\%$ of the number of iterations after burn-in.

Table 3 also shows means and standard deviations of the sampled model rank $\tau$ for the models with the highest posterior probability, indicating that estimation uncertainty (i.e., the posterior SD) increased for models with smaller posterior probabilities. Moreover, the proportion of replications is shown for which the sampled rank $\tau$ was identical to the model index ($\tau=\#$) and smaller than or equal to 10 ($\tau\leq 10$). According to these proportions, exact ranks were estimated precisely only for the two best models, whereas the set of the 10 models with highest posterior probabilities was relatively stable across posterior samples (with the exception of model 10). Importantly, the Markov approach provided mean estimated probabilities $\overline{P(\tau=\#)}$ and $\overline{P(\tau\leq 10)}$ that matched the corresponding empirical proportions across replications.

Note that these results regarding estimation uncertainty are in line with our expectations — if models have small posterior probabilities, they are also sampled infrequently, which in turn results in estimation uncertainty. To quantify this variability, the proposed Markov chain approach provides an estimate for the achieved precision to assess the quality of the results and to find an appropriate stopping rule for MCMC sampling.

6 Conclusion

We proposed a novel approach for estimating the precision of transdimensional MCMC output. Essentially, a first order Markov model is fitted to the observed model-indexing variable $z^{(t)}$ to quantify estimation uncertainty of the corresponding stationary distribution. We showed that the method accounts for autocorrelation in a given sequence $z^{(t)}$ and provides a good assessment of estimation uncertainty. Importantly, the method does not require output of multiple independent MCMC chains and thus reduces the computational costs for adaption and burn-in. Besides being useful for transdimensional MCMC output, the method provides an estimate of the precision and effective sample size of discrete parameters in MCMC samplers in general. Thereby, researchers can easily decide whether the obtained precision is sufficiently high for substantive applications of interest.

Acknowledgments

Daniel W. Heck was supported by the research training group Statistical Modeling in Psychology (GRK 2277), funded by the German Research Foundation (DFG).

7 Appendix: Estimating the Shape Parameters of a Dirichlet Distribution

In the following, we outline the fixed-point algorithm proposed by Minka (2000) to estimate the vector of shape parameters $\bm{\alpha}=(\alpha_{1},\dots,\alpha_{I})^{\top}$ of a Dirichlet distribution. Given a set of $R$ probability vectors $\bm{\pi}^{(r)}$ (in the proposed method, these are the derived samples of the posterior model probabilities), the likelihood function of the shape parameters $\bm{\alpha}$ is

[TABLE]

To maximize this likelihood function, Minka (2000) developed an efficient fixed-point algorithm and proved its convergence to the unique maximum likelihood estimate $\hat{\bm{\alpha}}$. The computational steps are outlined in Algorithm 2. At its core, the current estimates $\alpha_{i}$ are updated in line 8 by using the digamma function $\Psi$ and its inverse $\Psi^{-1}$. As remarked by Minka (2000), the algorithm converges very fast even for a large number of shape parameters $I$ (e.g., 80 milliseconds on an Intel® i7-7700HQ for $I=1,000$).