## 3 Implementation and Software
### Implementation of Distances and MDS Methods  
For the implementation in R of the distance measures ED, cED, MD1, cMD1, MD2, cMD2, MMD, and UMD, we used the code presented in [15]. For BMD1 and cBMD1, we properly modified the code of the MD1 and cMD1 distances. In particular, we used the _weightedCorr()_ function of the _wCorr_ package to calculate tetrachoric correlations and the _qnorm()_ function for the probit function. Note that the probit function cannot be calculated when \(p=0\) or \(p=1\). In these cases, we approximated \(p=0\) by \(p=1/2n\) and \(p=1\) by \(p=1\)-\(1/2n\), where \(n\) is the number of observations used in the calculation of \(p\). The MDS methods mMDS and nMDS were implemented using the _mds()_ function of the _smacof_ package. This function performs metric MDS when using the argument _type_=_\"ratio\"_ and nonmetric MDS when using _type_=_\"_ordinal\"_. PCoA [22] was implemented using the _cmdscale()_ function of R. Note that nMDS is not compatible with the BOOTres method since this combination cannot generate ellipses in most of the datasets applied. Note also that the _cmdscale()_ function may exhibit computational problems when the input dissimilarities are not Euclidean. These problems can be handled via the _add_ argument of this function, which can be used to compute a constant added to the nondiagonal dissimilarities such that the modified dissimilarities are Euclidean. Finally, in what concerns the _corMDS_, the MDS solution was determined by minimizing the stress function:  
\[\text{Stress}=1-\text{corr}(\mathbf{d},\mathbf{D}) \tag{6}\]  
where corr is the correlation coefficient between \(\mathbf{d},\mathbf{D},\mathbf{d}=\left(d_{11},\ d_{12},\ldots,d_{ij},\ \ldots,\ d_{(n-1)n}\right)\) is the vector with initial distances between all samples \(i\) and \(j\) (\(i<j\)), and \(\mathbf{D}=\left(D_{11},\ D_{12},\ldots,D_{ij},\ \ldots,\ D_{(n-1)n}\right)\) is the vector with the MDS Euclidean distances calculated by:  
\[\text{D}_{\text{ij}}=\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}} \tag{7}\]  
Here, \((x_{i},\ y_{i})\) and \((x_{j},\ y_{j})\) are the coordinates of the samples \(i\) and \(j\) in the MDS map. It is seen that based on this approach, a MDS solution arises from the determination of \(2n\) adjustable parameters \((x_{i},\ y_{i})\), where \(i\) ranges from \(1\) to \(n\), which minimize the stress function. The MDS solution is centered. For its determination, optimization techniques may be used. In this study, we used the _optim()_ function of the base library of R with the \"BFGS\" option, which applies a quasi-Newton method. Note that not all \(2n\) adjustable parameters are needed to minimize the stress function, but \(2n\)-\(3\) of them, since the minimum stress value is invariant if we select arbitrarily and keep constant three of the \(2n\) values of \(x_{i},y_{i}\), for example, the values of \(x_{n},y_{n}\), and \(y_{n-1}\). This is because the values of \(x_{n}\) and \(y_{n}\) determine the position of the MDS solution in the two-dimensional Euclidean space, whereas the value \(y_{n-1}\) determines its orientation at the point \((x_{n},y_{n})\). Therefore, for \(x_{n},y_{n}\), and \(y_{n-1}\), we may use a constant value, for example, \(x_{n}=y_{n}=y_{n-1}=100\), or we may use the values obtained from another MDS technique, say the mMDS or nMDS method. Both approaches give convergent results. Therefore, the corMDS can be directly applied to the cED, cMD1, cMD2, MMD, UMD, and cBMD1 distance measures, where the correction for small sample sizes may lead to negative distance values. In contrast, this is not possible for the conventional mMDS, nMDS, and PCoA methods, which demand a nonnegative input distance matrix. To overcome this problem, we have two options. Negative distance values may be transformed into zero [25], or a proper fixed value may be added to all pairwise distances [26]. In the present study, we have adopted the arguments presented by [PERSON] et al.