Classically, clustering has been divided into hierarchical and partitional methods. The main difference between the two is that the hierarchical method breaks the data up into hierarchical clusters, whereas the partitional method divides the data into mutually disjoint partitions. The pillar of the partitional algorithms is the k-Means algorithm [²⁶], introduced almost four decades ago. As a consequence, the k-Means paradigm has been extended to various versions, including the k-Modes algorithm [²⁵] for categorical data.

The k-Modes algorithm owes its existence to the ineffectiveness of the k-Means algorithm for handling categorical data. Ralambondrainy [²⁷] attempted to rectify this using a hybrid numeric–symbolic method based on the binary characters 0 and 1. However, this approach suffered from an unacceptable computational cost, particularly when the categorical attributes had many categories. Since then, a variety of k-Modes-type algorithms have been introduced, such as k-Modes with new dissimilarity measures [²¹,²²], k-Population [²³], and a new Fuzzy k-Modes [²⁰].

Partitional algorithms use an objective function in their optimization process, and the determination of this function was described as the P problem by Bobrowski and Bezdek [²⁸] and Salim and Ismail [²⁹]. When he proposed the k-Modes clustering algorithm, Huang [²⁵] split P into P₁ and P₂. P₁ denotes the minimization problem of obtaining values for the partition matrix w_li of 0 or 1 (for the hard clustering approach) or 0 to 1 (for the fuzzy clustering approach); see Eq. (1b) as an example. Furthermore, P₂ denotes the minimization problem of obtaining the value that occurs most often (or the mode of a categorical data set) to represent the center of a cluster (often called the centroid). The minimization of P₂ by obtaining the appropriate mode essentially causes the minimization of problem P₂, and vice versa. As an example of the optimization process for problem P in the Fuzzy k-Modes algorithm, we wish to solve Eq. (1) subject to Eqs. (1a), (1b), and (1c).

subject to:

And

where:

• w_li is a (k × n) partition matrix that denotes the degree of membership of object i in the lth cluster that contains a value of 0 to 1,

• k (≤ n) is a known number of clusters,

• Z is the centroid such that [Z₁, Z₂,…,Z_k] ∈ R^mk,

• α [1, ∞) is a weighting exponent,

• d(X_i, Z_l) is the distance measure between the object X_i and the centroid Z_l, as described in Eqs. (2) and (2a).

where:

Huang and Ng [²⁴] described the optimization process of P₁ and P₂ as follows:

• Problem P₁: Fix Z = Z^ and solve the reduced problem P(W,Z^) as in Eq. (3). This process obtains the minimized values of 0–1 of the partition matrix w_li.

• Problem P₂: Fix W = Ŵ and solve the reduced problem P(Ŵ, Z) as in Eq. (4) subject to Eq. (4a). This process obtains the most frequent attributes, or the modes, which give the centroids.

where:

and α ∈ [1, ∞) is a weighting exponent.

Due to the characteristics of Y-STR data, there are two optimization problems for existing partitional algorithms: non-unique centroids and local minima problems. These two problems are caused by the drawback of the modes mechanism of determining the centroids. Non-unique centroids would result in empty clusters, whereas the local minima problem leads to poorer clustering results. Both problems are a result of the obtained centroids, which are not sufficient to represent their classes.

Therefore, problems will occur for the following two cases:

i)The total number of objects in a dataset is small while the number of classes is large. To illustrate this case, consider the following example.

Example I: Figure 1 shows an artificial example of a dataset consisting of nine objects in three classes: Class A = {A₁, A₂, A₃}, Class B = {B₁, B₂, B₃}, and Class C = {C₁, C₂, C₃}. Each object is composed of three attributes, represented in lower case; e.g., for object A₁, the attributes are a₁, a₂, and a₃. The dataset is considered to have a higher degree of similarity between objects in intra-classes, while the number of objects is small and number of classes is large. Thus, the appropriate modes for representing the classes are: Class A – [a₁, a₂, a₃], lass B – [a₁, b₂, c₃], and Class C – [b₁, c₂, d₄]. However, attribute a₁ in DOMAIN (A₁), a₂ in DOMAIN (A₂), and c₃ in DOMAIN (A₃) are too dominant, and would therefore dominate the process of updating P₂. Figure 2 shows the possibility that each cluster is formed by the dominant attributes.

As a result, the mode that consists of [a₁, a₂, c₃] would be obtained twice. Thus, P₂ would not be minimized due to this non-unique centroid. Another possibility is that the two modes are different, but are not distinctive enough to represent their clusters, such as modes [a₁, a₂, a₃] or [a₁, a₂, b₃] for Cluster 2. As a consequence, this case would fall into a local minima problem.

ii)An extreme distribution of objects in a class. To illustrate this case, consider the following example.

Example II: Figure 3 shows a dataset consisting of eight objects in two classes: Class A = {A₁, A₂, A₃, A₄, A₅, A₆} and Class B = {B₁, B₂}. Each object consists of three attributes, again represented in lower case. The appropriate modes to represent the classes are: Class A – [a₁, a₂, b₃] and Class B – [a₁, b₂, c₃] or [a₁, b₂, d₃]. The distribution of objects in Class A is considerably larger than in Class B, covering approximately 75% of the total set of objects. This characteristic of the data is found to be problematic for P₂, particularly for the fuzzy approach. The problem is actually caused by the initial centroid selection. Figure 4 shows the objects in Class A would be equally distributed into clusters 1 and 2.

As a result, object A becomes dominant in both clusters, and so the obtained modes might be represented solely by objects in Class A, e.g., [a₁, a₂, a₃] and [a₁, a₂, b₃].

The above situations cause P not to be fully optimized, thus producing poor clustering results. Therefore, a new algorithm with a new concept of P₂ is proposed in order to overcome these problems and improve the clustering accuracy results of Y-STR data.

The mode mechanism for the center of a cluster (problem P₂) is not appropriate for handling the characteristics of Y-STR data, and therefore, it cannot be used as a mechanism to represent the center of a cluster (centroid). Instead, the center of Y-STR data should be the modal haplotypes, which are required to calculate the distance of Y-STR objects. The distance between a Y-STR object and its modal haplotype can be formalized as in Eq. (5) subject to Eq. (5a).

subject to:

where m is the number of markers.

The modal haplotype is controlled by groups of objects that are similar or almost similar in Y-STR data. The similar and almost similar objects have a lower distance, or a higher degree of membership values in a fuzzy sense. Thus, these two groups are considerably the most dominant objects required to find the Approximate Modal Haplotype. Consider four objects x₁, x₂, x₃, and x₄ and two clusters c₁ and c₂. The membership value for each object and its cluster are as shown in Table 1, whereby objects x₁ and x₃ have a 100% chance of being the most dominant object in cluster c₁, but only a 50% chance of being the dominant object in cluster c₂, and so on. A dominant weighting value of 1.0 is given to any dominant object and a weight of 0.5 is given to the remaining objects.

Let X ={X₁, X₂,…, X_n} be a set of n Y-STR objects and A ={A₁,A₂,…, A_m} be a set of markers (attributes) of a Y-STR object. Let H = {H₁, H₂,.,H_k} ∈ X be the set of Approximate Modal Haplotypes for k clusters. Suppose k is known a priori. Let H_l be the Approximate Modal Haplotype, represented as [h_l,1, h_l,2,…,h_l,m], and therefore, H_l,j = X_i,j for 1≤ j ≤ m and 1≤ i ≤ n. The objective of the algorithm is to partition the categorical objects X into k clusters. Thus, the H_l can be replaced by X_i until n provided they satisfy the condition described in Eq. (6).

Here, P(Á) is the cost function described in Eq. (7), which is subject to Eqs. (7a), (8), (8a), (8b), (9), (9a), (9b), and (9c).

subject to:

• W_li^∝ is a (k × n) partition matrix that denotes the degree of membership of Y-STR object i in the lth cluster that contains a value of 0 to 1 as described in Eq. (8), subject to Eqs. (8a) and (8b).

• subject to:

and

where,

• k (≤ n) is a known number of clusters.

• H is the Approximate Modal Haplotype (centroid) such that [H₁, H₂,…,H_k] ∈ X.

• α ∈ [1, ∞) is a weighting exponent and used to increase the precision of the membership degrees. Note that this alpha is typical based on 1.1 until 2.0 as introduced by Huang and Ng [²⁴].

• d_ystr(X_i,H_l) is the distance measure between the Y-STR object X_i and the Approximate Modal Haplotype H_l as described in Eq. (5) and subject to Eq.(5a).

• D_li is another (k × n) partition matrix which contains a dominant weighting value of 1.0 or 0.5, as explained above (See Table 1). The dominant weighting values are based on the value of W_li^∝ above. D_li is described in Eq. (9), subject to Eqs. (9a), (9b), and (9c).

subject to:

The basic idea of the k-AMH algorithm is to find k clusters in n objects by first randomly selecting an object to be the Approximate Modal Haplotype h for each cluster. The next step is to iteratively replace the objects x one-by-one towards the Approximate Modal Haplotype h. The replacement is based on Eq. (6) if the cost function as described in Eq. (7) and subject to (7a), (8), (8a), (8b), (9), (9a), (9b) and, (9c) is maximized. Thus, the differences between the k-AMH algorithm and the other k-Mode-type algorithms are as follows.

i. The objects (the data themselves) are used as the centroids instead of modes. Since the distance of Y-STR objects is measured by comparing the objects and their modal haplotypes, we need to approximately find the objects that can represent the modal haplotypes. In finding the final Approximate Modal Haplotype for a particular group (cluster), each object needs to be tested one-by-one and replaced on a maximization of a cost function.

ii. A maximization process of the cost function is required instead of minimizing it as in the k-mode-type algorithms.

A detailed description of the k-AMH algorithm is given below.

Step 1 – Select k initial objects randomly as Approximate Modal Haplotype (centroids). E.g. if k = 4, then choose randomly 4 objects as the initial Approximate Modal Haplotype.

Step 2 – Calculate distance d_ystr(X_i,H_l) according to Eq. (5) and subject to (5a).

Step 3 – Calculate partition matrix w_li^∝ according to Eq. (8), subject to Eqs. (8a) and (8b). Note that the w_li^∝ is based on the distance calculated in Step 2.

Step 4 – Assign a weighting dominant of 1.0 or 0.5 for partition matrix D_li according to Eqs. (9), (9a), (9b) and (9c).

Step 5 – Calculate cost function P(Á) based on W_li^∝D_li according to Eqs (7) and (7a).

Step 6 – Test for each initial modal haplotype by the other objects one-by-one. If current cost function is greater than previous cost function according to Eq. (6), then replace it.

Step 7 – Repeat Step 2 until Step 6 for each x and h

Step 8 – Once the final Approximate Modal Haplotypes are obtained for all clusters, assign the objects to their corresponding crisp clusters C_li according to Eq. (10).

Furthermore, the implementation of the steps above of the algorithm is formalized in the form of pseudo-code as follows.

INPUT: Dataset, X, the number of cluster, k, the number of dimensional, d and the fuzziness index,

OUTPUT: A set of clusters, k

01: Select H_l randomly from Xsuch that 1≤l≤ k

02: for each H_l an Approximate Modal Haplotype do

03: for each X_i do

04: Calculate P(À)  =   ∑ _l = 1^k ∑ _i = 1ⁿÀ_li

05: if P(À)  =   ∑ _l = 1^k ∑ _i = 1ⁿÀ_li is maximized, then

06: Replace H_l by X_i

07: end if end for

09: end for

10: Assign X_i to C_l for all l, 1≤ l ≤ k; 1≤i≤ n as Eq. (10)

11: Output Results

In optimizing the problem P, the k-AMH algorithm uses a maximization process instead of the minimization process imposed by the k-Mode-type algorithms. This process is formalized in the k-AMH algorithm as follows.

Step 1 - Choose an Approximate Modal Haplotype, H^(t)∈ X. Calculate P(Á); Set t=1

Step 2 - Choose X^(t+1) such that P(Á)^t+1 is maximized; Replace H¹ by X^(t+1)

Step 3 - Set t=t+1; Stop when t=n; otherwise go to Step 2.

*Note:n is the number of objects

The convergence of the algorithm is proven as P₁ and P₂ are maximized accordingly. The function P(Á) incorporates the P(W, H) function imposed by the Fuzzy k-Modes algorithm, where W is a partition matrix and H is the approximate modal haplotype that defines the center of a cluster. Thus, P₁ and P₂ are solved by Theorems 1 and 2, respectively.

Theorem 1 – Let Ĥ be fixed. P(W, Ĥ) is maximized if and only if

Proof

Let X= {X₁,X₂,.,X_n} be a set of n Y-STR categorical objects and H= {H₁,H₂,.,H_k} be a set of centroids (Approximate Modal Haplotypes) for k clusters. Suppose that P= {P₁,P₂,.,P_k} is a set of dissimilarity measures based on d_ystr(X_i,H_l), as described in Eqs. (5) and subject to (5a), ∀ i and l 1  ≤  i  ≤  n; 1  ≤  l  ≤  k

Definition 1 - For X_i  =  H_l and X_i  =  H_z, where z  ≠  l, the membership value for all i is

For any P that is obtained from d_ystr(X_i,H_l) where X_i  =  H_l, the maximum value of w_li^∝ is 1 and X_i  =  H_z, z  ≠  l the value of w_li^∝ is 0. Therefore, because H_l is fixed, w_li^∝ is maximized.

Definition 2 – For the case of H_i  ≠  X_iand X_i ≠ H_z,   ∀  z,  1  ≤  z  ≤  k, the membership value for all i is

Suppose that p_li ∈ P is the minimum value, we write as

Therefore,

where z ≠ l

Thus, ∑z=1kPlizi1∝−1<∑z=1kPtiPzi1∝−1 where

t ≠ l and ∀  z and t,  1  ≤  z  ≤  k;  1  ≤  t  ≤  k It follows that

where t ≠ l

Therefore, based on definitions 1 and 2, w_li^∝ is maximal. Because Ĥ is fixed, PW,H^ is maximized.

Theorem 2 – Let h_l ∈ X be the initial center of a cluster for 1 ≤ l ≤ k. h_l is replaced by x_i as the Approximate Modal Haplotype if and only if

Proof

Let D= {D₁,D₂,.,D_k} be a set of dominant weighting values. For any maximum value of w_li^∝ as proved by Theorem 1, we assign an optimum value of 1.0 as a dominant weighting value, otherwise 0.5 as described in Eq, (9) and subject to Eqs. (9a), (9b) and (9c). We write

Because w_li^∝ and D_li are non-negative, the product W_li^∝D_li must be maximal. It follows that the sum of all quantities ∑ _l = 1^k ∑ _i = 1ⁿÁ_li is also maximal. Hence, the result follows.

The Y-STR data were mostly obtained from a database called worldfamilies.net [³⁰]. The first, second, and third datasets represent Y-STR data for haplogroup applications, whereas the fourth, fifth, and sixth datasets represent Y-STR data for Y-surname applications. All datasets were filtered for standardization on 25 similar attributes (25 markers). The chosen markers include DYS393, DYS390, DYS19 (394), DYS391, DYS385a, DYS385b, DYS426, DYS388, DYS439, DYS389I, DYS392, DYS389II, DYS458, DYS459a, DYS459b, DYS455, DYS454, DYS447, DYS437, DYS448, DYS449, DYS464a, DYS464b, DYS464c, and DYS464b. These markers are more than sufficient for determining a genetic connection between two people. According to Fitzpatrick [³¹], 12 markers (Y-DNA12 test) are already sufficient to determine who does or does not have a relationship to the core group of a family.

All datasets were retrieved from the respective websites in April 2010, and can be described as follows:

1) The first dataset consists of 751 objects of the Y-STR haplogroup belonging to the Ireland yDNA project [³²]. The data contain only 5 haplogroups, namely E (24), G (20), L (200), J (32), and R (475). Thus, k = 5.

2) The second dataset consists of 267 objects of the Y-STR haplogroup obtained from the Finland DNA Project [³³]. The data are composed of only 4 haplogroups: L (92), J (6), N (141), and R (28). Thus, k = 4.

3) The third dataset consists of 263 objects obtained from the Y-haplogroup project [³⁴]. The data contain Groups G (37), N (68), and T (158). Thus, k = 3.

4) The fourth dataset consists of 236 objects combining four surnames: Donald [³⁵], Flannery [³⁶], Mumma [³⁷], and William [³⁸]. Thus, k = 4.

5) The fifth dataset consists of 112 objects belonging to the Philips DNA Project [³⁹]. The data consist of eight family groups: Group 2 (30), Group 4 (8), Group 5 (10), Group 8 (18), Group 10 (17), Group 16 (10), Group 17 (12), and Group 29 (7). Thus, k = 8.

6) The sixth dataset consists of 112 objects belonging to the Brown Surname Project [⁴⁰]. The data consist of 14 family groups: Group 2 (9), Group 10 (17), Group 15 (6), Group 18 (6), Group 20 (7), Group 23 (8), Group 26 (8), Group 28 (8), Group 34 (7), Group 44 (6), Group 35 (7), Group 46 (7), Group 49 (10), and Group 91 (6). Thus, k = 14.

The values in parentheses indicate the number of objects belonging to that particular group. Datasets 1–3 represent Y-STR haplogroups and datasets 4–6 represent Y-STR surnames.

Table 2 shows the clustering accuracy scores for all datasets (boldface indicates the highest clustering accuracy). Based on these results, the performance of the k-AMH algorithm was very promising. Out of six datasets, our algorithm obtained the highest clustering accuracy scores for datasets 1, 2, 4, 5, and 6. In fact, the algorithm also achieved the optimal clustering accuracy for two datasets (4 and 5). However, for dataset 3, the results show that the accuracy of the k-AMH algorithm was 0.01 lower than that of the k-Population algorithm. A statistical t-test was carried out for further verification. This indicated that t(101.39) = 0.65, and p = 0.51. Thus, there was no significant difference at the 5% level between the accuracy score of our k-AMH algorithm and the k-Population algorithm. This means that both algorithms displayed an equal performance for this dataset.

During the experiments, the k-AMH algorithm did not encounter any difficulties. However, the Fuzzy k-Modes and the New Fuzzy k-Modes algorithms faced problems with datasets 1, 5, and 6. For dataset 1, the problem was caused by the extreme number of objects in Class R (475), which covered about 63% of the total objects. Further, for datasets 5 and 6, the problem was caused by many similar objects in a larger number of classes. In particular, both algorithms faced the problem P₂ caused by the initial centroid selections. Note also that the results for both algorithms were based on the diverse method, an initial centroid selection proposed by Huang [²⁵].

For an overall comparison, Table 3 shows the results of all Y-STR datasets. It clearly indicates that the k-AMH algorithm obtained the highest accuracy score of 0.93. The closest score of 0.91 belongs to the k-Population algorithm. Furthermore, the k-AMH algorithm also recorded the best results in terms of standard deviation (0.07), the lower bound (0.93), the upper bound (0.94), and the minimum accuracy score (0.79).

For further verification, a one-way ANOVA test was carried out. This indicated that the assumption of homogeneity of variance was violated; therefore, the Welch F-ratio is reported. There was a significant variance in the clustering accuracy scores among the nine algorithms, in which F(8, 2230) = 378, p < 0.001, and ω² = 0.25. Thus, the Games–Howell procedure was used for a multiple comparison among the nine algorithms. Table 4 shows the result of this comparison with regard to the k-AMH algorithm against the other eight algorithms. At the 5% level of significance, it is clearly shown that the k-AMH algorithm (M = 0.93, 95% CI [0.93, 0.94]) differed from the other eight algorithms (all P values < 0.001). Thus, the performance of k-AMH algorithm exhibited a very significant difference compared to the other algorithms.

We now consider the time efficiency of the k-AMH algorithm. The computational cost of the algorithm depends on the nested loop for k(n-k), where k is the number of clusters and n is the number of data required to obtain the cost function, P(À). The function P(À) involves the number of attributes m in calculating the distances and the membership values for its partition matrix w_li. Thus, the overall time complexity is O(km(n-k)). However, the time efficiency of the k-AMH algorithm will not reach O(n²) because the value of k in the outer loop will not become equivalent to the value of n-k in the inner loop. See pseudo-code for a detailed implementation of these loops.

A scalability test was also carried out for the k-AMH algorithm. These experiments were based on a dataset called Connect [⁴²]. This dataset consisted of 65,000 data, 42 attributes, and three classes. Two scalability tests were conducted: (a) scalability against the number of objects, when the number of clusters was three, and (b) scalability against the number of clusters, when the number of objects was 65,000. The test was performed on a personal computer with an Intel® Core™ 2 DUO Processor with 2.93 GHz and 2.00 GB memory. Figure 5(a) and (b) illustrate the results of the tests. In conclusion, the runtime of the k-AMH algorithm increased linearly with the number of clusters and data.

*Note: p < 0.05.