On mathematical optimization for clustering categories in contingency tables

Many applications in data analysis study whether two categorical variables are independent using a function of the entries of their contingency table. Often, the categories of the variables, associated with the rows and columns of the table, are grouped, yielding a less granular representation of the categorical variables. The purpose of this is to attain reasonable sample sizes in the cells of the table and, more importantly, to incorporate expert knowledge on the allowable groupings. However, it is known that the conclusions on independence depend, in general, on the chosen granularity, as in the Simpson paradox. In this paper we propose a methodology to, for a given contingency table and a fixed granularity, find a clustered table with the highest χ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi ^2$$\end{document} statistic. Repeating this procedure for different values of the granularity, we can either identify an extreme grouping, namely the largest granularity for which the statistical dependence is still detected, or conclude that it does not exist and that the two variables are dependent regardless of the size of the clustered table. For this problem, we propose an assignment mathematical formulation and a set partitioning one. Our approach is flexible enough to include constraints on the desirable structure of the clusters, such as must-link or cannot-link constraints on the categories that can, or cannot, be merged together, and ensure reasonable sample sizes in the cells of the clustered table from which trustful statistical conclusions can be derived. We illustrate the usefulness of our methodology using a dataset of a medical study.


Introduction
Data science comprises a plethora of methods and strategies to extract useful knowledge from raw data. Among other disciplines, such as statistics, computer science or information technology, mathematical optimization plays a crucial role across its tasks (Bottou et al. 2018;Gambella et al. 2021;Olafsson et al. 2008). Much effort has gone into incorporating recent developments in optimization theory and software to tackle data science problems more effectively, such as in regression and classification (Bertsimas and King 2016;Bertsimas and Shioda 2007;Blanquero et al. 2020;Carrizosa and Romero Morales 2013;Carrizosa et al. 2021;Toriello and Vielma 2012), clustering strategies (Benati and García 2014;Hansen and Jaumard 1997;Hochbaum and Liu 2018;Park et al. 2000;Saglam et al. 2006), correspondence analysis (van de Velden et al. 2020), dimensionality reduction methods (Carrizosa and Guerrero 2014;Carrizosa et al. 2020;Cunningham and Ghahramani 2015), deep learning (Anderson et al. 2020;Fischetti and Jo 2018) or data visualization (Carrizosa et al. 2017a(Carrizosa et al. , 2018a(Carrizosa et al. , b, 2019. There are still many data science problems which do not take advantage of such advancements and ad-hoc strategies are still used, such as, the analysis of the independence of two categorical variables through a function of the entries of their contingency table. Let U and V be two categorical variables, which take on a finite number of values, u 1 , . . . , u r and v 1 , . . . , v c , respectively. Given a set of n entities for which these variables have been observed, a first summary of their distribution is provided by their contingency table, in which the frequency of the event (u i , v j ), o i j , is collected for i = 1, . . . , r and j = 1, . . . , c. Table 1 contains an example of a contingency table in which n = 98 observations (lowest right corner) are cross-classified according to variable U , which has categories u 1 and u 2 and variable V , which has three categories (v 1 , v 2 and v 3 ). Whereas the inner rectangle in the table contains the joint frequencies o i j , the last row (resp. column) contains the marginal frequencies o . j of V , j = 1, . . . , c (resp. o i. of U , i = 1, . . . , r ).
When the data is cross-classified as in Table 1, the statistical (in)dependence of two categorical variables is usually investigated using the classical χ 2 measure (Pearson 1900;Mirkin 2001), although different approaches exist in the literature (Goodman and  Kruskal 1979; Joe 1989). The χ 2 coefficient is an estimate of the deviation between the empirical probability distribution of the variables U and V and the probability distribution that we would have if the two variables were statistically independent, and it is given by where e i j = o i. o . j n . The minimum value of χ 2 in (1) is 0, which occurs if and only if the variables U and V are statistically independent. Therefore, the larger χ 2 , the stronger the evidence against independence. The χ 2 statistic in (1) approximates a Chi-squared distribution with (r −1)×(c−1) degrees of freedom. In Table 1, one gets χ 2 = 6.355, which provides evidence against the statistical independence of the two variables involved ( p−value= 0.04, which tests the null hypothesis of statistical independence) for a 5% significance level. These inferential properties can be derived as long as the observed joint frequencies o i j are large enough for all i = 1, . . . , r and j = 1, . . . , c. To ensure this, the categories of the variables, associated with the rows and columns of the table, are often grouped, yielding a less granular representation of the categorical variables. Clustering categories in rows and/or columns of a contingency table is also desirable to enhance interpretability and transparency (Baesens et al. 2003;Carrizosa et al. 2017bCarrizosa et al. , 2022Goodman and Flaxman 2017;Ustun and Rudin 2016), by easing the presentation of the table as well as the conclusions of the analysis from a statistical perspective. Furthermore, constrained clustering allows the analyst to incorporate knowledge about the problem under study and support meaningful decision making (Abin 2019;Śmieja and Wiercioch 2017). However, it is known that the conclusions on independence depend, in general, on the granularity chosen for each of the categorical variables. For instance, let us consider that variable V in the example above is encoded as v 1 = v 1 & v 3 and v 2 = v 2 . Thus, observations in v 1 and v 2 are now grouped together, yielding the contingency table in Table 2, for which χ 2 = 3.519. Thus, the clustered table has c = 2 columns instead of the c = 3 in the initial table and is the one yielding the largest χ 2 among all the tables of that granularity, namely with two columns and two rows. In this case, there is not a significant evidence at a 5% significance level to reject the statistical independence assumption between the variables U and V ( p−value= 0.06). Therefore, the Simpson's paradox (Blyth 1972) arises in this example since the less granular representation of the categorical variables in Table 2 supports a conclusion, namely statistical independence, different from that suggested by the variables before the grouping of categories, namely statistical dependence in Table 1 (Shmueli and Yahav 2017;Tsumoto 2009). Thus, we have identified a socalled extreme grouping, namely Table 1 has the largest granularity of the variables U and V for which statistical dependence between them is found. In this paper we propose a mathematical optimization model to cluster the rows/columns of a contingency table so that the χ 2 statistic is maximized for a fixed granularity of the variables. Solving this problem for different sizes allows us to either identify extreme groupings or conclude that they do not exist. Our approach is flexible enough to include constraints on the desirable structure of the clusters, such as must-link or cannot-link constraints on the categories that can, or cannot, be merged together, and ensure reasonable sample sizes in the cells of the clustered table from which trustful statistical conclusions can be derived. This constrained clustering approach allows us to incorporate background knowledge to support the analysis and extract meaningful conclusions (Abin 2019;Śmieja and Wiercioch 2017).
The problem of clustering the categories of a contingency table to find extreme groupings for a fixed granularity has not been studied as such in the literature. There are however related approaches that involve distributional assumptions or use adhoc heuristic procedures that are not flexible enough to include constraints on the clusters. Indeed, Greenacre (1988) proposed a greedy procedure, based on hierarchical clustering, which uses a χ 2 related distance function between row (resp. column) vectors. However, this approach does not guarantee that clustered tables with the highest dependence for a fixed granularity are necessarily found because only a reduced family of allowable groupings, namely a hierarchical structure, are considered. The classical k-means clustering algorithm has been also adapted to the particular case of contingency tables (Govaert 1995;Govaert and Nadif 2007) as well as geometrical approaches, such as the maximum-tangent-plane (Bock 2003), have been developed. Ciampi et al. (2005) propose using the coordinates obtained with correspondence analysis to find a clustering and Álvarez de Toledo et al. (2018) use a similarity measure between the categories to obtain a partition. In order to find homogeneous clusters in document-term matrices, Ailem et al. (2016) propose maximizing a graph modularity criterion and Labiod and Nadif (2011) a community detection one. Whereas these approaches are based on the optimization of a measure of association, some probabilistic approaches have been also studied. In this case, it is assumed that each element of the contingency table is generated according to a probability model, which is tried to be recovered from the data. In this context, Ailem et al. (2017a, b); Riverain and Nadif (2022) propose latent block models to identify a diagonal structure of homogeneous blocks in document-term matrices. Proceeding this way, blocks of zeroes are identified and clustered together, thus yielding joint frequencies in the clustered table which are (close to) zero. A unified framework about the optimization of measures of association and probabilistic approaches is studied by Govaert and Nadif (2010). The aforementioned approaches are unable to deal with the analysis of dependence between variables in sparse tables, namely tables for which some of the observed joint frequencies are equal or close to zero and thus statistical conclusions cannot be inferred. It is well known that the common practice of adding constants to small joint frequencies can disturb the possible statistical dependence structure underlying in a sparse table (Agresti and Yang 1987). If some categories were properly clustered together, this sparsity problem would be removed without damaging the underlying possible relationship between the variables. However, the existing methods cannot incorporate the corresponding constraints to ensure that this removal of sparsity is achieved.
The idea of clustering categories in contingency tables has been also applied to the discretization of continuous variables involved in supervised learning algorithms. To guide the search of the partition, a criterion which assesses the relationship between the intervals in which the continuous variable is split and the target values to be predicted is optimized. For instance, Kerber (1992) uses the χ 2 distance between adjacent intervals to merge them if they are similar enough according to a given threshold. Boulle (2004) proposes a greedy approach and uses the p−value associated to the χ 2 statistic of the clustered table to select the discretization. However, these methods fail when constraints have to be imposed to the discretization being sought, such as that each interval in the partition has to have a large enough number of observations or there are rules that have to be accomplished (e.g. minimum or maximum length of the intervals which form the discretization).
In this paper, we propose an assignment and a set partitioning mathematical optimization formulations to cluster rows and/or columns of contingency tables maximizing the χ 2 statistic in (1), as a measure of the strength of the dependence, for a fixed size of the clustered table. Solving this model for different sizes, we can decide whether the statistical dependence can be preserved with the chosen granularity of the variables. If this is the case, we reduce the size of the parameter and solve the maximization problem again. We do this until we find the extreme groupings, or conclude that they do not exist, namely for any size of the reduced table the dependence of the variables can be preserved. Our model can easily be enriched with constraints to incorporate user knowledge on the allowable groups of categories, or to successfully handle sparse tables. With the proposed formulations, even contingency tables as the ones in the numerical section can be tackled using off-the-shelf optimization solvers.
The remainder of the paper is structured as follows. Section 2 states the mathematical optimization model to cluster categories in a contingency table maximizing the χ 2 statistic and imposing structural properties in the clusters. An assignment and a set partitioning formulations for such model are presented in Sect. 3. Finally, Sect. 4 illustrates our methodology and Sect. 5 concludes the paper with some remarks and future research.

Problem definition
This section is devoted to presenting a mathematical optimization model to cluster the rows and/or columns of a given contingency table which maximizes the χ 2 statistic in (1) for a fixed granularity of the categorical variables whereas requirements on the clusters, that is conditions about allowable groups of categories or thresholds over the sample sizes in the cells of the clustered table, are also imposed.
Let T 0 be a contingency table representing the counts of outcomes of two categorical variables U and V , which both take a finite set of values (categories), u 1 , . . . , u r and v 1 , . . . , v c , respectively. Recall that given a sample of n entities, o i j denotes the joint observed frequency of the pair (u i , v j ), o i. the marginal frequency of u i and o . j the marginal frequency of v j , for i = 1, . . . , r and j = 1, . . . , c. In order to measure the strength of the association between the variables U and V , the χ 2 statistic as stated in (1) is used. Let χ 2 (T 0 ) be the value of (1) for data in table T 0 .
Given the contingency table T 0 , a clustered table T is obtained from it by merging the rows and/or columns of T 0 into a new set of categories (clusters). In other words, a set of k (k ≤ c) clusters of columns of T 0 , {ṽ 1 , . . . ,ṽ k }, is a partition of the set {v 1 , . . . , v c } into k groups such that, for l, l = 1, . . . , k: Similarly, row clusters can be also defined asũ 1 , . . . ,ũ s (s ≤ r ). The clustered contingency table T from T 0 has a less granular representation of its categorical variables U and V and has as joint frequenciesõ ml , which are obtained from the sum of the corresponding joint frequencies in T 0 , namelyõ ml = i: u i ∈ũ m j: v j ∈ṽ l o i j . In other words, a clustered table T accumulates the corresponding frequencies in T 0 . Let χ 2 (T ) be the value of (1) in table T .
Clustering the rows and/or columns of a contingency table reduces the value of the χ 2 statistic, that is χ 2 (T ) ≤ χ 2 (T 0 ) (see Govaert and Nadif (2018) for a detailed proof). Therefore, to see whether we can preserve the dependence structure between the variables U and V when their categories are clustered, we seek, for a fixed size, the clustered table T which maximizes χ 2 (T ). Repeating this procedure for different values of the granularity, we can either identify an extreme grouping or conclude that it does not exist, namely the two variables are dependent regardless of the size of the clustered table. Indeed, in the event of obtaining a clustered table T so that statistical dependence is assumed and its clustered table exhibits independence, then we say that T is an extreme grouping.
In order to obtain new categories in table T , namely the clusters, which are meaningful for the analyst, being able to incorporate prior knowledge about the groups of categories which are allowed or not to be merged would be helpful. In other words, not every possible combination of categories is allowed. Besides easing the interpretability, clustering can be used to deal with sparsity issues in the entries of T 0 looking for aggregations of columns and/or rows which accumulate at least a certain number of observations. Let T (T 0 ) be the set of all possible contingency tables resulting from allowable groups of rows and columns of T 0 .
The problem of Clustering a Contingency Table (CCT) described above is stated as the following combinatorial optimization problem: (CCT) seeks the table T ∈ T (T 0 ) that maximizes the strength of the association between the variables in the clustered table T measured through the χ 2 statistic in (1), and satisfies the structure imposed on the clusters through the definition of the feasible set T (T 0 ). (CCT) is a combinatorial problem for which an assignment (0-1 nonlinear) formulation is proposed in the next section. Other formulations are also possible, such as a set partitioning one, which is stated in Sect. 3.3.

An assignment formulation and its set partitioning counterpart
This section is devoted to developing a mathematical optimization formulation for the (CCT) model stated in Sect. 2. An assignment formulation is proposed in Sect. 3.1, in which the decisions to be made are whether a column of the observed contingency table T 0 is assigned to a cluster of categories in the clustered table T or not, yielding a 0 − 1 nonlinear optimization model. Section 3.2 is devoted to formally model some structures which could be demanded to the clustered table T to, for instance, get meaningful clusters by incorporating expert knowledge on the allowable groupings or reduce sparsity. Recall that these conditions naturally arise from the problem under study. Finally, a set partitioning formulation is proposed in Sect. 3.3. Clustering a contingency table T 0 can be done either row-wise (only the rows are clustered while the initial columns in T 0 are maintained), column-wise (only the columns are clustered while the initial rows in T 0 are maintained), or in both directions, this is, column and rows are both clustered into new categories. Whereas our approach is valid for any of these three options, the assignment formulation for (CCT) and its extensions are fully developed column-wise for the sake of clarity.

The assignment formulation for (CCT) with k clusters
Recall that {v 1 , . . . , v c } is the set of categories (columns) of variable V in the observed contingency table T 0 . The categories in T 0 are aimed to be clustered into k new categories, named asṽ 1 , . . . ,ṽ k , k ≤ c in such a way that the categories in T form a partition of the ones in T 0 and the χ 2 statistic of the clustered table, that is χ 2 (T ), is maximized.
Let y jl for all j ∈ {1, . . . , c} and l ∈ {1, . . . , k} be a binary decision variable defined as be the χ 2 statistic for the clustered table T , which is defined by the y-variables and (1) as: The problem of clustering the columns of a contingency table T 0 into k clusters maximizing the χ 2 statistic is stated as a 0-1 nonlinear optimization model, which consists of the maximization of a convex function subject to linear constraints, as follows: Constraint (3) ensures that each category in T 0 goes to just one of the new categories (clusters) in T and constraint (4) imposes that each cluster has at least one category. Finally, constraint (5) defines the binary nature of y-variables. Note that problem (2)-(5) can be enriched with constraints to break the symmetry associated with the clusters.

Modelling some clustering structures
The assignment formulation in (2)-(5) provides a flexible framework to incorporate additional requirements on the clusters in table T in a straightforward manner, namely as additional linear constraints. We describe in what follows some of the most natural cases, although more complex structures, such as the constrained discretization of continuous variables for supervised learning algorithms, can also be modeled using the y-variables in the assignment formulation stated above.
• Non-sparsity constraints: In contingency tables analysis usually happens that some of the observed joint frequencies are equal or close to zero, and thus statistical conclusions cannot be inferred from the distribution of the χ 2 statistic. In order to be able to apply statistical inference theory, sparsity problems in the observed table T 0 might be mitigated by clustering some of its columns by imposing a threshold over the number of observations in each row of the new column (cluster). In other words, the user might require that in each row of the columns in the new table T there are at least β observations. A common value for β is 5. Such condition can be added as a constraint to problem (2)-(5) as follows: c j=1 o i j y jl ≥ β, i = 1, . . . , r , l = 1, . . . , k.
• Cannot-link constraints: A cannot-link constraint is used to specify that two or more specified categories in T 0 cannot be associated with the same cluster in T . In its simplest case, namely two categories v j and v j in T 0 cannot be grouped together in T , the cannot-link constraint is modeled as Condition (7) can be easily generalized to accommodate groups of categories in T 0 which cannot belong to the same cluster.
A complementary set of conditions to cannot-link ones are the so-called must-link constraints, which are used to specify that two or more specified categories in T 0 must be assigned to the same cluster in T . Although these kind of conditions could be also easily modeled in a similar fashion, they can be imposed in a preprocessing step. • Relational constraints: There might be structural conditions among categories which are more complex than the ones given by cannot or must-link constraints.
That is the case of, for instance, the existence of a partial order relation ≺ between the categories implying that, if two categories belong to one cluster then all the categories in-between must belong to the same cluster too. In its simplest case, namely two categories v j and v j in T 0 such that v j precedes v j in the partial order, the so-called relational constraint is modeled as y jl + y j l ≤ y j l + 1, for v j ≺ v j ≺ v j and l = 1, . . . , k.
• Demand / capacity constraints: We may also require that each column l (clusters) in T contains at least a l categories of T 0 and/or no more than b l , thus establishing demand and/or capacity constraints, respectively, for l = 1, . . . , k. Such conditions can be added as constraints to problem (2)-(5) as follows: • 'et al.' clustering: Given a contingency table T 0 with c columns, the analyst might be interested in obtaining a clustered table T with k columns in which k − 1 of its categories are exactly k − 1 of the categories in T 0 and the k-th category is made up of the aggregation of the remaining c − k + 1 categories in T 0 , that is the 'et al. ' category. This structure is a particular case of (CCT) with k clusters, in which k −1 clusters are singletons and the k-th category in the new table comprises c − (k − 1) categories. In order to get such structure in T , constraint (4) in the formulation (2)-(5) for (CCT) must be replaced by c j=1 y jl = 1, l = 1, . . . , k − 1, (4a) c j=1 y jk = c − k + 1.

A set partitioning formulation for (CCT)
In this section, an alternative formulation is proposed for (CCT), which assumes that there is a list of permissible aggregations of columns in T 0 that can be used to build the clustered table T . Given such a list, a set partitioning formulation is proposed whose benefits with respect to an assignment one are twofold: first, its continuous relaxation is, in general, tighter (Freling et al. 2003), and second, the so-obtained formulation becomes 0-1 linear.
In order to state a set partitioning formulation for (CCT), let {S 1 , . . . , S K } be a family of K subsets of the categories in T 0 given by {v 1 , . . . , v c }. These subsets represent the list of allowable aggregations of columns in T 0 , that is the list of clusters which can be used to build the columns in the clustered table T . Let A be a c × K 0-1 matrix with entries a j p for all j ∈ {1, . . . , c} and p ∈ {1, . . . , K } defined by and let x p for all p ∈ {1, . . . , K } be a binary decision variable defined by Let χ 2 ({x p } p∈{1,...,K } ) be the χ 2 statistic for the clustered table stated as Then, the set partitioning formulation for problem (CCT) is stated as: s.t.
Whereas (14) ensures that each column of T 0 belongs to just one cluster in T , constraint (15) imposes the binary nature of the x-variables. We point out that (13)-(15) is a 0-1 linear optimization problem, which can easily accommodate the structures discussed in Sect. 3.2. For instance, the relational condition requiring that if columns v j and v j in T 0 belong to the same cluster in table T then column v j belongs also to that cluster is defined through one of the S-sets as {v j , v j , v j }.

Illustrative examples
In order to illustrate the methodology proposed in this paper, a contingency table T 0 obtained from a medical study by Kandoth et al. (2013) is considered. This table comprises information about n = 9786 biological samples and its joint frequencies correspond to the number cross-classified cases of the categorical variables cancer type (U ), which has r = 11 categories, and significantly mutated gene (V ), which has c = 127 categories. The set of genes are divided into 20 groups, which are defined according to biological features. This contingency table can be obtained from the supplementary material in Kandoth et al. (2013) and it is also included in Tables 8  and 9 in the Appendix. In Table 3 we show, for each of the groups g = 1, . . . , 20: its description according to Kandoth et al. (2013), a color to represent it in the upcoming results, the number of genes (categories) in each group (S g ), and the percentage of cells the contingency table within each group which are sparse (that is, for each i = 1, . . . , 11 the cardinality of {o i j : o i j < 5, j = 1, . . . , S g } divided by S g times 100). We point out the noticeable amount of joint frequencies which are below the usual threshold of 5, being the level of sparsity greater than 49% within all the groups. For the interpretation of references to color in Table 3, the reader is referred to the web version of this article.   In what follows, we present results obtained from clustering the contingency table of Kandoth et al. (2013) using two different clustering structures. The optimization models involved in the experiments have been solved using Bonmin (Bonami and Lee 2017) under Pyomo. Bonmin is an open-source numerical optimization procedure for solving general Mixed Integer Nonlinear Programs by means of Branch-and-Bound and Branch-and-Cut algorithms, thus avoiding an explicit complete enumeration of all the feasible solutions of the models stated in Sect. 3. In order to avoid being stuck  at local optima, each problem has been solved 10 times with a time limit of one hour in a PC Intel ® Core TM i7-7700, 16GB of RAM. As abovementioned, the contingency table T 0 in Kandoth et al. (2013) is highly sparse. It is well known that in such case the asymptotic distribution of the χ 2 statistic fails and thus statistical conclusions about the statistical dependence between U and V cannot be inferred (Agresti and Yang 1987). In order to overcome such limitation, a less granular representation of the genes can be considered so that the columns of table T 0 (genes) are clustered into broader categories (groups of genes) in such a way that the aggregated joint frequencies are larger than a threshold, and thus an eventual statistical dependence between U and V in the original table T 0 could be revealed. To do so, the optimization model defined by (2)-(6) is solved for β = 5. Tables 4 and 5 contain the assignment of the genes in T 0 to the clusters (new categories made up of groups of genes) in T for the number of clusters k varying from 2 to 10. Thus, a initial table T 0 with r = 11 and c = 127 is reduced to tables T of r = 11 and c = k, for k = 2, . . . , 10. First column of Tables 4 and 5 contains the values of the χ 2 statistic in T (and the associated p−value in parentheses). We point out that statistical dependence between U and V is detected when the granularity of V is fixed to k = 2, . . . , 10 and a significance level of α = 5%. Thus, we can conclude that, under the aforementioned conditions, there does not exist a extreme grouping since the null hypothesis of statistical independence is rejected for all k ≥ 2.
The genes in T 0 are split into 20 groups defined through biological features. A plausible requirement could be that the genes of some groups in T 0 must belong to the same clusters in T to avoid having genes of the same groups spread out across different clusters as in Tables 4 and 5. These must-link conditions can be imposed in a preprocessing step. In this case, we impose that genes in Groups 1, 2, 4 and 6, respectively, must belong to the same cluster. In addition, our preprocessing incorporates the requirement that Groups 4 and 6 belong to the same cluster. Conversely, some groups might be required to belong to different clusters. This structure is illustrated by imposing constraint (7) for j ∈ Group 1 and j ∈ Group 2. Tables 6 and 7 contain the assignment of the genes in T 0 to the clusters (new categories) in T for the number of clusters k varying from 2 to 10. As before, the first column of such tables contains the values of the χ 2 statistic in the clustered table T (and the associated p−value in parentheses). In this case, statistical dependence between U and V is also detected when the granularity of V is fixed to k = 2, . . . , 10, a significance level of α = 5% is considered as well as the group structures in the clustering process. Thus, we can conclude that, under the aforementioned conditions, there does not exist a extreme grouping since the null hypothesis of statistical independence is rejected for all k ≥ 2.
The contingency tables obtained form the clusterings shown in Tables 4, 5, 6 and 7 are depicted in the Supplementary Material. We refer the reader to the web version of this article for the interpretation of references to colors in Tables 4-7.

Conclusions
In this paper we have addressed the problem of clustering categories in contingency tables maximizing the χ 2 statistic (Mirkin 2001;Pearson 1900). Solving this clustering problem for different sizes, namely different granularities of the categorical variables under study, allows us to identify extreme groupings or, in other words, the way categories can be clustered into larger ones so that the dependence of the variables is no longer detected if the granularity of the variables is reduced. To do so, a combinatorial mathematical optimization model has been stated, which allows to accommodate structural properties of the clusters in the clustered table which naturally arise in the context of the dataset under study. An assignment formulation has been proposed for such model, namely (CCT), yielding a 0-1 nonlinear optimization problem. Requirements on the clusters, such as non-sparsity conditions, relational and cannot link constraints, have been stated as linear constraints. In addition, a set partitioning reformulation of (CCT) is also proposed. Our methodology is illustrated using a dataset in a medical study, which naturally demands the use of the tools proposed in this paper to handle the study of statistical dependence between its variables under structural conditions on the clusters.
The problem studied in this paper can be extended in a few directions. First, other criteria different from χ 2 could be considered to measure statistical dependence (Goodman and Kruskal 1979;Joe 1989). Second, exploring different criteria to group the categories in contingency tables different from statistical dependence, and which are defined through appropriate combinatorial optimization models, could be also explored as an extension to this paper. Some interesting examples could be to explore patterns in the observed joint frequencies to group the categories in a contingency table, or to identify those patterns in the coordinates given by Correspondence Analysis (Ciampi et al. 2005;Pledger and Arnold 2014;van de Velden et al. 2020). Third, tighter formulations for (CCT) could be explored, in combination with metaheuristic approaches such as the Variable Neighborhood Search (Mladenović and Hansen 1997) or the Large Neighborhood Search (Pisinger and Ropke 2010), to address larger tables. Finally, extensions of the proposed methodology for dealing with multi-way tables require further research (Agresti and Gottard 2007).