Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
64,634 characters · 19 sections · 39 citation commands
MaxEnt-Copula-Oct2020
\noindentKeywords: Maximum entropy; Self-adaptive copula model; LP-Fourier transform; Categorical data analysis; Copula-logistic regression; United Statistical learning. \linespread{1.27} {1.5ex}
Copulas are the `bridge’ between the univariate and the multivariate statistics world, with applications in a wide variety of science and engineering fields---from economics to finance to marketing to healthcare. Because of the ubiquity of copula in empirical research, it is becoming necessary to develop a general theory that can unify and simplify the copula learning process. In this paper, we present a new class specially-designed nonparametric maximum-entropy (MaxEnt) copula model that offers the following advantages: First, it yields a bonafide (smooth, non-negative, and integrates to $1$) copula density estimate with interpretable parameters that provide insights into the nature of the dependence between the random variables $(X, Y)$. Secondly, the method is data-type agnostic---which is to say that it automatically (self) adapts to mixed-data types (any combination of discrete, continuous, or even categorical). Thirdly, and most notably, our copula-theoretic framework subsumes and unifies a wide range of statistical learning methods using a common mathematical notation---unlocking deep, surprising connections and insights, which were previously unknown. In the development of our theory and algorithms, the LP-Fourier method of copula modeling (which was initiated by D20copula), plays an indispensable role.
We introduce two new classes of maximum-entropy (MaxEnt) copula density models. But before diving into technical details, it will be instructive to review some basic definitions and concepts related to copula.
\vskip.25em {\bf Sklar’s Copula Representation Theory} sklar1959. The joint cumulative distribution function (cdf) of any pair of random variables $(X,Y)$ \[F_{X,Y}(x,y)=\Pr(X\le x, Y\le y),~~\mathrm{for}~(x,y) \in \cR^2\] can be decomposed as a function of the marginal cdfs $F_X$ and $F_Y$ \beq F_{X,Y}(x,y)\,=\,\Cop_{X,Y}\big( F_X(x), F_Y(y)\big), \mathrm{for} (x,y) \in \cR^2 \eeq where $\Cop_{X,Y}$ denotes a copula distribution function with uniform marginals. To set the stage, we start with the continuous marginals case, which will be generalized later to allow mixed-$(X,Y)$. Taking derivative of Eq. (ref), we get \beq f_{X,Y}(x,y)\,=\,f_X(x) f_Y(y)\cop_{X,Y}\big( F_X(x), F_Y(y)\big), \mathrm{for} (x,y) \in \cR^2 \eeq which decouples the joint density into the marginals and the copula. One can rewrite Eq. (ref) to represent copula as a “normalized” joint density function \beq \cop_{X,Y}\big( F_X(x), F_Y(y)\big)\,:=\,{\rm dep}_{X,Y}(x,y)\,=\dfrac{f_{X,Y}(x,y)}{f_X(x) f_Y(y)}, \eeq which is also known as the dependence function, pioneered by Hoeff40. To make (ref) a proper density function (i.e., one that integrates to one) we perform quantile transformation by substituting $F_X(x)=u ~\text{and}~ F_Y(y)=v$: \[ \cop_{X,Y}(u,v)={\rm dep}_{X,Y}( F_X^{-1}(u), F_Y^{-1}(v))=\dfrac{f_{X,Y}\big( F_X^{-1}(u), F_Y^{-1}(v) \big)}{f_X( F_X^{-1}(u) ) f_Y(F_Y^{-1}(v) )}, ~~(u,v) \in [0,1]^2. \] We are now ready to extend this copula density concept to the mixed $(X,Y)$ case. \vskip.5em {\bf Pre-Copula: Conditional Comparison Density} D13a. Before we introduce the generalized copula density, we need to introduce a new concept---conditional comparison density (CCD). For a continuous $X$, CCD is defined as: \beq d(u;X,X|Y=y)\,=\,\dfrac{f_{X|Y}\big( F_X^{-1}(u) \mid y \big)}{f_X\big( F_X^{-1}(u)\big)}, 0<u<1 \eeq For $Y$ discrete, we represent it using probability mass function (pmf): \beq d(v;Y,Y|X=x)\,=\,\dfrac{p_{Y|X}\big( Q_Y(v)|x \big)}{p_Y\big( Q_Y(v)\big)}\,=\,\dfrac{\Pr(Y=Q_Y(v)|X=x)}{\Pr(Y=Q_Y(v))}, 0<v<1. \eeq where $Q_Y(v)$ is the quantile function of $Y$. It is easy to see that the CCDs (ref) and (ref) are proper densities in the sense that \[\int_0^1 d(u;X,X|Y=y)\dd u\,=\,\int_0^1 d(v;Y,Y|X=x) \dd v\,=\,1.\] \vskip.12em {\bf Generalized Copula Representation Theory} D20copula. For the mixed case, when $Y$ is discrete and $X$ is continuous the joint density of (ref) is defined by either side of the following identity: \[ \Pr(Y\mid X=x) f_X(x)\,=\,f_{X|Y}(x|y) \Pr(Y=y). \] This can be rewritten as the ratios of conditionals and their respective marginals: \beq Pre-Bayes' Rule: \dfrac{\Pr(Y=y\mid X=x)}{\Pr(Y=y)}\,=\,\dfrac{f_{X|Y}(x|y)}{f_X(x)} \eeq This formula (ref) can be interpreted as the slices of the mixed copula density, since \beq \cop_{X,Y}\big( F_X(x), F_Y(y) \big)\,=\,\dfrac{f_{X|Y}(x|y)}{f_X(x)} \,=\,\dfrac{\Pr(Y=y\mid X=x)}{\Pr(Y=y)}. \eeq Substituting $F_X(x)=u$ and $F_Y(y)=v$, we get the following definition of the generalized copula density in terms of conditional comparison density (CCD): \beq \cop_{X,Y}(u, v)\,=\,d\big(u;X,X|Y= Q(v; Y)\big)\,=\,d\big(v;Y, Y|X=Q(u;X)\big), 0<u,v<1.\eeq Bayes' theorem ensures the equality of two CCDs, with copula being the common value. Equipped with this fundamentals, we now develop the nonparametric theory of MaxEnt copula modeling.
An exponential Fourier series representation of copula density function is given. The reasons for entertaining an exponential model for copula is motivated from two different perspectives.
\vskip.35em {\bf The Problem of Unboundedness}. One peculiar aspect of copula density function is that it can be unbounded at the corners of the unit square. In fact, many common parametric copula families---Gaussian, Clayton, Gumbel, etc.---tend to infinity at the boundaries. So naturally the question arises: How to develop suitable approximation methods that can accommodate a broader class of copula density shapes, including the unbounded ones? The first key insight: logarithm of the copula density function is far more convenient to approximate (due to its well-behaved nature) than the original density itself. We thus express the logarithm of copula density $\log \cop_{X,Y}$ in the Fourier series---instead of doing canonical $L_2$ approximation, which expands $\cop_{X,Y}$ directly in an orthogonal series D20copula. Accordingly, for densities with rapidly changing tails, `log-Fourier' method leads to an improved estimate that is less wiggly and more parsimonious than the $L_2$-orthogonal series model. In addition, the resulting exponential form guarantees the non-negativity of the estimated density function. \vskip.1em
Choice of Orthonormal Basis. To expand log-copula density function, we choose the LP-family of polynomials (see Appendix (ref)), which are especially suited to approximate functions of mixed random variables. In particular, we approximate $\log \cop_{X,Y}$ by expanding it in the tensor-product of LP-bases $\{S_j \otimes S_k\}$, which are orthonormal with respect to the empirical-product measure $\{\wtF_X \otimes \wtF_Y\}$. LP-bases' appeal lies in its ability to approximate the quirky shapes of mixed-copula functions in a completely automated way; see Fig. (ref). Consequently, it provides a unified way to develop nonparametric smoothing algorithms that simultaneously hold for mixed data types. \vskip.3em
{\bf The Maximum Entropy Principle}. Another justification for choosing the exponential model comes from the principle of maximum entropy (MaxEnt), pioneered by E. T. jaynes1957. The maxent principle defines a unique probability distribution by maximizing the entropy $H(\cop)=-\int \cop_{X,Y} \log \cop_{X,Y}$ under the normalization constraint $\int \cop_{X,Y}=1$ and the following LP-co-moment conditions: \beq \Ex_{\cop_{{\bm \te}}}[S_j(U;X) S_k(V;Y)]\,=\,\LP_{jk}. \eeq
LP-co-means are orthogonal “moments” of copula, which can be estimated by \beq \tLP_{jk}\,=\,\Ex_{\widetilde \Cop}\big[ S_j(U;X) S_k(V;Y) \big]\,=\,\dfrac{1}{n} \sum_{i=1}^n S_j\big( \wtF_X(x_i) ;X\big) \, S_k\big( \wtF_Y(y_i) ;Y\big).\eeq Applying calculus of variations, one can show that the maxent constrained optimization problem leads to the exponential (ref) form. The usefulness of Jaynes' maximum entropy principle lies in providing a constructive mechanism to uniquely identify a probability distribution that is maximally non-committal (flattest possible) with regard to all unspecified information beyond the given constraints.
{\bf Estimation}. We fit a truncated exponential series estimator of copula density \[\cop_{\teb}(u,v;X,Y)~=~\dfrac{1}{Z_{\bm \te}} \exp\Big\{ \sum_{j=1}^{m_1}\sum_{k=1}^{m_2} \theta_{jk}S_j(u;X) S_k(v;Y)\Big\}.\] The task of finding the maximum likelihood estimates (MLE) of ${\bm \te}$ boils down to solving the following sets of equations for $j=1,\ldots,m_1$ and $k=1,\ldots,m_2$: \beq \frac{\partial \log \cop_{\teb}}{\partial \te_{jk}} \equiv \frac{\partial \log Z_{\bm \te}}{\partial \te_{jk}}\,-\,\dfrac{1}{n} \sum_{i=1}^n S_j\big( \wtF_X(x_i) ;X\big) \, S_k\big( \wtF_Y(y_i) ;Y\big)\,=\,0. \eeq Note that the derivative of the log-partition function is equal to the expectation of the LP-co-mean functions: \beq \dfrac{\partial \log Z_{\bm \te}}{\partial \te_{jk}}\,=\,\Ex_{\cop_{{\bm \te}}}[S_j(U;X) S_k(V;Y)]. \eeq Replacing (ref) and (ref) into (ref) implies that the MLE of MaxEnt model is same as the method of moments estimator satisfying the following moment conditions: \[ \iint_{[0,1]^2} S_j(u;X) S_k(v;Y) \cop_{\teb}(u,v;X,Y)\dd u \dd v\,=\,\tLP_{jk}.\] At this point, one can apply any convex optimization\footnote{since the second derivative of the log-partition function is a s positive semi-definite covariance matrix} routine (e.g., Newton's method, gradient descent, stochastic gradient descent, etc.) to solve for $\widehat{{\bm \te}}$.
\vskip.2em {\bf Asymptotic}. Let the sequence of $m_1$ and $m_2$ increase with sample size with an appropriate rate $\frac{(m_1 m_2 )^3}{n} \to 0$ as $\nti$. Then, under certain suitable regularity conditions, the exponential $\cop_{\widehat \teb}$ is a consistent estimate in the sense of Kullback-Leibler distance; see barron91 for more details.
\vskip.2em {\bf Determining Informative Constraints}. Jayne's maximum entropy principle assumes that a proper set of constraints (i.e., sufficient statistic functions) are given to the modeler, one that captures the phenomena under study. This assumption may be legitimate for studying thermodynamic experiments in statistical mechanics or for specifying prior distribution in Bayesian analysis, but certainly not for building empirical models.
Which comes first: a parametric model or sufficient statistics? After all, the identification of significant components (sufficient statistics) is a prerequisite for constructing a legitimate probability model from the data D11a2. Therefore the question of how to judiciously design and select the constraints from data, seems inescapable for nonparametrically learning maxent copula density function from data; also see Appendix (ref), which discusses the `two cultures' of maxent modeling. We address this issue as follows: (i) compute $\tLP_{jk}$ using the formula eq. (ref); (ii) sort them in descending order based on their magnitude (absolute value); (iii) compute the penalized ordered sum of squares \[{\rm PenSum}(q)~=~\text{Sum of squares of top $q$ LP comeans}~-~\dfrac{\gamma_n}{n}q.~~\] For AIC penalty choose $\gamma_n=2$, for BIC choose $\gamma_n=\log n$, etc. Further details can be found in D20copula. (iv) Find the $q$ that maximizes the ${\rm PenSum}(q)$. Store the selected indices $(j,k)$ in the set $\mathcal{I}$. (v) Carry out maxent optimization routine based only on the selected LP-sufficient statistics-based constraints: $$\Big\{ S_j(u;X)S_k(v;Y)\Big\},~~ (j,k) \in \mathcal{I}.~~~~~~~$$ This pruning strategy guards against overfitting. Finally, return the estimated reduced-order (with effective dimension $|\mathcal{I}|$) maxent copula model.
We provide a second parameterization of copula density.
{\bf Connection}. Two fundamental representations, namely the log-bilinear (ref) and loglinear (ref) copula models, share some interesting connections\footnote{See D20copula for a parallel result on the LP-orthogonal series copula model.}. To see that perform singular value decomposition (SVD) of the $\Theta$-matrix whose $(j,k)$th entry is $\te_{jk}$: \[\Theta\,=\,U\Omega V^T\hspace{-.5em},\] $u_{ij}$ and $v_{ij}$ are the elements of the singular vectors with singular values $\mu_1 \ge \mu_2 \ge \cdots \ge 0$. Then the spectral bases can be expressed as the linear combinations of the LP-polynomials: \bea \phi_k(u;X)\,=\,\sum\nolimits_{j}u_{jk}S_{j}(u;X) \\ \psi_k(u;Y)\,=\,\sum\nolimits_{l}v_{lk}S_{l}(v;Y). \eea Hence, the LP-spectral functions ((ref)-(ref)) satisfy the following orthonormality conditions: \beas \int\phi_k(u;X) \dd u\,=\,\int\psi_k(v;Y) \dd v\,=\,0 \\ \hskip1em\int \phi_j(u;X) \psi_k(u;Y) \dd u\,=\,\delta_{jk}, \,for j\neq k. \eeas
We demonstrate the flexibility of the LP-copula models using real data examples. \vskip1.15em
The scope of the general theory of the preceding section goes far beyond simply a tool for nonparametric copula approximation. In this section, we show how one can derive a large class of applied statistical methods in a unified manner by suitably reformulating them in terms of the LP-maxent copula model. In doing so, we also provide statistical interpretations of LP-maxent parameters under different data modeling tasks.
Categorical data analysis will be viewed through the lens of LP-copula modeling. Let $X$ and $Y$ denote two discrete categorical variables; $X$ with $I$ categories and $Y$ with $J$ categories. The data are summarized in an $I\times J$ contingency table $\F$: $f_{kl}$ is the observed cell count in row $k$ and column $l$ of $\F$ and $n=f_{++}$ is the total frequency. The row and column totals are denoted as $f_{k+}$ and $f_{+l}$. The observed joint $\Pr(X=k,Y=l)$ is denoted by $\tp_{kl}=f_{kl}/f_{++}$; the respective row and column marginals are given by $\tp_{k+}=f_{k+}/f_{++}$ and $\tp_{+l}=f_{+l}/f_{++}$.
\vskip.5em {\bf LP Log-linear Model}. We specialize our general copula model (ref) for two-way contingency tables. The discrete LP-copula for the $I\times J$ table is given by \beq \cop(F_X(k), F_Y(l)) \,=\, \exp\Big( \mu_0 + \sum_{j=1}^m \mu_{j} \phi_{jk}\psi_{jl}\Big), \eeq where we abbreviate the row and columns scores $\phi_j(F_X(k))=\phi_{jk}$ and $\psi_j(F_Y(l))=\psi_{jl}$ for $k=1,2,\ldots,I$ and $l=1,2,\ldots,J$. The number of components $m \le M=\min(I-1,J-1)$; we call the log-linear model (ref) `saturated' (or `dense') when we have $m=M$ components. The non-increasing sequence of model parameters $\mu_j$'s are called “intrinsic association parameters” that satisfy\footnote{Compare our equation (ref) with equation (34) of good96.} \beq \mu_j\,=\,\sum_{k=1}^I\sum_{l=1}^J \Big( \log p_{kl}\Big) p_{k+} p_{+l} \phi_{jk} \psi_{jl}, for j=1,\ldots,m. \eeq Note that the discrete LP-row and column scores, by design, satisfy (for $j \neq j'$): \bea \sum_{k=1}^I\phi_{jk} p_{k+} \,= \sum_{l=1}^J \psi_{jl}p_{+l} = 0 \\ \sum_{k=1}^I\phi^2_{jk} p_{k+} \,= \sum_{l=1}^J \psi^2_{jl}p_{+l} = 1 \eea \bea \sum_{k=1}^I \phi_{jk} \phi_{j'k} p_{k+} \,= \sum_{l=1}^J \phi_{jl} \psi_{j'l} p_{+l}\,=\,0. \eea \vskip.5em {\bf Interpretation}. It is clear from (ref) that the parameters $\mu_j$'s are fundamentally different from the standard Pearsonian-type correlation ${\rm Cor}(\phi_{j}, \psi_{j})=\Ex[\phi_{j}\psi_{j}]$, due to (ref)--(ref): \beq \rho_j\,=\,\sum_{k=1}^I\sum_{l=1}^J p_{kl}\phi_{jk} \psi_{jl}, for\,j=1,\ldots,m.\eeq The coefficients of the LP-MaxEnt-copula expansion for contingency tables carry a special interpretation in terms of log-odds-ratio. To see this we start by examining the $2\times 2$ case. \vskip.5em {\bf The $2\times 2$ Contingency Table}. Applying (ref) for two-by-two tables we have \beq \mu_1\,=\,\sum_{k=0}^1\sum_{l=0}^1 \Big( \log p_{kl}\Big) p_{k+} p_{+l} \phi_{1k} \psi_{1l}.\eeq Note that for dichotomous $X$ the LP-spectral basis $\phi_1(F_X(x))$ is equal to $T_1(x;F_X)$. Consequently, we have the following explicit formula for $\phi_1$ and $\psi_1$: \bea \phi_1(F_X(x))&=&\dfrac{x-p_{2+}}{\sqrt{p_{1+} p_{2+}}}\\ \psi_1(F_Y(y))&=&\dfrac{y-p_{+2}}{\sqrt{p_{+1} p_{+2}}} \eea Substituting this into (ref) yields the following important result.
To deduce (ref) from (ref), verify the following, utilizing the LP-basis formulae (ref)-(ref) \beas \phi_{11}-\phi_{10}:=\phi_1(F_X(1)) - \phi_1(F_X(0))=\dfrac{p_{1+}+p_{2+}}{\sqrt{p_{1+}p_{2+}}} = \big(p_{1+}p_{2+}\big)^{-1/2}\\ \psi_{11}-\psi_{10}:=\psi_1(F_Y(1)) - \psi_1(F_Y(0))=\dfrac{p_{+1}+p_{+2}}{\sqrt{p_{+1}p_{+2}}} = \big(p_{+1}p_{+2}\big)^{-1/2} . \eeas \vskip.4em {\bf Reproducing Goodman's Association Model}. Our discrete copula-based categorical data model (ref) expresses the logarithm of “dependence-ratios”\footnote{ good96 calls it “Pearson ratios.”} \beq \cop(F_X(k),F_Y(l))\,=\,\dfrac{p_{kl}}{p_{k+}p_{+l}}\eeq as a linear combination of LP-orthonormal row and column scores satisfying (ref)-(ref). The copula-dependence ratio (ref) measures the strength of association between the $k$-th row category and $l$-the column category. To make the connection even more explicit, rewrite (ref) for two-way contingency tables as follows: \beq \log p_{kl}\,=\,\mu_0 + \mu_k^{{\rm R}} + \mu_l^{{\rm C}} + \sum_{j=1}^m \mu_{j} \phi_{jk}\psi_{jl}, \eeq where $\mu_k^{{\rm R}}$ denotes the logarithm of row marginal $\log p_{k+}$ and $\mu_l^{{\rm C}}$ denotes the logarithm of column marginal $\log p_{+l}$. goodman1991 called this model (ref) a “weighted association model” where weights are marginal row and column proportions. He used the term “association model” (to distinguish it from correlation (ref) based model) as it studies the relationship between rows and columns using odds-ratio.
We describe a graphical exploratory tool---logratio biplot, which allows a quick visual understanding of the relationship between the categorical variables $X$ and $Y$. In the following, we describe the process of constructing logratio biplot from the LP-copula model (ref).
\vskip.3em {\bf Copula-based Algorithm}. Construct two scatter plots based on the top two dominant components of the LP-copula model: the first one is associated with the row categories, formed by the points $(\mu_1 \phi_{1k}, \mu_2 \phi_{2k})$ for $k=1,\ldots,I$; and the second one is associated with the column categories, formed by the points $(\mu_1 \psi_{1l}, \mu_2 \psi_{2l})$ for $l=1,\ldots,J$. Logratio biplot is a two-dimensional display obtained by overlaying these two scatter plots--the prefix `bi' refers to the fact that it shares a common set of axes for both the rows and columns categories.
\vskip.3em {\bf Interpretation}. Here we offer an intuitive explanation of the logratio biplot from the copula perspective. We start by recalling the definition of conditional comparison density (CCD; see eq. (ref)-(ref)), as the copula-slice. For fixed $X=k$, logratio-embedding coordinates $(\mu_j \phi_{jk})$ can be viewed as the LP-Fourier coefficients of the $d(F_Y(y); Y,Y|X=k)$, since \[ d(F_Y(y); Y,Y|X=k)\,=\, \exp\Big\{ \mu_0 + \sum_{j=1}^m \big(\mu_{j} \phi_{jk}\big)\psi_{jl}\Big\}, \] Similarly, the logratio coordinates $(\mu_j \psi_{jl})$ for fixed $Y=l$ can be interpreted as the LP-expansion coefficients of $d(F_X(x); X,X|Y=l)$. Hence, the logratio biplot can alternatively be viewed as follows: (i) estimate the discrete LP-copula density; (ii) Extract the copula slice $\whd(u; Y,Y|X=k)$ along with its LP-coefficients $(\mu_1 \phi_{1k}, \mu_2 \phi_{2k})$; (iii) similarly, get the estimated $\whd(v; X,X|Y=l)$---the copula slice at $Y=F_Y(l)$ along with its LP-coefficients $(\mu_1 \psi_{1l}, \mu_2 \psi_{2l})$; (iv) Hence, the logratio biplot (see Fig. (ref)(b)) measures the association between the row and column categories $X=k$ and $Y=l$ by measuring the similarity between the `shapes' of $\whd(u; Y,Y|X=k)$ and $\whd(v; X,X|Y=l)$ through their LP-Fourier coefficients.
\vskip.25em
It has been known for a long time that classical maximum likelihood-based log-linear models break down when applied to large sparse contingency tables with many zero cells; see fienberg20073c. Here we discuss a new maxent copula-based smooth method for fitting a parsimonious log-linear model to sparse contingency tables.
\vskip.35em
\vskip.2em {\bf A Parsimonious Model}. The estimated smooth log-linear LP-copula model for the Zelterman data is given by: \[ \scalebox{0.9}{\mbox{\ensuremath{\displaystyle \whcop_{X,Y}(u,v)\,=\, \exp\Big\{0.52 \,\phi_1(u;X)\,\psi_1(v;Y)+0.37 \,\phi_2(u;X)\,\psi_2(v;Y)+0.18 \,\phi_3(u;X)\,\psi_3(v;Y)\,-\,0.27\Big\},}}} \] \vskip.1em displayed in Fig. (ref). This shows a strong positive correlation between salary and number of years of experience. However, the most notable aspect is the effective model dimension, which can be viewed as the intrinsic degrees of freedom (df). Our LP-maxent approach distills a compressed representation with reduced numbers of parameters that yields a smooth estimates: only requiring $m=3$ components to capture the pattern in the data---a radical compression with negligible information loss! Contrast this with the dimension of the saturated loglinear model: $(28-1)\times (26-1)=675$ \textemdash a case of a severely overparameterized non-smooth model with inflated degrees of freedom, which leads to an inaccurate goodness-of-fit test for checking independence between rows and columns. More on this in Sec. (ref).
\vskip.3em {\bf Smoothing Ordered Contingency Tables}. The nonparametric maximum likelihood-based cell probability estimates $\tp_{X,Y}(k,l)=f_{kl}/n$ are very noisy and unreliable for sparse contingency tables. By sparse, we mean tables with a large number of cells relative to the number of observations.
Using Sklar's representation theorem, one can simply estimate the joint probability $\hp_{X,Y}(k,l)$ by multiplying the empirical product pmf $\tp_X(k)\tp_Y(l)$ with the smoothed LP-copula. In other words, the copula can be viewed as a data-adaptive bivariate discrete density-sharpening function that corrects the independent product-density to estimate the cell probabilities. \beq dKernel(k,l)\,=\,\whcop_{X,Y}\big( \wtF_X(k), \wtF_Y(l)\big), for k=1,\ldots,I; l=1,\ldots,J.\eeq where the discrete-kernel function satisfies \[ \dfrac{1}{n^2} \sum_{k=1}^I\sum_{l=1}^J \texttt{dKernel}(k,l) f_{k+} f_{+l}\,=\,1. \] This approach can be generalized for any bivariate discrete distribution; see next section.
The topic of nonparametric smoothing for multivariate discrete distributions has received far less attention than the continuous one. Two significant contributions in this direction include: aitchison1976 and simonoff1983. In what follows, we discuss a new LP-copula-based procedure for modeling correlated discrete random variables.
\vskip.35em
Algorithm. The main steps of our analysis are described below.
Step 1. Modeling marginal distributions. We start by looking at the marginal distributions of $X$ and $Y$. As seen in Fig. (ref), negative binomial distributions provide excellent fit. To fix the notation, by $G_{\mu,\phi}= {\rm NB}(y;\mu,\phi)$, we mean the following probability distribution: \[ {\rm NB}(y;\mu,\phi) = \binom{y + \phi - 1}{y} \, \left( \frac{\mu}{\mu+\phi} \right)^{\!y} \, \left( \frac{\phi}{\mu+\phi} \right)^{\!\phi} \!,~~y \in \mathbb{N},\] where $\Ex(X)=\mu$ and $\Var(X)=\mu+\frac{\mu^2}{\phi}$. Using the method of MLE, we get: $X \sim G_1 ={\rm NB}(x;\hat \mu=0.97,\hat \phi=3.60)$, and $Y \sim G_2 ={\rm NB}(y; \hat \mu=4.30,\hat \phi=1.27)$.
\vskip.3em Step 2. Generalized copula density. The probability of bivariate distribution at $(x,y)$ can be written as follows (generalizing Sklar's Theorem): \beq \Pr(X=x,Y=y):=p_{X,Y}(x,y)=g_1(x) g_2(y) \,\cop_{X,Y}\big(G_1(x), G_2(y)\big), \eeq where generalized log-copula density admits the following decomposition: \beq \log \big\{\cop_{X,Y}(G_1(x), G_2(y))\big\}\,=\, \sum_{j=1}^{m_1}\sum_{k=1}^{m_2} \theta_{jk}T_j(x;G_1) T_k(y;G_2) \,-\,\log Z_{\bm \te}\,. \eeq It is important to note that the set of LP-basis functions $\{T_j(x;G_1)\}$ and $\{T_k(y;G_2)\}$ are specially designed for the parametric marginals $G_1$ and $G_2$, obeying the following weighted orthonormality conditions: \beas \sum\nolimits_x g_1(x) T_j(x;G_1)=0,&and&\sum\nolimits_x g_1(x) T_j(x;G_1)T_k(x;G_1)=\delta_{jk}; \\ \sum\nolimits_x g_2(x)T_j(x;G_2)=0,&and&\sum\nolimits_x g_2(x) T_j(x;G_2)T_k(x;G_2)=\delta_{jk}. \eeas We call them gLP-basis, to distinguish them from the earlier empirical LP-polynomial systems $\{T_j(x;\wtF_X)\}$ and $\{T_k(y;\wtF_Y)\}$; see Appendix (ref).
Step 3. Exploratory goodness-of-fit. The estimated log-bilinear LP-copula is \beq \whcop_{X,Y}(u,v)\,=\,\exp\Big\{ 0.287S_1(u;G_1)S_1(v;G_2)-0.043 \Big\}. \eeq There are three important conclusions that can be drawn from this non-uniform copula density estimate: (i) Goodness of fit diagnostic: the independence model (product of parametric marginals) $g_\perp(x,y)=g_1(x) g_2(y)$ is not adequate for the data. (ii) Nature of discrepancy: the presence of significant $\hat\te_{11}=0.287$ in the model (ref) implies that the tentative independence model should be updated by incorporating the strong (positive) `linear' correlation between $X$ and $Y$. (iii) Nonparametric repair: how to update the initial $g_\perp(x,y)$ to construct a “better” model? Eq. (ref) gives the general updating rule, which simply says: copula provides the necessary bivariate-correction function to reduce the `gap' between the starting misspecified model $g_\perp(x,y)$ and the true unknown distribution $p_{X,Y}(x,y)$. (iv) In contrast to unsmoothed empirical multilinear copulas genest2013b, our method produces smoothed and compactly parametrizable $\whcop(u,v)$ for discrete data.
\vskip.3em Step 4. LP-smoothed probability estimation. The bottom panel Fig. (ref) shows the final smooth probability estimate $\hp_{X,Y}(x,y)$, computed by substituting (ref) into (ref). Also compare Tables (ref) and (ref) of Appendix (ref).
Mutual information (MI) is a fundamental quantity in Statistics and Machine Learning, with wide-ranging applications from neuroscience to physics to biology. For continuous random variables $(X,Y)$, mutual information is defined as \beq \MI(X,Y)\,=\,\iint f_{X,Y}(x,y) \log \dfrac{ f_{X,Y}(x,y) }{f_X(x) f_Y(y)} \dd x \dd y. \eeq Among non-parametric MI estimators, $k$-nearest-neighbor and kernel-density-based methods moon1995MI,kraskov2004estimating, zeng2018jackknife are undoubtedly the most popular ones. Here we are concerned with a slightly general problem of developing a flexible MI estimation algorithm that is: (D1) applicable for mixed\footnote{Reliably estimating MI for mixed case is notoriously challenging task gao2017estimating.} $(X,Y)$; (D2) robust in the presence of noise; and, (D3) invariant under monotone transformations\footnote{This is essential to make the analysis less sensitive to various types of data preprocessing, which is done routinely in applications like bioinformatics, astronomy, and neuroscience.}. To achieve this goal, we start by rewriting MI (ref) using copula: \beq \MI(X,Y)\,=\,\int_{[0,1]^2} \cop_{X,Y}(u,v) \log \cop_{X,Y}(u,v) \dd u \dd v. \eeq
The next theorem presents an elegant closed-form expression for MI in terms of LP-copula parameters, which allows a fast and efficient estimation algorithm. \vskip1em
{\bf Proof}.\, Express mutual information as: \[ \MI_{\teb}(X,Y)\,=\,\Ex_{X,Y}\big[ \log \cop_{\teb}\big]\,=\,\sum_j\sum_k \te_{jk} \Ex_{X,Y}\big[ S_j(U;X) S_k(V;Y)\big]\,-\,\log Z_{\teb}\,.\] The first equality follows from (ref) and the second one from (ref). Complete the proof by replacing by $\LP_{jk}$ by $\Ex[ S_j(U;X) S_k(V;Y)]$ by virtue of (ref). As a practical consequence, we have the following efficient and direct MI-estimator, satisfying D1-D3: \beq \widehat{\rm MI}_{\teb}(X,Y)\,=\,\mathop{\sum\sum}_{j,k>0} \hte_{jk}{\tLP}_{jk}\, -\, \log Z_{\widehat \teb}\,. \eeq \vskip.2em {\bf Bootstrap inference.} Bootstrap provides a convenient way to estimate the standard error of the estimate (ref). Perform bootstrap sampling, i.e., sample $n$ pairs of $(x_i,y_i)$ with replacement and compute $\widehat{\rm MI}$. Repeat the process, say, $B=500$ times to get the sampling distribution of the statistic. Finally, return the standard error of the bootstrap sampling distribution along with 95% percentile-confidence interval.
Continuous $(X,Y)$ example. Consider the kidney fitness data, discussed in Example (ref). The LP-copula-based (using $m=4$) method yields: $\hMI=0.230 \,\,(\pm 0.021)$. To understand how precise is the estimate, we have reported the bootstrap standard error in parentheses.
Given $n$ independent samples from an $I\times J$ contingency table, the $G^2$-test of goodness-of-fit, also known as the log-likelihood ratio test\footnote{In 1935, Samuel Wilks introduced log-likelihood ratio test as an alternative to Pearson’s chi-square test. In our notation, Pearson proposed $\int \cop^2_{X,Y}(u,v)$ and Wilks proposed $2\times \int \cop_{X,Y}(u,v)\log \cop_{X,Y}$---both are conceptually equivalent: measuring how much the copula density deviates from the uniformity.}, is defined as \beq G^2(X,Y)\,=\,2n\sum_{k=1}^I\sum_{l=1}^J \,\tp_{kl} \log \dfrac{\tp_{kl}}{\tp_{k+}\tp_{+l}}, \eeq which under the null hypothesis of independence has asymptotic $\chi^2_{(I-1)(J-1)}$ distribution. From (ref) one can immediately conclude the following.
The problem arises when we try to apply $G^2$-test for large sparse tables, and it is not hard to see why: the adequacy of asymptotic $\chi^2_{(I-1)(J-1)}$ distribution depends both on the sample size $n$ and the number of cells $p = IJ$. koehler1986goodness showed that the approximation completely breaks down when $n/p<5$, leading to erroneous statistical inference due to significant loss of power; see Appendix (ref). The following example demonstrates this.
\vskip.3em
We consider the two-sample feature selection problem where $Y$ is a binary response variable, and $X$ is a predictor variable that can be either discrete or continuous.
A few remarks on the interpretation of the above formula:
$\bullet$ Distributional effect-size: Bearing in mind Eqs. ((ref), (ref)), note that $d(v;X,X|Y=1)$ compares two densities: $f_{X|Y=1}(x)$ with $f_X(x)$, thereby capturing the distributional difference. This has advantages over traditional two-sample feature importance statistic (e.g., Student's t or Wilcoxon statistic) that can only measure differences in location or mean.
\vskip.3em $\bullet$ Explainability: The estimated $\whd(v;X,X|Y=1)$ involves three significant LP-components of $X$; the presence of the 1st order `linear' $S_1(v;X)$ indicates location-difference; the 2nd order `quadratic' $S_2(v;X)$ indicates scale-difference; and 4th order `quartic' $S_4(v;X)$ indicates the presence of tail-difference in the two RBC-distributions. In addition, the negative sign of the linear effect $\hte_{11}=-0.95$ implies reduced mean level of RBC in the ckd-population. Medically, this makes complete sense, since a dysfunctional kidney cannot produce enough Erythropoietin (EPO) hormone, which causes the RBC to drop.
We describe a new copula-based nonparametric logistic regression model. The key result is given by the following theorem, which provides a first-principle derivation of a robust nonlinear generalization of the classical linear logistic regression model.
{\bf Proof}. The proof consists of four main steps.
Step 1. To begin with, notice that for $Y$ binary and $X$ continuous, the general log-bilinear LP-copula density function (ref) reduces to the following form: \beq \cop_{\teb}(u,v;X,Y) = \dfrac{1}{Z_\te} \exp\Big\{ \sum_j \te_{j1} S_j(u;X) S_1(v;Y)\Big\}, \eeq since we can construct at most $2-1=1$ LP-basis function for binary $Y$.
Step 2. Apply copula-based Bayes Theorem (Eq. (ref)) and express the conditional comparison densities as follows: \beq d_1(x)\,\equiv\, d(F_X(x);X, X|Y=1)\,=\,\dfrac{\mu(x)}{\mu}, \eeq and also, \beq d_0(x)\,\equiv\, d(F_X(x);X, X|Y=0)\,=\,\dfrac{1-\mu(x)}{1-\mu}. \eeq Taking logarithm of the ratio of (ref) and (ref), we get the following important identity: \beq \log \left( \dfrac{\mu(x)}{1-\mu(x)}\right)\,=\,\log \left( \dfrac{\mu}{1-\mu}\right)\,+\,\log d_1(x)\,-\,\log d_0(x). \eeq
Step 3. From (ref), one can deduce the following orthonormal expansion of maxent-conditional copula slices $\log d_1$ and $\log d_0$: \bea \log d_1(x)&=& \sum_j\big( \te_{j1} T_1(1;F_Y) \big)\, T_j(x;F_X) \,-\,\log Z_\te \\ \log d_0(x)&=& \sum_j\big( \te_{j1} T_1(0;F_Y) \big) \, T_j(x;F_X) \,-\,\log Z_\te \eea
Step 4. Substituting (ref) and (ref) into (ref) we get: \[ \log \left( \dfrac{\mu(x)}{1-\mu(x)}\right)\,=\,\log \left( \dfrac{\mu}{1-\mu}\right)\,+\,\sum_j \Big\{ \frac{\te_{j1}}{\sqrt{\mu(1-\mu)}}\Big\} T_j(x;F_X), \] since for binary $Y$ we have (see Appendix (ref)): \[ T_1(1;F_Y)-T_1(0;F_Y)\,=\,\dfrac{1-\mu}{\sqrt{\mu(1-\mu)}}\,+\,\dfrac{\mu}{\sqrt{\mu(1-\mu)}}\,=\,\dfrac{1}{\sqrt{\mu(1-\mu)}}. \] Substitute $\al_0=\logit(\mu)$ and $\al_j=\frac{\te_{j1}}{\sqrt{\mu(1-\mu)}}$ to complete the proof. \qed
Generalize the univariate copula-logistic regression model (ref) to the high-dimensional case as follows: \beq \logit\big( \mu(x) \big)\,=\, \al_0 + \sum_{j=1}^p h_j(x_j). \eeq Nonparametrically approximate the unknown smooth $h_j$'s by LP-polynomial series of $X_j$ \beq h_j(x_j)\,=\,\sum_{k=1}^m \al_{jk} T_k(x_j;F_{X_j}), for\,j=1,\ldots,p. \eeq
\vskip.1em
This paper makes the following contributions: (i) we introduce modern statistical theory and principles for maximum entropy copula density estimation that is self-adaptive for the mixed(X,Y)---described in Section (ref). (ii) Our general copula-based formulation provides a unifying framework of data analysis from which one can systematically distill a number of fundamental statistical methods by revealing some completely unexpected connections between them. The importance of our theory in applied and theoretical statistics is highlighted in Section (ref), taking examples from different sub-fields of statistics: Log-linear analysis of categorical data, logratio biplot, smoothing large sparse contingency tables, mutual information, smooth-$G^2$ statistic, feature selection, and copula-based logistic regression. We hope that this new perspective on copula modeling will offer more effective ways of developing united statistical algorithms for mixed-($X,Y$). \vskip.4em
This paper is dedicated to the birth centenary of {\bf E. T. Jaynes} (1922--1998), the originator of the maximum entropy principle.
I also like to dedicate this paper to the memory of {\bf Leo Goodman} (1928--2020)---a transformative legend of categorical data analysis, who passed away on December 22, 2020, at the age of 92 due to COVID-19.
This paper is inspired in part by the author's intention to demonstrate how these two modeling philosophies can be connected and united in some ways. This is achieved by employing a new nonparametric representation theory of generalized copula density.
\vskip1em