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

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.


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(M i | x) for a set of models M i , i = 1, . . ., I given the data x (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, Hayakawa, and Yip, 2010), generalized linear models (Forster, Gill, and Overstall, 2012), factor models (Lopes and West, 2004), and mixtures models (Frühwirth-Schnatter, 2001), and are widely used in substantive applications such as selection of phylogenetic trees (Opgen-Rhein, Fahrmeir, and Strimmer, 2005), gravitational wave detection in physics (Karnesis, 2014), or cognitive models in psychology (Lodewyckx et al., 2011).
Crucially, transdimensional MCMC methods always include a discrete parameter Z with values in 1, . . ., I indexing the competing models.At iteration t = 1, . . ., T , posterior samples are obtained for the indicator 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 indicator 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 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 stationary distribution π = (π 1 , . . ., π I ).Due to the design of the sampler, this distribution is identical to the posterior model probabilities of interest, π i = p(M i | x) and, given uniform model priors p(M i ) = 1/I, also proportional to the marginal likelihoods p(x | 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, πi = 1/T t I(z (t) = i), where 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 (Cowles and Carlin, 1996).In contrast, however, the estimate π = (π 1 , . . ., πI ) based on the sequence of discrete samples z (t) is usually reported without quantifying estimation uncertainty.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 π should be more reliable in the first case compared to the second case, in which the autocorrelation reduces the amount of information available about π.Moreover, the standard error SE(π i ) = πi (1 − πi )/T that assumes independent sampling will obviously underestimate the true variability of the estimate π (Green, 1995;Sisson, 2005).To obtain a measure of precision, Green (1995) proposed running several independent MCMC chains c = 1, . . ., C and computing the standard error of the estimate π(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 π, which depends on the autocorrelation of the sequence of discrete samples z (t) , is of major importance.In case of model selection, it must be ensured that the estimated posterior probabilities p(M i | x) are sufficiently precise for drawing substantive conclusions.This issue is especially important when estimating a ratio of marginal probabilities, that is, the Bayes factor  1961).More generally, it is often of interest to compute the effective sample size, that is, the number of independent samples that would provide the same amount of information as the given MCMC output.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).Software for convergence diagnostics such as the R package coda (Plummer et al., 2006) estimate the effective sample size for continuous parameters by computing the spectral density at zero (Heidelberger and Welch, 1981).
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 indicator variable.As a solution, we propose to fit a discrete Markov model to the MCMC output z (t) to assess the precision of the estimated stationary distribution π.Note that our approach differs from diagnostics previously proposed to monitor convergence of transdimensional MCMC methods.These diagnostics usually compare the variance of z (t) and other continuous parameters across and within chains and models (Brooks and Giudici, 2000;Castelloe and Zimmerman, 2002), similar to the widely used potential scale reduction factor (Gelman and Rubin, 1992).
Other methods estimate the convergence rate of transdimensional MCMC chains, rely on Kolmogorov-Smirnov or χ 2 tests (Brooks, Giudici, and Philippe, 2003), or use distancebased diagnostics (Sisson and Fan, 2007).However, none of these methods estimates the precision of the point estimate π.

Posterior Distribution of Model Probabilities
To estimate the fixed, but unknown stationary distribution π (e.g., the posterior model probabilities), we explicitly model the sequence of observed samples as a random variable Z (t) (t = 1, . . ., T ) by assuming that it has emerged from a discrete, homogeneous Markov chain M Markov .For this purpose, we only consider the marginal distribution of the discrete indicator variable and define the transition matrix P with the switching probabilities p ij = P (Z (t+1) = j | Z (t) = i) for all i, j = 1, . . ., I. Accordingly, P is a probability matrix with non-negative entries and rows summing to one, i.e.I j=1 p ij = 1.For this Markov model, the resulting probability distribution at iteration t is given by multiplying the transposed initial distribution π 0 by the transition matrix t times, P (Z (t) = i) = [π 0 P t ] i .Within the model M Markov , the transition matrix P is a free parameter that is to be estimated based on the sampled sequence z (t) .
Due to the construction of the transdimensional MCMC sampler, Z (t) has π as a stationary distribution.Therefore, the transition matrix P in the Markov model M Markov satisfies the condition which implies that the probability vector π is the left eigenvector (normalized to sum to one; Anderson and Goodman, 1957) of the matrix P with eigenvalue one.This eigenvector exists if the Markov chain is finite and irreducible (Stewart, 2009, Ch. 9).Hence, given an estimate for P , we can directly compute π as the corresponding eigenvector.
However, we are less interested in a point estimate π of the stationary distribution but rather in the precision of this estimate.For this purpose, the sequence z (t) is summarized by the sufficient statistic N (Anderson and Goodman, 1957), the observed matrix of transition frequencies, where n ij is the number of switches from z (t) = i to z (t+1) = j.Conditional on the MCMC output, the posterior distribution of P can be approximated by drawing r = 1, . . ., R samples according to By computing the (normalized) eigenvectors with eigenvalue one (Eq.1), posterior samples of the stationary distribution π are directly obtained.As a prior for the transition matrix P , we assume independent Dirichlet distributions with parameter ≥ 0 for each of the rows, Therefore, independent posterior samples P (r) can efficiently be drawn from the conjugate posterior distribution, With regard to the prior parameter , small values should be chosen to reduce its influence on the estimation of P .Note that this choice becomes less influential as the number of MCMC samples increases (especially if the row sums of N are large).Here, we use the prior = 1/I, which has an impact equivalent to one observation for each row of the observed transition matrix N .This prior has proved to be numerically robust in the two examples in Sections 4 and 5, and resulted in point estimates close to the default i.i.d. estimates.
Alternatively, the improper prior = 0 can be used, which minimizes the impact of the prior on the estimated stationary distribution.Note that this 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 N are filled with zeros.With = 0, the relevant eigenvector of the posterior P | N is thus identical to that of a reduced matrix P | N that includes only the transitions for the subset of models sampled in z (t) .Moreover, if the set of competing models is large, N is likely to be a sparse matrix because many transitions will never occur.When choosing = 0, P (r) will also be a sparse matrix when sampling from the conjugate Dirichlet distribution D(n i1 , . . ., n iI ), which facilitates an efficient computation of the eigenvectors π (r) .However, in our simulations, this improper Dirichlet prior proved to be numerically unstable and resulted in more variable point estimates than the proper prior = 1/I.
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 ij can be set to zero selectively.However, such adjustments will be dependent on the chosen MCMC sampling scheme.Therefore, we propose the weakly informative prior = 1/I as a default, which provides a good compromise of being very general and numerically robust, while having a small effect on the posterior.

Precision
Based on the posterior samples π (r) , it is straightforward to estimate the stationary distribution by the posterior mean π (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 π i .To quantify the precision of the estimate π, one can report posterior standard deviations or credibility intervals for the components πi .Note that these component-wise summary statistics are most useful if the number of models I is relatively small.
For very large numbers of sampled models, the assessment of estimation uncertainty can be focused on the subset of models with the highest posterior model probabilities.
Besides summarizing the estimated posterior model probabilities, estimation uncertainty for the k best-performing models can also be assessed by computing ranks for each of the posterior samples π (r) .Then, the variability of these model rankings 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 (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 ⊂ {1, . . ., I}, the posterior probability for each subset of models is directly obtained by summing the posterior samples π (r) i for all i ∈ M s .Note that it is invalid to aggregate across model subsets 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).

Effective Sample Size
Besides quantifying estimation uncertainty, the posterior samples π (r) can be used to compute the effective sample size for the transdimensional MCMC output.For this purpose, we consider the benchmark model M iid under the ideal scenario of drawing independent samples z(t) from the categorical distribution with probabilities π (which is equivalent to sampling from a multinomial distribution).For this model, we also assume a Dirichlet prior, but this time directly on the stationary distribution, π ∼ D(γ, . . ., γ) with a fixed parameter γ ≥ 0. Since the prior is conjugate, the posterior for the estimated distribution is given by based on the observed frequencies ñi = I(z (t) = i) = j n ij .By considering only the proportion that each model is visited, the transition frequencies are rendered irrelevant in this i.i.d.approach, thereby ignoring possible dependencies in the sample sequence z (t) .
Given the output of a transdimensional MCMC chain, we can now compare the empirical posterior of π derived from the model M Markov against the theoretically expected posterior π under the model M iid to estimate the sample size T = i ñi .For this purpose, a Dirichlet distribution with parameters α 1 , . . ., α I is fitted to the samples π (r) , which can be accomplished by an efficient fixed-point iteration scheme (Minka, 2000).Second, a comparison of the estimated Dirichlet parameters with the conjugate posterior in Eq. 5 yields α i ≈ ñi + γ.Therefore, after subtracting the prior sample size I 2 of the I × I transition matrix P (Eq.3), the effective sample size under the assumption of independent sampling from a multinomial distribution can be estimated as To minimize the impact of the prior, one can assume the improper prior γ = 0 as a default.
Importantly, this estimate takes the discreteness of the indicator variable Z into account and does not change under permutations of arbitrary, numerical values of the model indices.

Remarks
If output from multiple independent chains c = 1, . . ., C is available, the transition frequency matrices N (1) , . . ., N (C) can simply be summed before applying the method.This follows directly from Bayesian updating of the stationary distribution π.Essentially, each chain provides independent evidence for the posterior, which is reflected by using the sums ij for the conjugate Dirichlet prior in Eq. 4. Note that this feature can be used to compare the efficiency of many short versus few long MCMC chains.
The proposed method bears a resemblance to the convergence diagnostic of Brooks, Giudici, and Philippe (2003).To assess the convergence rate of a transdimensional MCMC chain, Brooks, Giudici, and Philippe (2003) computed the second largest eigenvalue of the maximum-likelihood estimate of the transition matrix, pij = n ij / k n ik .Since P is a probability matrix, this eigenvalue must be smaller than one and provides a measure of the dependency of the samples z (t) .
Note that the simplifying assumptions underlying our approach are not guaranteed to hold.Whereas samples of the full model space (z (t) , θ (t) ) necessarily follow a Markov process by construction, this does not imply that the samples z (t) follow a Markov chain marginally (Brooks, Giudici, and Roberts, 2003;Lodewyckx et al., 2011).Intuitively, this is due to the fact that the transition probabilities depend on the exact locations 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.
The proposed method can be applied irrespective of specific transdimensional MCMC implementations and requires only the sampled sequence z (t) of the discrete parameter or the matrix N with the observed frequency of transitions.In the R package MCMCprecision (Heck et al., 2017), 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 (Eddelbüttel and Sanderson, 2014).On a notebook with an Intel R i5-3320M processing unit, drawing R = 1, 000 samples from the posterior distribution for 10 (100) sampled models requires approximately 70 milliseconds (10 seconds).

Illustration: Effect of Autocorrelation
Before applying the proposed method to actual MCMC output, we first illustrate its use in an idealized setting.To investigate the effect of autocorrelation in the case of discrete parameters, we generate sequences z (t) of length T = 1, 000 from the Markov model M Markov for a given stationary distribution π = (.85,.13,.02) .To induce autocorrelation, we define a mixture process for each iteration t.With probability β, the discrete indicator variable will be identical to the current model, z t+1 = z t .In contrast, with probability 1 − β, the value z t+1 is sampled from the given stationary distribution π.Thereby, increasing values of β result in larger autocorrelation of the sequence z (t) .
For varying levels of β = 0, 0.1, . . ., 0.8, we sampled 500 replications, applied the proposed method and computed the precision of the estimate πMarkov .Besides the posterior SD, we were interested in the coverage probability of the data-generating value π being in the 90% credibility interval.As a benchmark, we also computed the corresponding posterior SD under the (false) assumption that the samples z (t) were independently drawn by fitting the model M iid with the Dirichlet prior parameter γ = 0.In contrast, the model M iid does assume independence a priori.Thereby, the posterior uncertainty is independent of β.As a result of this, the 90% credibility interval is less likely to include the data-generating value π as shown in Panel B, whereas the Markov model provides an accurate description of the estimation uncertainty.

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, Forster, and Ntzoufras (2000) to highlight the implementation of transdimensional MCMC in BUGS (see also Dellaportas Forster, and Ntzoufras, 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 Bernoulli distribution B of the survival frequencies x jl and a linear model on the logit-transformed survival probabilities p jl , where n jl are the total number of patients in condition jl and β 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, 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, Forster, and Ntzoufras (2000) and assume a centered Gaussian prior with variance σ 2 = 8 for each regression parameter, β k ∼ N (0, 8).Moreover, the model probabilities were set to be uniform, p(M i ) = 1/5.Note that efficiency can be increased by selecting prior probabilities that result in approximately uniform posterior probabilities Lodewyckx et al., 2011).Moreover, 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 γ ik ∈ {0, 1} for each regression coefficient β k in model M i .The parameter γ i determines which regression coefficients are included by removing some of the additive terms of the linear model in Equation 8. Details about the full and conditional posterior distributions are provided by Dellaportas, Forster, and Ntzoufras (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 θ = (z, β 1 , . . ., β I ), where β i is the vector of regression parameters of model M i .Thereby, this approach samples a total of 12 regression parameters along with the indicator variable z.
Note that the method of Carlin and Chib (1995) uses pseudo-priors p( that do not influence the statistical inference about p(x | M i ) and p(β i | x, M i ).However, these pseudo-priors determine the conditional proposal probabilities p(z | x, β 1 , . . ., β 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(β i | 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 use β ik | M i ∼ 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 of the posterior model probabilities using one Markov chain with 21,000 iterations (including 1,000 burn-in samples).The vertical black lines show the correct reference values, approximated by eight independent chains with one million samples each.As expected, the assumption that z (t) are sampled independently results in overconfidence in the point estimates of the CC95 approach.For all models, the corresponding posterior distributions miss the correct value and do not  identify this estimation uncertainty.In contrast, the proposed Markov approach results in a posterior distribution that covers the target values with sufficiently high probability.
Moreover, the novel estimation method reveals that the KM98 implementation has a higher precision compared to the (intentionally not optimized) CC95 approach.
To test the validity of the proposed method more rigorously, we replicated the previous analysis 100 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 from the Markov and the i.i.d.approach were similar with mean absolute differences smaller than 0.02% and 0.49% for the KM98 and CC95 implementations, respectively.To judge whether the estimated precision (i.e., the mean posterior standard deviations SD i and SD M ) is valid, Table 2 shows the empirical SD of the estimates πMarkov across replications.The results clearly show that the assumption of independent samples z (t) leads to an overconfident assessment of the precision for the estimated model probabilities, whereas the proposed Markov approach provides a good estimate of the actual estimation uncertainty.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 likely to be 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 indicator variable z (t) (i.e., when implicitly assuming independent samples).
The higher efficiency of the KM98 approach becomes even clearer when assessing the mean effective sample size across replications, which was estimated to be 4,514 compared to only 163 for the CC95 method.Note that commonly used estimators of effective sample size (e.g., Plummer et al., 2006) depend on the exact numerical labels of the modelindicator variable Z.To illustrate this, Figure 4 shows the estimated sample size for all 120 permutations of the indices (1, . . ., 5) for a fixed sequence z (t) , computed by the spectral decomposition available in the R package coda (Plummer et al., 2006).This estimate varied  (Plummer et al., 2006) for all permutations of the model indicator labels of the MCMC output z (t) (based on 20,000 samples of the method by Kuo and Mallick, 1998).
considerably depending on the arbitrary labeling of the models.In contrast, the proposed Markov approach results in a well-defined, invariant estimate by explicitly accounting for the discreteness of Z.
Finally, we show that the posterior samples π (t) of the model M 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 20,000 samples, Figure 5 shows the resulting posterior distribution of the Bayes factor B A+B,AB in favor for the absence of an interaction.Similar to the posterior model probabilities, the i.i.d.approach results in overconfidence regarding the estimate and most of the probability mass excludes the correct value 8.50 (approximated with a precision of SD = 0.014).In contrast, the Markov approach corrects for dependencies in the samples z (t) and includes the correct value.The same pattern emerged across 100 replications, that is, the mean estimated SD of the Bayes factor matched the corresponding empirical SD of the Bayes factor estimates (KM98: 0.40 vs. 0.47; CC95: 6.93 vs. 6.24).When using transdimensional MCMC, Bayes factors cannot be expected to be reliably estimated if models  and Mallick, 1998).
are never or very infrequently sampled (e.g., Model 1 in Table 2).For instance, the Bayes factor B A,B ≈ 44.4 was estimated very imprecisely even in the KM98 approach (mean SD = 7.5; empirical SD = 14.8).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 α and β 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, . . ., 2 6 ) by assuming a Poisson distribution with mean µ j and log µ j = φ + x j β, where φ is the intercept, β the vector of regression parameters, and x j the (transposed) design vector, which selects the elements of β 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, 2014a).
To select between all 7.8 million possible hierarchical log-linear models (Dellaportas and Forster, 1999), we use the reversible jump algorithm proposed by Forster, Gill, and Overstall (2012), which is implemented in the R package conting (Overstall and King, 2014b).Assuming a unit information prior (Ntzoufras, Dellaportas, and Forster, 2003), we sampled 100,000 iterations, discarded 10,000 as burn-in, and applied the proposed Markov chain method by drawing 2,000 samples for the posterior model probabilities.To assess whether the estimated uncertainty accurately quantifies sampling variability, we ran 100 replications initialized with randomly chosen models.
Table 3 shows the 10 models with the highest posterior probabilities.The relatively large posterior standard deviations of the estimated posterior model probabilities indicate that the samples z (t) are autocorrelated to a substantial degree, despite the large number of iterations.This is also reflected by the effective sample size, which was estimated to be T eff = 4, 484 on average (SD = 186), approximately 5% of the number of iterations after burn-in.
Table 3 also shows means and standard deviations of the sampled rank R 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 posterior samples is shown for which the sampled rank R was identical to the rank across all replications (R = #) and smaller than or equal to 10 (R ≤ 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 mean estimated probabilities P (R = #) and P (R ≤ 10) 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.

Conclusion
We proposed a novel approach for estimating the precision of transdimensional MCMC output.Essentially, a Markov model is fitted to the observed model-indicator variable z (t) to obtain posterior samples of the corresponding stationary distribution.We showed that the method corrects 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.

Figure 1 :
Figure 1: Illustration of the MCMC iterations of the model-index parameter z (t) .Whereas samples are independent multinomial samples in Panel A, samples are drawn from a Markov chain with positive autocorrelation in Panel B.

Figure 2 :
Figure 2: Estimation uncertainty for the stationary distribution π.(A) The Markovmethod (black dots) correctly indicates that estimation error of the posterior model probabilities increases as autocorrelation increases.When assuming i.i.d.samples (gray triangles), the estimated precision does not depend on the autocorrelation.(B) Proportion of 500 replications for which the 90% CI intervals include the data-generating stationary distribution π.

Figure 2
Figure2shows the result of this simulation.In Panel A, the estimation uncertainty increases for larger values of β, thereby taking the increasing autocorrelation into account.

Figure 3 :
Figure 3: Posterior distribution of the posterior model probabilities π in the logistic regression example based on the Markov and the i.i.d.model.Vertical black lines show the target values (CC95 = Carlin and Chib, 1995; KM98 = Kuo and Mallick, 1998).

Figure 4 :
Figure4: Effective sample size as estimated by the spectral density at zero(Plummer et al., 2006) for all permutations of the model indicator labels of the MCMC output z (t) (based on 20,000 samples of the method byKuo and Mallick, 1998).

Figure 5 :
Figure 5: Posterior distribution for the Bayes factor in favor of Model A+B vs. AB.The vertical black line shows the target value (CC95 = Carlin and Chib, 1995; KM98 = Kuo and Mallick, 1998).

Table 2 :
Estimated posterior model probabilities in percent.Posterior model probability estimates π are shown in percent.Mean(π) and SD(π) were computed across 100 replications.As a measure for the estimated precision, means of the posterior SD are shown (SD i assumes independent sampling; SD M assumes a Markov chain model). Note:

Table 3 :
Models with the highest posterior probability for the 2 6 contingency table.Posterior model probabilities π are shown in percent.Mean(π) and SD(π) were computed across 100 replications.All models include the six main effects, A: smoking, B: strenuous mental work, C: strenuous physical work, D: systolic blood pressure, E: ratio of α and β lipoproteins, F: family anamnesis of coronary heart disease, and the first-order interactions AC, AD, AE, BC, and DE.