Paper deep dive
Geometric mean-based pairwise comparison method with the reference values -- statistical approach
Konrad Kułakowski, Jacek Szybowski
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%
Last extracted: 7/14/2026, 8:35:31 AM
Summary
This paper introduces a statistical framework for the pairwise comparison method in decision-making, focusing on the geometric mean-based Heuristic Rating Estimation (HRE). It formulates incomplete geometric HRE as a linear regression problem, treating pairwise comparisons as observations and logarithmic weights as parameters. The study proposes statistical indicators to quantify weight vector quality, model inconsistency as statistical error, and handle ranking uncertainty via tie-clustering, extending traditional AHP and GMM approaches with probabilistic rigor.
Entities (7)
Relation Signals (6)
Konrad Kułakowski → affiliatedwith → AGH University of Krakow
confidence 99% · Konrad Kułakowski konrad.kulakowski@agh.edu.pl AGH University of Krakow, WEAiIB
Jacek Szybowski → affiliatedwith → AGH University of Krakow
confidence 99% · Jacek Szybowski szybowsk@agh.edu.pl AGH University of Krakow, WMS
Analytic Hierarchy Process (AHP) → subtypeof → Pairwise Comparison Method
confidence 97% · The best-known example of this method is the Analytic Hierarchy Process (AHP).
Incomplete Geometric HRE → formulatedas → Linear Regression
confidence 96% · we present a statistical interpretation of the geometric HRE method and show that estimating unknown alternative weights can be formulated as a linear regression problem.
Heuristic Rating Estimation (HRE) → incorporates → Reference Values
confidence 95% · HRE method assumes that there are two sets of alternatives, the reference ones with weights a priori known, and the non-reference ones weights
Linear Regression → models → Inconsistency
confidence 94% · providing an intuitive interpretation of inconsistency as statistical error.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:For many years, the pairwise comparison method has been widely used for decision-making involving experts. The best-known example of this method is the Analytic Hierarchy Process (AHP). In this now classic approach, the weights of alternatives are calculated using the principal eigenvector of the comparison matrix. In this paper, we present a statistical view of the pairwise comparison method, using reference values and the geometric mean to calculate alternative priorities. Thanks to this approach, we can simultaneously capture both the phenomenon of inconsistency in pairwise comparisons and the preference distance between alternatives. In this paper, we define indicators that measure the quality of the obtained weight vector, which, thanks to the statistical approach, have a clear and intuitive interpretation.
Tags
Links
- Source: https://arxiv.org/abs/2607.10038v1
- Canonical: https://arxiv.org/abs/2607.10038v1
Trouble viewing inline? Open PDF directly →
Full Text
91,938 characters extracted from source content.
Expand or collapse full text
Geometric mean-based pairwise comparison method with the reference values - statistical approach Konrad Kułakowski konrad.kulakowski@agh.edu.pl Jacek Szybowski szybowsk@agh.edu.pl AGH University of Krakow, WEAiIB, al. Mickiewicza 30, 30-059 Poland AGH University of Krakow, WMS, al. Mickiewicza 30, 30-059 Poland Abstract For many years, the pairwise comparison method has been widely used for decision-making involving experts. The best-known example of this method is the Analytic Hierarchy Process (AHP). In this now classic approach, the weights of alternatives are calculated using the principal eigenvector of the comparison matrix. In this paper, we present a statistical view of the pairwise comparison method, using reference values and the geometric mean to calculate alternative priorities. Thanks to this approach, we can simultaneously capture both the phenomenon of inconsistency in pairwise comparisons and the preference distance between alternatives. In this paper, we define indicators that measure the quality of the obtained weight vector, which, thanks to the statistical approach, have a clear and intuitive interpretation. keywords: pairwise comparisons , decision-making , heuristic rating estimation , AHP , GMM , inconsistency 1 Introduction Pairwise comparisons are one of the oldest approaches to decision-making. Their first well-documented application is generally attributed to Ramon Llull, although the underlying idea is likely much older. Early applications by Llull [11] and later thinkers, including Nicholas of Cusa, the Marquis de Condorcet, and Arthur Copeland [12, 5], were based on qualitative comparisons that allowed only for the identification of a preferred alternative. Quantitative pairwise comparisons emerged with the work of Louis Thurstone [45] and gained widespread recognition thanks to the Analytic Hierarchy Process (AHP) proposed by Thomas L. Saaty [41]. Currently, pairwise comparisons are used not only in AHP but also in methods such as the Best-Worst Method (BWM) [7, 37], PROMETHEE, ELECTRE [17], and MACBETH [2]. Despite the long history of the pairwise comparison method, it remains an active area of research, particularly in measuring inconsistency [6, 4, 9, 21], modeling uncertainty [46, 29], handling incomplete comparisons [20] or studying the properties of the method [43, 19]. Some decision-making methods use reference values or have variants that allow for such use. Examples of such methods include TOPSIS [33], VIKOR [32], the Helwig method [38], BWM [34, 37], COMET [42], DARP, BIPOLAR [39], or the Parsimonious Analytic Hierarchy Process [1]. For quantitative pairwise comparison methods, one such method is HRE (heuristic rating estimation) [27, 23, 22, 26]. In the latter, it is possible to specify virtually any number of alternatives with reference values, depending on data availability. The result of a single pairwise comparison can be interpreted as a single measurement of a decision-maker’s (or expert’s) preference. This understanding of comparison makes it possible to use statistical methods to assess the quality of a set of pairwise comparisons. This line of research—which involves pairwise comparisons but does not use reference data—appears in a number of scholarly works. In particular Carriere and Finster [8] developed a statistical framework for the ratio model of paired comparisons, treating observed comparison ratios as random variables generated by underlying preference parameters. The authors investigated the estimation and inferential properties of the model, deriving asymptotic results and statistical procedures for assessing preference intensities. In their subsequent paper, Genest and Rivest [16] interpreted pairwise judgments as noisy observations of underlying preference ratios and formulated the problem within a statistical estimation framework. By introducing a log-linear model for comparison data, they demonstrated that priority estimation can be viewed as a least-squares problem and showed the close relationship between this approach and the geometric mean method. De Jong [15] examined Saaty’s priority scaling method from the perspective of statistical estimation theory. The paper models pairwise comparison judgments as observations affected by random errors and investigates the properties of the resulting priority estimates. Ramsay developed a maximum likelihood model for multidimensional scaling, treating the observed distances as noisy measurements arising from the underlying geometric configuration of the objects. He presented the theoretical foundations and estimation procedures [35], and in his subsequent research [36], he examined the performance of maximum likelihood estimators for small samples. Lin, Gang, and Ergu [28] proposed a statistical approach to assessing the consistency of pairwise comparison matrices. Rather than relying solely on deterministic inconsistency indices, the authors modeled comparison judgments within a statistical framework and derived measures that account for random variations in decision-maker responses. Luo et al. [30] investigated statistical procedures for assessing the multiplicative consistency of fuzzy preference relations. Using Monte Carlo simulations, the authors evaluated the effectiveness and robustness of several hypothesis-testing approaches under different levels of inconsistency and uncertainty. Previous research on the statistical interpretation of the quantitative pairwise comparison method has focused on approaches that do not use reference data. Therefore, in the study, we present a statistical interpretation of the geometric HRE method and show that estimating unknown alternative weights can be formulated as a linear regression problem. In this formulation, pairwise comparisons are treated as observations, while logarithms of unknown weights are interpreted as model parameters. This perspective makes it possible to derive not only the priority vector itself, but also standard statistical characteristics of the estimation process, including the residual variance, covariance matrix of estimators, standard errors, and confidence intervals for the obtained weights. The main contribution of this paper is twofold. First, we show that incomplete geometric HRE naturally fits within the least-squares regression framework, providing an intuitive interpretation of inconsistency as statistical error. Second, we propose probabilistic measures for assessing the reliability of the resulting ranking. In particular, we introduce indicators based on the probability that the order of two alternatives implied by the estimated weights is preserved in the underlying preference model. These indicators complement classical inconsistency measures by taking into account not only the dispersion of judgments but also the preference distance between alternatives. This is important because even highly inconsistent data can still yield a reliable ordering of clearly separated alternatives. In contrast, alternatives with very similar weights may remain difficult to distinguish despite relatively low inconsistency. The proposed framework also supports the identification of alternatives whose ranking positions are statistically uncertain. For such cases, we introduce a tie-threshold mechanism that allows alternatives with insufficiently reliable ordering relations to be grouped into tie clusters. It leads to a more cautious and interpretable representation of the decision result. Finally, we discuss how the proposed approach can be extended to group decision-making, where pairwise comparisons provided by several experts are treated as a common set of observations, possibly with expert-specific error variances. In this way, the present study extends recent developments in HRE with reference values [27] by providing a coherent statistical foundation for uncertainty assessment, ranking reliability, and quality evaluation of weight vectors derived from incomplete pairwise comparison data. The remainder of the paper is organized as follows. Section 2 introduces the basic notions of pairwise comparisons and recalls the heuristic rating estimation method in both its arithmetic and geometric variants. Section 3 shows that incomplete geometric HRE can be formulated as a linear regression problem. Section 4 develops the statistical interpretation of the model, including variance estimation, confidence intervals, and probabilities of preserving ranking relations. Section 5 introduces quality indicators for the obtained weight vector and proposes a tie-clustering procedure for alternatives with uncertain ordering. Section 6 discusses extending the proposed framework to group decision-making. Finally, Section 7 summarizes the main findings and conclusions. 2 Preliminaries 2.1 Pairwise comparisons The pairwise comparisons (PC) method is a procedure used to convert a set of comparative judgments into a ranking of alternatives. Let A=a1,…,anA=\a_1,…,a_n\ be a set of alternatives, and let C=[cij]C=[c_ij] represent the corresponding set of comparisons expressed as an n×n× n matrix, where cij∈ℝ+c_ij _+ for i,j=1,…,ni,j=1,…,n. Each element cijc_ij reflects the outcome of an expert’s comparison between alternatives aia_i and aja_j. Sometimes it may happen that there is no comparison. For example, an expert may refuse to compare two specific decision options due to personal beliefs. If the result of a particular comparison is unavailable, it is denoted by cij=cji=?c_ij=c_ji=?. Definition 1. We will call a matrix C reciprocal if cij=1/cjic_ij=1/c_ji for all i,j=1,…,ni,j=1,…,n, except in cases where cij=?c_ij=?. In practice, in most cases, the pairwise comparison matrices considered are reciprocal. This approach saves time spent on evaluation and is consistent with the belief that there is no need to ask the same question twice. In the literature, there are examples of models in which the order of the alternatives being compared can significantly affect the assessment and the property of reciprocity is not preserved. The result of the priority setting procedure is a weight vector. Definition 2. The weight vector for the set of alternatives A will be given as a mapping: w:A→ℝ+w:A _+. If for certain 0<i0<i, j≤nj≤ n, w(ai)>w(aj)w(a_i)>w(a_j) occurs, then alternative aia_i is considered more preferable than aja_j. This fact is denoted by ai≻aja_i a_j. The function establishes the priority ordering of the alternatives: those that are more preferred are assigned larger weights. Similarly, if for two alternatives w(ai)<w(aj)w(a_i)<w(a_j), we will denote the fact that aia_i is less preferred than aja_j as ai≺aja_i a_j. Finally, if w(ai)=w(aj)w(a_i)=w(a_j), this means that both alternatives are equally preferred, which can be written as ai∼aja_i a_j. The weight vector w usually takes the form w=(w(a1),w(a2),…,w(an))Tw= (w(a_1),w(a_2),…,w(a_n) )^T. If there is no need to distinguish additional attributes related to the mapping itself, the notation can be shortened to w=(w1,w2,…,wn)Tw= (w_1,w_2,…,w_n )^T. Since the value wiw_i corresponds to the significance (importance, priority) of the i-th alternative, it is natural to expect that a direct comparison cijc_ij of two alternatives aia_i and aja_j will translate into the ratio wi/wjw_i/w_j. In practice, the equality cij=wi/wjc_ij=w_i/w_j rarely occurs due to inconsistency of the model. Definition 3. An n×n× n PC matrix C=[cij]C=[c_ij] is called inconsistent if there exist indices i,k,j=1,…,ni,k,j=1,…,n such that cij≠cikckjc_ij≠ c_ikc_kj. There are many methods for calculating the weight vector based on a pairwise comparison matrix [10, 25]. The two most widely used methods are the Eigenvector Method (EVM), originally proposed by Saaty, the creator of AHP [40], and Geometric Mean Method (GMM), formulated by Crawford and Williams [13, 14]. According to the first method, the weight vector is a suitably scaled eigenvector of C corresponding to the principal eigenvalue (spectral radius) of this matrix. Therefore, assuming that wˇ w is the solution to the equation Cwˇ=λmaxwˇC w= _max w, the weight vector takes the form111Where for v∈ℝv holds ‖v‖1=∑i=1n|vi| \|v \|_1=Σ^n_i=1 |v_i |. w=wˇ/‖wˇ‖1w= w/ \| w \|_1. According to this proposal, the weight of each alternative is determined as the weighted arithmetic mean of the weights of all the other alternatives. The GMM approach proposes replacing the arithmetic mean with the geometric mean. Hence, the weight assigned to the i-th alternative is calculated as the appropriately normalized geometric mean of the i-th row of the comparison matrix. Therefore, assuming that the mean vector w~ w is given as w~=((∏j=1nc1j)1/n(∏j=1nc2j)1/n⋮(∏j=1ncnj)1/n), w= ( array[]c (Π^n_j=1c_1j )^1/n\\ (Π^n_j=1c_2j )^1/n\\ \\ (Π^n_j=1c_nj )^1/n array ), the weight vector in the GMM approach is w=w~/‖w~‖1w= w/ \| w \|_1. A strong argument in favor of using GMM is its optimality in the context of the least squares logarithm (LLSM) criterion [13]. Thus, the results obtained using GMM are identical to those obtained using the LLSM approach [25]. In the original papers in which the EVM and GMM procedures were proposed, the authors assumed that the set of comparisons was complete, i.e., that every cij∈ℝ+c_ij _+. However, variants of these methods have been developed: incomplete EVM (IEVM) [18] and incomplete GMM (IGMM) [24]. In addition to the two methods mentioned above, many others have also been defined [25, 44]. A pairwise comparison matrix can be represented as a directed graph, where the vertices correspond to alternatives and the edges correspond to comparisons [25]. The direction of the edges can indicate the preferred alternative. Formally, this can be defined as follows. Definition 4. A directed graph TC=(V,E,L)T_C=(V,E,L) is called a graph of pairwise comparisons matrix C if V=A=a1,…,anV=A=\a_1,…,a_n\ is the set of vertices, E=eij:cij≠?∧cij<1E=\e_ij:c_ij≠? c_ij<1\ is the set of edges, and L:E→ℝ+L:E _+ is the labeling function such that L(eij)=cijL(e_ij)=c_ij. The discussion below also introduces the concept of the Laplace matrix for a graph G. This concept is not directly related to the pairwise comparison matrix; however, assuming that for a matrix C=[cij]C=[c_ij] its graph is TCT_C, it can be derived that the corresponding Laplace matrix L(C)=[lij]L(C)=[l_ij] takes the form lij=sii=j−1i≠j∧cij≠?,l_ij= casess_i&i=j\\ -1&i≠ j c_ij≠? cases, where sis_i denotes the number of defined comparisons plus 11 in the i-th row of matrix C. 2.2 Heuristic rating estimation Heuristic rating estimation (HRE) method assumes that there are two sets of alternatives, the reference ones with weights a priori known, and the non-reference ones weights, i.e., those whose weights need to be estimated. Let the indices of the reference alternatives be elements of the set IK⊂ℕI_K , and let the indices of the non-reference alternatives be elements of the set IU⊂ℕI_U . To make it easier to follow the theoretical discussion, we will number the non-referential alternatives from 11 to k and the referential alternatives from k+1k+1 to n. Therefore, unless otherwise specified, we will assume that the set of non-reference alternatives AUA_U is given by AU=a1,…,akA_U=\a_1,…,a_k\, and similarly, the set of reference alternatives is given by AK=ak+1,…,anA_K=\a_k+1,…,a_n\. The values of the reference alternatives, i.e., the weights wk+1,…,wn∈ℝ+w_k+1,…,w_n ^+, constitute the data of the decision model. The existence of reference alternatives means that comparisons between them are also known and are not subject to expert evaluation. HRE defines two methods for calculating the weight vector that are analogous to the EVM and GMM approaches [22, 26]. Both of them require solving a system of linear equations. These systems differ from each other. Below there is a brief description of both solutions in an extended version that takes into account the incompleteness of a pairwise comparison matrix. 2.2.1 Arithmetic approach In the first (incomplete) arithmetic approach [22, 27] to calculate the weight vector (AHRE), it is necessary to solve a system of equations of the form C¯w=b Cw=b where C¯=(1−1n−t1−1d1,2⋯−1n−t1−1d1,k−1n−t2−1d2,11⋯−1n−t2−1d2,k⋮−1n−tk−1−1dk−1,1⋯⋱−1n−tk−1−1dk−1,k−1n−tk−1dk,1⋯−1n−tk−1dk,k−11), C= ( array[]c1&- 1n-t_1-1d_1,2&·s&- 1n-t_1-1d_1,k\\ - 1n-t_2-1d_2,1&1&·s&- 1n-t_2-1d_2,k\\ & & & \\ \,\,- 1n-t_k-1-1d_k-1,1&·s& &- 1n-t_k-1-1d_k-1,k\\ - 1n-t_k-1d_k,1&·s&- 1n-t_k-1d_k,k-1&1 array ), (1) where dij=cijifcij≠?0ifcij=?,d_ij= casesc_ij&if\,\,c_ij≠?\\ 0&if\,\,c_ij=? cases, and tit_i denotes the number of missing elements in the i-th row of matrix C. The constant term vector is given as: b=(1n−t1−1c1,k+1wk+1+…+1n−t1−1c1,nwn1n−t2−1c2,k+1wk+1+…+1n−t2−1c2,nwn⋮1n−tk−1ck,k+1wk+1+…+1n−tk−1ck,nwn).b= ( array[]c\,\,\,\, 1n-t_1-1c_1,k+1w_k+1+…+ 1n-t_1-1c_1,nw_n\\ \,\,\,\, 1n-t_2-1c_2,k+1w_k+1+…+ 1n-t_2-1c_2,nw_n\\ \\ \,\,\,\, 1n-t_k-1c_k,k+1w_k+1+…+ 1n-t_k-1c_k,nw_n array ). In the vector received in w∈ℝkw ^k, the value wiw_i is equal to the arithmetic mean of all other values assigned to alternatives weighted by comparisons, i.e. wi=1/(n−1)∑j=1,i≠jncijwjw_i=1/(n-1)Σ^n_j=1,i≠ jc_ijw_j. Vector w∈ℝkw ^k can be extended to w∈ℝnw ^n by adding the reference values. After completion and normalization, we obtain a standard weight vector covering all alternatives. 2.2.2 Geometric approach In the case of the (incomplete) geometric approach [26, 27], we solve the system of equations C^w^=b C w=b (2) where C^=((n−t1−1)q1,2⋯q1,k⋮⋱⋮⋱⋮qk,1qk,2⋯(n−tk−1)), C= ( array[]c(n-t_1-1)&q_1,2&·s&q_1,k\\ & && \\ && & \\ q_k,1&q_k,2&·s&(n-t_k-1) array ), (3) and qij=−1ifcij≠?0ifcij=?.q_ij= cases-1&if\,\,c_ij≠?\\ 0&if\,\,c_ij=? cases. Vectors take the form of w^=[w^1,…,w^k]T w= [ w_1,…, w_k ]^T, b=[b1,…,bk]Tb=[b_1,…,b_k]^T, where every bi=∑j=1,j≠i,cij≠?nlncij+∑j=k+1,cij≠?nlnw(aj)b_i=Σ^n_j=1,j≠ i,c_ij≠? c_ij+Σ^n_j=k+1,c_ij≠? w(a_j). The solution of (2) is the logarithmic weight vector of the original problem given by the pairwise comparison matrix C=[cij]C=[c_ij] and the set of reference values wk+1,…,wnw_k+1,…,w_n. Therefore, in order to calculate the desired weights, we use exponential transformation, i.e. w=expw^=(ew^1,…,ew^k)T.w= w= (e w_1,…,e w_k )^T. (4) If each alternative can be indirectly or directly compared with at least one reference alternative, then equation (2) has a real solution [27]. Let us additionally denote the set of pairs of indices as O(C)=(i,j):cij≠?∧(i,j∈IU∨(i∈IU∧j∈IK)∨(i∈IK∧j∈IU)).O(C)= \(i,j):c_ij≠? (i,j∈ I_U (i∈ I_U j∈ I_K ) (i∈ I_K j∈ I_U ) ) \. The pairs in set O(C)O(C) correspond to those comparisons that required determination by an expert. Additionally, assuming reciprocity of matrix C, it is convenient to limit the number of comparisons to those “above the diagonal.” So let us denote O<(C)=(i,j):i<j∧(i,j)∈O(C).O_<(C)= \(i,j):i<j (i,j)∈ O(C) \. The solution obtained is optimal [27], i.e., it minimizes the quadratic logarithmic error function ℰ:ℝk→ℝE:R^k given as: ℰ(w1,…,wk) (w_1,…,w_k) =∑(i,j)∈O<(C)n(lncij−lnwiwj)2 =Σ^n_ subarrayc(i,j)∈ O_<(C) subarray ( c_ij- w_iw_j )^2 (5) The calculated vector w∈ℝkw ^k can be extended with reference values to form a complete vector corresponding to the weights of all alternatives. After completion, the weight vector normalization process can be performed to standardize its length as defined by ∥⋅∥1 \|· \|_1. That is, for a given wˇ=(ew^1,…,ew^k⏟estimated values,wk+1,…,wn⏟reference values)T, w= ( \,\,estimated values e w_1,…,e w_k, reference values\,\, w_k+1,…,w_n )^T, (6) after normalization, the final weight vector takes the form w=wˇ/‖wˇ‖1w= w/\| w\|_1. 3 Incomplete geometric HRE as a linear regression problem Due to the fact that the incomplete geometric HRE is optimal, i.e., it minimizes the value of the function ((5)). The calculation of the vector w∈ℝkw ^k boils down to solving a linear regression problem. The classic linear regression model [31, p. 2] is defined by the equation y=β0+β1x+ϵ,y= _0+ _1x+ε, where x is called the predictor or regressor variable, y is called the response variable, and ϵε is a statistical error, i.e., a random variable corresponding to the inaccuracy of the model’s fit to the data. Estimating the solution involves finding the best fit of the model to the data. Assuming that we have n data pairs of the form (yi,xi)(y_i,x_i) for i=1,…,ni=1,…,n, the least squares criterion is as follows: S(β0,β1)=∑i=1n(yi−β0−β1xi)2.S( _0, _1)=Σ^n_i=1 (y_i- _0- _1x_i )^2. Hence, the estimators β1 _1 and β2 _2 of β1 _1 and β2 _2 must satisfy the equations: ∂S∂β0|β^0,β^1=−2∑i=1n(yi−β^1xi−β^0)=0, . ∂ S∂ _0 |_ β_0, β_1=-2Σ^n_i=1 (y_i- β_1x_i- β_0 )=0, ∂S∂β1|β^0,β^1=−2∑i=1n(yi−β^1xi−β^0)xi=0. . ∂ S∂ _1 |_ β_0, β_1=-2Σ^n_i=1 (y_i- β_1x_i- β_0 )x_i=0. The above equations pave the way for the analytical calculation of the values β^0 β_0 and β^1 β_1 [31, p. 15]. The difference between the observed value yiy_i and the estimated value y^i=β^0+β^1xi y_i= β_0+ β_1x_i is the residual (estimation error) ϵi=yi−yi _i=y_i- y_i for i=1,…,ni=1,…,n. For the purposes of the model based on pairwise comparisons, let us assume that the data of the (possibly incomplete) matrix C=[cij]C=[c_ij] are independent observations. Let us denote lncij=yij c_ij=y_ij and lnwi=θi w_i= _i. With these notations, the linear regression model takes the form: yij=θi−θj+εij.y_ij= _i- _j+ _ij. (7) Denoting xij=ei−ejx_ij=e_i-e_j where eie_i is the i-th vector of the canonical basis in ℝnR^n, i.e., xijT=(0,…,1i,…,−1j,…,0)Tx^T_ij=(0,…,1_i,…,-1_j,…,0)^T, we obtain: yij=xij⊤θ+εij,y_ij=x _ijθ+ _ij, where the set of pairs R=(yij,xijT):(i,j)∈O<(C)R= \(y_ij,x^T_ij):(i,j)∈ O_<(C) \ is the data set needed to estimate the vector θ. The least squares criterion takes the form: S(θ)=∑(i,j)∈O<(C)(yij−xij⊤θ)2.S(θ)= _(i,j)∈ O_<(C) (y_ij-x _ijθ )^2. The estimator θ θ calculated as a result of minimizing the above criterion must satisfy the following equations: ∂S∂θi|θ^=0,for i=1,…,k. . ∂ S∂ _i |_ θ=0,\,\,\,for \,\,\,i=1,…,k. Let us organize the set R by adopting a certain (but arbitrary chosen) order R=(yij,1,xij,1T),(yij,2,xij,2T),…,(yij,r,xij,rT)R=\(y_ij,1,x^T_ij,1),(y_ij,2,x^T_ij,2),…,(y_ij,r,x^T_ij,r)\. Of course, the number of pairs in R is |R|=r|R|=r. Let vector y of length r and matrix X of dimensions r×nr× n have the form: y=(yij,1yij,2⋮yij,r),X=(xij,1Txij,2T⋮xij,rT).y= ( array[]cy_ij,1\\ y_ij,2\\ \\ y_ij,r array ),\,\,\,X= ( array[]cx^T_ij,1\\ x^T_ij,2\\ \\ x^T_ij,r array ). Then, let XUX_U be the matrix obtained from X by deleting the columns corresponding to the reference values, i.e., the columns with indices from IKI_K. Following the accepted convention, where AU=a1,…,akA_U= \a_1,…,a_k \, XUX_U is the matrix formed from the first k columns of the matrix X. Similarly, let XRX_R be the matrix obtained from X by deleting the columns indexed by the elements of IUI_U. Under the assumed convention, XKX_K is formed from the columns of X located at positions k+1k+1 to n. The XKX_K matrix is responsible for reference data and will be used in the model to modify the response variable, while XUX_U models variables whose estimators have yet to be determined. Let us denote y~=y−XKθK y=y-X_K _K, where θK _K is a vector consisting of logarithmic values of reference alternatives in the same order as they appear in XKX_K. Then the linear regression model takes the form: y~=XUθU+ε, y=X_U _U+ , (8) and the minimization condition of the least squares criterion is: XU⊤XUθ^U−XU⊤~=0,X _UX_U θ_U-X _U y=0, (9) where θ^U θ_U is the estimator of the value θU _U. After expanding the above equation we get XU⊤XUθ^U=XU⊤y−XU⊤XKθK.X _UX_U θ_U=X _Uy-X _UX_K _K. Therefore, to determine the value of the estimator, it is sufficient to calculate: θ^U=(XU⊤XU)−1(XU⊤y−XU⊤XKθK). θ_U= (X _UX_U )^-1 (X _Uy-X _UX_K _K ). (10) By denoting LU=XU⊤XUL_U=X _UX_U, we obtain an even more concise form of the equation: θ^U=LU−1XU⊤(y−XKθK). θ_U=L^-1_UX _U (y-X_K _K ). (11) The above equation allows for the effective calculation of the estimator θ^U θ_U, whose exponentially transformed and normalized values will become the weights of non-reference alternatives, i.e., (expθ^U)⊤=(w~1,…,w~k)⊤ ( θ_U ) = ( w_1,…, w_k ) . The matrix LUL_U is the Laplace matrix for the graph G(C)G(C) after deleting the rows and columns corresponding to the indices from the set IKI_K. It can be shown that the above equation (11) has a solution, i.e., LUL_U is non-singular, as long as each non-reference alternative is indirectly or directly compared with the reference alternative [27]. Example 5. Let us consider a model in which the set of non-reference alternatives is given as AU=a1,a2,a4,a6A_U=\a_1,a_2,a_4,a_6\ and the set of reference alternatives is AK=a3,a5,a7A_K=\a_3,a_5,a_7\ where the weights of the reference alternatives are w~(a3)=3,w~(a5)=5 w(a_3)=3, w(a_5)=5, and w~(a7)=7 w(a_7)=7, respectively. Thus, IU=1,2,4,6I_U=\1,2,4,6\ and IK=3,5,7I_K=\3,5,7\. Furthermore, let us assume that the pairwise comparison matrix C=[cij]C=[c_ij] looks as follows: C=(10.161?0.1590.2490.155?6.2161?0.340?0.569???11.06335?376.2572.9340.94010.7350.520?4?531.35912.181576.4431.757?1.9190.45811.033??73?750.9671).C= ( array[]c1&0.161&?&0.159&0.249&0.155&?\\ 6.216&1&?&0.340&?&0.569&?\\ ?&?&1&1.063& 35&?& 37\\ 6.257&2.934&0.940&1&0.735&0.520&?\\ 4&?& 53&1.359&1&2.181& 57\\ 6.443&1.757&?&1.919&0.458&1&1.033\\ ?&?& 73&?& 75&0.967&1 array ). The values of cijc_ij where i,j∈IKi,j∈ I_K are known a priori, or are determined by an expert i.e. cij∈O(C)c_ij∈ O(C), or remain undefined i.e. (cij=?c_ij=? ). The set of index pairs determining the observations is O<(C)=(1,2)O_<(C)=\(1,2), (1,4)(1,4), (1,5)(1,5), (1,6)(1,6), (2,4)(2,4), (2,6)(2,6), (3,4)(3,4), (4,5)(4,5), (4,6)(4,6), (5,6)(5,6), (6,7)(6,7)\. Since |O<(C)|=11|O_<(C)|=11, the design matrix X consists of 1111 rows, each corresponding to one observation (one pairwise comparison under consideration), and 77 columns (the total number of alternatives). X=(1−100000100−10001000−10010000−10010−100001000−10001−10000001−10000010−1000001−10000001−1)=(x12Tx14Tx15Tx16Tx24Tx26Tx34Tx45Tx46Tx56Tx67T).X= ( array[]c1&-1&0&0&0&0&0\\ 1&0&0&-1&0&0&0\\ 1&0&0&0&-1&0&0\\ 1&0&0&0&0&-1&0\\ 0&1&0&-1&0&0&0\\ 0&1&0&0&0&-1&0\\ 0&0&1&-1&0&0&0\\ 0&0&0&1&-1&0&0\\ 0&0&0&1&0&-1&0\\ 0&0&0&0&1&-1&0\\ 0&0&0&0&0&1&-1 array )= ( array[]cx^T_12\\ x^T_14\\ x^T_15\\ x^T_16\\ x^T_24\\ x^T_26\\ x^T_34\\ x^T_45\\ x^T_46\\ x^T_56\\ x^T_67 array ). Let us denote the matrix XU=[:,IU]X_U=[:,I_U] as the matrix X with the columns corresponding to reference alternatives removed (11, 22, 44 and 66 are left), and let XK=[:,IK]X_K=[:,I_K] be the matrix X with the columns corresponding to non-reference alternatives removed (33, 55, and 77 left). Thus, XU=(1−10010−101000100−101−10010−100−100010001−1000−10001),XK=(0000000−100000000001000−1000001000−1).X_U= ( array[]c1&-1&0&0\\ 1&0&-1&0\\ 1&0&0&0\\ 1&0&0&-1\\ 0&1&-1&0\\ 0&1&0&-1\\ 0&0&-1&0\\ 0&0&1&0\\ 0&0&1&-1\\ 0&0&0&-1\\ 0&0&0&1 array ),\,\,\,X_K= ( array[]c0&0&0\\ 0&0&0\\ 0&-1&0\\ 0&0&0\\ 0&0&0\\ 0&0&0\\ 1&0&0\\ 0&-1&0\\ 0&0&0\\ 0&1&0\\ 0&0&-1 array ). Matrix LU=XUTXUL_U=X^T_UX_U (being the matrix C C, eq. (3)) and corresponding LU−1L^-1_U are as follows: LU=(4−1−1−1−13−1−1−1−15−1−1−1−15),LU−1=(51331321321331371352652621352623785392135265392378).L_U= ( array[]c4&-1&-1&-1\\ -1&3&-1&-1\\ -1&-1&5&-1\\ -1&-1&-1&5 array ),\,\,\,L^-1_U= ( array[]c 513& 313& 213& 213\\ 313& 713& 526& 526\\ 213& 526& 2378& 539\\ 213& 526& 539& 2378 array ). The remaining components of the equation used to determine the estimator, i.e., y and θK _K, are as follows: y=(lnc1,2lnc1,3lnc1,4lnc1,6lnc2,4lnc2,6lnc3,4lnc4,5lnc4,6lnc5,6lnc5,7)=(−1.82719−1.83379−1.38636−1.86307−1.0765−0.5636440.061674−0.307264−0.652260.7798820.0332809),θK=(lnw(a3)lnw(a5)lnw(a7))=(ln3ln5ln7).y= ( array[]c c_1,2\\ c_1,3\\ c_1,4\\ c_1,6\\ c_2,4\\ c_2,6\\ c_3,4\\ c_4,5\\ c_4,6\\ c_5,6\\ c_5,7 array )= ( array[]c-1.82719\\ -1.83379\\ -1.38636\\ -1.86307\\ -1.0765\\ -0.563644\\ 0.061674\\ -0.307264\\ -0.65226\\ 0.779882\\ 0.0332809 array ),\,\,\, _K= ( array[]c w(a_3)\\ w(a_5)\\ w(a_7) array )= ( array[]c 3\\ 5\\ 7 array ). Combining all components of the solution into a single equation θ^U=LU−1XU⊤(y−XRθR) θ_U=L^-1_UX _U (y-X_R _R ), we obtain: θ^U=(−0.3826150.8937331.330841.54594). θ_U= ( array[]c-0.382615\\ 0.893733\\ 1.33084\\ 1.54594 array ). Hence, the unnormalized weight vector w~ w for the non-reference alternatives is expθ^U=(e0.682,e2.444,e3.784,e4.693)T θ_U=(e^0.682,e^2.444,e^3.784,e^4.693)^T i.e. w~1=e0.682 w_1=e^0.682, w~2=e2.444 w_2=e^2.444, w~4=e3.784 w_4=e^3.784 and w~6=e4.693 w_6=e^4.693 (and of course w~3=3 w_3=3, w~5=5 w_5=5 and w~7=7 w_7=7). Finally, after normalization, taking into account the reference alternatives, we obtain w=w~1Tw~=(0.0256,0.092,0.112,0.142,0.187,0.176,0.263)T.w= w1^T w= (0.0256,0.092,0.112,0.142,0.187,0.176,0.263 )^T. (12) 4 A statistical perspective 4.1 Standard error and variance of estimation In the pairwise comparison method, the estimation error of the weight vector is most often considered as δij _ij where cij=wiwjδij.c_ij= w_iw_j _ij. In the literature, it is usually assumed that δij _ij has a log-normal distribution [28, 16, 8, 3, 15, 36, 35]. Therefore, after logarithmization, the distribution ϵij=lnδij _ij= _ij becomes a normal distribution with mean 0 and variance σ2σ^2, i.e., ηij∼(0,σ2) _ij (0,σ^2). In our model, we have yij=θi−θj+ϵijy_ij= _i- _j+ _ij, thus ϵ^ij=yij−(θ^i−θ^j) ε_ij=y_ij- ( θ_i- θ_j ). The sum of the squares of the residuals (errors) is given as SSRes=∑(i,j)∈O<(C)ϵ^ij 2.S_SRes= _(i,j)∈ O_<(C) ε^\,2_ij. (13) Let us denote the number of observations as |O<(C)|=r|O_<(C)|=r. Since the number of estimated parameters is k=|AU|k= |A_U | (the number of non-referential alternatives), the number of degrees of freedom is r−kr-k. Thus, the unbiased estimator of variance σ^2 σ^2 (σ σ is the standard error of regression) is given [31, p. 80 - 81] as: σ^2=SSResr−k. σ^2= S_SResr-k. (14) Let us go back to the equation for a moment. ((8)). The equation θ^U θ_U of the estimator using the variable y~ y is θ^U=LU−1XU⊤y~. θ_U=L^-1_UX _U y. By transforming the above equation, we obtain: θ^U−θU=LU−1XU⊤ε. θ_U- _U=L^-1_UX _U . From the definition of covariance matrix Cov(θ^U)=Cov(θ^U−[θ^U])Cov( θ_U)=Cov( θ_U-E[ θ_U]), and since the estimator θ^U θ_U is unbiased, i.e., E[θ^U]=θUE[ θ_U]= _U, then Cov(θ^U)=Cov(θ^U−θU)Cov( θ_U)=Cov( θ_U- _U). Thus, Cov(θ^U)=Cov(LU−1XU⊤ε).Cov( θ_U)=Cov (L^-1_UX _U ). Now, from the covariance property of linear mapping (covariance matrix property), we have that if Z=AεZ=A , then Cov(Z)=ACov(ε)A⊤Cov(Z)=ACov( )A . Taking the expression LU−1XU⊤εL^-1_UX _U as matrix A, we obtain Cov(θ^U)=LU−1XU⊤Cov(ε)XULU−1.Cov( θ_U)=L^-1_UX _UCov( )X_UL^-1_U. Because Cov(ε)=σ2ICov( )=σ^2I hence, assuming the estimated value σ^2 σ^2 Cov(θ^U)=σ^2LU−1XU⊤XU⏟LULU−1=σ^2LU−1.Cov( θ_U)= σ^2L^-1_U L_U X _UX_UL^-1_U= σ^2L^-1_U. (15) The covariance matrix allows us to determine the variance values for individual elements of θ i.e. Var(θ^i)=σ^2[LU−1]iiVar( θ_i)= σ^2 [L^-1_U ]_i, and the value of their mutual covariance Cov(θ^i,θ^j)=σ^2[LU−1]ijCov( θ_i, θ_j)= σ^2 [L^-1_U ]_ij. 4.2 Quantitative interpretation and confidence intervals The weights for individual alternatives can be interpreted quantitatively or qualitatively. In the quantitative approach, it is not so much the position of a given alternative in the ranking that is important, but the value of the weight itself. It is easy to imagine that such a weight, multiplied by 100%100\%, represents the percentage share of a given alternative in the reward, or the number of loyalty program points we can use later. In such a situation, one may ask about the confidence interval within which the weight of the alternative of interest lies. The limits of this interval allow us to estimate our potential loss or gain that may result from an incorrect estimate. If the confidence interval at a given level of certainty is wide, i.e., the potential loss may be large, this may be a reason to question such a weight vector. To determine the confidence interval, we need to calculate the standard deviation (standard error): SE(θ^i)=Var(θ^i)=σ^[LU−1]iiSE( θ_i)= Var( θ_i)= σ [L^-1_U ]_i. Therefore, we can state that in a model with normally distributed errors, θi _i belongs to the confidence interval defined by the formula: θi∈[θ^i−t1−α/2,r−kSE(θ^i),θ^i+t1−α/2,r−kSE(θ^i)], _i∈ [ θ_i-t_1-α/2,r-k\,SE( θ_i),\; θ_i+t_1-α/2,r-k\,SE( θ_i) ], where t1−α/2,νt_1-α/2,ν - is the value quantile from the Student’s t-distribution. In the above formula, α is the probability that θi _i does not belong to the specified interval, and r−kr-k is the number of degrees of freedom of the model. The value of t1−α/2,νt_1-α/2,ν can be taken from mathematical tables or calculated numerically222Formally, t1−α/2,νt_1-α/2,ν is a number that, for the Student’s t-distribution given as Ftν(t)=∫−∞tΓ(ν+12)νπΓ(ν2)(1+u2ν)−ν+12uF_t_ν(t)= ^t_-∞ \! ( ν+12 ) νπ\, \! ( ν2 ) (1+ u^2ν )^- ν+12du, satisfies the equation Ftν(t1−α/2,ν)=1−α/2.F_t_ν(t_1-α/2,ν)=1-α/2. Fortunately, most tools have a built-in function for calculating the value of t1−α/2,νt_1-α/2,ν.. The confidence interval for alternative weights in the initial model with a multiplicative pairwise comparisons matrix is obtained through exponential transformation. Thus, wi=exp(θ^i±t1−α/2,r−kSE(θ^i)).w_i= ( θ_i± t_1-α/2,r-k\,SE( θ_i) ). Example 6. For the model considered in the example (5) the sum of the squares of the logarithmic errors ((13)) is SSResS_SRes == (−1.827−(0.682−2.444))2+…=(-1.827-(0.682-2.444))^2+…= 2.084952.08495. Since |O<(C)|=r=11|O_<(C)|=r=11 and k=4k=4, then σ^2=2.08495/(11−4)=0.2978 σ^2=2.08495/(11-4)=0.2978. The estimated covariance matrix Cov(θ^U)Cov( θ_U) takes the form: Cov(θ^U)=σ^2LU−1=(0.11450.06870.04580.04580.06870.16040.05730.05730.04580.05730.08780.03820.04580.05730.03820.0878).Cov( θ_U)= σ^2L^-1_U= ( array[]c0.1145&0.0687&0.0458&0.0458\\ 0.0687&0.1604&0.0573&0.0573\\ 0.0458&0.0573&0.0878&0.0382\\ 0.0458&0.0573&0.0382&0.0878 array ). (16) Based on the above matrix, we determine Var(θ^1)=0.1145Var( θ_1)=0.1145, Var(θ^2)=0.1604Var( θ_2)=0.1604, and Var(θ^4)=Var(θ^6)=0.0878Var( θ_4)=Var( θ_6)=0.0878 and, accordingly, the standard deviation SE(θ^1)=0.3384SE( θ_1)=0.3384, SE(θ^2)=0.4SE( θ_2)=0.4, SE(θ^4)=0.2963SE( θ_4)=0.2963, and SE(θ^4)=0.2963SE( θ_4)=0.2963. Let us assume a popular five percent confidence level, i.e., α=0.05α=0.05. Therefore, after calculation (checking in the tables), we obtain t1−0.05/2,7=t0.975,7=2.3646t_1-0.05/2,7=t_0.975,7=2.3646. This allows us to determine confidence intervals at a confidence level of 95%95\% for the θ θ. After performing the calculations, we obtain θ1∈[−1.183,0.4177],θ2∈[−0.0532,1.841], _1∈ [-1.183,0.4177 ],\,\,\,\, _2∈ [-0.0532,1.841 ], θ4∈[0.6301,2.0316],θ6∈[0.8452,2.2467]. _4∈ [0.6301,2.0316 ],\,\,\,\,\, _6∈ [0.8452,2.2467 ]. To obtain confidence intervals for a non-logarithmic weight vector, simply use an exponential transformation. After transformation, we obtain w~1∈[0.3125,1.489] w_1∈ [0.3125,1.489 ], w~2∈[0.9707,6.1548] w_2∈ [0.9707,6.1548 ], w~4∈[1.911,7.495] w_4∈ [1.911,7.495 ], and w~6∈[2.3691,9.2937] w_6∈ [2.3691,9.2937 ]. These values correspond to the weight vector before normalization. For normalization, these intervals need to be divided by 1Tw~1^T w, i.e., by the sum of the elements w~ w ((12)). Thus, after normalization, we ultimately obtain w1∈[0.0117,0.0559],w2∈[0.0365,0.2313], w_1∈ [0.0117,0.0559 ],\,\,\,\,w_2∈ [0.0365,0.2313 ], w4∈[0.0718,0.2817],w6∈[0.089,0.3493]. w_4∈ [0.0718,0.2817 ],\,\,\,\,w_6∈ [0.089,0.3493 ]. By lowering the confidence threshold (i.e., increasing the value of α), we can narrow the confidence intervals accordingly. 4.3 Qualitative interpretation and probability of rank reversal Although the weight vector resulting from the priority deriving method is quantitative in nature, its outcome is often interpreted qualitatively. That is, the order of alternatives, which results from ranking them according to their decreasing weight values, is taken into account. With this approach, we are more interested in the order of the alternatives than in the specific weight values. It is often not all the positions in the ranking that matter, but only the top few, and it is only their order that may be the subject of dispute among the stakeholders. Adopting a qualitative perspective entails asking a series of questions regarding the validity of the resulting ranking. More specifically, these questions concern the probability, given the decision-making data, that the estimated ranking is correct. In its simplest form, this question boils down to determining the probability that one alternative is more (or less) preferred than another. 4.3.1 A pair of alternatives Let us therefore calculate, for two estimators θ^i θ_i and θ^j θ_j in question, the probability that the order relation between the two random variables they represent (the true but unknown weights) is given. Since the estimators θ^i θ_i and θ^j θ_j are correlated, this correlation, i.e., the matrix Cov(θ^U)Cov( θ_U), must be taken into account in the calculations. To this end, let us consider the difference δ=θi−θjδ= _i- _j and, in practice, its estimated value δ^=θ^i−θ^j δ= θ_i- θ_j. Since for any random variables ξ and ζ, we have Var(ξ−ζ)=Var(ξ)+Var(ζ)−2Cov(ξ,ζ)Var(ξ-ζ)=Var(ξ)+Var(ζ)-2\,Cov(ξ,ζ), and recalling that Cov(θ^i,θ^j)=σ2vijCov( θ_i, θ_j)=σ^2v_ij where333It should be noted that the successive rows and columns in the matrix LU−1L^-1_U correspond to successive non-reference alternatives. Therefore, if we index the reference and non-reference alternatives together, as is the case here, we must remember to remap the indices appropriately when referring to the matrix LU−1L^-1_U. In particular, this means that only those elements of [LU−1]p,q [L^-1_U ]_p,q are well-defined for which p,q∈IUp,q∈ I_U. For example, if IU=1,2,4,6I_U=\1,2,4,6\, then the element [LU−1]6,6 [L^-1_U ]_6,6 actually lies at the intersection of the fourth column and the fourth row of the matrix LU−1L^-1_U. vij=[LU−1]ijv_ij= [L^-1_U ]_ij, we obtain variance estimate: Var(δ^)=σ^2(vii+vjj−2vij).Var( δ)= σ^2 (v_i+v_j-2v_ij ). Noting that δ∼(δ^,σ^2(vii+vjj−2vij))δ ( δ, σ^2(v_i+v_j-2v_ij)) the question of the probability that θi<θj _i< _j can be reduced to the question of whether δ<0δ<0. Since if ξ∼(μ,σ2)ξ (μ,σ^2), then the variable ζ=ξ−μσζ= ξ-μσ has a distribution of ζ∼(0,1)ζ (0,1). Thus, after performing this standardization we obtain P(θi<θj∣data)=P(δ<0∣data)≈Ftr−k(0−δ^σ^vii+vjj−2vij),P( _i< _j )=P(δ<0 )≈ F_t_r-k ( 0- δ σ v_i+v_j-2v_ij ), that is, P(θi<θj∣data)≈Ftr−k(θ^j−θ^iσ^vii+vjj−2vij),P( _i< _j )≈ F_t_r-k ( θ_j- θ_i σ v_i+v_j-2v_ij ), where F is the cumulative distribution function of the t-student distribution. For large samples, i.e., when the sample size is sufficiently large, F can be replaced by the cumulative distribution function of the normal distribution Φ (for larger values, the distribution functions are nearly identical). Since our starting point was the two estimators θ^i θ_i and θ^j θ_j, this means that the above formula applies to two non-reference alternatives. That is, given the monotonicity of the exp transformation, we obtain that P(ai≺aj∣data)≈Ftr−k(θ^j−θ^iσ^vii+vjj−2vij),P(a_i a_j )≈ F_t_r-k ( θ_j- θ_i σ v_i+v_j-2v_ij ), (17) where ai,aj∈AUa_i,a_j∈ A_U. If one of the alternatives under consideration is the reference alternative, i.e., aj∈AKa_j∈ A_K, meaning that θj=θ^j=lnw(aj) _j= θ_j= w(a_j), this implies that there is no relationship between the variables under consideration. Thus, given the data, θi≈(θ^i,σ^2vii) _i ( θ_i,\ σ^2v_i). Thus, after standardization, the formula takes the form P(θi<θj∣data)≈Ftr−k(θj−θ^iσ^vii).P( _i< _j )≈ F_t_r-k ( _j- θ_i σ v_i ). Thus P(ai≺aj∣data)≈Ftr−k(θj−θ^iσ^vii),P(a_i a_j )≈ F_t_r-k ( _j- θ_i σ v_i ), (18) where ai∈AUa_i∈ A_U and aj∈AKa_j∈ A_K. If, on the other hand, ai∈AKa_i∈ A_K and aj∈AUa_j∈ A_U, then P(ai≺aj∣data)≈Ftr−k(θ^j−θiσ^vjj).P(a_i a_j )≈ F_t_r-k ( θ_j- _i σ v_j ). (19) Ultimately, for the two reference alternatives P(ai≺aj)=0θi≥θj1θi<θj,P(a_i a_j)= cases0& _i≥ _j\\ 1& _i< _j cases, (20) where ai,aj∈AKa_i,a_j∈ A_K. The above formulas (17), (18), (19) and (20) allow us to estimate the probability of a specific order for any pair of alternatives. 4.3.2 Three alternatives This reasoning can be extended to a larger number of alternatives. For example, for three non-reference alternatives aia_i, aja_j, and aka_k, we can calculate the probability that ai≺ak≺aja_i a_k a_j. That is, for example, to verify the reliability with which the three “top spots” in the ranking were correctly identified. To calculate the probability of a given sequence occurring P(ai≺ak≺aj|data)=P(θi<θk<θj|data)P(a_i a_k a_j|\, data)=P( _i< _k< _j\,|\, data) we need to calculate the probability that θk−θi>0 _k- _i>0 and θj−θk>0 _j- _k>0. Let us define the difference vector δ as follows: δ=[δ1δ2]=[θk−θiθj−θk],δ= bmatrix _1\\ _2 bmatrix= bmatrix _k- _i\\ _j- _k bmatrix, which can be written as δ=Aθδ=Aθ, where the matrix A is given by: A=[(ek−ei)⊤(ej−ek)⊤].A= bmatrix(e_k-e_i) \\ (e_j-e_k) bmatrix. For the purposes of our calculations, we will, of course, use estimates of δ, i.e. δ^=[δ^1δ^2]=[θ^k−θ^iθ^j−θ^k]. δ= bmatrix δ_1\\ δ_2 bmatrix= bmatrix θ_k- θ_i\\ θ_j- θ_k bmatrix. (21) From (15) we get that Cov(δ^)=ACov(θ^)A⊤=σ^2ALU−1A⊤Cov( δ)=A\,Cov( θ)\,A = σ^2AL^-1_UA . Let us denote Σδ:=ALU−1A⊤ _δ:=AL^-1_UA . Then Cov(δ^)=σ^2Σδ.Cov( δ)= σ^2 _δ. The matrix Σδ _δ takes the form: Σδ=[(ek−ei)⊤LU−1(ek−ei)(ek−ei)⊤LU−1(ej−ek)(ej−ek)⊤LU−1(ek−ei)(ej−ek)⊤LU−1(ej−ek)], _δ= bmatrix(e_k-e_i) L^-1_U(e_k-e_i)&(e_k-e_i) L^-1_U(e_j-e_k)\\[4.0pt] (e_j-e_k) L^-1_U(e_k-e_i)&(e_j-e_k) L^-1_U(e_j-e_k) bmatrix, (22) i.e. Σδ=[vkk+vii−2vikvkj−vkk−vij+vikvkj−vkk−vij+vikvjj+vkk−2vjk]. _δ= bmatrixv_k+v_i-2v_ik&v_kj-v_k-v_ij+v_ik\\[4.0pt] v_kj-v_k-v_ij+v_ik&v_j+v_k-2v_jk bmatrix. where vij=[LU−1]ijv_ij= [L^-1_U ]_ij. Hence, the distribution of the difference vector δ is (δ^,σ^2Σδ).N ( δ,\ σ^2 _δ ). Let us define the standardized variables Z1=δ1−δ^1σ^Σδ,11,Z2=δ2−δ^2σ^Σδ,22.Z_1= _1- δ_1 σ _δ,11, Z_2= _2- δ_2 σ _δ,22. The variables (Z1,Z2)(Z_1,Z_2) follow a two-dimensional standard normal distribution with a correlation of ρ=Σδ,12/Σδ,11Σδ,22ρ= _δ,12/ _δ,11 _δ,22. From the given inequalities δ1>0 _1>0 and δ2>0 _2>0, we obtain: δ1−δ^1σ^Σδ,11>−δ^1σ^Σδ,11,δ2−δ^2σ^Σδ,22>−δ^2σ^Σδ,22, _1- δ_1 σ _δ,11> - δ_1 σ _δ,11, _2- δ_2 σ _δ,22> - δ_2 σ _δ,22, and then Z1>−δ^1σ^Σδ,11,Z2>−δ^2σ^Σδ,22.Z_1>- δ_1 σ _δ,11, Z_2>- δ_2 σ _δ,22. So, we get P(δ1>0,δ2>0)=P(Z1>−δ^1σ^Σδ,11,Z2>−δ^2σ^Σδ,22).P( _1>0, _2>0)=P\! (Z_1>- δ_1 σ _δ,11,\;Z_2>- δ_2 σ _δ,22 ). Given the symmetry of the Student’s t-distribution, i.e., that P(Z>a)=1−Ftr−k(a)=Ftr−k(−a)P(Z>a)=1-F_t_r-k(a)=F_t_r-k(-a), and bearing in mind that P(ai≺ak≺aj|data)=P(θi<θk<θj|data)=P(δ1>0,δ2>0)P(a_i a_k a_j|\, data)=P( _i< _k< _j\,|\, data)=P( _1>0, _2>0) we obtain: P(ai≺ak≺aj|data)=Ftr−k(2)(δ^1σ^Σδ,11,δ^2σ^Σδ,22,ρ),P(a_i a_k a_j|\, data)=F^(2)_t_r-k ( δ_1 σ _δ,11, δ_2 σ _δ,22,ρ ), where ρ=Σδ,12/Σδ,11Σδ,22ρ= _δ,12/ _δ,11 _δ,22, ai,ak,aj∈AUa_i,a_k,a_j∈ A_U and Ftr−k(2)F^(2)_t_r-k is the bivariate t-Student’s cumulative distribution function. Of course, when one of the alternatives under consideration is the reference alternative, the uncertainty structure changes. One of the variables in the differences ceases to be random. As a result, the variances change, and the correlation changes. For example, assuming that aka_k is a reference alternative, i.e., ak∈AKa_k∈ A_K, we have that Var(δ1)=Var(θi)Var( _1)=Var( _i), Var(δ2)=Var(θj)Var( _2)=Var( _j) that is, the variance of the differences depends on only one variable. As a result Σδ,11=vii _δ,11=v_i and Σδ,22=vjj. _δ,22=v_j. Similarly, the covariance Cov(δ1,δ2)=Cov(−θi,θj)=−Cov(θi,θj)Cov( _1, _2)=Cov(- _i, _j)=-\,Cov( _i, _j), i.e. Σδ,12=−vij _δ,12=-v_ij. Therefore, the correlation, that is, the covariance divided by the standard deviations, is expressed by the formula: ρ=−vij/viivjj.ρ=-v_ij/ v_iv_j. This leads to the conclusion: P(ai≺ak≺aj|data)=P(θi<θk<θj∣data)=Ftr−k(2)(δ^1σ^vii,δ^2σ^vjj,ρ),P(a_i a_k a_j|\, data)=P( _i< _k< _j )=F^(2)_t_r-k ( δ_1 σ v_i, δ_2 σ v_j,ρ ), where ai,aj∈AUa_i,a_j∈ A_U, ak∈AKa_k∈ A_K, δ^1=θk−θ^i δ_1= _k- θ_i, and δ^2=θ^j−θk δ_2= θ_j- _k. When we are dealing with two reference alternatives, the differences δ1 _1 and δ2 _2 become functions of a single random variable. For example, let us consider a situation in which we want to estimate the probability that aka_k “has been caught” by two reference variables, aia_i and aja_j. From the model’s assumptions, we have θk≈(θ^k,σ^2vkk). _k ( θ_k, σ^2v_k). For the sake of standardization, let us put Z=(θk−θ^k)/σ^vkk∼(0,1)Z=( _k- θ_k)/ σ v_k (0,1). Then P(θi<θk<θj)=P(θi−θ^kσ^vkk<Z<θj−θ^kσ^vkk),P( _i< _k< _j)=P\! ( _i- θ_k σ v_k<Z< _j- θ_k σ v_k ), i.e. P(ai≺ak≺aj|data)=Ftr−k(θj−θ^kσ^vkk)−Ftr−k(θi−θ^kσ^vkk),P(a_i a_k a_j|\, data)=F_t_r-k ( _j- θ_k σ v_k )-F_t_r-k ( _i- θ_k σ v_k ), where ai,aj∈AKa_i,a_j∈ A_K and ak∈AUa_k∈ A_U. The reasoning presented here can be scaled to any number of alternatives by appropriately increasing the vector δ, the matrix Σδ _δ, and extending ρ to the correlation matrix between the individual random variables. Example 7. For the model considered in the previous examples (5 and 6), given the LU−1L^-1_U matrix we can easily calculate the probability, for example, that a1≺a2a_1 a_2. From 17, we therefore have P(a1≺a2|data)≈Ft11−4(θ^2−θ^1σ^v1,1+v2,2−2v1,2)=P(a_1 a_2|\, data)≈ F_t_11-4 ( θ_2- θ_1 σ v_1,1+v_2,2-2v_1,2 )= =Ft7(0.8937−(−0.3826)0.5457513+713−2⋅313)=0.9946.=F_t_7 ( 0.8937- (-0.3826 )0.5457 513+ 713-2· 313 )=0.9946. Since the actual weights calculated for alternatives a1a_1 and a2a_2 satisfy the condition w1<w2w_1<w_2, the value of P(a1≺a2)P(a_1 a_2) being close to 11 suggests that, based on the data (the collected pairwise comparisons), we can be fairly confident in this result. The situation is quite different for alternatives a5a_5 and a6a_6. Although w6=0.176<0.186=w5w_6=0.176<0.186=w_5, P(a6≺a5)=0.5818P(a_6 a_5)=0.5818. This relatively low probability value may indicate that, for this particular pair of alternatives, the calculated order may differ from the actual result. For three alternatives, it is also possible to calculate the probability of a given order. Let us consider three non-reference alternatives a1a_1, a2a_2, and a4a_4. To calculate, for example, P(a1≺a2≺a4)P(a_1 a_2 a_4), following 21 and 22, we first obtain: δ^=[δ^1δ^2]=[θ^k−θ^iθ^j−θ^k]=[1.2760.437],andΣδ=(613−726−7263578). δ= bmatrix δ_1\\ δ_2 bmatrix= bmatrix θ_k- θ_i\\ θ_j- θ_k bmatrix= bmatrix1.276\\ 0.437 bmatrix,\,\, and\,\, _δ= ( array[]c 613&- 726\\ - 726& 3578 array ). With these values, we can easily calculate the correlation coefficient σ and determine the probability: P(a1≺a2≺a4)=Ft11−4(2)(δ^1σ^Σδ,11,δ^2σ^Σδ,22,ρ)=P(a_1 a_2 a_4)=F^(2)_t_11-4 ( δ_1 σ _δ,11, δ_2 σ _δ,22,ρ )= Ft7(2)(3.442,1.195,−0.5916)=0.8593.F^(2)_t_7 (3.442,1.195,-0.5916 )=0.8593. This way, we can calculate the probability of each arrangement of the three alternatives. For example P(a2≺a1≺a4)=0.0051P(a_2 a_1 a_4)=0.0051 and P(a1≺a4≺a2)=0.1348P(a_1 a_4 a_2)=0.1348 etc. 5 Quality of the weight vector 5.1 Quality indicators The probability that, for two alternatives aia_i and aja_j, the first precedes the second – i.e., P(ai≺aj∣data)P(a_i a_j ) – contains information about both the inconsistency of the pairwise comparison matrix and the preference distance between these two alternatives. We can expect (we will confirm this experimentally later) that the greater the inconsistency of the pairwise comparison matrix C, the smaller the expected value of P(ai≺aj)P(a_i a_j), and the greater the preference distance between aia_i and aja_j (the more the weights wiw_i and wjw_j differ), the higher the probability P(ai≺aj)P(a_i a_j) is. In particular, in practice it may happen that for a pairwise comparison matrix with relatively high inconsistency (CR(C)≫0.1CR(C) 0.1), for a pair of alternatives that are preferentially distant from each other, their mutual ordering relationship can be determined with a high degree of certainty, whereas in a matrix with acceptable inconsistency (CR(C)<0.1CR(C)<0.1), for two alternatives with very similar weights, their mutual ordering relationship may prove uncertain. Therefore, the order probabilities for pairs of alternatives, which are relatively easy to calculate, can be used to determine the quality of the resulting weight vector. To introduce quality measures for a weight vector, let us define the following concepts. Let Rank(C)=(i1,i2,…,in)Rank(C)=(i_1,i_2,…,i_n) where wiq≤wiq+1w_i_q≤ w_i_q+1 will be referred to as the alternative ranking for the extended weight vector (6) given as w=(w1,w2,…,wn)Tw= (w_1,w_2,…,w_n )^T. In other words, the Rank(C)Rank(C) contains a list of indices of alternatives ranked from the least to the most preferred. Furthermore, let PSPPSP be a list of probabilities of maintaining order for pairs of alternatives, i.e., PSP(C)=dfP(ai⪯aj):wi≤wj∧i≠j. PSP(C) df= \P (a_i a_j ):w_i≤ w_j i≠ j \. Similarly, let us denote the restriction of the set PSP(C) PSP(C) to pairs of indices from the set G as PSPG(C)=dfP(ai⪯aj):wi≤wj∧i≠j∧(i,j)∈G. PSP_G(C) df= \P (a_i a_j ):w_i≤ w_j i≠ j (i,j)∈ G \. Thus, PSPIU×IU(C) PSP_I_U× I_U(C) is the set of probabilities of the order of pairs formed from non-reference alternatives. Definition 8. Let lcPOI(C,G)=dfminPSPG(C) lcPOI(C,G) df= PSP_G(C) is said to be the least certain pairwise order index. The value of lcPOIU(C)=dflcPOI(C,IU×IU) lcPOI_U(C) df= lcPOI(C,I_U× I_U) indicates the least certain preference relationship of two alternatives for which a weight vector has been computed. That is, the probability that, in light of the collected data (the results of pairwise comparisons), the indicated preference relationship actually exists. The value 1−lcPOIU(C)1- lcPOI_U(C) can be interpreted as the probability of a reversal in the ranking. That is, the value 1−lcPOIU(C)1- lcPOI_U(C) means the probability that the order of the alternatives resulting from the calculated weight vector is, in fact, different. Additionally, let us introduce an indicator of the average order probability for pairs of alternatives. Definition 9. Let alPOI(C,G)=dfmeanPSPG(C) alPOI(C,G) df=mean\,\, PSP_G(C) is said to be the average likelihood pairwise order index. In other words, the value of, for example, alPOIU(C)=dfalPOI(C,IU×IU) alPOI_U(C) df= alPOI(C,I_U× I_U) is the probability that any randomly selected pair of non-reference alternatives is, in fact, in the same order relation as that implied by the computed weight vector. The lcPOI and alPOI indicators can also be used in the context of reference alternatives. In this case, we will want to determine the values of lcPOIK(C)=dflcPOI(C,IU×IK∪IK×IU) lcPOI_K(C) df= lcPOI(C,I_U× I_K∪ I_K× I_U) and alPOIK(C)=dfalPOI(C,IU×IK∪IK×IU) alPOI_K(C) df= alPOI(C,I_U× I_K∪ I_K× I_U), respectively. When considering the set IU×IK∪IK×IUI_U× I_K∪ I_K× I_U, we are interested in the probability that the preference relations between the reference alternatives and the non-reference alternatives are preserved. In addition to their probabilistic interpretation, the calculated values can also be understood as confidence indicators regarding the extent to which the calculated ranks fit the ranks of reference alternatives. Finally, let us define lcPOIUK(C)=dflcPOI(C,IU×IK∪IK×IU∪IU×IU) lcPOI_UK(C) df= lcPOI(C,I_U× I_K∪ I_K× I_U∪ I_U× I_U) and, accordingly, alPOIUK(C)=dfalPOI(C,IU×IK∪IK×IU∪IU×IU) alPOI_UK(C) df= alPOI(C,I_U× I_K∪ I_K× I_U∪ I_U× I_U). The last two indices take into account all pairs of alternatives in which at least one is non-reference. 5.2 The tie threshold and tie clustering The indicators defined above allow us to assess the reliability of preference relationships for a set of alternatives. In a situation where this assessment is unfavorable — i.e., for example, the lcPOI value for a selected pair is low — the question arises regarding the reliability of the resulting weight vector. One solution may be to attempt to modify the expert assessment, increasing the consistency of the data so that the collected data more unambiguously indicates the winner in comparisons of pairs (ai,aj)(a_i,a_j) with an unsatisfactory probability P(ai≺aj)P(a_i a_j). However, this may be difficult for procedural reasons and costly due to the need to pay for additional expert (or experts’) working time. Another approach is to group all those alternatives whose relative order may be disputed into clusters. Let us consider a set of alternatives A=a1,…,anA= \a_1,...,a_n \ with weights w1,…,wnw_1,...,w_n. The weight values determine the order of the alternatives, i.e. wi1≤wi2≤…≤winw_i_1≤ w_i_2≤…≤ w_i_n implies ai1⪯ai2⪯…⪯aina_i_1 a_i_2 … a_i_n. Definition 10. The tied pairs will be all those pairs (aik,aik+1)(a_i_k,a_i_k+1) for which wik≤wik+1w_i_k≤ w_i_k+1 and P(aik≺aik+1)<δP(a_i_k a_i_k+1)<δ. The coefficient δ will be referred to as the tie threshold. Tie clusters should group together those alternatives for which the probability of ordinal concordance is less than the tie threshold δ. Thus, let us define a partition of the set of alternatives. Definition 11. A set Q=Q1,…,QrQ= \Q_1,…,Q_r \ is called a tie partition of A if for any two QgQ_g and QhQ_h, Qg∩Qh=∅Q_g∩ Q_h= , ⋃Q=A Q=A, and for all aik,aik+1∈Qs∈Qa_i_k,a_i_k+1∈ Q_s∈ Q such that wik≤wik+1w_i_k≤ w_i_k+1 holds P(aik≺aik+1)<δP(a_i_k a_i_k+1)<δ. It may happen that there are multiple valid tie partitions for A. For example, for three alternatives aik≺aik+1≺aik+2a_i_k a_i_k+1 a_i_k+2, it may be that P(aik≺aik+1)<δP(a_i_k a_i_k+1)<δ, P(aik+1≺aik+2)<δP(a_i_k+1 a_i_k+2)<δ but P(aik≺aik+2)>δP(a_i_k a_i_k+2)>δ. Therefore, both, the partition in which aik,aik+1\a_i_k,a_i_k+1\ are in one cluster, and the partition in which aik+1,aik+2\a_i_k+1,a_i_k+2\ form a cluster, are valid. In such a case, the disputed alternative aik+1a_i_k+1 should rather tie with the one of the two alternatives aika_i_k and aik+2a_i_k+2 for which the probability of ordinal concordance is lower. That is, if P(aik<aik+1)<P(aik+1<aik+2)P(a_i_k<a_i_k+1)<P(a_i_k+1<a_i_k+2), then the suggested partition should look as follows Q=…,aik,aik+1,aik+2,…,…Q=\…,\a_i_k,a_i_k+1\,\a_i_k+2,…\,…\. After calculating the partition weights for set A, the weights of the alternatives can be recalculated as follows: wki=∑kj∈Qdwkj|Qd|,for eachwki∈Qd.w_k_i= _k_j∈ Q_dw_k_j |Q_d |,\,\,for each\,\,w_k_i∈ Q_d. (23) As a result, all mutually tied alternatives will be assigned the same average weight. The weights of the alternatives, which will be the only elements in the cluster, will remain unchanged. The following algorithm can be used to calculate the proposed partition of the set of alternatives. 1. Sort the alternatives according to the calculated weights and determine the order ai1⪯ai2⪯…⪯aina_i_1 a_i_2 … a_i_n. 2. Determine probability of order concordance P(aip≺aiq)P(a_i_p a_i_q) for all pairs of alternatives for which wip<wiqw_i_p<w_i_q. 3. Create a partition Q=Qi1,…,QinQ=\Q_i_1,…,Q_i_n\ composed of n singletons such that aig∈Qiga_i_g∈ Q_i_g. 4. Create a priority queue L of pairs of ranking-adjacent alternatives (aik,aik+1)(a_i_k,a_i_k+1) sorted in ascending order by the probability P(aik≺aik+1)P(a_i_k a_i_k+1). 5. Extract the first pair p=(aik,aik+1)p=(a_i_k,a_i_k+1) from L, where aik∈Qika_i_k∈ Q_i_k and aik+1∈Qik+1a_i_k+1∈ Q_i_k+1. 6. If, for every pair (ag,ah)(a_g,a_h) where ag∈Qika_g∈ Q_i_k and ah∈Qik+1a_h∈ Q_i_k+1, we have P(aig≺aih)<δP(a_i_g a_i_h)<δ, then join sets QikQ_i_k and Qik+1Q_i_k+1. 7. Repeat steps 5–6 until the L queue is empty 8. For the calculated Q, update the weights of the alternatives according to formula (23) 9. If necessary, the modified weight vector is renormalized. Tie clustering The above algorithm, for an arbitrarily defined tie threshold, produces a partition in which alternatives with relatively low probabilities of matching the specified ranking form tie subsets. Example 12. Let us consider the model discussed in the previous examples (5), (6) and (7). Based on the formulas (17), (18), (19) and (20) we can construct a matrix =[P(ai≺aj)]P=[P(a_i a_j)] specifying the probability of the order in each pair. For the data in the example, we have =(00.99460.99840.99930.99970.99970.99990.005400.68770.86460.94150.94120.9830.00160.312300.77051.0.91251.0.00070.13540.229500.81080.74160.96170.00030.058500.189200.41821.0.00030.05880.08750.25840.581800.89040.00010.01700.038300.10960).P= ( array[]c0&0.9946&0.9984&0.9993&0.9997&0.9997&0.9999\\ 0.0054&0&0.6877&0.8646&0.9415&0.9412&0.983\\ 0.0016&0.3123&0&0.7705&1.&0.9125&1.\\ 0.0007&0.1354&0.2295&0&0.8108&0.7416&0.9617\\ 0.0003&0.0585&0&0.1892&0&0.4182&1.\\ 0.0003&0.0588&0.0875&0.2584&0.5818&0&0.8904\\ 0.0001&0.017&0&0.0383&0&0.1096&0 array ). Based on this set of comparisons, we can determine PCP(C)=0.9946, 0.99840.9984, 0.99930.9993, 0.99970.9997, 0.99970.9997, 0.99990.9999, 0.68770.6877, 0.86460.8646, 0.94150.9415, 0.94120.9412, 0.9830.983, 0.77050.7705, 11, 0.91250.9125, 11, 0.81080.8108, 0.74160.7416, 0.96170.9617, 11, 0.58180.5818, 0.89040.8904\. Based on this, we calculate the following indicators: lcPOIU(C)=0.74161 lcPOI_U(C)=0.74161, lcPOIK(C)=0.5817 lcPOI_K(C)=0.5817, alPOIU(C)=0.9235 alPOI_U(C)=0.9235, and alPOIK(C)=0.8781 alPOI_K(C)=0.8781. One immediately notices the relatively low value of the index lcPOIK(C)=0.5817 lcPOI_K(C)=0.5817, suggesting that in at least one case, the preference relation derived from the calculated weight vector is not very strong. The probability of error (i.e., making the rank reversal) is high and amounts to 1−0.5817=0.41831-0.5817=0.4183. A closer look at the matrix P reveals that the preference values for a6a_6 and a5a_5 are problematic. The average value for the reference and non-reference alternatives is remarkable higher (better): alPOIK(C)=0.8781 alPOI_K(C)=0.8781, suggesting that the actual order of a randomly selected pair of such alternatives matches the calculated weights in 87.887.8 out of 100100 cases. For the reference alternatives, both indices are clearly better: lcPOIU(C)=0.74161 lcPOI_U(C)=0.74161 and alPOIU(C)=0.9235 alPOI_U(C)=0.9235. The risk of rank reversal between non-reference alternatives has decreased to 1−0.7416=0.25841-0.7416=0.2584, while the average probability that the actual ranking will match the result derived from the weight values is alPOIU(C)=0.9235 alPOI_U(C)=0.9235. Calculating the matrix P also allows us to propose ties between the most contested neighboring alternatives. According to the algorithm defined above, for the data in the example, the sorted set of pairs of alternatives looks as follows: L=(a6,a5)withP(a6≺a5)=0.5818(a2,a3)withP(a2≺a3)=0.6877(a4,a6)withP(a4≺a6)=0.7416(a3,a4)withP(a3≺a4)=0.7705(a6,a7)withP(a6≺a7)=0.8904(a1,a2)withP(a1≺a2)=0.9946L= cases(a_6,a_5)&with\,\,P(a_6 a_5)=0.5818\\ (a_2,a_3)&with\,\,P(a_2 a_3)=0.6877\\ (a_4,a_6)&with\,\,P(a_4 a_6)=0.7416\\ (a_3,a_4)&with\,\,P(a_3 a_4)=0.7705\\ (a_6,a_7)&with\,\,P(a_6 a_7)=0.8904\\ (a_1,a_2)&with\,\,P(a_1 a_2)=0.9946 cases for Q initially equal to Q(0)=a1,a2,a3,a4,a5,a6,a7.Q^(0)= \ \a_1 \, \a_2 \, \a_3 \, \a_4 \, \a_5 \, \a_6 \, \a_7 \ \. Assuming δ=0.75δ=0.75, the result of the first iteration (steps 5 and 6 of the algorithm) is the union of the sets containing a6a_6 and a5a_5, i.e. Q(1)=a1,a2,a3,a4,a5,a6,a7,Q^(1)= \ \a_1 \, \a_2 \, \a_3 \, \a_4 \, \a_5,a_6 \, \a_7 \ \, In the next iteration, merge the sets for alternatives a2a_2 and a3a_3, i.e., Q(2)=a1,a2,a3,a4,a5,a6,a7.Q^(2)= \ \a_1 \, \a_2,a_3 \, \a_4 \, \a_5,a_6 \, \a_7 \ \. In the third iteration, although P(a4≺a6)=0.7416<δP(a_4 a_6)=0.7416<δ but P(a4≺a5)=0.8108>δP(a_4 a_5)=0.8108>δ so the sets a4 \a_4 \ and a5,a6 \a_5,a_6 \ are not joined. The remaining iterations also do not change the structure of the set Q. At the end of the procedure, a set of weights is determined such that the weights of alternatives a1a_1, a4a_4, and a7a_7 remain unchanged, i.e., w1=0.0256w_1=0.0256, w4=0.142w_4=0.142, and w7=0.263w_7=0.263 while the weights of a2,a3,a5a_2,a_3,a_5 and a6a_6 are determined as follows: w2=w3=0.092+0.1122=0.102,w_2=w_3= 0.092+0.1122=0.102, w5=w6=0.187+0.1762=0.1815.w_5=w_6= 0.187+0.1762=0.1815. Finally, after normalization, the weight vector in the example under consideration takes the form: w=(0.0257,0.1022,0.1022,0.1423,0.1819,0.1819,0.2636)T.w= (0.0257,0.1022,0.1022,0.1423,0.1819,0.1819,0.2636 )^T. The two pairs of alternatives whose ranking positions were the least certain were ranked exactly in second (alternatives a5a_5 and a6a_6) and fourth (alternatives a2a_2 and a3a_3) place. 6 Towards Group Decision Making The model presented in Section (3) is based on a single pairwise comparison matrix. That is, as in Example (5) the observations under consideration—i.e., the pairwise comparisons taken into account when calculating the weight vector (i.e., the set of observations O<(C)O_<(C)) — come from a single pairwise comparison matrix C. However, this does not have to be the case. In the practice of group decision-making using the pairwise comparison method, many experts may participate in the decision-making process, and each of them is required to provide their own set of pairwise comparisons in the form of a matrix. Furthermore, these matrices may differ in terms of the structure of missing comparisons; that is, what one expert was able to estimate for another may pose a serious problem. By treating pairwise comparisons as distinct observations, we can freely increase their number. Furthermore, it is not necessary to require that the structure of missing comparisons in matrices from different experts be identical. In this situation, the interpretation of the LUL_U matrix changes. It is no longer the Laplacian of the graph of C restricted to the rows and columns corresponding to non-reference alternatives, but rather the Laplacian of the multigraph formed by summing all the matrices after restricting them to the rows and columns corresponding to non-reference alternatives. Suppose we have m experts who provided m PC matrices C(1),…,C(m)C^(1),…,C^(m). For each of them, we construct a projection matrix XU(s)X^(s)_U. Then LU=∑s=1m(XU(s))TXU(s).L_U=Σ^m_s=1 (X^(s)_U )^TX^(s)_U. As a result, we obtain a matrix in which (LU)ii (L_U )_i is the total number (across all C(t)C^(t)) of observations in which alternative i occurs, and (LU)ij (L_U )_ij is the number of C(t)C^(t) matrices in which the pair (i,j)(i,j) exists and has been included in the model as an observation. As before, LUL_U is positive-definite if there exists a sequence of comparisons in the set of comparisons from all experts that connects each non-reference alternative with at least one reference alternative. That is, let G(1),G(2),…,G(m)G^(1),G^(2),…,G^(m) be the pairwise comparison graphs corresponding to the matrices C(1),…,C(m)C^(1),…,C^(m). Then, let G⋆=⋃s=1mG(s).G = ^m_s=1G^(s). Assuming G⋆=(V,E⋆)G =(V,E ) the matrix LUL_U is invertible if, for every vertex vi∈Vv_i∈ V where i∈IUi∈ I_U, there exists at least one vertex uj∈Vu_j∈ V such that j∈IKj∈ I_K and there is a path between them. In theory, therefore, it may happen that no expert matrix is complete enough to calculate the weight vector based on it, but after aggregating the pairwise comparisons into a single list of observations, it becomes possible to calculate the weight vector. Conversely, if at least one of the experts provides a sufficiently complete pairwise comparison matrix to calculate the weight vector, it will be possible to calculate the weight vector even after aggregation. In practice, it is convenient to require that each expert provides a matrix that is sufficiently complete to induce a nonzero number of degrees of freedom. This requirement allows for a subsequent attempt to determine the individual variance for a given expert and to distinguish the impact of individual experts on the final aggregated weight vector. In the context of the statistical analysis in Section (4) increasing the number of experts, and thus the number of matrices, leads to an increase in the number of degrees of freedom in the equation (14), which, in turn, may result in greater confidence in the result, estimated as a probability (Section 4.3) and narrower confidence intervals for the weights of the alternatives (Section 4.2). In practice, therefore, compared to the original single-expert model, the multi-expert model would require expanding the corresponding matrices and vectors XUX_U, XKX_K, and y in the expressions (8, 9, 10, 11) so that they account for the observations from successive experts arranged in sequence (while preserving the order), as well as the sum of all observations in the definitions of values such as SSResS_SRes or σ^2 σ^2 (13, 14). However, in the model under consideration for a single expert, we assumed that errors of their observations are characterized by a certain amount of normal random noise, as reflected by the variance (14) estimated using σ^2 σ^2 . In the opinion of many experts, the assumption that all observations from different matrices will have the same variance is quite strong. It can be met when the experts have similar subject-matter expertise, conduct their assessments under similar conditions, and so on. If these conditions are not met, we suggest to consider that the observation errors of each expert have their own variance, i.e., ε(t)∼N(0,σt2) ^(t) N(0,σ^2_t), for t=1,…,mt=1,…,m. This assumption leads to a weighted regression model [31, p. 190] in which the equation (10) takes the form of θ^U=(XU⊤WXU)−1(XU⊤Wy−XU⊤WXRθR), θ_U= (X _UWX_U )^-1 (X _UWy-X _UWX_R _R ), where W=diag(1σ12In1,…,1σm2Inm).W= diag ( 1σ^2_1I_n_1,…, 1σ^2_mI_n_m ). In this approach each expert’s final ratings are weighted by 1/σt21/σ^2_t for t=1,…,mt=1,…,m. In this way, an expert who is more confident in their judgments (i.e., whose opinions exhibit less variance) will have a greater influence on the aggregated result. 7 Summary In this paper, we present a statistical framework for the pairwise comparisons with reference values method, based on the geometric mean and the logarithmic least-squares method. The proposed approach treats incomplete geometric HRE as a linear regression problem, where expert comparisons are regarded as observations and the logarithms of the alternative weights as estimated parameters. As a result, it enables not only the determination of the weight vector for non-reference alternatives but also a quantitative assessment of the uncertainty associated with the obtained results. The main contribution of this work is the introduction of statistical tools for evaluating the quality of the resulting ranking, including an estimator of the error variance, a covariance matrix of the estimators, confidence intervals for alternative weights, and probability estimates for preserving order relations among alternatives. Based on these concepts, we define the quality indicators lcPOI and alPOI, which measure, respectively, the least certain and the average reliability of ordering relations in the ranking. Unlike traditional inconsistency indices, these indicators simultaneously account for both the inconsistency of the comparison data and the preference distances between alternatives. We also propose a mechanism for identifying alternatives whose ranking positions are uncertain and for grouping them into tie clusters. This approach allows replacing an overly precise yet statistically weakly supported ranking with a more conservative, reliable representation of the results. Furthermore, we show that the proposed framework can be naturally extended to group decision-making settings by treating comparisons provided by multiple experts as a common set of observations, potentially characterized by different error variances. In result, the work provides a coherent statistical foundation for the geometric HRE method and enhances the interpretability of weights and rankings derived from incomplete pairwise comparison data. Acknowledgments The research has been supported by the National Science Centre, Poland within the grant VIRGO 2024/55/B/HS4/00860. References Abastante et al. [2019] Abastante, F., Corrente, S., Greco, S., Ishizaka, A., Lami, I.M., 2019. A new parsimonious ahp methodology: Assigning priorities to many objects by comparing pairwise few reference objects. Expert Systems with Applications 127, 109–120. doi:10.1016/j.eswa.2019.02.036. Bana e Costa et al. [2005] Bana e Costa, C., De Corte, J.M., Vansnick, J., 2005. On the mathematical foundation of MACBETH, in: Figueira, J., Greco, S., Ehrgott, M. (Eds.), Multiple Criteria Decision Analysis: State of the Art Surveys. Springer Verlag, Boston, Dordrecht, London, p. 409–443. Basak [1991] Basak, I., 1991. Inference in pairwise comparison experiments based on ratio scales. Journal of Mathematical Psychology 35, 80–91. doi:https://doi.org/10.1016/0022-2496(91)90035-R. Bortot et al. [2023] Bortot, S., Brunelli, M., Fedrizzi, M., Pereira, A.R.M., 2023. A novel perspective on the inconsistency indices of reciprocal relations and pairwise comparison matrices. Fuzzy Sets and Systems 454, 74–99. doi:https://doi.org/10.1016/j.fss.2022.04.020. decision sciences. Brandt et al. [2016] Brandt, F., Conitzer, V., Endriss, U., Lang, J., Procaccia, A.D. (Eds.), 2016. Handbook of Computational Social Choice. Cambridge University Press, 32 Avenue of the Americas, New York, NY 10013-2473, USA. Brunelli and Corrente [2025] Brunelli, M., Corrente, S., 2025. Do inconsistency indices measure inconsistency of preferences? Journal of Multi-Criteria Decision Analysis 32, e70026. doi:https://doi.org/10.1002/mcda.70026. e70026 MCDA-25-0066. Brunelli and Rezaei [2019] Brunelli, M., Rezaei, J., 2019. A multiplicative best–worst method for multi-criteria decision making. Operations Research Letters 47, 12–15. URL: https://doi.org/10.1016/j.orl.2018.11.008, doi:10.1016/j.orl.2018.11.008. Carriere and Finster [1992] Carriere, J., Finster, M., 1992. Statistical theory for the ratio model of paired comparisons. Journal of Mathematical Psychology 36, 450–460. doi:10.1016/0022-2496(92)90031-2. Cavallo [2020] Cavallo, B., 2020. Functional relations and Spearman correlation between consistency indices. Journal of the Operational Research Society 71, 301–311. doi:10.1080/01605682.2018.1516178. Choo and Wedley [2004] Choo, E.U., Wedley, W.C., 2004. A common framework for deriving preference values from pairwise comparison matrices. Computers and Operations Research 31, 893 – 908. doi:10.1016/S0305-0548(03)00042-X. Colomer [2011] Colomer, J.M., 2011. Ramon Llull: from ‘Ars electionis’ to social choice theory. Social Choice and Welfare 40, 317–328. Condorcet [1785] Condorcet, M., 1785. Essay on the Application of Analysis to the Probability of Majority Decisions. Paris: Imprimerie Royale. Crawford and Williams [1985] Crawford, G., Williams, C., 1985. The Analysis of Subjective Judgment Matrices. Technical Report R-2572-1-AF. The Rand Corporation. Crawford [1987] Crawford, G.B., 1987. The geometric mean procedure for estimating the scale of a judgement matrix. Mathematical Modelling 9, 327 – 334. doi:http://dx.doi.org/10.1016/0270-0255(87)90489-1. De Jong [1984] De Jong, P., 1984. A statistical approach to Saaty’s scaling method for priorities. Journal of Mathematical Psychology 28, 467–478. doi:10.1016/0022-2496(84)90013-0. Genest and Rivest [1994] Genest, C., Rivest, L.P., 1994. A Statistical Look at Saaty’s Method of Estimating Pairwise Preferences Expressed on a Ratio Scale. Journal of Mathematical Psychology 38, 477–496. doi:https://doi.org/10.1006/jmps.1994.1034. Greco et al. [2016] Greco, S., Ehrgott, M., Figueira, J.R. (Eds.), 2016. Multiple Criteria Decision Analysis: State of the Art Surveys. International Series in Operations Research & Management Science, Springer, New York, NY. Harker [1987] Harker, P.T., 1987. Alternative modes of questioning in the analytic hierarchy process. Mathematical Modelling 9, 353 – 360. doi:https://doi.org/10.1016/0270-0255(87)90492-1. Janicki and Koczkodaj [1996] Janicki, R., Koczkodaj, W.W., 1996. A weak order approach to group ranking. Comput. Math. Appl. 32, 51–59. doi:10.1016/0898-1221(96)00102-2. Koczkodaj et al. [2013] Koczkodaj, W.W., Herman, M.W., Orlowski, M., 2013. Managing Null Entries in Pairwise Comparisons. Knowledge and Information Systems 1, 119–125. Koczkodaj and Urban [2018] Koczkodaj, W.W., Urban, R., 2018. Axiomatization of inconsistency indicators for pairwise comparisons. International Journal of Approximate Reasoning 94, 18–29. Kułakowski [2014] Kułakowski, K., 2014. Heuristic Rating Estimation Approach to The Pairwise Comparisons Method. Fundamenta Informaticae 133, 367–386. URL: https://doi.org/10.3233/FI-2014-1081, doi:10.3233/FI-2014-1081. Kułakowski [2016] Kułakowski, K., 2016. Notes on the existence of a solution in the pairwise comparisons method using the heuristic rating estimation approach. Annals of Mathematics and Artificial Intelligence 77, 105–121. URL: http://dx.doi.org/10.1007/s10472-015-9474-6, doi:10.1007/s10472-015-9474-6. Kułakowski [2020a] Kułakowski, K., 2020a. On the geometric mean method for incomplete pairwise comparisons. Mathematics 8. Kułakowski [2020b] Kułakowski, K., 2020b. Understanding the Analytic Hierarchy Process. Chapman and Hall / CRC Press, 6000 Broken Sound Parkway, Boca Raton, FL, 33487, USA. doi:10.1201/b21817. Kułakowski et al. [2015] Kułakowski, K., Grobler-Dębska, K., Wąs, J., 2015. Heuristic rating estimation: geometric approach. Journal of Global Optimization 62, 529–543. doi:10.1007/s10898-014-0253-4. Kułakowski et al. [2026] Kułakowski, K., Kędzior, A., Szybowski, J., Mazurek, J., 2026. Mean-based incomplete pairwise comparisons method with the reference values. arxiv.org. URL: https://arxiv.org/abs/2207.10783, arXiv:2207.10783. Lin et al. [2014] Lin, C., Kou, G., Ergu, D., 2014. A statistical approach to measure the consistency level of the pairwise comparison matrix. Journal of the Operational Research Society 65, 1380–1386. URL: http://dx.doi.org/10.1057/jors.2013.92, doi:10.1057/jors.2013.92. Liu et al. [2021] Liu, F., Qiu, M.Y., G., Z.W., 2021. An uncertainty-induced axiomatic foundation of the analytic hierarchy process and its implication. Expert Systems with Applications 183, 115427. doi:https://doi.org/10.1016/j.eswa.2021.115427. Luo et al. [2024] Luo, D., Zhang, C., Su, W., Zeng, S., Balezentis, T., 2024. Statistical tests for multiplicative consistency of fuzzy preference relations: A monte carlo simulation. Information Sciences 664, 120333. doi:https://doi.org/10.1016/j.ins.2024.120333. Montgomery et al. [2012] Montgomery, D.C., Peck, E.A., Vining, G.G., 2012. Introduction to Linear Regression Analysis. Wiley Series in Probability and Statistics. 5 ed., Wiley-Blackwell, Hoboken, NJ. Opricovic and Tzeng [2004] Opricovic, S., Tzeng, G., 2004. Compromise solution by mcdm methods: A comparative analysis of vikor and topsis. European Journal of Operational Research 156, 445–455. doi:10.1016/s0377-2217(03)00020-1. Panda and Jagadev [2018] Panda, M., Jagadev, A.K., 2018. Topsis in multi-criteria decision making: A survey, in: 2018 2nd International Conference on Data Science and Business Analytics (ICDSBA), IEEE. p. 51–54. URL: http://dx.doi.org/10.1109/ICDSBA.2018.00017, doi:10.1109/icdsba.2018.00017. Peykani et al. [2026] Peykani, P., Emrouznejad, A., Nouri, M., 2026. Best-worst multi-criteria decision-making method: A review of the literature. Socio-Economic Planning Sciences 104, 102345. doi:https://doi.org/10.1016/j.seps.2025.102345. Ramsay [1977] Ramsay, J.O., 1977. Maximum likelihood estimation in multidimensional scaling. Psychometrika 42, 241–266. doi:10.1007/BF02294052. Ramsay [1980] Ramsay, J.O., 1980. Some Small Sample Results for Maximum Likelihood Estimation in Multidimensional Scaling. Psychometrika 45, 139–144. doi:10.1007/BF02293604. Rezaei [2026] Rezaei, J., 2026. Best-worst method: A decade of evolution and future prospects. Omega 143, 103546. doi:https://doi.org/10.1016/j.omega.2026.103546. Roszkowska [2024] Roszkowska, E., 2024. A comprehensive exploration of hellwig’s taxonomic measure of development and its modifications—a systematic review of algorithms and applications. Applied Sciences 14, 10029. doi:10.3390/app142110029. Roszkowska et al. [2020] Roszkowska, E., Filipowicz-Chomko, M., Wachowicz, T., 2020. Using individual and common reference points to measure the performance of alternatives in multiple criteria evaluation. Operations Research and Decisions 30. doi:10.37190/ord200305. Saaty [1977] Saaty, T.L., 1977. A scaling method for priorities in hierarchical structures. Journal of Mathematical Psychology 15, 234 – 281. doi:10.1016/0022-2496(77)90033-5. Saaty and Vargas [2013] Saaty, T.L., Vargas, L.G., 2013. Decision Making with the Analytic Network Process. volume 195 of Economic, Political, Social and Technological Applications with Benefits, Opportunities, Costs and Risks. Springer Science and Business Media, Boston, MA. Sałabun [2014] Sałabun, W., 2014. The characteristic objects method a new distancebased approach to multicriteria decisionmaking problems. Journal of Multi-Criteria Decision Analysis 22, 37–50. doi:10.1002/mcda.1525. Shiraishi and Obata [2025] Shiraishi, S., Obata, T., 2025. Calculating maximum eigenvalues in pairwise comparison matrices for the analytic hierarchy process. Operations Research Forum 6, 10. doi:10.1007/s43069-024-00412-x. Srdjevic and Srdjevic [2023] Srdjevic, B., Srdjevic, Z., 2023. Prioritisation in the analytic hierarchy process for real and generated comparison matrices. Expert Systems with Applications 225, 120015. doi:https://doi.org/10.1016/j.eswa.2023.120015. Thurstone [1927] Thurstone, L.L., 1927. The Method of Paired Comparisons for Social Values. Journal of Abnormal and Social Psychology , 384–400. Wang [2025] Wang, Z.J., 2025. Closed-form solution-based fuzzy utility vectors acquired from trapezoidal fuzzy pairwise comparison matrices using logarithmic quadratic programming for improving fuzzy ahp decision-making systems. Journal of Computational and Applied Mathematics 468, 116647. doi:https://doi.org/10.1016/j.cam.2025.116647.