Assessing Copula Models for Mixed Continuous-Ordinal Variables

For a pair of continuous variables, normal score plots are often used to visualize the copula for the variables and check for deviations from Gaussian dependence. With the ordinal variable $X$, $\{u_{iX}^{+}\}$ is a set of discrete values. To visualize the relationship between $X$ and $Y$, we use an appropriate latent continuous variable associated with $X$.

For the normal score plots, we assume that the ordinal variable $X$ is generated from a latent standard normal variable $Z$ with estimated cutpoints $\zeta_{j}=\Phi^{-1}\left((n_{1}+\cdots+n_{j})/{n}\right)$ for $j=1,\ldots,k-1$, where $\Phi$ and $\Phi^{-1}$ are the CDF and quantile function of the standard normal distribution. Let $\zeta_{0}=-\infty$ and $\zeta_{k}=\infty$. We would like the latent variable $Z$ to be generated in such a way that the correlation between $Z$ and $Y$ is close to the correlation between $X$ and $Y$ in order for the visualization to preserve the strength of correlation from the original variables.

One measure of the association between an ordinal and a continuous variable is the polyserial correlation (see Olsson et al., (1982) and Section 2.12.7 of Joe, (2014)). The polyserial correlation between $X$ and $Y$ is defined as

where $\phi$ is the probability density function (PDF) of the standard normal distribution.

Let $\mathcal{S}_{j}=\{i:x_{i}=j\}$, $j=1,\ldots,k$, be the set of indices of $X$ taking value $j$. We propose to generate the normal scores of the latent variable $Z$ for the ordinal variable $X$ in the following steps.

1. 1.

   Compute the polyserial correlation $\rho_{N}$ between $X$ and $Y$. For $j=1,\ldots,k$, perform steps 2 to 5 on $\mathcal{S}_{j}$:

2. 2.

   For each $i\in\mathcal{S}_{j}$, independently generate $\omega_{i}\sim U(0,1)$. Let

   $$u_{iW}=\Phi\left(\sqrt{1-\rho_{N}^{2}}\Phi^{-1}(\omega_{i})+\rho_{N}\Phi^{-1}(u_{iY})\right),$$

   assuming a bivariate Gaussian copula with correlation $\rho_{N}$. Note that the generated $u_{iW}$ are in the range of $[0,1]$. Let $\boldsymbol{u}_{W,j}=\{u_{iW}:i\in\mathcal{S}_{j}\}$ be a vector of length $n_{j}$.

3. 3.

   For each $i\in\mathcal{S}_{j}$, independently generate $\psi_{i}\sim U(n^{-1}\sum_{t=1}^{j-1}n_{t},n^{-1}\sum_{t=1}^{j}n_{t})$ for $j>1$, or $\psi_{i}\sim U(0,n_{1}/n)$ for $j=1$. Let $\boldsymbol{\psi}_{j}=\{\psi_{i}:i\in\mathcal{S}_{j}\}$ be a vector of length $n_{j}$.

4. 4.

   For each $u_{iW}$ with $i\in\mathcal{S}_{j}$, find its rank in $\boldsymbol{u}_{W,j}$. Replace the value of $u_{iW}$ by $\psi_{\ell}$ that has the same rank in $\boldsymbol{\psi}_{j}$, and denote this value by $u_{iZ}$. Let $\boldsymbol{u}_{Z,j}=\{u_{iZ}:i\in\mathcal{S}_{j}\}$.

5. 5.

   Let $z_{i}=\Phi^{-1}(u_{iZ})$ for each $i\in\mathcal{S}_{j}$.

The $\{u_{iW}\}$ generated in step 2 does not have a uniform distribution. The additional steps 3 and 4 match the ranks of $u_{iW}$ with a random vector generated from the uniform distribution. This procedure ensures that the generated normal scores $\{z_{i}\}$ fall into the desired bins separated by the cutpoints $\boldsymbol{\zeta}$ and marginally follow the standard normal distribution. Within each bin, $\{z_{i}\}$ and $\{\Phi^{-1}(u_{iY})\}$ also have similar correlation to $\rho_{N}$. The two sets of normal scores, $\{z_{i}:i=1,\ldots,n\}$ and $\{\Phi^{-1}(u_{iY}):i=1,\ldots,n\}$, can then be plotted against each other to decide on suitable parametric copula families. If a bivariate copula model can fit $\{(u_{iZ},u_{iY})\}$ well, then it can also potentially provide adequate fits to $\{(x_{i},y_{i})\}$.

Let $F_{n,X}(j)={\widehat{F}}_{X}(j)$ for $j=1,\dots,k$. As $n\to\infty$, $F_{n,X}(j)=n^{-1}\sum_{t=1}^{j}n_{t}\to F_{X}(j)$ for $j=1,\dots,k$. Therefore, each element $\psi_{i}$ of $\boldsymbol{\psi}_{j}$ in Step 3 of the proposed algorithm above follows $U\left(F_{n,X}(j-1),F_{n,X}(j)\right)$ and satisfies the condition (1) as $n\to\infty$. Similarly, each element of $\boldsymbol{u}_{Z,j}$ also satisfies the condition (1) as $n\to\infty$, since $\boldsymbol{u}_{Z,j}$ is a permutation of $\boldsymbol{\psi}_{j}$. However, the correlation for $\{(u_{iZ},u_{iY})\}$ is stronger than the correlation for $\{(\psi_{i},u_{iY})\}$. The bivariate empirical copulas of the two sets $\{(\psi_{i},u_{iY}):i=1,\dots,n\}$ and $\{(u_{iZ},u_{iY}):i=1,\dots,n\}$ evaluated at $(u_{x},u_{y})$ are the same for empirical CDF values $u_{x}=F_{n,X}(j)$ and arbitrary $u_{y}\in(0,1)$.

With $F_{XY}(x,y)=C_{XY}(F_{X}(x),F_{Y}(y))$ for a copula, the conditional distribution of the continuous variable $Y$ given $X=x$, where $x$ is a possible value of the ordinal variable $X$, is

If ${\widehat{C}}_{XY}=C(\cdot;{\widehat{\boldsymbol{\theta}}})$ is an estimate from a parametric family, let

be the estimated conditional distribution. Given a quantile level $q$, the conditional quantile ${\widehat{F}}^{-1}_{Y|X}(q|x)$ is obtained as ${\widehat{F}}^{-1}_{Y}(v)$ where $v$ is the root to the equation

To assess the fit of ${\widehat{C}}_{XY}$, conditional Q-Q plots can be generated based on ${\widehat{F}}^{-1}_{Y|X}(\cdot|j)$ for $j=1,\ldots,k$. For the $j$th conditional Q-Q plot, the quantiles ${\widehat{F}}_{Y|X}^{-1}(\cdot|j)$ are plotted against the quantiles of the empirical distribution of $\{y_{i}:x_{i}=j\}$, that is, the model-based quantiles ${\widehat{F}}_{Y|X}^{-1}\left((m-0.5)/{n_{j}}\big{|}j\right)$ for $m=1,\ldots,n_{j}$ are plotted against the sorted values of $\{y_{i}:x_{i}=j\}$. If the points in each conditional Q-Q plot lie closely along the $45^{\circ}$ diagonal line, it indicates that the copula estimate fits the sample well.

The use of these conditional Q-Q plots for copula diagnostics is illustrated in Section 6 for simulated datasets and Section 7 for a real dataset.

There are two classes of models for a mix of continuous and ordinal variables depending on the order of conditioning.

For the first class, a multinomial distribution is assumed for the ordinal variable and the conditional distribution of the continuous variable given the ordinal variable can be either Gaussian with different mean and variance parameters, or a more general location-scale family (Little and Schluchter, (1985) and Krzanowski, (1993)). For the second class, the continuous variable is transformed to have a parametric distribution such as $N(0,1)$ and the ordinal variable conditioned on the continuous variable is an ordinal probit or logit model; references are Cox, (1972) and Krzanowski, (1993). We introduce the notation for these two classes of models with an ordinal variable $X$ and a continuous variable $Y$.

For the first class of models, a location-scale mixture model is considered. The conditional distribution $F_{Y|X}(\cdot|x)$ is a general location-scale model with

where $p$ is a density function on the real line. The ordinal variable $X$ has a probability distribution representing the mixture components. If $Y$ can only take positive values, it can be transformed via a logarithm to get a location-scale model.

For the second class of models, the conditional distribution $F_{X|Y}(\cdot|y)$ is

where $F$ is a standard normal or logistic CDF, $a$ is a slope parameter, and $b$ is an offset depending on the ordinal category. The continuous variable $Y$ is assumed to have a unimodal distribution. If $Y\sim N(0,1)$, simpler calculation results can be obtained when $F$ is the standard normal CDF.

Let $H_{X,Y},G_{X,Y}$ be two bivariate distributions with densities $h_{X,Y},g_{X,Y}$ on the respective relevant measure spaces. For ordinal $X$ with values in $\{1,\ldots,k\}$ and absolutely continuous $Y$, the product measure comes from a counting measure and the Lebesgue measure on the real line. The non-negative KL divergence of $h_{X,Y}$ from $g_{X,Y}$ is

Let $C(\cdot;\theta)$ be a bivariate copula family with $C_{1|2}(u|v;\theta)={\partial C(u,v;\theta)/\partial v}$. If $H_{X,Y}(x,y;\theta)=C(F_{X}(x),F_{Y}(y);\theta)$ with conditional probability mass function $h_{X|Y}(x|y;\theta)=C_{1|2}(F_{X}(x)|F_{Y}(y);\theta)-C_{1|2}(F_{X}(x^{-})|F_{Y}(y);\theta)$, then the joint density is $h_{X,Y}(x,y;\theta)=f_{Y}(y)\,h_{X|Y}(x|y;\theta)$ with $f_{Y}(y)=\mathrm{d}F_{Y}(y)/\mathrm{d}y$. In the setting of a mixture model that specifies $f_{X}$ and $g_{Y|X}$, $f_{Y}(y)=\sum_{i=1}^{k}f_{X}(i)\,g_{Y|X}(y|i)$.

The copula with parameter estimate

is the member in the family that is the closest to $g_{X,Y}$. The value $KL(h_{X,Y}(\cdot;{\hat{\theta}}),g_{X,Y})$ is considered as the KL divergence of the family $\{h_{X,Y}(\cdot;\theta)\}$ from $g_{X,Y}$.

In the results below, we consider different parametric copula families denoted by $C^{(m)}(\cdot;\theta^{(m)})$ with corresponding densities $h^{(m)}_{X,Y}(\cdot;\theta^{(m)})$. Given $g_{X,Y}$, this leads to $KL(h^{(m)}_{X,Y}(\cdot;{\hat{\theta}}^{(m)}),g_{X,Y})$ for model $m$. A copula family that has a smaller value of KL divergence is considered to have a better approximation of $g_{X,Y}$.

Next are steps to find a sequence of $C^{(m)}$ that leads to smaller values of the KL divergence of a family from a given $g_{X,Y}$. Because many simple parametric copula families only have positive dependence, we assume that $X$ and $Y$ have been oriented to have positive dependence. We also assume that $Y$ has a distribution that is close to unimodal (otherwise one might use a non-copula-based model for $(X,Y)$). Unless an additional reference is given, properties of the listed copula families can be found in Chapter 4 of Joe, (2014).

1. 1.

   Compute the KL divergence of the Gaussian copula family as a baseline.

2. 2.

   For copula family candidates with more dependence in one joint tail, compute the KL divergence of the Gumbel and survival Gumbel copula families with asymmetric dependence in the joint upper and lower tails. If one of these two families leads to smaller KL divergence, then try copula families with even more tail asymmetry in the same direction.

3. 3.

   For copula family candidates with less dependence in both joint tails (compared to Gaussian), compute the KL divergence of the Frank and Plackett copula families with reflection symmetric tail quadrant independence. If one of these two copula families leads to smaller KL divergence, then compute the KL divergence for the BB8 and BB10 ¹¹1Note that Kadhem and Nikoloulopoulos, (2021) show that there are parameters that can lead to non-convex contours of copula densities with N(0,1) margins for the BB10 copula. families with reflection tail asymmetry and their survival counterparts.

4. 4.

   For copula family candidates with more dependence in both joint tails, compute the KL divergence for the t copula family.

5. 5.

   If permutation asymmetry is possible in $g_{X,Y}$ (e.g., with heterogeneous components in a mixture model), compute the KL divergence for some copula families with permutation asymmetry; examples are skew normal (Yoshiba, (2018)) and asymmetric Gumbel (Section 4.15 of Joe, (2014)) copula families.

The concept of tail order (Hua and Joe, (2011)), covering copula families with different strengths of dependence in the joint lower and upper tails, is used in the above steps. The main distinctions are intermediate tail dependence (such as Gaussian copula with positive dependence), strong tail dependence (more dependence in the joint tail than Gaussian), and tail quadrant independence (less dependence in the joint tail than Gaussian). The concept of tail order may be less important when one variable is ordinal, because there are no observations in the joint upper or lower tail. However, copula families with tail asymmetries are still important to provide more flexibility when finding approximations to a given $g_{X,Y}$.

In this section, we consider a variety of concrete examples of the probability models in Section 4.1 which lead to $g_{X,Y}$ in (3). In all of these examples, the marginal density $f_{Y}$ of $g_{X,Y}$ is close to unimodal. Otherwise, $\text{cor}(X,Y)$ could be large and “clusters" might be seen in bivariate scatterplots. After examining a large number of cases, we find that a KL divergence value less than 0.003 usually indicates a good approximation from a copula family, and a KL divergence value greater than 0.01 usually indicates a poor approximation when using the conditional Q-Q diagnostic plots in Section 3.2.

The tables in Section 4.3.1 summarize some representative examples to illustrate what happens for (a) mixture models with different separation of location parameters and common versus varying scale parameters, and (b) conditional probit or logit models. Section 4.3.2 contains conclusions drawn from these examples.

Suppose there are $M$ parametric bivariate copula models to consider as candidate families. For a parametric bivariate copula model $C^{(m)}$ with parameter vector $\boldsymbol{\theta}^{(m)}$ and $C_{1|2}^{(m)}(u|v;\boldsymbol{\theta}^{(m)})=\partial C^{(m)}(u,v;\boldsymbol{\theta}^{(m)})/\partial v$, its log-likelihood function is

The maximum likelihood estimator is $\widehat{\boldsymbol{\theta}}^{(m)}_{n}=\operatorname*{arg\,max}_{\boldsymbol{\theta}^{(m)}}\mathcal{L}_{m}(\boldsymbol{\theta}^{(m)})$ for a sample of size $n$. As $n\to\infty$, $\widehat{\boldsymbol{\theta}}^{(m)}_{n}$ converges in probability to the value $\widetilde{\boldsymbol{\theta}}^{(m)}$ that minimizes the KL divergence for family $C^{(m)}(\cdot;\boldsymbol{\theta}^{(m)})$. Parametric copula models are then compared using model selection criteria such as AIC and BIC. Models with smaller values of AIC or BIC are usually considered to fit the data better.

Since (2) only involves the copula CDF but not the copula density, we consider the empirical beta copula (Segers et al., (2017)) as a nonparametric alternative fitted to pairs $(u_{iZ},u_{iY})$, $i=1,\ldots,n$. The empirical beta copula density performs less well than the KDE copula estimator in Nagler, (2018), even after some averaging over distinct subsamples. However, the empirical beta copula CDF has the advantage of being a proper copula, while the KDE approach only leads to a distribution that is approximately a copula.

Let $\boldsymbol{u}_{Z}=\left\{u_{iZ}:i=1,\dots,n\right\}$ and $\boldsymbol{u}_{Y}=\left\{u_{iY}:i=1,\dots,n\right\}$. For $\{(u_{iZ},u_{iY})\}$, the bivariate empirical beta copula CDF is given by

where $F_{B}(\cdot;\alpha,\beta)$ is the CDF of the $\text{Beta}(\alpha,\beta)$ distribution, $R_{iZ}^{(n)}$ is the rank of $u_{iZ}$ in $\boldsymbol{u}_{Z}$, and $R_{iY}^{(n)}$ is the rank of $u_{iY}$ in $\boldsymbol{u}_{Y}$. Note that (4) is a continuous differentiable function on $[0,1]^{2}$. The $r$th order statistic $U_{(r:n)}$ in a random sample of size $n$ generated from $U(0,1)$ follows a $\text{Beta}(r,n+1-r)$ distribution, which leads to the two beta distributions in (4).

Assuming a consistent estimator for $F_{Y}$, the consistency of the empirical beta copula with the latent vector $\boldsymbol{u}_{Z}$ is shown in Section 5.3 below.