Stochastic approximation Monte Carlo importance sampling for approximating exact conditional probabilities

Abstract

Importance sampling and Markov chain Monte Carlo methods have been used in exact inference for contingency tables for a long time, however, their performances are not always very satisfactory. In this paper, we propose a stochastic approximation Monte Carlo importance sampling (SAMCIS) method for tackling this problem. SAMCIS is a combination of adaptive Markov chain Monte Carlo and importance sampling, which employs the stochastic approximation Monte Carlo algorithm (Liang et al., J. Am. Stat. Assoc., 102(477):305–320, 2007) to draw samples from an enlarged reference set with a known Markov basis. Compared to the existing importance sampling and Markov chain Monte Carlo methods, SAMCIS has a few advantages, such as fast convergence, ergodicity, and the ability to achieve a desired proportion of valid tables. The numerical results indicate that SAMCIS can outperform the existing importance sampling and Markov chain Monte Carlo methods: It can produce much more accurate estimates in much shorter CPU time than the existing methods, especially for the tables with high degrees of freedom.

Appendices

Appendix A: Markov basis

Let M denote a Poisson log-linear model, as described in Sect. 2, for a contingency table. Let T(M) denote the set of all contingency tables having the marginal totals of the sufficient statistics of model M. Let δ denote a data swapping move (Dalenius and Reiss 1982) of M, in which cell entries are moved from one cell to the other while the marginals of the sufficient statistics of M are left unchanged.

Definition A.1

A Markov basis $\mathcal{M}$ of the model M is a finite collection of data swapping moves that connect any two tables in T(M). In other words, for any two tables x,y∈T(M), there exists a sequence of data swapping moves $\boldsymbol{\delta }_{1}, \ldots, \boldsymbol{\delta}_{s} \in\mathcal{M}$ such that

$$\boldsymbol{y}-\boldsymbol{x}=\sum_{k=1}^s \boldsymbol{\delta}_k, \quad\mbox{and} \quad\boldsymbol{x}+\sum _{k=1}^{s'} \boldsymbol { \delta}_k \in T(M), $$

for 1≤s′≤s.

The Markov basis can be computed by finding a Gröbner basis of a well-specified polynomial ideal, which has been implemented in the package 4ti2 (4ti2 team 2006). Although a Markov basis always exists, computing the Markov basis can be very time-consuming or even practically impossible for certain types of models. For example, De Loera and Onn (2006) showed that for the three-way contingency tables with fixed two margins, the Markov basis can contain entries of arbitrarily large size.

Appendix B: Contingency tables of examples

Table 6 Data for the sexual fun study example (source: Hout et al. 1987)

Table 7 Data for the health care study example (source: Long et al. 2010). Notation: A=Level of care with 1 = independent, 2 = assisted living, and 3 = continuous care; B = religious preference with 1 = Protestant, 2 = Catholic, 3 = Jewish, and 4 = Other; and C = marital status with 1 = single, 2 = married, 3 = widowed, 4 = divorced, and 5 = separated

Table 8 Data for the abortion opinion study example (source: Christensen 1997)

Table 9 Whittaker’s survey data on women’s economic activity (source: Whittaker 1990)