1 Introduction
The use of asymptotic tests relied on large-sample approximations to the sample dis- tributions of test statistics is unreliable when a contingency table is large and sparse. Exact tests evaluating statistical significance using p-values are preferred over asymp- totic tests based on large sample approximations, such as χ 2 approximations. However, computing exact p-values is harder compared with computing asymptotic p-values, since the former requires enumeration of all contingency tables conditionally on the set This research is supported by the Japan Society for the Promotion of Science (JSPS), Grant- in-Aid for Scientific Research (C), No 20500263. Masahiro Kuroda Department of Socio-Information, Okayama University of Science, 1-1 Ridaicho, Okayama 700- 000. JAPAN E-mail: kuroda@soci.ous.ac.jp Hiroki Hashiguchi Graduate School of Science and Engineering, Saitama University, 255 Shimo-Okubo, Sakura, Saitama 338-8570, JAPAN Shigakazu Nakagawa Department of Computer Science and Mathematics, Kurashiki University of Science and the Arts, 2640 Nishinoura, Tsurajima-cho, Kurashiki-shi, Okayama 712-8505, JAPAN of marginals. For a detailed review of exact inference for contingency tables, refer to Agresti [1][2]. When the exact p-value can not be calculated due to the infeasibility of completely enumerating all contingency tables, MCMC is a valuable tool for evaluating p-values. Diaconis and Sturmfels [10] proposed an algebraic algorithm for constructing MCMC for sampling from conditional distributions. The algorithm generates random tables by the use of Markov bases computed using Gröbner bases. Further, MCMC with a Markov basis generates an irreducible Markov chain. Dobra [11] presented explicit formulae that identify Markov bases for decomposable graphical models and showed an algorithm for generating Markov bases dynamically. Forster et al. [13] used Gibbs sampling to perform exact tests for goodness of f t of the all-two-way interaction model for a 2 d contingency table. Caff o and Booth [4] suggested a different MCMC approach that computes p-values by making use of importance sampling with rounded normal deviates. Although Gibbs sampling and the MCMC method using importance sampling enable the computation of exact p-values in various log-linear models for multi-way tables, the Markov chains generated may be reducible. Our aim is to generate a Markov basis characterized by the Markov properties of Lauritzen [17] associated with a given conditional independence model for a four- way contingency table and develop a MCMC method for generating random tables by using the Markov basis. In Section 2, we describe conditional independence models for four-way contingency tables and show the Markov properties that form conditional independence structures for these models. In Section 3, we define the Markov bases for conditional independence models. In Section 4, we present an algorithm for generating a random table by employing the Markov basis. In Section 5, we give MCMC with the Markov basis for generating random tables of a Markov chain. The estimates of exact p-values can be calculated from random tables generated by MCMC. Numerical experiments in Section 6 examine the performance of the proposed MCMC method in comparison with χ 2 approximation. In Section 7, we present our concluding remarks.
2 Conditional independence models for a four-way contingency table
Consider an I × J × K × L contingency table n = { n abcd } with non-negative integer cell entries formed from categorical variables A , B , C and D . To denote conditional independence models for n , we use notations pertaining to their minimal sufficient statistics. For example, [ ABC ][ BCD ] denotes the model with conditional independence between A and D given ( B , C ). Using the notation of Dawid [9], its relationship can be expressed as A ⊥ D ( B , D ). The minimal sufficient statistics for [ ABC ][ BCD ] are marginals n ABC = { n abc+} and n BCD = { n+bcd}, where
For a four-way contingency table, there exist four conditional independence models, as described in Darroch et al. [8]. Table 1 shows these conditional independence models and their interpretations. In Table 2, the conditional independence relationships for these models in the form of Markov properties are presented. Figure 1 illustrates the graphical representations for the conditional independence models in Table 1.
Lauritzen [17] and Sundberg [18] showed that the exact conditional distributions of n for M 1 to M 4 have hypergeometric distributions conditionally on the sets of marginal constraints: Then, the MCMC method draws random tables from the above exact conditional distributions and is applied to estimate exact p -values from the generated random tables.
3 Markov bases for conditional independence models
Let T(M) denote the set of all tables having the marginal totals of sufficient statistics of model M. To generate table n ′ ∈ T ( M ), we use the data swapping technique of Dalenius and Reiss [7] such that cell entries are moved from one cell to the other, while the fixed marginals of sufficient statistics of M are left unchanged. A data swap associated with n ′ ∈ T ( M ) is an array f = { fabcd } with fabcd ∈ {0, ±1, ±2, . . } for all a, b, c, d. Definition 1 A move f = { abcd } for model M is a data swap that preserves the marginals of sufficient statistics of M . Then, we have n + f ∈ T ( M ) if and only if n abcd + f abcd ≥ 0 for all a, b, c, d. From the marginal constraints for M 1 to M 4 , we have
Definition 2  A Markov basis M is a finite collection of moves that preserve the fixed marginals of the sufficient statistics of model M . For any two tables n,n ′ ∈ T ( M ), there exists a sequence of moves f(1), . . ,f(R) E M such that
for 1 ≤ R ≤ R . Diaconis and Sturmfels [10] showed that there exists a Markov basis for T ( M ), which allows the construction of a Markov chain on T(M). A move f = { fabcd } in which two entries are equal to 1, two entries are equal to −1 and the remaining entries are 0 is called a primitive move . Diaconis and Sturmfels [10] described a Markov basis that is the set of primitive moves for the independence model for a two-way contingency table. Dobra [11] also showed that the set of primitive moves is a Markov basis for the class of decomposable graphical models.
4 Algorithm for generating random tables conditionally on the fixed
marginals
Consider an I × J × K contingency table formed from variables X, Y and Z. Let n = { n xyz} be the table of counts with non-negative integer cell entries. We assume that X and Z are conditionally independent given Y , i.e., X ⊥ Z | Y . The fixed marginals for the conditional independence model [ XY ][ Y Z ] are n xy = { n xy +} and n z = { n +y z }. By employing Theorem 3.1 of Diaconis and Sturmfels [10], the Markov basis for the model can be defined by the set of primitive moves with the following marginal constraints: The algorithm for generating a random table n ′ ∈ T ([ XY ][ Y Z ]) is described below. Algorithm: Random table generation for [XY][YZ] Step 1: Select an index yo in {1, . . , J}. Step 2:  Choose indices x 1 , x 2 , z 1 , z 2 such that 1 ≤ x 1 < x 2 ≤ I and 1 ≤ z 1 < z 2 ≤ K . Step 3:  Define move f = { fxyz } by
Step 4:  Generate n' from
n' = n + εf,
where ε is 1 or −1 with probability 1/2 each. All conditional independence structures for M 1 to M 4 can be transformed into the form X ⊥ Z|Y , as described below. Thus, Random table generation for [XY][Y Z] provides an essential algorithm for producing random tables in T(M 1 ) to T(M 4 ) and enables the generation of random tables for each model. The set {f(r)} in a Markov basis for each of M 1 to M 4 is given by the set of primitive moves conditionally on the marginal constraints specified by Equations (5) to (14). In the following subsections, we show how the Random table generation algorithm is applied to produce a move f in a Markov basis for each of M 1 to M 4 . In Appendix, we describe the procedure for selecting indices in Steps 1 and 2 of the Random table generation algorithm.
4.1 Model M 1 : [ABC][BCD]
We identify pair ( B,C ) with new variable Y that has the categories defined by all possible combinations of values of B and C , and set X = A and Z = D . By denoting M 1 as [XY][YZ], we can generate a random table in T(M 1 ) by applying the Random table generation algorithm.
4.2 Model M2: [AB][BCD]
Similar to model M 1 above, we define new variable Y with K L combinations of cat- egories of C and D . We denote M 2 as [XY][YZ] by setting X = A and Y = B and apply the Random table generation algorithm to produce a random table in T ( M 2 ).
4.3 Model M3: [AB][BC][CD]
As indicated in Table 2, there exist two Markov properties associated with M 3: (i) M 31 : A ⊥ (C, D)|B and (ii) M 32 : (A, B) ⊥ D|C . We set
Then, we uniformly generate a move with probability proportional to the number of moves associated with M 31 and M 32. Weights w 31 and w 32 represent the number of moves for M 31 and M 32 , respectively, and are calculated by Thus, we have p 31 = w 31/(w 31 + w 32) and p 32 = 1 − p 31. We select M31 or M32 with probability p3 = (p31,p32) and apply the Random table generation algorithm to generate a random table in T ( M 3 ).
4.4 Model M 4 : [AB][BC][BD]
For M4, we have three Markov properties: (i) M41 : A ⊥ C |(B,D), (ii) M42 : A ⊥ D|( B, C ) and (iii) M43 : C ⊥ D|( A, B ). We set
To generate a move uniformly, we select one model from M 41 , M 42 and M 43 with probability p4 = (p41,p42,p43). Weights for these models are calculated by
Thus, we have p 41 = w 41/(w 41 + w 42 + w 43), p 42 = w 42/(w 41 + w 42 + w 43) and P43 = 1 — (p41 + p42). We generate a move in a Markov basis for the model chosen by probability p4 and obtain a random table in T ( M 4 ) by using the Random table generation algorithm.
5 MCMC for the computation of p-values
In the two sections below, we present the Metropolis-Hastings (M-H) algorithm of Hastings [16] for generating Markov chains, followed by our approach to computing the corresponding p-values. 5.1 Metropolis-Hastings algorithm for generating a Markov chain Diaconis and Sturmfels [10] showed the M-H algorithm for generating a Markov chain using a Markov basis and then proved that the Markov chain is an irreducible, aperiodic Markov chain with stationary distribution π. The M-H algorithm for generating random tables of a Markov chain on T ([ XY ][ Y Z ]) is presented below.
Algorithm: M-H algorithm for [XY][YZ]
Step 1: Initialize the iteration counter r = 1 and set n(0) as the initial contingency table. Step 2: Generate candidate table n ′ = n(r − 1) + εf by using the Random table gener- ation algorithm. 6 Step 3: If all cell entries of n ′ are non-negative integers, accept n ′ as the next table n(r) with probability
otherwise retain n(r-1) and n(r) = n(r-1) Step 4: Increment counter r and return to Step 2. We apply the M-H algorithm presented above to generate {n ( r )} of a Markov chain on each of T ( M 1 ) to T ( M 4 ). When n ′ is drawn from hypergeometric distributions (1) to (4), the ratio π(n ′ )/π(n ′ − of Equation (15) involves only four cell counts associated with indices yo, x1, x2, z1, z 2 selected in Steps 1 and 2 of the Random table generation algorithm. For example, this algorithm generates table n ′ ∈ T(M 1 ) in which the four cell counts in n(r−1) are modified as Then, the ratio π(ń)/π(n
5.2 Computation of p-values
As the measure of the goodness of f t of model M , we use the Pearson chi-squared statistic where m = { mabcd } is the set of f fitted values under M . It follows that the p-value for M can be computed from
where I{ } denotes an indicator function. However, the computation of the exact p- value in Equation (16) is not feasible when complete enumeration of contingency tables in T ( M ) is very large or more complicated. The MCMC method is applicable to cir- cumvent this enumeration problem. Let {n(r)}0<r<(S+R) be random tables in T(M) generated by the M-H algorithm described above. After discarding the first S tables as burn-in, the MCMC p -value can be obtained by

6 Numerical experiments
We provide two numerical examples: the first is that exact results obtained by complete enumeration are available for comparison, and the second is that MCMC is the only feasible method for exact inference. The computation is performed by using the R language. Example 1 Table 3 is a study of nonmetastatic osteosarcoma by Goorin et al. [15]. The response shown in the table indicates whether the subject achieved a three-year disease-free interval. For M 1 and M 2 , the exact p-values are calculated from complete enumeration of T( M 1 ) and T( M 2 ) using MIMWIN of Edwards [12]. When exact p -values are unknown, we use the Monte Carlo estimation method from the importance sampling (IS) algorithm of Booth and Butler [3] instead of exact computations. The Monte Carlo p-values for M 3 and M 4 are obtained from 1,000,000 random tables by using the R package exactLoglinTest of Caffo [5], an implementation of the IS algorithm. We compare the p -values from the χ 2 approximation and the proposed MCMC approach with the exact or Monte Carlo p-values. The MCMC method described in Section 4 generates 1,000,000 random tables with 10,000 tables as burn-in. To evaluate the accuracy of a MCMC p-value, we compute its standard error using the batch means method of Geyer [14]. When dividing random tables {n(r)}s+1<r<S+R,R2 into R1 batches each of size R2, the MCMC p-value for the k-th batch is calculated as The variance estimate of ˆ ρ is obtained from and the batch means estimate of the MCMC standard error is calculated from √ Var/R 1 . In our experiments, we set R 1 = 100 and R 2 = 10,000. The values in parenthesis are the standard errors of the MCMC p -values. Results are shown in Table 4. As s shown in the table, MCMC p-values are in good agreement with the exact and Monte Carlo p-values, while the asymptotic p-values using the chi-squared distribution differ greatly. Figure 2 shows histograms of values 80 {χ 2 (n ( r ), m)} of the Pearson chi-squared statistics for M 1 to M 4 . For models other than M 3 , we see substantial discrepancies between the asymptotic and MCMC estimated exact distributions. They indicate that the asymptotic χ 2 approximation is unreliable for sparse data. Example 2  Table 5 shows abortion opinion data from Christensen [6]. Observations are classified by four factors: Race ( A ), Sex ( B ), Opinion ( C ) and Age ( D ). For C , three different opinions are possible: a "Yes" answer supports legalized abortion; a "No" answer opposes legalized abortion; an "Undecided" answer is undecided. Because the contingency table contains many small expected values and the sample size is large, the computation of exact p -values obtained by enumeration is not feasible. Table 6 presents the asymptotic and MCMC p -values for M1 to M4 . The MCMC p -values for each model are obtained from 1,000,000 simulated tables after a burn-in of 10,000 tables. Standard error values of p are shown in parentheses. The small MCMC p -value for M 4 may lack accuracy, because a MCMC estimate from 1,000,000 tables guarantees at most three-digit precision.
7 Concluding remarks
In this paper, we proposed the use of MCMC with a Markov basis for estimating exact p-values for the conditional independence models for a four-way contingency table. A Markov basis is required to generate random tables with fixed marginals specified by a conditional independence model. We defined the Markov basis characterized by the Markov properties associated with a given conditional independence model. Then, we developed an algorithm for generating primitive moves in the Markov basis and presented a MCMC method for producing a Markov chain by using the Markov basis. Results from our first experiment demonstrated that the proposed MCMC method finds the estimates of exact p -values in close agreement with exact p -values. Our second experiment illustrated that the MCMC method is feasible for exact inference when the asymptotic χ 2 approximation is unreliable and the exact computation of p -values by complete enumeration is not feasible. In the future, we intend to develop a MCMC method for decomposable graphi- cal models in multi-way contingency tables by extending our MCMC method with a Markov basis characterized by Markov properties. Acknowledgements The authors would like to thank the editor and a referee for their valu- able comments and helpful suggestions.
Appendix
We describe the procedure for selecting indices in Steps 1 and 2 of the Random table generation algorithm for each model. Algorithm: Random table generation for M 1 Step 1: Select yo = (bo,co) in {1,...,J} × {1,...,K}. Step 2: Choose x 1 , x 2 in {1,..., I} and z 1 , z 2 in {1,..., L}. Algorithm: Random table generation for M2 Step 1: Select yo in {1,..., J}. Step 2: Choose x 1 , x 2 in { 1 , . . , I } and z 1 = (c 1 , d 1 ), z 2 = (c 2 , d 2 ) in { 1 , . . , K } x {1, . . . , L}. Algorithm: Random table generation for M3 Determine M31 or M32 with probability p3. If M31 is selected, set (X,Y,Z) = (A,B,(C,D)), else (X,Y,Z) = ((A,B),C,D). Step 1: Select yo: – for M31, yo in {1, . . . , J}, − for M32, yo in {1,...,K}. Step 2: Choose x 1 , x 2 , z 1 , z 2: − for M 31, x 1, x 2 in {1, . . . , I} and z 1 = (c 1,d 1 ), z 2 = (c 2,d 2 ) in {1, . . . ,K} x {1, . . . , L}, − for M32, x1 = (α1, b1), x2 = (α2, b2) in {1, . . . , I} × {1, . . , J} and z1, z2 in {1, . . . , L}. Algorithm: Random table generation for M4 Select one model among M41, M42 and M43 with probability p4. Set
( A, ( B, D ), C ), when choosing M41, ( A , ( B , C ), D ), when choosing M42, (C, (A, B), D), when choosing M43. Step 1: Select yo: — for M41, yo = (bo,do) in {1,...,J} × {1,...,L}, − for M42, yo = (b0,c0) in {1,...,J} × {1,...,K}, − for M43, yo = (ao,bo) in {1,...,I} × {1,...,J}. Step 2: Choose x 1 , x 2 , z 1 , z 2: − for M 41 , x 1 , x 2 in {1, . . . , I} and z 1 , z 2 in {1, . . . , K}, − for M 42 , x 1 , x 2 in {1, . . . , I} and z 1 , z 2 in {1, . . . , L} , − for M43, x 1 , x 2 in {1, . . . , K } and z 1 , z 2 in {1, . . . , L } .