Variability Analysis for Nominal Variables I have a very large datasets (billions of observations) made of multiple nominal categorical variables (nominal, not ordinal), and I want to outline the set of variables that accounts for the most variability in the overall dataset. Multiple Correspondence Analysis could be used for that purpose, but it relies on the use of an indicator matrix (aka complete disjunctive table), which can become extremely large when dealing with large sets of records and large numbers of categorical dimensions that have large cardinalities, because its size is the product of the number of records by the product of the cardinalities.
Furthermore, MCA will outline which variables account for the most variability, but it won’t help identify which combinations of variables account for it the most. This is due to the fact that when using MCA, data is represented in a low-dimensional Euclidian space, while interesting patterns might emerge when representing data in spaces of higher dimensionality, and the optimal dimensionality is usually not known in advance.
In such a context, I am trying to figure out if Pivot Tables could be used as an alternative to Indicator Matrices or Burt Tables. According to this method, a set of n discrete variables are selected. For simplification purposes, all variables are considered as nominal, even though some of them might actually be ordinal.
From there, we compute multiple sets of pivot tables, taking all combinations of k dimensions as pivot axes, with k ranging from 1 to n, and using COUNT as aggregation function. Consequently, for a given value of k, the number of computed pivots will be equal to the binomial coefficient, or n! / (k! × (n - k)!). In total, 2ⁿ pivots will be computed.
For every pivot, we then compute sets of average counts (averages of counts). For a pivot defined with k pivot axes, we will compute averages across k-1 dimensions, and we will do that for every combinations of k-1 values of the k-1 dimensions. For a given combination of k-1 dimensions, which is defined by the dimension being excluded from the set, the number of averages will be equal to the product of the cardinalities of the dimensions that are included in the combination.
For example, if k=2 and the first dimension (displayed in columns) is made of c individual values and the second dimensions (displayed in rows) is made of r individual values, there would be c averages on a column-by-column basis, and r averages on a row-by-row basis.
These averages will then be distributed across k buckets, one for every pivot axis, with each average going to the bucket related to the dimension not included in the k-1 dimensions used to compute the average. Then, we will compute for every bucket the relative standard deviation of the averages contained in the bucket. Finally, we will compute the average of all relative standard deviations across all buckets. This will give us a normalized variability coefficient for every pivot, which (we hope) quantifies the degree to which the pivot’s dimensions account for the variability of the overall dataset.
My questions are the following:


*

*Are there better ways of achieving similar results?

*Has such a method been used before (links please)?

*Is such a method mathematically correct?


Note: This article provides a more complete description of the method.
 A: There is no single best, textbook answer to your question. The literature offers many heuristics and workarounds, each with a large number of possible options. Your creativity in making some of these recommendations work in drilling down to a final solution will determine the success of this initiative.
In my opinion, it is computationally prohibitive to explore the combinatorics of these multinomial factors on the full set of "billions" of records. I would consider making your life easier in any preliminary, exploratory phase and take a 5% random sample or k-fold set of samples, reserving the full set of billions of records for the close-to-final analyses. 
Regardless of whether or not you sample, there are many workarounds to the sheer technical challenges associated with analyzing massive numbers of candidate predictors in the presence of massive amounts of information. Here's a paper that is the most recent and best summary of many of these workarounds that I'm aware of:  A Survey of Statistical Methods and Computing for Big Data by Wang, Chen, Schifano and Win (2015). 
In addition to the purely mechanical challenges, there are a multitude of published solutions to variable selection with massive numbers of candidate predictors. In point of fact, the number of possible variable selection methods is practically limitless. This literature began with Breiman's paper on "random forests," but Breiman was concerned only with predictive accuracy and minimizing the MSE. Other considerations of model performance are just as important, parameter stability being one since that can directly impact truly out-of-sample model performance.  In my opinion, Peter Buhlmann at ETH Zurich is one of the best in this field, e.g., see his paper Stability Selection with Meinshausen (2009). 
Many software packages offer logit model modules for contingency table analysis, which is a functional form appropriate for your type of data (it's not the only one, of course). Typically included in those modules is an option that allows you to automatically take up to p-way interactions or combinations based on p, the number of model inputs. 
In a first, exploratory phase, combining this logit modeling approach with one of the "random forest" type workarounds reviewed by Wang, et al., would create potentially millions of "mini-models" where standardized metrics (chi-square values?) of the performance of each combination of attributes are retained. In this way, a distribution of the conditional performance across many randomly drawn subsamples of data and attributes for a specific combination of attributes can be developed. Running these potentially millions of "mini-models" on a massively parallel platform might require a few hours of CPU, depending on the machine, etc. 
If you've understood things up to this point, the next step involves aggregating (averages? medians?) and ranking the scale invariant metric(s) associated with each parameter for a specific combination of attributes across the "mini-models" should give you a good line-of-sight into the rank-order relationships across the full range of possible combinations. 
Variable selection based on the standardized metric(s) by identifying a cut-off or threshold for inclusion or exclusion into a next, closer-to-final phase of modeling and analysis. Again, there are no rules of thumb available for making these decisions: the criteria would entirely be up to you. The fact is, and even if you're doing this for a dissertation, as long as you can motivate your approach and protocol based on a careful reading of the literature, at this stage of the development of workarounds and heuristics for modeling massive amounts of information, there is enough diversity of opinion out there that you can't be faulted or found wrong. Whichever variable selection method you've chosen, the entire process amounts to a kind of Darwinian, "survival-of-the-fittest" approach to attribute attrition.
This is klugey, I agree, but today's kluge is tomorrow's best practice.
