ϕ-Divergence in Contingency Table Analysis

The ϕ-divergence association models for two-way contingency tables is a family of models that includes the association and correlation models as special cases. We present this family of models, discussing its features and demonstrating the role of ϕ-divergence in building this family. The most parsimonious member of this family, the model of ϕ-scaled uniform local association, is considered in detail. It is implemented and representative examples are commented on.


Introduction
Contingency tables and their analysis are of special importance for various diverse fields, like medical sciences, psychology, education, demography and social sciences. In these fields, categorical variables and their cross-classification play a predominant role, since characteristics of interest are often categorical, nominal or most frequently ordinal. For example, diagnostic ratings, strength of opinions or preferences, educational attainments and socioeconomic characteristics are expressed in ordinal scales. The origins of contingency table analysis (CTA) lie back in 1900 with the well-known contributions by Pearson and Yule, while, for the history of CTA before 1900, we refer to the interesting paper by Stigler [1]. The interest lies mainly in identifying and describing structures of underlying association in terms of appropriate models or measures. Divergence measures have been employed in the CTA mainly for hypothesis testing (model fit) and estimation, leading to general families of test statistics and estimators. The family of φ-divergence test statistics contain the classical likelihood ratio and Pearson test statistics as special cases while the maximum likelihood estimators (MLEs) belong to the family of the minimum φ-divergence estimators (MφEs). Families of φ-divergence based test statistics as well as MφEs for various standard models in CTA have a long history. However, their consideration and discussion is out of our scope. For a detailed overview and related references, we refer to the comprehensive book of Pardo [2]. For log-linear models, see also [3].
Here, we aim at highlighting a different structural role of φ-divergence in contingency tables modelling, namely that of linking phenomenological different models, forming thus a family of models and providing a basis for their comparison, understanding and unified treatment. Through this approach, new insight is gained for the standard association and correlation models (see [4,5]) while further alternatives are considered. We restrict our discussion on two-dimensional contingency tables, but the models and approaches discussed are directly extendable to tables of higher dimension.
The organization of the paper is as follows. Preliminaries on log-linear models, divergence measures, association and correlation models for two-way tables are provided in Section 2. In the sequel, the general family of φ-divergence based association models (AMs), which includes the classical association and correlation models as special cases, is reviewed in Section 3. The most parsimonious φ-divergence based association model, that of φ-scaled uniform local association, and its role in conditional testing of independence is considered and discussed in Section 4. For this family of models, the effect of the specific φ-function used is illustrated by analysing representative examples in Section 5. Some final comments are provided in Section 6.

Preliminaries
Consider an I × J contingency table n = (n ij ) with rows and columns classification variables X and Y, respectively, where n ij is the observed frequency in cell (i, j). The total sample size n = ∑ i,j n ij is fixed and the random table N is multinomial distributed N ∼ M(n, π), with probability table be the table of expected cell frequencies. Since the mapping (n, π) → m is one-to-one on ∆ I J , m ∈ {m = (m ij ) : m ij > 0, ∑ i,j m ij = n} and models (hypotheses) for π can equivalently be expressed in terms of m. Furthermore, let π r = (π 1+ , . . . , π I+ ) T and π c = (π +1 , . . . , π +J ) T be the row and column marginal probabilities vectors, respectively, and p = (p ij ) the table of sample proportions with p ij = n ij /n.
The classical independence hypothesis for the classification variables X and Y (π = π r π T c = π I ) corresponds to the log-linear model of independence (I), defined in terms of expected cell frequencies as where λ X i and λ Y j are the i-th row and j-th column main effects, respectively, while λ is the intercept. If the independence model is rejected, the interaction between X and Y is significant and the only alternative in the standard log-linear models set-up is the saturated model which imposes no structure on π. Identifiability constraints need to be imposed on the main effect and interaction parameters of these models, like λ X 1 = λ Y 1 = 0 and λ XY 1j = λ XY i1 = 0, for all i, j. An important generalized measure for measuring the divergence between two probability distributions is the φ-divergence. Let π = (π ij ), q = (q ij ) ∈ ∆ I J be two discrete finite bivariate probability distributions. Then, the φ-divergence between q and π (or Csiszar's measure of information in q about π), is given by where φ is a real-valued strictly convex function on [0, ∞) with φ(1) = φ (1) = 0, 0φ(0/0) = 0, 0φ(y/0) = lim x→∞ φ(x)/x (cf. [2]). Setting φ(x) = x log x, (3) is reduced to the Kullback-Leibler (KL) divergence For λ → −1 and λ → −0, (5) converges to the I KL (π, q) and I KL (q, π), respectively, while λ = 1 corresponds to Pearson's divergence.

Association Models
In case of ordinal classification variables, the association models (AMs) impose a special structure on the underlying association and thus provide non-saturated models of dependence. AMs are based on scores µ = (µ 1 , . . . , µ I ) and ν = (ν 1 , . . . , ν J ) assigned to the rows and columns of the ordinal classification variables, respectively. They are defined by the expression where the row and column scores are standardized subject to weights w 1 = (w 11 , . . . , w 1I ) T and w 2 = (w 21 , . . . , w 2J ) T , respectively, i.e., it holds Usually, the uniform (w 1 = 1 I , w 2 = 1 J , where 1 k is a k × 1 vector of 1s) or the marginal weights (w 1 = π r , w 2 = π c ) are used.
If µ and ν are both known and ordered, then (6) has just one parameter more than independence, parameter ζ, and is known as the Linear-by-Linear (LL) AM. In case the vector µ is unknown, (6) is the Row effect (R) AM, while the Column effect (C) AM is defined analogously. Finally, when the row and the column scores are all unknown parameters to be estimated, (6) is the multiplicative Row-Column (RC) AM. Scores that are unknown need not necessarily to be ordered. Note that models LL, R and C are log-linear while the RC is not. The degrees of freedom (d f ) of these . The special LL model for which the row and column scores are equidistant for successive categories is known as the Uniform (U) AM.
In case the RC model is not of adequate fit, multiplicative row-column AMs of higher order can be considered. Such a model of M-th order is defined as Vectors µ m = (µ 1m , . . . , µ Im ) and ν m = (ν 1m , . . . , ν Jm ) are the corresponding row and column eigenvectors, m = 1, . . . , M, which are orthonormalized with respect to the weights w 1 and w 2 , i.e., the following constraints are satisfied: AMs have been mainly developed by Goodman (see [4,5] and references therein) and are thus often referred to as Goodman's AMs. For a detailed presentation of the association models, their inference, properties, interpretation, the role of the weights used and associated literature, we refer to the book of Kateri [7] (Chapter 6).
For ease in understanding but also for interpretation purposes, it is convenient to think in terms of the local associations of the table and define thus AMs through the local odds ratios (LORs) Recall that the (I − 1) × (J − 1) table of LORs θ = (θ ij ), jointly with the marginal probabilities vectors π r and π c , specify uniquely the corresponding I × J probability table π. Thus, given π r and π c , a model on π can equivalently be expressed in terms of θ. Hence, model (8) is alternatively defined as for i = 1, . . . , I − 1 , j = 1, . . . , J − 1. In this set-up, the diffrences of successive row and column scores are constant for the U model, equal to since scores for successive categories are equidistant. Hence, the U model is equivalently defined as and is the model under which all local odds ratios are equal across the table, justifying its 'uniform association' characterization.

Correlation Models
A popular, mainly descriptive method for exploring the pattern of association in contingency tables is correspondence analysis (CA). The detailed discussion of CA is beyond our scope. For this, we refer to the book of Greenacre [8]. Correspondence analysis is a reparameterized version of the canonical correlation model of order M with 1 ≤ M ≤ M * . The row and column scores (x m = (x 1m , . . . , x Im ) and y m = (y 1m , . . . , y Jm ), m = 1, . . . , M) satisfy constraints (9) subject to the marginal weights. Usually, it is assumed M = 2 and the row and column scores (coordinates) are graphically displayed as points in two-dimensional plots.
The similarities between (8), expressed in terms of π in its multiplicative form, and (14) are obvious. Thus, motivated by the inferential approaches for AMs, MLEs have been considered for model (14), leading to the row-column correlation model of order M, while, for M = 1, special correlation models of U, R or C type have also been discussed by Goodman [4,5].

φ-Divergence Based Association Models
AMs and correlation models were developed competitively and opposed to each other (cf. [5]), until Gilula et al. [9] linked them in an inspiring manner under an information theoretical approach. They proved that, under certain (common) conditions, both of them are the closest model to independence. Their difference lies on the divergence used for measuring their closeness to independence. AMs are the closest in terms of the KL divergence and correlation models in terms of the Pearson's divergence. This result motivated subsequent research and led to the definition of general classes of dependence models by substituting the KL and Pearson divergences through generalized families of divergences.
Under the conditions of [9], namely for given marginal distributions (π r and π c ), given scores (µ and ν) and given their correlation ρ = corr(µ, ν), Kateri and Papaioannou [10] derived that the joint distribution π that is closest to independence in terms of the φ-divergence is of the form where F −1 is the inverse function of F(x) = φ (x). The scores µ and ν satisfy the constraints (7) with marginal weights. Under these constraints, it can easily be verified that ρ = corr(µ, ν) = ∑ i,j µ i ν j π ij . Additionally, the identifiability constraints are imposed on parameters α = (α 1 , . . . , α I )) and β = (β 1 , . . . , β J )). The parameter ζ is measuring association. It can be verified that (15), due to (7) with marginal weights and (16), leads to and that ζ = 0 if and only if the independence model holds (π = π I ). Furthermore, under model (15), the correlation ρ between the row and column scores is increasing in ζ and the φ-divergence measure I C φ (π, π I ), for given φ-function, is increasing in |ζ|.
Model (15), with known row and column scores, is the φ-divergence based extension of the LL model and is denoted by LL φ . If the scores are additionally equidistant for successive categories, (15) becomes the φ-divergence based U model, U φ , while the classes of models R φ , C φ and RC φ are defined analogously. The standard LL, R, C or RC models correspond to φ(x) = x log x. The analogue correlation models, defined by (14) for M = 1, are derived for φ(x) = (1 − x) 2 , setting µ = x 1 , ν = y 1 and ζ = ρ 1 . For the standard association and correlation models, (17) simplifies to ζ(π, µ, ν) = ∑ i,j π i+ π +j µ i ν j log(π ij ) (18) and respectively. For the power divergence (for φ(x) = x λ+1 −x λ(λ+1) , λ = −1, 0), model (15) becomes considered by Rom and Sarkar [11]. Expression (20) defines parametric classes of AMs, controlled by the parameter λ, which are denoted by LL λ , U λ , R λ , C λ or RC λ , according to the assumption made for the row and column scores. The RC(M) model, 1 ≤ M ≤ M * , is analogously generalized to RC φ (M), the class of φ-divergence AMs of order M, given by where the scores µ m and ν m satisfy restrictions (9) with marginal weights. The standard RC(M) model (8) is derived for φ(x) = x log x, while, for φ(x) = (1 − x) 2 , model (21) becomes the correlation model (14) with µ m = x m , ν m = y m and ζ m = ρ m , m = 1, . . . , M. Models (15) and (21) can alternatively be expressed as and is a scaled measure of local dependence, defined for i = 1, . . . , I − 1, j = 1, . . . , J − 1 as For φ(x) = x log x, (24) becomes the well-known log local odds ratio log(θ ij ), modelled in (13) and from now on denoted as θ Remark 1. Kateri and Papaioannou [10] studied properties of the class of φ-divergence based AMs. Additionally, they developed a test based on a φ-divergence statistic for testing the GOF of φ-divergence AMs and studied its efficiency. The choice of the φ-function in the test statistic is independent of the φ-function used for the model definition. Thus, it serves as a φ-divergence based GOF test for the traditional well-known association or correlation models.
In order to understand the role and nature of the scaled measures of local association (24), one can examine a simple 2 × 2 contingency table. In this case, (10) becomes the well-known odds ratio θ (= θ 11 ). Its φ-divergence scaled generalization, (24) for I = J = 2, has been explored by Espendiller and Kateri [12]. For these φ-scaled association measures for 2 × 2 tables, asymptotic tests of significance and confidence intervals (CIs) are constructed based on the fact that they are asymptotically normal distributed. An interesting feature of the family of φ-scaled association measures for 2 × 2 tables is that when the odds ratio θ cannot be finitely estimated (due to sampling zeros), some members of this family provide finite estimates, while, for a subset of them, their variance is also finite. Extensive evaluation studies verified earlier statements in the literature about the low coverage of the log odds ratio CIs when the association is very high (log θ > 4). In such cases, focusing on the power divergence scaled odds ratios θ (λ) , θ (1/3) for λ = 1/3 is to be preferred, since the corresponding CI is of better coverage when approaching the borders of the parameter space and is in general less conservative than the classical log odds ratio CI [12]. Here, the role of the scale effect on measuring the local dependence will be further clarified in the examples discussed in Section 5.

Remark 2.
The idea of viewing a model as a departure from a parsimonious reference model, with the property of being the closest to this reference model under certain conditions in terms of the KL divergence, can be adopted for other types of models as well, such as the quasi symmetry (QS) model and the logistic regression, lightening thus a different interpretational aspect of these models. Substitution of the KL divergence by the φ-divergence, leads further, in an analogous manner to AMs, to the derivation of generalized φ-divergence based classes of QS models [13], ordinal QS [14] and logistic regression models [15]. For example, in case of the QS, the role of the scaled measures (24) take the analogously defined scaled deviations from the model of complete symmetry.

Uniform Local Association
We shall focus on the case of uniform local association and the corresponding φ divergence scaled model U φ , defined by (15) or (22) with row and column scores satisfying (12). Model U φ is equivalently expressed as and forms a family of models. Compare to (13) for the U model defined in terms of the usual log odds ratio. The associated probability table under model U φ is uniquely specified by the one-to-one map (π r , π c , θ (φ) ) → π U φ , not given in a closed-form expression. The MLEs of the marginal probabilities are the corresponding marginal sample proportionŝ π r = p r = (p 1+ , . . . , p I+ ) T ,π c = p c = (p +1 , . . . , p +J ) T , while the MLE of θ (φ) ,θ (φ) , is not available in explicit form. For a given φ-function, model U φ belongs to the family of homogeneous linear predictor (HLP) models [16]. It is straightforward to verify that it satisfies the two conditions of Definition 3 in [16]. In practice, it can be fitted using Lang's mph R-package. The standard U model, denoted in the sequel as U 0 , has the equivalent HLP model expression where X and β are the model's design matrix and parameter vector (scalar for U 0 ), respectively. I−1,J−1 ) T . For more details on inference for model U 0 through HLP models, we refer to [7] (Sections 5.6 and 6.6.4). In Section 6.6.4 of [7], the approach is implemented in R for the example of Table 1 (see Section 5), while an R-function for constructing the design matrix C for two-way tables of any size is provided in the web appendix of the book. This approach is easily adjusted for model U φ , by replacing log(m v ) in (27) by F(m vs ), where m vs is the I J × 1 vector with entries m ij m i+ m +j /n = π ij π i+ π +j , expanded by rows.
Under U φ , all 2 × 2 subtables, formed by any successive rows and any successive columns, share the same φ-scaled local association θ (φ) . It is of practical interest to have estimators of this common local association, alternative to the MLEs, that are provided in explicit forms. One such estimator, based on the sample version of (17), is given bỹ Another option is which is the mean of the sample φ-scaled local association measures (24). For the power-divergence based models U λ , estimators (28) and (29) are denoted byθ (λ) andθ (λ) , respectively. For the U correlation model, derived for λ = 1 and denoted thus by U 1 , estimator (28) takes the form where r is the sample correlation between the row and column scores (compare to (19)).
Under the U φ model, θ (φ) is the single association parameter, measuring the strength of the local association that is uniform across the table. Furthermore, θ (φ) = 0(⇔ ζ = 0) if and only if the independence (I) model holds. Since model I is nested in U φ , the following test of independence, conditional on U φ , can be considered which is asymptotically X 2 1 distributed. This test is well-known for the standard U 0 model (see [17] and Section 6.3 in [7]).

Remark 3.
For model U 0 , Tomizawa [18] proposed a measure of divergence from uniform association, based on the KL-divergence, taking values in the interval [0, 1) and being equal to 0 if and only if U 0 holds. He constructed also asymptotic confidence interval for this measure, provided U 0 does not hold. Conde and Salicrú [19] extended his work by considering such a measure based on the φ-divergence and developed asymptotic inference for it, and also for the case that U 0 holds. Their approach and measures should not be confused with the approach followed here. They aim at detecting departures from U 0 (in favor of a non-uniform association structure) while here we focus on measuring the strength of uniform association, provided that U 0 (or U φ ) holds.
Tomizawa [18] as well as Conde and Salicrú [19] based the estimation of their measures on the following closed-form estimator of θ (0) under U 0 Remark 4. For square contingency tables with commensurable classification variables, analogous to the measure of departure from U 0 (see Remark 3), Tomizawa et al. [20] introduced a measure of departure from complete symmetry relying on the power divergence and Tomizawa [21] a measure of departure from marginal homogeneity. Kateri and Papaioannou [22] extended these measures to corresponding φ-divergence based measures and proposed further φ-divergence based measures of departure from the QS and triangular symmetry models. The work of Menéndez et al. [23,24] is also related.

Illustrations
We revisit a classical data set of Grizzle [25], provided in Table 2, which is adequately modelled by the U 0 model (see [18,19]). Our second data set in Table 1 corresponds to a study in [26] and provides strong evidence in favor of U 0 (see [7]). The maximum likelihood estimates of the expected cell frequencies under U 0 are also given in parentheses in Tables 1 and 2. The G 2 test statistics for the U λ models fitted on these data, for λ → 0 and λ = 1/3, 2/3, 1 are provided in Table 3, along with correspondingθ (λ) ,θ (λ) and θ (λ) values, i.e., the estimates of the θ (λ) s discussed in Section 4. For U 0 , the θ (0) * estimates (32) for Tables 1 and 2 are equal to 0.2890 and 0.7972, respectively.
We observe that all considered U λ models are of very similar (acceptable) fit for the first example while they differ enormously for the second one. For the data of Table 1, the fit of U 0 is impressive, that of U 1/3 acceptable while U 2/3 and U 1 are of very bad fit. A reverse situation appears for another data set, given in Table 4. In this case model, U 1 can be accepted while U 0 can not (see Table 5). The maximum likelihood estimates of the expected cell frequencies under U 1 are provided in Table 4 in parentheses. Table 3. Goodness of fit of the U λ models for the data of Tables 1 and 2 along with estimates of the common λ-scaled local association θ (λ) under U λ . Table 1 λ  Table 5. Goodness of fit of the U λ models and maximum likelihood estimates of the common λ-scaled local association θ (λ) under U λ for the data of Table 4. Observe (in Table 3) that the closed form estimates for the λ-scaled local associations are close to the corresponding maximum likelihood estimates in case the assumed model is of adequate fit while they diverge for models of bad fit.

Conclusions
We revealed the role of φ-divergence in modelling association in two-way contingency tables and illustrated it for the special case of uniform association in ordinal contingency tables. Targeting at pointing out the potential of this modelling approach and the generated families of models, we avoided presenting technical details, properties and inferential results for these models, which can be found in the initial sources cited.
Crucial quantities are the θ (φ) ij s, the φ-scaled measures of local association. The generalized family of φ-divergence based AMs enriches the modelling options in CTA, since the pattern of underlying association structure in a table may be simplified and thus described by a more parsimonious model when considering a different scale. A crucial issue, as also pointed out by one of the reviewers, is how to decide on the scale. So far, such a decision is based on trials of various alternative options. A formal approach selecting the scale is missing. In case of a parametric family, like the one based on the power divergence, the problem can be tackled by considering λ as an unknown parameter and estimating it from the data. Such an approach has been followed in the logistic regression set-up by Kateri and Agresti [15].
It is important to realize that, due to the scale difference, θ (φ) ij s are not directly comparable for different φ-function (or λ-values in case of the power divergence). Thus, comparisons across different φs (or λs) are possible only in terms of the corresponding expected cell frequencies or a common measure of local association evaluated on them. AMs can also be considered for other types of generalized odds ratios, like, for example, the global odds ratios. The extension of such models through the φ-divergence and their study is the subject of a work in progress. Inference for closed form estimators of the common θ (φ) of the U φ model and comparisons among them is the content of a paper under preparation.
The conditional test of independence (31) can be based not only on U φ but on LL φ models as well. Another 1 d f test of independence for ordinal classification variables is the linear trend test of Mantel (see [27]). It considers the testing problem H 0 : ρ = 0 vs. H 1 : ρ = 0, where ρ is the correlation between ordered scores assigned to the categories of the classification variables of the table. It is thus applicable only when the underlying association exhibits a linear trend for the assigned scores. The test of Mantel uses the test statistic M 2 = (n − 1)r 2 , which is under H 0 asymptotically X 2 1 distributed. The way this test is linked to the above-mentioned conditional tests, in view also of (19), is interesting to be investigated further.
Throughout this paper, we assumed a multinomial sampling scheme. For the models considered here, the other two classical sampling schemes for contingency tables (independent Poisson and product multinomial) are inferentially equivalent. Furthermore, for ease of presentation, we restricted here to two-way tables. The proposed models extend straightforwardly to multi-way tables. For twoor higher-dimensional tables, the subset of models that are linear in their parameters (i.e., RC-and RC(M)-type terms are excluded) belong to the family of homogeneous linear predictor (HLP) models [16] and can thus be fitted using the R-package mph.