Goodman and Kruskal’s Gamma Coefficient for Ordinalized Bivariate Normal Distributions

We consider a bivariate normal distribution with linear correlation ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} whose random components are discretized according to two assigned sets of thresholds. On the resulting bivariate ordinal random variable, one can compute Goodman and Kruskal’s gamma coefficient, γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}, which is a common measure of ordinal association. Given the known analytical monotonic relationship between Pearson’s ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} and Kendall’s rank correlation τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} for the bivariate normal distribution, and since in the continuous case, Kendall’s τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} coincides with Goodman and Kruskal’s γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}, the change of this association measure before and after discretization is worth studying. We consider several experimental settings obtained by varying the two sets of thresholds, or, equivalently, the marginal distributions of the final ordinal variables. This study, confirming previous findings, shows how the gamma coefficient is always larger in absolute value than Kendall’s rank correlation; this discrepancy lessens when the number of categories increases or, given the same number of categories, when using equally probable categories. Based on these results, a proposal is suggested to build a bivariate ordinal variable with assigned margins and Goodman and Kruskal’s γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} by ordinalizing a bivariate normal distribution. Illustrative examples employing artificial and real data are provided.


Introduction
The use of the multivariate normal distribution as a latent construct for modelling observed correlated or associated discrete or ordinal variables can be dated back to the seminal book by Lazarsfeld and Henry (1968) and to the later work by Muthen (1983), in the context of structural equation modeling. Latent variable modelling has since then gradually become an integral part of mainstream statistics and is currently used for a multitude of applications in different subject areas (Beaujean 2014).
It is also an indisputable fact that the ability to simulate artificial data resembling the main features of some observed dataset or following the specifications of a study design is necessary when comparing and investigating the behaviour of statistical procedures and exploring their robustness; such features or specifications can often be conveniently summarized by the empirical marginal distributions and pairwise measures of correlation or association.
In this work, we will focus our attention on the relationship between a concordance measure of the bivariate normal latent variable (Kendall's rank correlation, directly related to the more popular Pearson's correlation) and a corresponding association measure of the ordinalized variable (Goodman and Kruskal's gamma). On the one hand, we will analyze how the correlation of the latent variable and the marginal distributions of the final ordinal variables affect their correlation/association: we can call it the "direct problem". Specifically, we will investigate the effect of Correspondence should be made to Alessandro Barbiero, Department of Economics, Management, and Quantitative Methods, Università degli Studi di Milano, Via Conservatorio, 7, 20122 Milan, Italy. Email: alessandro.barbiero@unimi.it 906 PSYCHOMETRIKA the number of categories and the probability distribution (uniform/symmetrical/skewed) on the distortion of the association measure. On the other hand, we will provide a procedure that allows construction and simulation of a random vector of discrete variables with assigned marginal distributions and association by discretizing a bivariate normal random vector with a pairwise correlation that is able to induce the desired association among the ordinal variables (the "inverse problem", i.e., finding the parameters of the continuous model that produce the (parameters of the) observed ordinal variables).
The paper is structured as follows. In the next section, we will recall the definition of Goodman and Kruskal's gamma coefficient and summarize its main properties, also in comparison with alternative measures of ordinal association. In Sect. 3, we will analyze the change in magnitude of association before and after ordinalization of a bivariate normal variable, measured by Kendall's tau and Goodman and Kruskal's gamma, respectively, under a wide array of experimental conditions. Section 4 suggests and examines a procedure for building a bivariate ordinal variable with assigned marginal distributions and association. Section 5 proposes a possible application to inference of the algorithm of Sect. 4; Sect. 6 presents an illustrative example with real data. The last section is devoted to some final remarks and possible research prospects.

Measures of Ordinal Association: Goodman and Kruskal's Gamma
Measures of ordinal association between two variables use the property that the categories of ordinal variables have a natural order; but the way in which these measures use this information may differ considerably. Some measures are based on the intuitive notion that ordinal association should have an interpretation analogous to that of association for metric variables, say, Pearson's correlation; in this group, we find Spearman's ρ, Kendall's τ , polychoric correlation, Somer's d, and Goodman-Kruskal's γ . These coefficients allow us to make statements of the general form "if scores on X increase, then most probably, scores on Y will increase" (Kampen and Swyngedouw 2000). Other measures of ordinal association have a background in information theory and have interpretations in terms of stochastic entropy (Bryson and Phillips 1975;Laird 1979;Gilula et al. 1988). Henceforth, we will focus on the former family of measures.
Let us consider a pair of ordinal variables (X, Y ), with H and K ordered categories, respectively, and introduce the concept of concordance/discordance by considering two independent realizations (x i , y i ) and (x j , y j ). (x i , y i ) and (x j , y j ) will be said to be concordant if x i < x j and y i < y j , or if x i > x j and y i > y j . Conversely, (x i , y i ) and (x j , y j ) will be said to be discordant if x i < x j and y i > y j , or if x i > x j and y i < y j . Thus, one can define the probability of concordance as and similarly the probability of discordance as The gamma coefficient belongs to a larger family of ordinal correlation measures (see, e.g., Woods 2009, for an exhaustive account); it is defined as the following ratio (see, e.g., Agresti 2010, pp.186-187): Π c and Π d can be conveniently expressed in terms of the joint probabilities p i j = P(X = x i , Y = y j ): For a bivariate continuous random variable, γ can be still computed through Eq. (1): in this case, since the probability of concordance and the probability of discordance sum up to 1 (there is no probability of tied values for either X or Y ), then γ = Π c − Π d = 2Π c − 1, and it coincides with Kendall's rank correlation τ (Kendall 1945), simply defined as the difference between the probability of concordance and the probability of discordance. We recall that for two continuous random variables X and Y , Kendall's τ depends only on their unique copula C (representing the dependence structure of the bivariate random vector (X, Y ); see McNeil et al., 2005, chapter 5), specifically, and not on the marginal distributions of X and Y . τ ranges between −1 and +1; it attains the upper bound if and only if X and Y are comonotonic, whereas it attains the lower bound if and only if X and Y are countermonotonic. Like Pearson's correlation and other correlation measures such as Spearman's rho and the aforementioned Kendall's tau,+1]. The values −1, 0, and +1 are attained when Π c = 0, Π c = Π d , Π d = 0, respectively. There are however some potential problems with the gamma coefficient, which were immediately recognized by Goodman and Kruskal (1954): • γ is unstable over various "cutting points", that is to say, γ tends to increase as the categories of a contingency table are collapsed, since γ gives no consideration to tied pairs, as can be seen from Eq. (1), and the number of tied pairs increases as the table is collapsed; • γ also usually yields greater association values than other measures of ordinal association, as it does not consider any of the tied pairs; • Finally, γ is a weakly monotonic measure of ordinal association, i.e., it reaches +1 under a variety of cell frequency configurations, not only in case of strict perfect association (Berry et al. 2018). For example, the bivariate distribution of Table 1, where all the joint probabilities are zero except those labelled with ×, presents a value of γ equal to +1, although the relationship between X and Y is not perfectly monotonic (Kendall's τ and Spearman's ρ would be strictly smaller than 1).  With the objective of overcoming these pitfalls, some modifications to the gamma coefficient or alternative measures of ordinal association have been proposed (Rousson 2007;Kvålseth 2017Kvålseth , 2018. A remark on the possible extension of the gamma coefficient to higher dimensions can be made. Let us consider the d-variate ordinal random vector X X X = (X 1 , X 2 , . . . , X d ) , with d > 2; then, a γ association matrix remains defined as ...,d; j=1,...,d with γ i j = γ (X i , X j ). This matrix however is not necessarily positive semidefinite, as happens for Spearman's or Kendall's rank correlation matrices, which are both valid correlation matrices (see McNeil et al. 2005, p.207). An easy counterexample is provided as follows. Consider the trivariate distribution of Table 2. The values of γ for the bivariate distributions (X, Y ), (X, Z ), and (Y, Z ) are γ xy = 1/3, γ xz = 1, and γ yz = −7/9, respectively. The corresponding γ association matrix is: its three eigenvalues are 2.119, 1.322, and −0.441; since one of them is negative, the matrix is clearly not positive semidefinite.
With n pairs of observations constituting a multinomial sample, the sample analog of γ is given byγ = (C − D)/(C + D), with C the number of concordant pairs and D the number of discordant pairs. It is possible to prove that √ n(γ − γ ) is asymptotically normal with mean zero and variance given by i j ) (see, e.g., Agresti 2010). In practice, replacing p i j , c , and d with their sample values in σ 2 yields the ML estimatê σ 2 of σ 2 . The term S E =σ / √ n is an estimated standard error forγ , and Wald confidence interval for γ is (γ ± z 1−α/2 S E). Rosenthal (1966); Gans and Robertson (1981) showed thatγ has a tendency to converge slowly to normality and to have distributional irregularity, bias, and skewness problems, especially when the true absolute value of γ is large. O'Gorman and Woolson (1988); Carr et al. (1989) pointed out that better convergence occurs using the Fisher-type transformξ = 1 2 log[(1+γ )/(1− γ )], whose asymptotic variance equals the asymptotic variance ofγ multiplied by (1 − γ 2 ) −2 . A confidence interval can be constructed for ξ and then inverted to one for γ , by using the inverse transformationγ = (e 2ξ − 1)/(e 2ξ + 1) .

Analysis of the Relationship Between ρ and γ for Ordinalized Bivariate Normal Distribution
Discretization of continuous variables is commonly encountered in practice. Based on observed nominal age, income, temperature, and depression score, one can derive ordinal variables such as young-middle-old age, low-medium-high income, cold-cool-average-hot temperature, and no-mild-moderate-severe depression. Discretization is usually avoided by statisticians for intuitive and valid reasons, the most prominent of which is the associated power and information loss. Some problems related to identification of statistical models obtained by discretization have been recently raised by Grønneberg and Foldnes (2019). However, simplicity, better interpretability and comprehension of the effects of interest, and superiority of some categorical data measures such as odds ratio have been argued by proponents of discretization (Liu et al. 2002).
The objective of this section is the determination of association magnitude changes when the univariate components of a bivariate normal distribution are ordinalized, i.e., the range of each continuous component is divided into contiguous intervals, which are assigned an ordered category or a consecutive positive integer. A similar and extensive analysis has been conducted by Demirtas and Vardar-Acar (2017) to investigate the effects of discretization of continuous variables on the magnitude of linear correlation. If the underlying continuous distribution is bivariate normal, it can be proved that discretization always preserves the sign of linear correlation and more importantly leads to a reduction in magnitude. This result, which has been empirically observed in many simulation studies (see, e.g., Bollen and Barb 1981) and was claimed to hold only "in large samples" by Demirtas and Vardar-Acar (2017), is just a consequence of a previous, more general theoretical result, named Lancaster's theorem (Lancaster 1957), which states that the correlation of a bivariate normal cannot increase whatever transformations are applied to its univariate components (see also Mari and Kotz 2001, p.155, where it is reported as an "extremal property" of the bivariate normal distribution).
For a bivariate normal distribution with correlation coefficient ρ, the following relationship holds between ρ and Kendall's rank correlation τ : Equation (2) holds also for most elliptical distributions, e.g., for bivariate Student's t, and for most bivariate distributions whose dependence structure is described by an elliptical copula (McNeil et al. 2005). If both univariate margins are discretized/ordinalized through pre-specified thresholds, one can compute some measure of ordinal association based on the resulting bivariate ordinal variable, such as Goodman and Kruskal's γ : we know (see the previous section) that for continuous bivariate random variables, such as the bivariate normal, the definitions of γ and τ coincide. Thus, it makes sense to analyze the relationship between τ and γ as a measure of the change in association before and after ordinalization. Let us start from a very simple case, i.e., dichotomization of the margins of a standard bivariate normal random variable (Z 1 , Z 2 ), with zero marginal means, unit marginal variances, and linear correlation ρ. Consider the bivariate ordinal variable (X, Y ) obtained as follows: with x 1 < x 2 and y 1 < y 2 being two ordered categories. Therefore, the marginal probabilities for X and Y are P(X = x 1 ) = P(X = x 2 ) = 1/2 and analogously P(Y = y 1 ) = P(Y = y 2 ) = 1/2. Since we know that for a bivariate normal distribution P( , McNeil et al. 2005), then the joint probability of (X, Y ) is that displayed in the following table: Based on it, by using Eq.
(1), one can compute the gamma coefficient as a function of ρ: or, recalling (2), in terms of τ : whose graph is plotted in Fig. 1. This graph clearly shows how Goodman-Kruskal's γ for the dichotomized random variable is a strictly increasing and odd function of its analogue Kendall's τ for the bivariate normal distribution and, in absolute value, is always larger than or equal to τ , the equality holding when τ = 0 (i.e., when the two normal components are independent), when τ equals 1 (perfectly positively correlated normal random variables), or when τ equals −1 (perfectly negatively correlated normal random variables). It is worth noting that in this context, γ is equal to +1 (−1) only in case of a perfect monotonic (countermonotonic) relationship between X and Y . By choosing thresholds different from zero for the random components Z 1 and Z 2 , we can numerically obtain the value of γ corresponding to any value of τ . Similarly, one can discretize either or both continuous random components into more than two categories, applying different sets of thresholds (or, equivalently, assigning different marginal distributions to the final ordinal Graph of the gamma coefficient for a dichotomized bivariate normal random variable with Kendall's rank correlation τ , see Eq. (5). The dashed line is the 45 degrees line passing through the origin variables). In general, let us consider a bivariate standard normal distribution (Z 1 , Z 2 ) whose two continuous components are discretized into two ordinal variables X and Y according to the following scheme: where θ 1 < θ 2 < · · · < θ H = ∞ and η 1 < η 2 < · · · < η K = ∞ constitute two sets of thresholds. If we denote with F(i, j) the bivariate joint cumulative probability P( where τ is the joint cdf of a bivariate normal with standard components and rank correlation τ , and F −1 1 (F −1 2 ) is the generalized inverse of F 1 (F 2 ). The joint probabilities p i j = P(X = x i , Y = y j ) are then obtained as Kruskal's γ for an ordinalized bivariate normal random variable with Kendall's rank correlation τ ; the categories of the two ordinal variables are set both equal to an integer K and have uniform probabilities 1/K . The dashed line is the 45 degrees line passing through the origin and based on these p i j , the gamma coefficient can be computed using Eq. (1). In the following study, we will consider three different possible types of ordinal marginal distributions, namely, (i) uniform (i.e., equiprobable categories) (ii) symmetrical (non-uniform), based on normal scores, and (iii) asymmetrical (triangular). We will examine what kind of relationship τ and γ have under these three macro-settings by considering the same distribution for the two margins and varying the number of categories. Relevant code developed in the R environment (R Core Team 2020) is freely available as supplementary material. Figure 2 displays the graphs of the (τ, γ ) curve when the two marginal distributions have the same number of equiprobable categories H = K = 3, 4, 5, 10, 20. Table 3 reports the values of γ for τ ranging from −1 to +1, with step length of 0.1, and several values of the common number of categories of the two identical margins (H = K = 2, 3, 5, 7, 10, 20, 50, 100). Note that the (τ, γ ) curve, for a given H , is strictly increasing and always passes through the points (−1, −1), (0, 0), and (1, 1), and, as one can expect, by increasing H (i.e., when the ordinal distributions "resemble" a continuous one), it tends to get close to the 45-degree line passing through the origin. Note that γ is an odd function of τ , i.e., γ (τ ) = −γ (−τ ). For a given value of τ > 0, γ is a decreasing function of H ; for a given value of τ < 0, γ is an increasing function of K ; as K tends toward ∞, γ seems to converge, though slowly, to τ . It is worth observing, however, that even for a moderately large number of categories, there is a substantial difference between the values of τ and γ before and after ordinalization. When H = 100, there is still a relative difference of up to 2% between the values of τ and γ .
We now move to non-uniform symmetrical ordinal distributions; in particular, we consider a probability mass function that somewhat resembles the probability density function of the normal variable. For this aim, the (−4, +4) interval, whose probability for a standard normal variable is almost equal to 1, is divided into H equal-width adjacent intervals, (−4 + 8(i − 1)/H, −4 + 8i/H ), i = 1, . . . , H . A similar construction was employed by Becker   (1989) (who named it the "uniform cut-points (UCP) method" and employed it to study similarities between uniform association models and ordinalized bivariate normal distribution) and Ferrari and Barbiero (2012) (to study the effects of discretization of a bivariate normal distribution on the correlation coefficient). To the i-th interval, the i-th ordered category is associated, whose probability is thus equal to (−4 + 8i/H ) − (−4 + 8(i − 1)/H ), i = 1, . . . , H ; the residual probabilities of (−∞, −4) and (4, +∞) are assigned to the first and last categories, respectively. Thus, the probability distribution is always symmetrical about the central category (if H is odd) or categories (if H is even). The probabilities obtained according to this scheme are graphically displayed, for K = 5 and K = 10, in Fig. 3. Table 4 reports the values of Goodman and Kruskal's γ for τ ranging from −1 to +1, with a step length of 0.1, and several values of the common number of categories of the two identical symmetrical margins. Note that as in the case of equiprobable categories, γ is an odd function of τ . For a given value of τ > 0 and for K > 3, γ is a decreasing function of K ; for a given value of τ < 0 and for K > 3, γ is an increasing function of K (for K = 2, we obtain again a uniform distribution with two equally probable categories). As K tends toward ∞, γ converges to τ , but much more slowly than in the case of uniform margins. When K = 100, the relative difference between τ and γ in is not negligible at all, attaining 4.5%. This fact can be easily explained by just looking at the barplot in the right side of Fig. 3: among the K categories, only the central ones have appreciable probabilities, so the "actual" number of categories is (much) smaller than K , and then a larger value of K is required to obtain the same difference between τ and γ that is achieved in case of uniform distributions.       In the end, we consider strongly asymmetrical distributions, namely triangular distributions with H categories, where the probability of the i-th category is proportional to i (through some positive constant α). Since under this assumption, H i=1 p i = H i=1 αi = α H (H + 1)/2, then p i = 2i/[H (H + 1)]. Choosing this distribution as the marginal distribution of both X and Y , Table 5 reports the corresponding values of γ for τ ranging from −1 to +1, with a step length of 0.1, and for several values of the common number of categories. One can notice that, differently from the two cases analyzed before, due to the asymmetrical nature of the margins, γ is no longer an odd function of τ , γ (τ ) being generally different from −γ (−τ ).
For a given value of τ > 0, γ is a decreasing function of K ; for a given value of −0.9 < τ < 0, γ is an increasing function of K . When τ = −0.9, we have a curious behavior of γ as a function of K , which is displayed in Fig. 4.
It is worth underlining how the values of |γ | are significantly larger than the corresponding |τ |, especially when the number of categories is small. This feature is more apparent for negative values of τ . For example, when τ = −0.5 and H = 5, the value of γ is −0.69721. Even when discretizing the two continuous components through an extremely large number of ordered categories, such as 100, the relative difference between the values of τ and γ is still non-negligible (not smaller than 2.6% when −0.9 ≤ τ ≤ 0.9).

Building a Bivariate Ordinal Variable with Prescribed Margins and Gamma Coefficient
Researchers can be interested in building and simulating samples from a bivariate (multivariate) ordinal distribution with prescribed marginal distributions and (pairwise) levels of association, Kruskal's γ for an ordinalized bivariate normal random variable with Kendall's rank correlation τ = −0.9 as a function of the common number of categories (from 2 to 10) of the two ordinal triangular variables expressed in terms of the gamma coefficient. Such a concern may arise when one has to show the appropriateness and study the performance or robustness of a novel statistical technique, since it is not always possible to do so by using analytic arguments, except in elementary cases; on the contrary, generating data replicates that mimic the real data's characteristics of interest allows one to study the performance of the statistical method in any given setting. For example, researchers may be interested in assigning and preserving in their simulation study the marginal structures as well as ordinal associations (Demirtas 2006). In other circumstances, their concern is matching marginal means, variances, skewnesses and kurtoses, and intercorrelations (Vale and Maurelli 1983). van der Ark and van Aert (2015), with the aim of investigating the performance of different types of confidence intervals for γ , constructed several bivariate ordinal distributions by assigning each a fixed value of γ and imposing constraints on the margins, considering uniform or skewed distributions. However, when assigning the margins and association (correlation) value for two random variables, a uniqueness issue arises for the corresponding bivariate distribution. In fact, one of the fallacies of Pearson's correlation is that if we consider a bivariate random vector (X, Y ) with assigned marginal distributions F 1 and F 2 and a feasible value of ρ, then F 1 , F 2 , and ρ do not in general determine the joint distribution F univocally; on the contrary, there may exist several (even infinite) joint distributions whose margins are F 1 and F 2 and whose linear correlation is ρ (see, e.g., McNeil et al. 2005, p. 202). This issue is not overcome even by other dependence measures such as Spearman's ρ and Kendall's τ , nor most likely by Kruskal and Goodman's γ . Nevertheless, if we restrict our focus to the family of joint distributions obtained by discretizing a specific class of continuous random vectors, we can find that there is a unique joint distribution satisfying the match of margins and association value.
Let us consider a bivariate ordinal variable with joint probabilities p i j , i = 1, . . . , H , j = 1, . . . , K , and marginal probailities p i· = K j=1 p i j and p · j = H i=1 p i j . Let us start with the simplest case: for a bivariate ordinal distribution with dichotomous margins, H = K = 2, the problem consists of finding the solution to the following equation system: and then From the first equation of the system above, we derive that p 22 is the unique real root (the one with minus sign) of the second-order equation ; the other joint probabilities can then be easily computed in cascade from the remaining equations. Thus, the problem has a unique feasible solution. For example, if we assign the margins p 1· = 0.5, p ·1 = 0.5, and set γ = 3/5, we obtain the following second-degree equation: 3x 2 − 4x + 1, whose roots are p 22,1 = 1 and p 22,2 = 1/3; taking the latter leads to the solution: Note that this set of probabilities corresponds to the joint distribution of the bivariate ordinal variable obtained by discretizing a bivariate normal (3) with ρ = 1/2 (or, equivalently, τ = 1/3; see Eq. (4)).
For larger numbers of categories, there will be infinite solutions, which are not easy to derive, due to the nonlinear nature of the problem: the gamma coefficient is not linear in the p i j 's, so the equation system is nonlinear, as we have already experienced with a 2 × 2 contingency table. This means that there are several bivariate ordinal distributions p i j whose margins are the assigned p i· and p · j and whose gamma coefficient is equal to a prespecified value γ ∈ [−1, +1]. If one is interested in building just one bivariate ordinal variable with assigned margins and gamma, he/she can restrict his/her attention to the class of bivariate ordinal variables with the assigned margins derived by ordinalization of a standard bivariate normal distribution; then he/she can exploit the results of the previous section, by accommodating the value of the parameter τ in order to match the assigned value γ between the assigned ordinal margins. Due to the monotonic relationship between γ and τ , whose form depends on the selected margins, and the fact that γ can span its entire natural range [−1, +1], which we proved empirically, we are confident that there will be a unique bivariate ordinalized distribution satisfying the requested requisites. This means that among all the bivariate distributions with assigned margins and γ , we select the unique one whose underlying continuous latent model is the bivariate normal.
Writing the relationship between the two association measures as γ = g(τ ; F 1 , F 2 ), we just need to find the (unique) root τ of the equation γ − g(τ ; F 1 , F 2 ) = 0, γ being an assigned number in [−1, +1] and F 1 , F 2 being assigned ordinal distributions. Recovering the correct τ is a task that can be carried out by using an iterative procedure, which requires setting an initial value. Since we have empirically found that γ after discretization is -in absolute value -larger than the rank correlation τ of the bivariate normal random variable, one can use τ (0) = γ as a starting value for the unknown τ . One can then determine the corresponding value of the bivariate ordinal random variable γ (0) and then iteratively adjust the value of τ according to some updating rule till γ (t) converges to the assigned value. The updating process shall take into account that when τ is zero, γ is zero, too, and should prevent the updated value of τ (t) from escaping the interval [−1, +1]; a proposal is suggested in the following algorithm: 1. Set τ (0) ← 0, γ (0) ← 0; let > 0 be an arbitrarily small number 2. Set t ← 1 and τ (t) ← γ 3. Compute F(i, j; τ (t) ) by using (6) 4. Compute p (i, j; τ (t) ) by using (7) 5. Compute γ (t) for p (i, j; τ (t) ) by using (1) 6. If |γ (t) − γ | < stop; else set t ← t + 1, go back to 3.
Given the monotonic relationship between τ and γ for the bivariate normal model and its ordinalized counterpart, the above algorithm should be able to recover the value of τ inducing the target γ in a few steps, for any choice of γ and of F 1 and F 2 . The algorithm is implemented in the R environment and is freely available as supplementary material.
Example 1. Consider the following margins for X : p 1· = p 2· = p 3· = p 4· = 0.25, and for Y : p ·1 = 0.1, p ·2 = 0.2 p ·3 = 0.3, p ·4 = 0.4, and assign to γ the value 0.5. After only five iterations, the algorithm illustrated above recovers the joint distribution ensuring the target γ and the assigned margins, by ordinalization of a bivariate standard normal variable, which is displayed in Table 6. This joint distribution is not the unique one satisfying the requirements on the margins and gamma. A different one can be obtained as a convex combination of the cograduation (denoted by the letter "M") and countergraduation ("W") tables (Table 7): The joint probabilities p * i j = λp M i j + (1 − λ) p W i j , with λ = 0.7293, preserve the margins and ensure the target γ . It is worth underlining that the nonlinearity of γ does not allow its direct derivation for a convex combination of two joint probability distributions. In fact, the value of γ for a convex combination of two joint distributions is not equal to the same convex combination of the corresponding γ 's; in the case of combination of cograduation and countergraduation tables, it is not equal to λ · 1 + (1 − λ) · (−1) = 2λ − 1.

Application to Inference
If a bivariate ordinal sample of size n is available, one can assume it is an i.i.d. sample from an ordinalized bivariate normal distribution, or, more generally, from an ordinalized bivariate continuous distribution whose unique copula is the Gaussian one (see, in this sense, Grønneberg and Foldnes 2019; Foldnes and Grønneberg 2019, for a more detailed account on identifiability issues). One has then to estimate the value of the unknown τ (or ρ) and those of the unknown thresholds that define the marginal distributions of X and Y ; this can be done by numerically maximizing the joint log-likelihood of the observed sample with respect to all the unknown parameters simultaneously (see Olsson 1979, for details).
As an alternative to this full maximum likelihood approach, one can resort to the following method of moments for estimating the thresholds and τ : 1. Compute the empirical cumulative distribution functionsF 1 andF 2 based on the bivariate sample, and the thresholds by using the inverse cdf (see Sect. 3). Compute the sample value of γ ,γ M , based on the bivariate sample. 2. Compute the value of τ ,τ M , of the underlying bivariate normal distribution, inducinĝ γ M givenF 1 andF 2 , by resorting to the iterative procedure of Sect. 4. That is, since we can write γ = g(τ ; F 1 , F 2 ) and τ = g −1 (γ ; F 1 , F 2 ),τ M = g −1 (γ M ;F 1 ,F 2 ). Given the way the estimates are derived, this method can be classified as a two-stage method of moments. The unknown parameters are in fact estimated in two sequential steps: first, the thresholds for Z 1 and Z 2 are estimated independently based on the empirical cdf of X and Y , respectively, and the sample value of γ is calculated as well; then, the dependence parameter τ is estimated based on the quantities computed at the first stage.
A numerical example with artificial data is provided in the R code provided as supplementary material, where the two methods are implemented and compared.

Real Data
In a now classic study of mental health in Manhattan, New York, Srole and Fischer (1978) explore the relationship, among others, between mental impairment (Y ) and parents' socioeconomic status (X ). Table 8, from that study, has been used extensively to illustrate the utility and application of models for ordered categorical data. The sample value of Goodman-Kruskal's γ is 0.15429. By assuming a bivariate standard normal distribution underlying the bivariate ordinal data, we can derive the maximum likelihood estimates (MLEs) for Kendall's tau and the thresholds for the two margins (see Table 9). Note that the MLE of τ is 0.10762, quite a bit smaller than the corresponding estimate of γ , confirming the empirical results of the study in Sect. 3. However, the estimate is significantly greater than zero, with an associated standard error of 0.01718, which provides strong evidence that mental health status and parents' socioeconomic status are positively associated (i.e., higher socioeconomic status is associated with better mental health), even if the strength of that association is quite small. Note also that all the MLEs of the thresholds are highly significant, except for θ 3 , which is the central threshold for the variable parents' socioeconomic status.
If one applies the two-stage method of moments in Sect. 5, the estimate of τ is equal to 0.10665, just slightly different from the maximum likelihood estimate.
Table 10 displays the expected joint frequencies n (o) i j under the bivariate ordinalized distribution whose parameters are set equal to the corresponding MLEs. Note that under this model, the marginal frequency distributions are not equal to the analogous ones of Table 8 (even if they are very close to each other). This is because the full MLEs of the thresholds are not equal to the corresponding marginal MLEs (which are used as starting values for the optimization routine). Comparing the observed and theoretical contingency tables, one can notice some slight discrepancies between homologous frequencies. In order to (approximately) evaluate the goodness-of-fit of the suggested ordinalized bivariate normal model to the data at issue, we computed the usual chi-squared statistic based on the theoretical and observed frequencies (all theoretical joint frequencies are greater than 5, so there is no need for pooling cells); its value is 8.84. Under the null hypothesis that the data come from the ordinalized bivariate normal distribution with parameters equal to their MLEs, this statistic approximately follows a chi-square distribution with a number of degrees of freedom equal to 24 − 9 − 1 = 14, where 24 is the number of pooled frequencies and 9 is the number of estimated parameters. The corresponding p-value is 0.841, leading us to comfortably accept the null hypothesis. The value of the log-likelihood ratio statistic, given by i j /n i j ), is equal to 8.959 ( p-value 0.834) and thus leads to the same conclusion.
Other more sophisticated models may fit the sample data better. For some examples of alternatives (namely, association models) fitting the Midtown Manhattan Mental study data, see for example Becker (2014). In particular, a uniform association model can be used: the expected cell frequencies can be calculated by using the vcdExtra R package (Friendly 2017) (see the supplementary material). Note that these expected frequencies are quite close to their analogs under the ordinalized bivariate normal model, confirming the findings of Becker (1989) and as stressed by Kateri (2014). To prove this, one can compute the Kullback-Leibler distances between the two theoretical joint distributions ( p (o) i j and p (u) i j , for the ordinalized normal and uniform association models, respectively): The values are computed as K 1 = 0.0001721 and K 2 = 0.0001724, which are very small, proving that the two distributions are very close to each other.

Conclusions and Further Research
We focused on the bivariate standard normal variable and studied the change in association before and after ordinalization, measured by Kendall's τ and Goodman and Kruskal's γ . This analysis has been facilitated (i) by the wide availability of software implementations of the bivariate normal cdf (despite its non-closed expression), which is required when computing the joint probabilities of the bivariate ordinal distribution resulting from discretization, and (ii) by the analytic relationship between Pearson's correlation ρ and Kendall's τ for the bivariate normal distribution. We empirically investigated the relationship between τ and γ by considering several specific configurations for the two final marginal distributions (uniform, unimodal or bimodal symmetric, triangular). The study confirmed a somewhat expected result, i.e., the association measure tends to inflate after discretization, and to a larger extent when the number of ordered categories is small. For a same number of categories, we also highlighted the effect of the marginal probabilities in the change of association. Based on these results, we also elaborated a scheme for building and simulating samples from a bivariate ordinal distribution with assigned margins and value of association.
Other bivariate continuous distributions or bivariate copulas can be explored when the interest is in analyzing the effects of ordinalization on Goodman and Kruskal's γ or in constructing a bivariate ordinal variable with assigned marginal distributions and association; one should prefer the parametric copula families satisfying conditions (i) and (ii) mentioned above. Moreover, one should focus on so-called "comprehensive" bivariate copulas (Nelsen 2006, p.15), i.e., copulas able to model continuously the whole range of dependence from the lower to the upper Fréchet bounds passing through the product copula. As for (i), there is a wide variety of parametric copulas with a closed-form expression of their joint cdf. As for (ii), an analytic expression linking the copula parameter to its Kendall's correlation is not always available. For example, for the Frank family of copulas, C(u 1 , u 2 ; θ) = − 1 θ ln 1 + (e −θu 1 − 1)(e −θu 2 − 1) we have the following relationship between Kendall's τ and θ : τ (θ) = 1 − 4[1 − D 1 (θ )]/θ , with D 1 (x) = 1 x x 0 t e t −1 dt, which can be computed numerically. Note that letting the parameter θ go to 0, the Frank copula boils down to the product copula, with τ (0) = 0; when θ → −∞, C reduces to the countermonotonicity copula; when θ → +∞, C reduces to the comonotonicity copula. The Plackett family of copulas has the following form: C(u 1 , u 2 ; θ) = 1 + (θ − 1)(u 1 + u 2 ) − [1 + (θ − 1)(u 1 + u 2 )] 2 − 4θ(θ − 1)u 1 u 2 2(θ − 1) , with θ ∈ (0, +∞) \ {1}. When θ → 1, it reduces to the product copula, whereas for θ → 0, it tends to the countermonotonicity copula, and for θ → ∞ to the comonotonicity copula. For this family, there does not appear to be a closed-form expression for Kendall's τ (Nelsen 2006, p.171). Another point to inspect is whether Kendall's tau for a bivariate continuous distribution is always smaller -in absolute value -than Goodman and Kruskal's gamma computed on any Relationship between τ and γ for a bivariate Student's t random variable before and after ordinalization (marginal probabilities for the two identical ordinalized variables are 1/3 and 2/3). In both graphs, the dotted line is the bisector of the first and third orthants. In the right graph, it can be better appreciated the behavior of γ when τ is closer to zero ordinalized version thereof. We have empirically shown that this happens for the bivariate normal distribution, in a certain sense reversing what happens with respect to the correlation coefficient, whose absolute value is always diminished by discretization of both components (Lancaster's theorem). Actually, it happens that if we discretize the two components of a bivariate Student's t distribution when ρ = τ = 0, the corresponding value of γ between the two ordinalized variables is generally different from zero, so γ does not generally inherit the sign of τ , and the inequality |γ | ≥ |τ | is no longer true. For example, if we consider two identical marginal distributions with probabilities 1/3 and 2/3 for the two ordered categories, then the bivariate Student's t with 3 degrees of freedom and with uncorrelated components is discretized into the following bivariate ordinal variable: X, Y y 1 y 2 x 1 0.1153 0.2180 x 2 0.2180 0.4487 whose value of γ is about 0.0424. In the graphs of Fig. 5, the relationship between τ and γ is displayed for the ordinalized bivariate Student's t. The reader can also refer to Jin and Yang-Wallentin (2017) for some potential alternatives to the bivariate normal distribution assumption as un underlying stochastic model for ordinal data, where the authors study robustness against misspecification of the underlying distribution with respect to the polychoric correlation estimation.
As an aside, the paper illustrated how bivariate ordinal models can be built for simulation studies, as an alternative to existing procedures employing categorical marginal models. Extension of the proposed procedures to the multivariate case is not straightforward. Analogous issues in determining the existence of a multivariate binary/ordinal variable with assigned margins and Pearson's correlation matrix have been examined by Cario and Nelson (1997), Chaganty and Joe (2006), and Barbiero and Ferrari (2017). As we noticed, for a multivariate ordinal variable, the gamma association matrix is not necessarily positive semidefinite, so once all the margins are assigned, it is not straightforward to establish whether the assigned values of pairwise gamma coefficients lead to a feasible association matrix or not. However, if one selects a valid correlation matrix, this is a feasible Kendall's tau correlation matrix for any choice of the (continuous) margins, and thus, per the arguments related to the change in magnitude after discretization, it should also be a feasible gamma association matrix for any choice of the ordinal margins. The points raised in this conclusive section can be addressed by future research.

Supplementary material
Computer code developed in the R environment is available at https://tinyurl.com/PMET-D-19-00173.