Multiple Ordinal Correlation Based on Kendall’s Tau Measure: A Proposal

: The joint analysis of various ordinal variables is necessary in many experimental studies within research ﬁelds such as sociology and psychology. Therefore, the necessary measures of multiple ordinal dependence must be easy to interpret and facilitate the interpretation of multivariate models that ﬁt ordinal data. The main objective of this article is to propose a multiple ordinal correlation measure based on a bivariate correlation measure: Kendall’s tau. A sample version of the measure is proposed for its estimation. Furthermore, a conﬁdence interval and a multiple ordinal independence test are proposed. The measure is applied to various simulations, covering a wide range of multiple ordinal dependency scenarios, in order to illustrate the adequacy of the measure and the proposed inferential techniques. Finally, the measure is applied to a real-world study based on a social survey of the levels of life satisfaction and the happiness index of a population.


Introduction
The presence of ordinal variables is common in studies and experiences based on structured questionnaires and evaluations by judges and expert observers. Thus, ordinal data frequently arise in many fields, especially in the social and biomedical sciences. However, it is also common for the statistical analysis of these variables to be carried out by applying techniques oriented to the analysis of nominal variables (if the number of ordered categories is small) or techniques oriented to the analysis of numerical or quantitative variables (if the number of ordered categories is large). One can suppose that this is due to the difficulty of interpreting ordinal statistical models, such as ordinal regression models, and especially the difficulty of transferring statistical and mathematical conclusions to conclusions in their specific field of study (economic, social, biomedical, etc.). Good works of presentation, development, and interpretation of techniques for modeling ordinal variables are presented by Agresti [1] and Fahremeir et al. [2].
This difficulty of interpretation does not arise when bivariate ordinal correlation measures are applied. The two most popular measures of association for ordinal random variables are Kendall's tau [3] and Spearman's rho [4]. However, to the best of our knowledge, it is necessary to propose new multiple ordinal correlation measures that are easily interpretable for researchers from different scientific areas. That is, measures describing the ordinal correlation of a target ordinal variable against a set of explanatory variables (whether they are ordinals or not). Agresti [5] proposes a generalization of Kendall's tau defined in terms of the proportional reduction in prediction errors obtained by predicting the ordering of pairs of observations on the objective variable based on orderings of the pairs on the explanatory variables. However, it does not have as simple an interpretation as bivariate Kendall's tau. Nagelkerke [6] proposes a general definition of a likelihoodbased coefficient of determination, from which correlation or association measures can be constructed. Specifically, Muñoz-Pichardo et al. [7] apply this idea to regression models for multivariate count data. In the case of continuous numerical variables, the sample multiple linear correlation coefficient (R 2 ) is available, which can be interpreted as a measure of linear statistical dependence and a measure of fit of a model.
In the present work, a multiple ordinal correlation measure is proposed that uses the concept of the Kendall's tau coefficient as a basis and is constructed following a similar path to the multiple linear correlation coefficient, considering the adjusted values through a model of ordinal regression. That is, it is not based on the model, nor on the prediction of the ordering through the joint ordering of the explanatory variables. The proposed measure is based on the prediction of the ordering through the ordering of the values fitted by the model. This model-based strategy of construction of statistical dependency measures can be applied to many models and types of variables. Once the measure is defined in what follows, then a methodology is proposed to develop its statistical inference, based on resampling methods to determine bias, confidence interval, and hypothesis testing. To illustrate the application of this methodology, simulations of various scenarios are carried out, from a scenario with statistical ordinal independence to a scenario of almost perfect ordinal dependence.
To illustrate the adequacy of the proposed measure and its easy interpretation, it is applied in the field of social studies. In recent years, there is an increasing need to evaluate the society with a focus on the population and with such essential objectives as happiness and satisfaction with life. Within this framework of social and scientific interest, a survey was carried out in the geographical area of Andalusia (southern Spain). Its main objective was to obtain a personal assessment of the level of happiness, the level of life satisfaction in general, and certain specific aspects of life (family, education, work, etc.). Hence, the study determined a happiness index and a level of satisfaction with life in general, in addition to a life satisfaction profile defined using ordinal measures on various topics or issues of personal life. The modeling of life satisfaction and the degree of happiness in relation to other specific aspects of life has been studied by other authors (for example, refer to the work of Chow [8] on life satisfaction among university students). Multiple ordinal dependence of the happiness index against the satisfaction profile and the degree of general satisfaction with the satisfaction profile will allow us to illustrate the adequacy of the ordinal correlation measure proposed in the present work.

Statistical Model: Ordinal Multiple Regression
A large number of parametric ordinal models can be found in the literature. A possible classification of these models determines three categories [9]: cumulative models, sequential models, and adjacent-category models. Bürkner and Vuorre [9] indicate that the adjacent-category models are literally chosen for their mathematical convenience rather than any quality of interpretation, because it is difficult to envisage a natural process leading to these models. Fahrmeir et al. [2] argue that, due to its construction, the sequential model is useful when the categories of the response variable are obtained sequentially, that is, a category can only be attained provided that the previous borderline category has been attained. Therefore, we have chosen to apply a cumulative model for the analysis of the experience presented in the introduction and developed in Section 5. This model is one that provides the best results in a data set associated with this experience. Furthermore, the methodological proposal to construct the multiple association measure can be adapted to the other ordinal regression models.
The cumulative model of ordinal regression [1] consists of relating the logit of the cumulative probability of the target variable with the linear predictor of the explanatory variables. Thus, given the objective variable Y that takes the values 1, 2, . . . , K, the advantages (odds) and its logarithmic transformation (logits) are considered: Hence, Agresti [1] calls this a cumulative link model. Thus, the model represents each logit as a function of the linear predictor of the explanatory or predictive variables. Each cumulative logit has its own intercept: the β 0,j are increasing in j, because Pr[Y ≤ j |x ] increases in j for x fixed. Furthermore, the model assumes identical effects of the covariates for each logit. If we denote by γ j = exp β 0,j , the conditional probability distribution function that describes the model is given by: The model described by the conditional probability distribution is denoted as OrdMod(Y X 1 , . . . , X p ) . We can easily obtain the probability function: In the case where the explanatory variables are not significant in the model, that is, the objective variable is not associated with the explanatory variables, the resulting model, which can be considered as the base model, has the following cumulative probabilities: Consequently, it must be verified that the constants γ j , generators of the probabilities of the base model, must be ordered in the form: 0 ≤ γ 1 ≤ . . . ≤ γ K−1 . In the ordinal regression model, the generating constants of the conditional ordinal distribution Y |x are transformed to a function of the covariates: Thus, for a covariant X 1 binary factor (values 0,1), if η 1 > 0 then the presence of that factor leads to a contraction (e −η 1 ) that affects all the generating constants. Similarly, if η 1 < 0 then the presence leads to an expansion of these constants. To facilitate an interpretation, let us consider the odds ratios (OR) associated with this covariate: otherwise, Consequently, if η 1 > 0 (equivalently, associated odds ratio less than 1), then that is, the binary factor X 1 is a covariant whose "presence" leads to improve or increase the objective variable Y. In the opposite direction, the case η 1 < 0 can be interpreted.
Conversely, if X 1 is a quantitative or ordinal variable, with η 1 > 0, the increase of one unit in it leads to the increase of the objective variable and vice versa. The estimators obtained through maximum likelihood, with approximate numerical methods, allow obtaining the fitted model. Specifically, the collection of parameters that must be estimated is β 0,1 . . . β 0,K−1 , η 1 . . . η p . For the study of the significance of the parameters, the asymptotic likelihood ratio test (LRT) can be applied, as well as the global significance test of the model. Furthermore, as measures of goodness of fit, measures based on the likelihood of the model can be considered: the residual deviance statistic and the Akaike Information Criterion (AIC) [10]. The detailed development of the inference on this model is collected in the works of Agresti [1], McCullagh [11], and Fahremeir et al. [2].
With such estimators, the model allows "predicting" the values of the target ordinal variable through the covariates:ŷ This makes it possible to analyze the residuals to evaluate the adjustment and the predictive capacity of the model OrdMod Y X 1 , . . . , X p . Following the approach of Agresti (2010), we consider as a classification or prediction rule the response category that has the highest estimated probability or mode of the conditional distribution Y|X = x . This approach may suffer from some problems such as: Cases in which one outcome is much more likely than the others, this can result in always or nearly always predicting that category. This problem leads to the absence of multiple ordinal association.
Cases in which the conditional distribution of the ordinal response variable is bimodal. For the purpose of defining the measure, it has been chosen to randomize between modal categories. A detailed study of the problem may be the subject of future work.
Other prediction rules have been proposed for this model (see Agresti [12] and Marshall [13]), but we believe that prediction in the modal category is the most appropriate rule for defining the measure of association.

Definition of the Model-Based Multiple Kendall's Tau
To complete the study of the ordinal regression model, it is necessary to obtain a measure of the multiple dependency or multiple association between the target ordinal variable and the collection of explanatory variables included in the model, that is, a measure similar to the multiple correlation coefficient in the multiple regression model under normal conditions. The multiple correlation coefficient between a quantitative random variable U and a quantitative random p-vector V is defined as the maximum correlation between U and any linear combination of V, β V [14]: Another equivalent way of defining this coefficient is, under normal conditions, the linear correlation coefficient between U and the random variable U * = E[U|V] , with support identical to U and a function of V: That is, given that, Consequently, two relevant aspects can be highlighted: first, the multiple correlation coefficient is a measure of the linear dependence between one random variable and a certain collection of random variables; second, the linear dependence is measured with respect to the best linear combination that approximates U according to the multiple regression model [9,10]. As a result, the multiple correlation coefficient is the linear correlation coefficient between U and its prediction β V based on the model. It is generally considered squared since the sign of the correlation is not understandable in the multiple case.
From a sampling perspective, the estimator of the sample multiple correlation coefficient is obtained as the sample linear correlation coefficient between the observed values of U and their values adjusted for the modelÛ =β V withβ as the maximum likelihood estimator of the regression coefficients or parameters of the regression model.
In parallel with the analysis of multiple ordinal dependence, we propose a new measure based on the ordinal regression model, OrdMod(Y X 1 , . . . , X p ) , and Kendall's bivariate ordinal dependence measure or Kendall's τ coefficient.
The non-parametric correlation coefficient (or measure of association) known as Kendall's tau [3] is defined as follows: let (U 1 , V 1 ) and (U 2 , V 2 ) be independent twodimensional random vectors with the same distribution as (U, V), then That is, the Kendall's τ coefficient is the probability of agreement minus the probability of disagreement between the variables. This coefficient can also be defined (see [15]) as the expected value of the random variable Its estimation, through a sample {(u i , v i ), i = 1 . . . n} is obtained asτ = 2S n(n−1) , named Kendall's sample tau, where S is the number of concordant pairs minus the number of discordant pairs.
Thus, a measure of multiple ordinal dependence between an ordinal variable Y and a random p-vector X = X 1 , . . . , X p can be defined as Kendall's τ coefficient between Y and the prediction with the best linear combination of X with respect to the OrdMod(Y X 1 , . . . , X p ) model. In other words, it can be called the model-based multiple Kendall's tau between Y and X and is denoted by with mode[Y|X = x] the modal value of conditioned distribution, or equivalently Thus, Y * = mode[Y|X] is a function of the random vector X, as in the regression model under normal conditions, In the case of the ordinal independence of Y with respect to the random vector X, the variable Y * will be constant, Y * = mode(Y), consequently, by the definition of Kendall's tau coefficient, In general, the converse is not true, but this disadvantage is shared by the multiple linear correlation coefficient. Appendix A illustrates this measure through a simple example.
If τ M (Y, X) = 1, the relationship is direct and perfect in the sense that, that is, perfect concordance or agreement between the target ordinal variable and the ordinal variable adjusted to the model OrdMod(Y X 1 , . . . , X p ) . In practice, there will be scenarios where perfect concordance and independence are absent. In these scenarios, increasing degrees of concordance between Y and Y * are reflected by increasing positive values of τ M (Y, X).
Since the ordering of an ordinal variable is arbitrary, we can wonder if the definition is coherent with respect to a change of variable from Y to K + 1 − Y. The answer is positive since the τ M 's Kendall measures the ordinal correlation between the target variable and the model prediction variable. We can also speculate about the problem of ties that occur often with ordinal data. Since the τ M measure is defined as a bivariate measure, this problem can be approached identically to the bivariate case (see Agresti [12]).

Inference: Estimation, Confidence Interval, and Test
Proceeding analogously to the multiple regression model, an estimate of the aforementioned measure via a sample {(y i , With this sample measure of multiple ordinal dependence,τ M (Y, X), it is possible to construct a test on the ordinal non-association of Y with respect to the vector X, that is, a test on the null hypothesis As noted by Bergsma and Dassios [16], a test of independence based on i.i.d. data can be obtained by application of the permutation test (or randomization test) to an estimator of τ M . For moderately large sample sizes, it is not appropriate to use an exact permutation test. As with bootstrapping, a permutation test constructs the sampling distribution of statistics of the test, under the null hypothesis, by resampling the observed data through permutations. Specifically, it can be applied to the Monte Carlo approximation or resampling test, that is, a test of randomly selected permutations. In this case, starting from the initial sample {(y i , Through this, the B permutated samples are obtained of the form: The sample multiple Kendall's tau is calculated for each of the permutated samples, thus obtaining a collection of values τ M(k) ; k = 1, . . . , B which allows obtaining a sampling distribution of the measure under the null hypothesis of ordinal independence between Y and X. This sampling distribution can be used to perform the test, since it is an approximate distribution to the distribution of the set of all n! permutations. Specifically, an approximation to the permutation-achieved significance level (ASL) is given by Finally, the confidence interval (CI) must be obtained to complete the inference on the parameter τ M (Y, X). Given the difficulty of obtaining an explicit expression and in accordance with the work of Ruscio [17], the conclusion is that "coverage of bootstrap CIs was usually as good or better than coverage of analytic CIs". Furthermore, with the application of the bootstrap methodology, estimates of the bias and the estimation error of the estimatorτ M (Y, X) can be obtained. According to Canty [18], under suitable regularity conditions, given a theoretical distribution F, an associated functional or parameter θ(F) and an estimatorθ = θ F defined on the empirical distribution functionF of a random sample, we can approximate the bias distribution of the estimator (θ − θ) through the empirical distribution ofθ * −θ, whereθ * = θ F * is the estimator defined on the bootstrap distribution. That is, considering R bootstrap resamples of the original sample, the estimated bias is and, consequently, its standard error of estimation is v 1/2 θ . The boot package [18,19] has implemented various methods of obtaining bootstrap CIs: basic bootstrap CI (also known as reverse percentile bootstrap CI), percentile bootstrap CI, studentized bootstrap CI, and bias-corrected and accelerated bootstrap CI (BCa). We believe that this last method may be the most suitable for our objective. As noted by Jung et al. [20], to overcome the problems of overcoverage in the bootstrap confidence intervals based on percentiles [21,22], the BCa method corrects both the bias and the skewness of the bootstrap estimates of the parameters by incorporating a bias correction factor and an acceleration factor [18,19]. In the present study, R = 2000 bootstrap resamples were both generated and analyzed because this value is recommended in the bootstrap literature [23].

Simulation Study
In order to illustrate the adequacy of the proposed multiple ordinal correlation measure, the situations described below are randomly generated, beginning with the simple case (p = 1), with a strong association, either positive or negative. Subsequently, the multiple case is analyzed (p = 3), from a strong multiple association to ordinal independence. In all simulations, the sample size generated is n = 500.
All computations have been performed in R [24]. To avoid the overfitting of the model and the impossibility of carrying out the calculations to obtain the estimators, an error is randomly included in 2.5% of cases, in which the value of Y is generated randomly according to a uniform discrete distribution with support {1, 2, . . . , 10}. That is, the errors are generated independently of the explanatory variables (see Appendix B). This generation process is also applied in the following simulations. • Simulation 2: negative association with a single variable. Identical to Simulation 1, with coefficient η 1 = −1 < 0.
• Simulation 3: association with a three-dimensional vector (2.5% error). Three explanatory variables X 1 , X 2 , X 3 that are independently distributed according to N(0, 1) are generated and also an ordinal variable Y is generated with support {1, 2, . . . , 10} following the pattern of the conditional distribution described in the model OrdMod(Y|X 1 , X 2 , X 3 ) with negative and positive coefficients η = (0.8, −1.0, 0.2) . This is intended to include positively (X 1 , X 3 ) and negatively (X 2 ) associated variables, with both strong (X 1 , X 2 ) and weak (X 3 ) intensity. The base distribution of Y is considered to be identical to the base distribution included in the previous simulations, with the specified error. • Simulation 4: association with a three-dimensional vector (5% error). Simulation procedure is identical to that collected in Simulation 3, with error in 5% of the cases. • Simulation 5: association with a three-dimensional vector (10% error). Simulation procedure identical to that in Simulation 3, with error in 10% of cases. • Simulation 6: association with a three-dimensional vector (error 25%). Simulation procedure identical to that in Simulation 3, with error in 25% of cases. • Simulation 7: association with a three-dimensional vector (50% error). Simulation procedure identical to that in Simulation 3, with error in 50% of cases. • Simulation 8: association with a three-dimensional vector (75% error). Simulation procedure identical to that in Simulation 3, with error in 75% of cases. • Simulation 9: association with a three-dimensional vector (90% error). Simulation procedure identical to that in Simulation 3, with error in 90% of cases. • Simulation 10: null association or multiple ordinal independence. Three explanatory variables X 1 , X 2 , X 3 independently distributed in accordance with N(0, 1) are generated and an ordinal variable Y is also generated according to a discrete uniform distribution with support {1, 2, . . . , 10}. That is, it can be considered to be generated following the ordinal model, including a 100% error.
The results of the simulations are shown in Table 1. For each simulation, the parameters used are indicated (error and coefficients of the linear predictor), the bivariate Kendall's tau coefficients are obtained with their corresponding level of significance (calculated through the function "cor.test ()" of R-Program [24]) and, finally, the multiple ordinal correlation coefficients and their significance obtained according to the permutations test with number of replications N rep = 2000 are collected. For fitting of the models, the function polr () of the MASS package [25] is applied.
The proposed measure adequately reflects the multiple ordinal correlation. In the scenario defined for Simulation 3, the data have been generated according to the multiple ordinal model, except for a disturbance in 2.5% of the data generated, obtaining an estimate of the multiple Kendall's tau measureτ M (Y, X) = 0.9501, which reflects the almost functional relationship of the target variable versus the explanatory variables.
From Simulation 4 to the final simulation, there is an increase in the percentage of data generated that does not respect the model, with a gradual decrease in the measurement (seen in Simulation 4τ M (Y, X) = 0.9397 and Simulation 9τ M (Y, X) = 0.0429). Finally, in the last scenarios (Simulations 9 and 10), the data of the objective variable have been generated quasy-independent or independently of the predictor variables, obtaining estimates of the measure of multiple ordinal correlation that are practically null (τ M (Y, X) = 0.0429, 0.0297), thus being able to accept the nullity of τ M (Y, X) according to the permutations test (p-value ≈ 0.2715, 0.3095). In Scenarios 1 and 2, the data have been generated according to the OrdMod(Y|X 1 ) model with positive and negative coefficients (η 1 ), respectively. The bivariate Kendall's tau coefficient does not coincide with the multiple Kendall's tau coefficient. This fact is logical since the former is a measure of the ability to predict the order of Y through the order of X 1 , τ(Y, X), while the latter is a measure of the ability to predict the order via the order of the variable fitted by the model, τ(Y, Y * ). Furthermore, given that the data have been generated with respect to this model, in both scenarios, the multiple measure is greater than the bivariate measure in absolute terms. Table 2 shows the BCa bootstrap confidence intervals for each of the simulations. They have been obtained from the bootstrap library [15,16], generating R = 2000 bootstrap resamples. In the estimates of bias and standard error of estimation, the accuracy of the estimation procedure can be observed: the absolute value of the bias is always less than 0.01, except in Simulation 9, with bias = 0.0208; the values of the std.error are all less than 0.05. To carry out these calculations, problems have arisen in Simulations 1, 2 and 3, generated with only 2.5% error. That is, only 2.5% of the sample observations differ from the exact value generated by the model. Consequently, in a large number of bootstrap resamples, all the observations respect the model exactly. This situation means that the fit of the ordinal regression model for these resamples cannot be estimated due to overfitting in the model, or that the results in these simulations are not very stable. However, with exercising caution in these cases, the general conclusion to be drawn is that the estimation procedure, except for singular cases, is unbiased and very accurate.

Application to a Social Study: Happiness Index and Life Satisfaction Level
As presented in the introduction, the data for the study consist of an anonymized data set obtained through a survey that was carried out in Andalusia (southern Spain) [26]. One of the objectives of the survey was the analysis of the level of happiness and the degree of life satisfaction of the population, as well as obtaining an approximation to the main problems that affect and concern that population. In this survey, the concepts of happiness and satisfaction with life, that were subjectively considered by the study participants who were interviewed, were assumed. As the concepts are subjective, consequently, they are not completely identical for all individuals, in the same way that there is no broad consensus in the scientific and philosophical fields. Although the concepts have different nuances for each individual, in a practical sense they coincide fundamentally for everyone, especially when restricted to a population under the same sociocultural and temporal patterns, which in our case is the population of Andalusia.
However, since the measures of happiness and life satisfaction are obtained through a survey, both are measures of the hic et nunc (here and now), filled with a wealth of sensations that provide spontaneity and possible volatility caused by poor reflection on the part of the individuals. Specifically, the people interviewed were asked about their overall satisfaction with life and about other dimensions related to living conditions, work, leisure, and social relationships. The respondents had to respond on a Likert-type scale of 0 to 10, where 0 is defined as completely dissatisfied and 10 is completely satisfied. Thus, it is possible to work with a life satisfaction profile through ordinal measurements of nine areas or aspects of personal life (see Table 3). Similarly, for the level of happiness the Likert-type scale from 0 to 10 is used, in which 0 is defined as completely unhappy and 10 is completely happy. Table 3 presents the variables considered.
The number of interviewees who answered all questions associated with the previous indices is n = 445. A statistical summary is included in Table 3.
The happiness index and the degrees of satisfaction can be considered as ordinal variables. Therefore, for the analysis of the bivariate correlation, Kendall's tau measure is considered. Given that the two target variables are Sgen and Ihap, the bivariate ordinal correlations are obtained using the remaining variables and applying the corrplot library [27]. The results are presented in Table 4. The degrees of satisfaction regarding the various aspects considered are positively correlated with the two objective variables. Furthermore, as one might suspect, the two target variables are correlated with each other, with sample Kendall's tau equal to 0.568.   Two ordinal regression models are carried out to determine the dimensions related to the living conditions that have the most influence over the level of happiness and the degree of satisfaction with life in general.

Sfam Sjob Sleis Sfrie Shea Sphys Seduc Sbeha Srela Sgen Ihap
First, Ihap is considered as the objective variable and the degrees of satisfaction with specific aspects of life are considered as explanatory variables: Sfam, Sjob, Sleis, Sfrie, Shea, Sphys, Seduc, Sbeha, and Srela. To estimate the model, the polr () function from the MASS package [25] of the R-Program is applied, thus obtaining the results shown in Table 5. Thus, the significant covariates in the model are the degrees of satisfaction with family, work, leisure, the way of being, and romantic relationships in addition to the intercepts. All the coefficients associated with these predictors are positive, that is, they are positively associated with the happiness index, particularly the degrees of satisfaction Sfam and Sbeha. A deeper and more extensive interpretation can be made regarding the field of sociology, but it is not the object of interest in the present study.
The proposed multiple ordinal correlation measure, the multiple Kendall's tau, attains the value of 0.501, with significance p < 0.001, this latter being obtained through the permutation test with 2000 permutations generated. That is, the ordinal correlation was improved through the adjustment of the model, with respect to the bivariate correlations (see Tables 4 and 7). Finally, the BCa method was applied to obtain a confidence interval. Table 7 shows the results obtained. The estimation of the bias and the standard error of estimation provides new evidence for the accuracy and unbiasedness of the estimation method.
The statistical interpretation of the adjusted model, the estimates of its parameters, and the multiple ordinal correlation measure can lead to an adequate sociological interpretation of the results toward determining the aspects of people's lives that most favorably influence their happiness and evaluating the capacity of the vital aspects considered in explaining the level of happiness. Similarly, the present study was carried out with consideration of the degree of satisfaction with life in general, Sgen, as an objective variable. According to the results of the model (see Table 6), the most relevant aspects are the degrees of satisfaction with family relationships, work, leisure and free time, physical appearance, the way of being and acting in society, and sentimental relationships. However, personal health and the level of professional training did not influence the degree of satisfaction with life. The multiple ordinal correlation measure reaches the value 0.599, exceeding the bivariate ordinal correlations (see Tables 4 and 7). Conclusions that are similar to the previous case can be obtained with the bias and the standard error of estimation [1].
Further, as with the previous model, the statistical interpretation of all these results allows obtaining a relevant sociological interpretation.

Discussion and Conclusions
In many economic, social, and biomedical studies, data sets containing ordinal variables are considered as a relevant part of the studies. However, statistical techniques not suitable for the ordinal scale are applied. This fact may seem strange given the set of ordinal regression models proposed in the body of literature. One of the reasons that their limited use may be justified is the difficulty of interpreting the models and their parameters. For this reason, it may be necessary to propose new parameters and functions that facilitate their interpretation and the transference of statistical conclusions to conclusions within the scientific scope of the study in a particular field.
The main contribution of this article is a multiple ordinal correlation measure and the inferential techniques necessary for its estimation. Important features of this contribution are considered to be the following. First, the measurement is based on a very intuitive idea, to measure the bivariate ordinal correlation between the target variable and the values fitted by an ordinal regression model. Kendall's tau coefficient is considered as a bivariate ordinal correlation measure and the cumulative model of ordinal regression is considered as a model. This intuitive idea, which is used in the definition of the multiple linear correlation coefficient, leads to easy interpretability of the correlation measure. Second, its easy interpretation can facilitate its use and the application of ordinal regression models in many scientific fields. It can also be considered as a measure of the adequacy of the model. Third, the idea of constructing a measure based on the model broadens the path, already opened by the multiple linear correlation coefficient, to construct other measures of multiple correlation and measures of multiple association at different scales (binary, nominal, ordinal, quantitative discrete, quantitative continuous). Fourth, the proposed inference methods are also easy to understand and apply. In particular, their point estimate is a sample version of the population measure. This estimate is almost unbiased according to the bias estimated through the bootstrap methodology, although its unbiasedness has not been mathematically demonstrated, to the best of our knowledge. Furthermore, this bootstrap methodology aids in obtaining confidence intervals with a computational procedure that is not complex and tedious. Finally, the significance test of the measure, which can be considered as a multiple ordinal independence test, is obtained through a permutation test that does not require excessive computation time.
Some aspects that are not in the present study need to be deepened and they should be the object of study in a future line of research. In the field of measurement inference, the aim is to obtain confidence intervals and significance tests, exact or asymptotic, which are not based on resampling techniques. However, the measurement construction methodology can be applied to other ordinal regression models (cumulative models, sequential models, and adjacent-category models), as well as with other ordinal correlation measures such as the Spearman rank correlation coefficient [4]. The methodology can also be extended to models based on binary and nominal variables through measures of bivariate association, as proposed in the body of literature. This methodology to construct the multiple association measure is easily applicable to the other models (for binary, categorical, ordinal, and continuous data). It is enough to combine two components: statistical model and bivariate association measure. Such components should be appropriately in accordance with the scale of the data and the objective of the study. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Simple example to illustrate the definition of the multiple Kendall's tau measure. Let X 1 and X 2 be two independent and identically distributed variables with support {1,2,3} and probabilities {1/3, 1/3, 1/3}, respectively. Let Y be an ordinal random variable with support {1,2,3,4,5}. The following parameters of the model OrdMod(Y|X 1 , X 2 ) are considered: (A1) Table A1 shows the conditioned probabilities {Pr[Y = j |x ], j = 1, . . . , 5} according to the model. In addition, the modal value of each conditional distribution is included, that is, the value of the random variable Y * = mode[Y|X = x]. Table A1. A simple example: probability distributions associated with the model (1) .  (1) OrdMod(Y|X 1 , X 2 ) with β 01 = −1.7, β 02 = −0.5, β 03 = 0.5, β 01 = 1.7, η 1 = −0.5, η 2 = 0.5. Table A2 shows the joint probability distribution of (Y, Y * ) and the marginal distribution of Y * , that is,  The population multiple Kendall's τ M (Y, X) can be obtained as the bivariate Kendall's τ of the joint distribution of (Y, Y * ), that is, the population version of tau-b measure (see Agresti [12]): with Π c and Π d concordance and discordance probabilities, respectively. For the joint probability distribution included in Table A2, we obtain τ M (Y, X) = 0.3205. Finally, in order to illustrate the estimation of this measure through the sample version, we randomly generated a sample of size 60 according to the probability distribution associated with the population model OrdMod(Y|X 1 , X 2 ) with the parameters listed above. The data set is included in Table A3. Fitting the model with the polr function of the MASS library package [25] provides the estimates:  Finally, the crosstab between observed values and adjusted values (Table A4) allows us to obtain the sample version of tau-b measure (see Agresti [12]) as an estimate of the multiple ordinal correlation measure:τ M (Y, X) = 0.3243.

Appendix B
Simulation generation procedure (see Section 4). Considering a vector of explanatory variables X = X 1 . . . . , X p and the model OrdMod Y X 1 , . . . , X p .