A robust method for clustering football players with mixed attributes

A robust fuzzy clustering model for mixed data is proposed. For each variable, or attribute, a proper dissimilarity measure is computed and the clustering procedure combines the dissimilarity matrices with weights objectively computed during the optimization process. The weights reflect the relevance of each attribute type in the clustering results. A simulation study and an empirical application to football players data are presented that show the effectiveness of the proposed clustering algorithm in finding clusters that would be hidden unless a multi-attributes approach were used.


Introduction and literature review
Data in sports are being collected and analyzed, with the integration of physical and digital sources, increasing the knowledge of professional sports for all parties involved. Statistical methodology and data-driven analytics in sports can drive decision-making in different fields: marketing strategies, performances of players or teams, forecasting of revenues, health. Depending on the type of sport, on the nature of the data at hand, and on the objectives of the analysis, a variety of statistical learning and operations research methods has been proposed. In order to analyse the massive and the different kinds of sport data and their complex structure and then to capture their extensive information many advanced statistical methodologies, strategies of analyses and data-driven procedures have to be considered in the analysis pro-cess. In this way, managing this kind of data with advanced theoretical tools we obtain an informational gain and then the knowledge that represents the basis of any decision-making process in sport.
The playing characteristics of football players are important both from a technical and economic point of view. From the technical point of view they allow to evaluate the playing characteristics that lead the player and the team to achieve winning results; from the economic point of view they allow to establish the value of a player in the transfer market. Clustering of football players on the basis of playing characteristics, position and performance variables is relevant for clubs, either to drive team formation and selection of players, or for determining the value of a football player in the transfer window period (Behravan and Razavi 2021;Shelly et al. 2020;Narizuka and Yamazaki 2019).
In the literature, many empirical studies and methodological proposals based on data science and data-driven approach have been carried out on many sports disciplines to analyze the large mass of sport data both in the field of performance and in the medical, social or economic fields (Table 1). For instance, Palacios-Huerta (2004) analyses the effects of rules of the game of football using an econometric methodology for dating structural breaks in tests with non-standard asymptotic distributions. Dawson et al. (2007) present a statistical analysis of patterns in the incidence of disciplinary sanction (yellow and black cards) that were taken against football players in the English Premier League over the period 1996. Goossens et al. (2012 compare league formats with respect to match importance in Belgian football. Yang et al. (2014) evaluate the efficiency of National Basketball Association (NBA) teams under a two-stage DEA (Data Envelopment Analysis) framework. Applying the additive efficiency approach, they decompose overall team efficiency into first-stage wage efficiency and second-stage on-court efficiency and find out the individual endogenous weights for each stage. Koopman and Lit (2015) propose a dynamic bivariate Poisson model for analysing and forecasting match results in the English Premier League. Nikolaidis (2015) builds a basketball game strategy through statistical analysis of data. In particular, the aim of his paper is on the one hand to present some indicative, simple ideas for the statistical analysis of basketball data, and on the other hand to show that any basketball team can improve significantly its decision-making process if it chooses to be statistically supported. Andrienko et al. (2017) propose a visual analysis of pressure in football. Carpita et al. (2019) explore and model team performances of the Kaggle European Soccer database. Galariotis et al. (2018) propose a twostage method for the concurrent evaluation of the business, financial and sports performance of football clubs analysing the case of France. Geenens and Cuddihy (2018) review the Wald confidence interval for a proportion, suggest new non-parametric confidence intervals for conditional probability functions, revisit the problems of bias and bandwidth selection when building confidence intervals in non-parametric regression and provide a novel bootstrapbased solution to them. The new intervals are used when analysing game outcome data for the UEFA (Union of European Football Associations) Champions and Europa Leagues from 2009-2010to 2014-2015. McHale and Relton (2018 identify key players in soccer teams using network analysis and pass difficulty. Metulini et al. (2018) model the dynamic pattern of surface area in basketball and its effects on team performance. Zuccolotto et al. (2018) use big data analytics for modeling scoring probability in basketball in order to study the effect of shooting under high-pressure conditions. Goes et al. (2018) propose a data-driven model to measure pass effectiveness in professional soccer matches. Van Bulck et al. (2019) consider a tabu search based approach by scheduling a non-professional indoor football league. Adhikari et al. (2020) propose a methodology for cricket player selection based on an efficiency data envelopment analysis, semi-variance approach, and Shannon-entropy. Cea et al. (2020) analyze the procedure used by FIFA up until 2018 to rank national football teams and define by random draw the groups for the initial phase of the World Cup finals. They calibrate a pblackictive model to form a reference ranking to evaluate the performance of a series of simple changes to that procedure. These proposed modifications are guided by a qualitative and statistical analysis of the FIFA ranking. Successively they analyze the use of this ranking to determine the groups for the World Cup finals. Dadeliene et al. (2020) analyse the effects of high intensity training on physical and functional capacities of elite kayakers by using the principal component analysis.
Papers specifically on clustering of sports data are Gates et al. (2017), Behravan and Razavi (2021), Fortuna et al. (2018), Lu and Tan (2003), Narizuka and Yamazaki (2019), Narizuka and Yamazaki (2020), Shelly et al. (2020), Ulas (2021), most of which with applications to football data. Gates et al. (2017) propose a unsupervised classification method that defines subgroups of individuals that have similar dynamic models. They apply this method on functional MRI from a sample of former American football players. In Behravan and Razavi (2021) a two phase method is proposed. In the first phase, the dataset is clustered using an automatic clustering method called APSO-clustering. In the second phase, a hybrid regression method which is a combination of particle swarm optimization (PSO) and support vector regression (SVR), is used to build a prediction model for each clustersâ data points. In Fortuna et al. (2018), focusing on top football players data, a comparison of functional k-means and functional hierarchical clustering for detecting specific patterns of google queries over time is presented. In Lu and Tan (2003) an unsupervised clustering of dominant scenes in sports video is presented, in which data are prepocessed by Principal Components and Linear Discriminant Analysis. Narizuka and Yamazaki (2019) develop a clustering algorithm to extract transition patterns of the formation of particular team during the game. Narizuka and Yamazaki (2020) perform a network bipartite graph and subgroup (cluster) analyses to clarify the injured player's experience and the cause of injury on longitudinal rugby data. Shelly et al. (2020) use K-means Clustering to Create Training Groups for Elite American Football Student-athletes Based on Game Demands. In Ulas (2021) NBA team's characteristics and similarities were assessed firstly with Machine Learning techniques (K-means and Hierarchical clustering) and secondly with Ordinary Linear Regression (OLS) to investigate the factors that affect the NBA team values.
Finally, we remark the interesting special issue on 'Statistical Modelling for Sports Analytic's by Groll et al. (2018).
The presented literature has shown the importance of partitioning and clustering of football players on the basis of performance, position and other variables. The proposed clustering model aims at targeting some relevant issues: (i) the variables of interest in sport are of different types (mixed data), e.g, quantitative, nominal, time series; (ii) these variables don't play the same role in measuring the within cluster similarity; (iii) robustness is a desirable property for a clustering method. The proposed model takes into account the three points. A mixed distance for the different attributes is considered; weights to distances related to different attribute types giving relevance to the variable types capable to increase the within cluster similarity among the units are objectively provided by the model; a noise cluster represented by a noise prototype is introduced to achieve robustness with respect to outliers.
The proposal is novel for the methodology used, a robust PAM Fuzzy clustering algorithm based on a weighted mixed distance, and for its application to positional and performance football players data.
The paper is structured as follows. In Sect. 2 the Robust Fuzzy C-Medoids Clustering for Mixed Data model (FCMd-MD-NC) is proposed. In Sect. 3 a simulation study is carried out to illustrate the performance of the proposed clustering model. Section 4 reports the results of the application of the model to clustering of football players, to show the substantive features of FCMd-MD-NC. Section 5 concludes the paper and provides directions for future work.

Robust fuzzy C-medoids clustering for mixed data model (FCMd-MD-NC model)
Let X = {X 1 , . . . , X P } be a set of P variables, or attributes, observed on n units, in which the P variables are of different types (mixed data), e.g, quantitative, nominal, time series, sequences of qualitative data, imprecisely observed data, textual data. More precisely, the set X contains S types of variables, with p s variables for each attribute type, with s = 1, . . . , S; 1 < S ≤ P; 1 ≤ p s < P; S s=1 p s = P.
Without loss of generality, assume that variables are arranged so that the first p 1 variables are of the same type (for instance, quantitative), the second p 2 variables are also of the same type, different from that of the first p 1 variables (for instance, qualitative), and so on, so that .., X p 1 +...+ p s } is the set of variables of the s-th type. Finally, X is is the set of values observed for the i-th unit on the p s variables of the s-th type.
Depending on the nature of the attribute, X is could be a vector, a matrix, or could have a more complicated structure. For instance, in the case of quantitative variables, X is ≡ x is is the vector of p s values observed on the i-th unit. In the case of time series of length T , X is ≡ X is is a T × p s matrix whose columns are represented by the p s time series observed on the i-th unit, and the rows are the values observed at time t (t = 1, . . . , T ). In the case of ordeblack sequences of qualitative items X is is a set of p s sequences (see D'Urso and Massari 2013).
The distance between units i and i computed according to the nature of the s-th variable type-on this, see Remark 2 below-can be formalized as: (1) Then is the overall weighted squared distance considering the S attribute types. As observed by Everitt (1988), the weights of the squared distance are in a quadratic form. The role of the weights will be discussed at large in Remark 3.
As an example, suppose that X = {X 1 , X 2 } where X 1 is a set of two quantitative variables, while X 2 is a set of two qualitative variables. Then, S = 2, p 1 = p 2 = 2, P = 4, and are the matrices of the pairwise distances-say, Euclidean distance for X 1 and overlapping distance for X 2 , respectively. Then Once the formal notation and the overall distance have been described, in the following the clustering algorithm can be illustrated. Following the PAM approach in a fuzzy framework, let X s ≡ { X 1s , . . . , X cs , . . . , X (C−1)s } be a subset of X s with cardinality C − 1, and X cs ∈ X s the values observed for the c-th elements of X s . Then, X s ≡ { X 1s , . . . , X cs , . . . , X (C−1)s } is a subset of X with cardinality C − 1.
This model achieves its robustness with respect to outliers by introducing a noise cluster, provided there is a way in which all the noise points could be dumped into that single cluster. By following Davé (1991), "Noise prototype is a universal entity such that it is always at the same distance from every point in the data-set." Provided the noise cluster distance is specified, objects closer to the noise cluster than to other objects would get classified into the noise cluster. In this proposal the noise cluster is represented by a noise prototype, i.e. a noise medoid, which is always at the same distance from all units. Let there be C − 1 good clusters and let the C-th cluster be the noise cluster. Let X C be the noise prototype (i.e. noise medoid). It is assumed that the distance measure of unit i from medoid C is equal to δ 2 , i = 1, . . . , n.
Formally, the proposed clustering model, called Fuzzy C-Medoids Clustering of Mixed Data model with Noise Cluster (FCMd-MD-NC model) is characterized in the following way: where: • u ic indicates the membership degree of the i-th objects to the c-th cluster; • m > 1 is a weighting exponent that controls the fuzziness of the obtained partition; • X cs is the s-th component of th c-th medoid, related to the s-th variable type; • C is the noise cluster; • s d ic = d(X is , X cs ) denotes the distance between the i-th observation and the c-th medoid, according to the s-th variable type; for comparison's sake across attribute types, the S distances s d ic are normalized to vary in the range [0, 1]; . . , C − 1 and represents the overall weighted squared distance between unit i and the medoid c based on all variable types; d 2 ic = δ 2 for c = C and represents the distance of each unit from the noise cluster; • w s is the weight associated to the s-th attribute type, and, hence, to the s-th distance (s = 1, . . . , S); The proposed model considers separately the distances for the different attributes and uses a suitable weighting system computed within the model for each distance component. Then the weights w s constitute specific parameters to be estimated within the clustering procedure.
Notice that, the model (3) represents an extension of Davé's model (1991) for fuzzy data with medoid prototypes and weighted mixed distance matrix. The distance from the noise cluster depends on the average distance among units and medoids The value of ρ may range between 0.05 and 0.5. In any case, the results do not seem very sensitive to the value of the multiplier ρ (Davé 1991). Due to the presence of δ 2 , units that are close to good clusters are correctly classified in a good cluster while the noise units that are away from good clusters are classified in the noise cluster.

Proposition 1
The solutions of (3) are: For c = C (4) becomes: The proof is in Appendix.

Remark 1 (Algorithm and computational issues)
1. The fuzzy clustering algorithm that minimizes (3) is built by adopting an estimation strategy based on the Fu and Albus heuristic algorithm (Fu and Albus 1977). Indeed, the alternating optimization estimation procedure cannot be adopted because the necessary conditions cannot be derived by differentiating the objective function in (3) with respect to the medoids. The fuzzy clustering procedure is illustrated in Algorithm 1.

Algorithm 1 Robust Fuzzy C-Medoids Clustering for Mixed Data (FCMd-MD-NC) algorithm
1: Fix C, max.iter and ρ, and generate randomly the degree matrix U ; 2: Set iter = 0; 3: Compute δ 2 ; 4: Pick initial medoids: Compute u i (i = 1, . . . , n) by using (4); 8: Compute w by using (5); 9: Select the new medoids:X cs , c = 1, . . . , C − 1, s = 1, . . . , S: 10: for c = 1 to C do 11: return ⇒X cs = X qs 13: end for 14: iter with the latter issue, it is possible to cope with the former two. First, the PAM approach requires that the distance matrix is computed only once at the beginning of the clustering process, and not at each iteration, thus decreasing dramatically the computing time required. Secondly, the search for the optimal medoids can be accelerated by "linearising" the clustering process, as in Krishnapuram et al. (2001), so that the complexity is linear in the number of units. 3. The degree of fuzziness of the resulting clusters is determined by m. The parameter can be pre-estimated by considering the usual fuzzy cluster-validity indices (see D'Urso and Maharaj 2009). However, since the medoid always has a membership of one in the cluster, raising its membership to the power of m has no effect on the medoid, while all other memberships decrease to 0. Thus, when m is high, the mobility of the medoids from iteration to iteration may be lost. For this reason, a value of m between 1 and 1.5 is recommended (Krishnapuram et al. 2001).

Remark 2 (Distances and dissimilarities)
One crucial decision in the clustering process for mixed data is the choice of suitable distance, or dissimilarity, measure for each attribute type. The choice is mainly heuristic, based on the data at hand and on the peculiar properties of each distance measure.
An admittedly non-exhaustive list of possible distance measures for several attribute types is reported in Table 2 (D'Urso and Massari 2019).
It should be highlighted that any kind of dissimilarity measure can be used in the proposed method. As in the standard non-hierarchical clustering algorithm e.g., k-means, k-medoids, the choice of the distance measure adopted in the clustering procedure is exogenous, so in the proposed method the choice of the distance measures for each attribute types is fixed beforehand. Any subset of variables can be managed with any of the dissimilarity measures presented in Table 2, and contribute to the "mixed" distance matrix in (2).
Remark 3 (Weighting system) The weights of the different attribute types in the clustering process are objectively provided by the model as shown in (5) and in the Appendix. An attribute type which displays a good separation into different groups should play a more significant role in clustering of data objects, against all other attribute types (Yeung and Wang 2002;Ahmad and Dey 2007). Indeed, the weight w s measures the within clusters similarity for the variables of the s-th type. Thus, the optimization procedure gives more relevance to the variable types capable to increase the within cluster similarity among the units.
Remark 4 (Determining the optimal number of clusters) A widely used cluster validity criterion for selecting C is the Xie-Beni criterion (Xie and Beni 1991), which can be suitably adapted for FCMd-MD-NC as follows: where C represents the set of possible values of C (C < n), and d(.) is the overall weighted distance (2). The numerator of I X B represents the total within-cluster distance. The ratio J /n measures the compactness of the fuzzy partition. The smaller this ratio, the more compact a partition with a given number of clusters. The minimum distance between centroids at the denominator of I X B is a measure of the separation between clusters. The greater this distance, the more separate a data partition with a given number of clusters. Therefore, letting the number of Canberra distance (Everitt et al. 2011) Mahalanobis distance (Everitt et al. 2011) Time series Auto-regressive based distance (Corduas and Piccolo 2008) Wavelet-based distance (Maharaj et al. 2010) Dynamic Time Warping (Berndt and Clifford 1994) Categorical data Jaccard distance (Everitt et al. 2011) Simple Matching Coefficient (Sokal 1958) Eskin dissimilarity index (Eskin et al. 2002) Geographical location data Geodesic distance (Karney 2013) Categorical-geographical data Geco distance (Hennig and Hausdorf 2006) Ordered sequences Hamming distance (Hamming 1950) Sequence Alignment Methods (Levenshtein 1966;Kruskal 1983) Fuzzy data External weighted distance Yang and Ko (1996) Internal weighted distance (D'Urso and Giordani 2006) Interval-valued data Distance for interval-valued data (D'Urso and Giordani 2004) Symbolic data Dissimilarity measure for symbolic data (Gowda and Diday 1991) clusters vary over the set C , the optimal number of clusters is identified in correspondence with the lower value of I X B .

Fuzzy profiling of the clusters
Results of cluster analysis can be summarized in the profiling phase where internal and external variables-i.e., variables involved and not involved in the cluster algorithm, respectively-are used to characterise and interpret the clusters (Everitt et al. 2011;Hair et al. 1998). In the case of fuzzy clustering algorithms, the (n × C) membership degrees matrix U = {u ic : i = 1, . . . , n, c = 1, . . . , C} can be used to properly weigh the observations on profiling variables and further describe the final clusters (D'Urso et al. , 2016. Let X = {x 1 , . . . , x n } be a quantitative variable observed on the sample. The weighted average of X in the c-th cluster is: As it can be seen, the greater is the membership degree of unit i to cluster c, the greater is the contribution of observation x i to the weighted average. Similarly, let Y = {y 1 , . . . , y n } be a categorical variable with L (L ≥ 2) categories. Let l be the generic category, and y il the observation in the i-th unit, which is equal to 1 if the category is observed on the i-th unit and 0 otherwise. The weighted proportion of the l-th category in the c-th cluster is: The concept of weighted averages and weighted proportions can be easily extended to other attribute types.

Simulation study
The aim of this simulation study is to highlight the capability of the FCMd-MD-NC model of correctly clustering objects in the presence of outliers. To this aim two distance matrices related to two groups of variables were generated, and outliers were added.
A dataset of n = 90 objects was simulated, with two numeric continuous variables, X 1 , X 2 and three numeric discrete variables, X 3 , X 4 , X 5 (S = 2). In particular, X 1 and X 2 are both generated from the Uniform distribution. (Different numbers of objects were considered with similar results). X 3 , X 4 and X 5 are discrete variables, with two, three and four values, respectively. Then, the set of variables is: Objects are grouped into three well separated and equal sized clusters according to both continuous and discrete variables.
The three clusters were obtained as follows (Figs. 1, 2): • cluster 1: X 1 and X 2 with Uniform density in the intervals A number of outliers equal to 6 (6.6%) and 12 (13.3%) was generated. Different simulation scenarios have been considered: 1. outliers for the continuous and the discrete variables; 2. outliers for the two continuous variables; 3. outliers for the three discrete variables.
The outliers of the continuous variables were generated according to Normal distributions; of the discrete variables according to discrete distributions. The euclidean distance was used to generate the two distance matrices related to the two groups of variables. We expected that, given the weighting structure, FMDd-MD-NC should be able to correctly classify the objects, despite the presence of outliers.
The correctness of clustering is evaluated by means of the Fuzzy Rand Index (FRI) to compare the obtained fuzzy partition with the reference crisp partition (30 objects in each cluster). The simulation study involved 100 replications (different numbers of replications were considered with similar results). Two values of the fuzziness parameter m were considered, 1.3 and 1.5. For each replication the weights computed for the two attributes and the FRI were collected. The results are reported in Table 3.
As for the correctness of the clustering, FCMd-MD-NC shows the expected robustness to outliers. Values of FRI are well above 0.90 in all scenarios and in all replications, thus indicating that the obtained fuzzy partitions are very close to the theoretical cluster partition. As for the attribute type weights, the two attribute types are weighted as expected according to their clustering structure. In all 100 replications FCMd-MD-NC attributed approximately equal weights to the two attributes in Scenario 1, greater weight to the continuous variables in Scenario 2 as the outliers of the numerical variables are assigned to the noise cluster, greater weight to the discrete variables in Scenario 3, as the outliers of the discrete variables are assigned to the noise cluster.
For comparative assessment the simulations have been run without the presence of the noise cluster. As it can be seen from Table 3, the results obtained with the model FCMd-MD-NC outperform the results obtained with the Fuzzy C-Medoids Clustering for Mixed Data (FCMd-MD) model. In all 100 replications FCMd-MD attributed approximately equal weights to the two attributes in Scenario 1, greater weight to the discrete variables in Scenario 2 as the outliers of the continuous variables are not assigned to the noise cluster, greater weight to the continuous variables in Scenario 3, as the outliers of the discrete variables are not assigned to the noise cluster.
The flipped behaviour of the two models in weighting the two distance matrices -of the continuos variables and of the discrete variables -deserves a comment. In the model FMDd-MD without noise cluster the outliers are assigned to the three clusters with great membership, and according to (20) they contribute to a great distance of the objects from the medoids; in the model FCMd-MD-NC with noise cluster the outliers are assigned to the three clusters with small membership being assigned to the noise cluster, and according to (20) they contribute to a small distance of the objects from the medoids, resulting in a greater weight.

Application: clustering of football players
The aim of this application is the clustering of football players based on attributes of different types. Professional football clubs invest a lot of resources in the recruitment of players. The clustering of football players based on performance data may be useful for prototyping successful players and for providing insights to football managers when assessing players.
Data for this application are drawn from whoscored. 1 The units are the players that have played in the "Serie A" tournement at least 200 min in the season 2018/19. The choice of 200 min follows from the choice of evaluating performances (goal, shots) per game (90 min). Less than 200 min would have led to re-proportioning the number of goals of a player, who scored a goal having played only 10 min, at 9 per game. The considered players were 397 from the initial number of 544. The teams of the season were: AC Milan, Atalanta, Bologna, Cagliari, Chievo, Empoli, Fiorentina, Frosinone, Genoa, Inter, Juventus, Lazio, Napoli, Parma Calcio 1913, Roma, Sampdoria, Sassuolo, SPAL 2013, Torino, Udinese. The most represented (greater than 3.0%) nationalities were Italy (39.3%), Argentina (6.3%), Brazil (5.8%), France (3.5%), Spain and Croatia (3.0%).
The considered variables are "Performance variables", "Success Variables" and "Position variables". The variables are grouped into four groups (S = 4) as described in Table 4 according to their type and meaning (Akhanli and Hennig 2017), in order to use an appropriate distance for each type and to obtain endogeneously the weight of each group in the definition of the clusters.
Performance ( Table 4; alongside with the weights computed in the clustering process for the different attributes types as in (5).

Data preprocessing
The Performance Upper level variables are represented per 90 min. The histograms of the performance variables are presented in Fig. 4.
The Success variables are represented as percentages. The Success Upper level variables Assists, Passes Aerials and Dribbles are percentages of accurate/successful actions over the total number of actions.
The Success Lower level compositions of Goals by body parts, situations and zones are percentages of success over the number of shots in the related sub-part (goal body part head over shots body part head).
The sub-parts of the Performance Upper level variable are transformed into percentages. If a player has scored 10 Goals per 90 min, of which 3 in situation out of box, 5 penalty area and 2 six yard box, the compound variables offensive Goals zones shows the values (30%, 50%, 20%) in the categories out of box, penalty area, six yard box.
The eleven position variables are binary variables.    The Performance Upper level count variables are standardised by average absolute deviation, whereas lower level compositions are standardised by the pooled average absolute deviation from all categories belonging to the same composition of lower level variables.
The Manhattan distance has been used for the Performance Upper level variables and for the Success Upper/Lower level variables.
The Performance Lower level variables are compositional data in the sense of Aitchison (Aitchison 1986), who set up an axiomatic theory for the analysis of compositional data. According to Akhanli and Hennig (2017), for the compositional percentage data the simple Manhattan distance is used as more appropriate. The distance measure used for the Position variables has been proposed in Hennig and Hausdorf (2006) and incorporates geographic distances. The proposed "geco" ("geographic distance and congruence") coefficient considers both the number of positions (geographic locations) shared by two players and the geographic relations of the occupied positions. Let A 1 , A 2 ⊆ R be the vectors of positions (presence-absence) of players a and b in the eleven positions of the football field (Fig. 3). Assume that there is a distance d R defined on R; d R is the geographic distance between positions (geographic locations) (euclidean distance in Table 5). For example, assuming segments of length 1 between adjacent positions on the same line, the geographic distance between positions DMC and MR is (1 2 + 1 2 ) = √ 2; between positions DMC and AMR (2 2 + 1 2 ) = √ 5 . The "geco" distance between players a in region A 1 and player b in region A 2 is defined as : where and A i denotes the number of elements in the geographical region of the i − th object. Then, d G is the mean of the average geographic distance of all units of A 1 to the respective closest unit in A 2 and the average geographic distance of all units of A 2 to the respective closest unit in A 1 . From the definitions it follows d G (A, A) The values are obtained by using Euclidean geometry based on Fig. 3

Results
Different values of the fuzziness parameter m (m = 1.3, 1.5, 2.0) and of the number of clusters C (C = 3, 4, 5, 6) have been considered. According to the Xie-Beni criterion, the optimal number of clusters is three, with a value of m = 1.3. The numerosity of the clusters is 109, 127, 112, 49 (noise cluster). The medoids of the three clusters are players 278, 288 and 91. The attributes' weights obtained during the optimization process are reported in the last column of Table 4. The Performance variables and the Success variables provide the greatest contribution to the clustering process.
Reported in Table 6 are the three medoids' characteristics, separately for each attribute type (a subset for the Success Upper level variables). The numeric Lower level compositional variables are not presented due to the small weight assigned in the clustering procedure.
The first cluster is represented by a player showing a discrete number of actions (Performance Upper level), successful (Goal/Shots, Assists accuracy, Aerials won and Dribbles successful) (Success Upper/Lower level).
The second cluster is represented by a player showing a limited number of actions (Performance Upper level), unsuccessful (Goal/Shots=0.0) (Success Upper/Lower level).
The third cluster is represented by a player showing a high number of actions per 90 min (Performance Upper level), very successful (Goal/Shots, Passes accuracy, and Dribbles successful) (Success Upper/Lower level).
The three clusters are similar with respect to the positions of the players. It is worth noting that the medoid players in cluster 1 and 3 have played in the Forward position. We observe that the player Ronaldo, who shows one of the highest values of Goal/Shots, is assigned to the third cluster.
As expected, the noise cluster shows heterogeneity of players with respect to the variables, as expected (otherwise there would be one more cluster).
In Figs. 5, 6, 7, and 8 the composition of the clusters with respect to some of the segmentation variables is reported. Percentages are calculated taking into account the membership degrees to each cluster, as explained in Sect. 2.2. The composition of the clusters is consistent with the three medoids.
We observe in Fig. 5 the prevalence of Goals in cluster 3. We observe in Fig. 6 that in the first cluster the successful percentage of shots that gives rise to a goal is played with body part left foot, situation setpiece and zone six yard box; in the second cluster with body part right foot, situation setpiece and zone six yard box; in the third cluster with body part right foot, situation open play and zone six yard box.
We observe in Fig. 7 that in the three clusters the greater percentage of Goals is played in the category of body part right foot, in situation open play and in zone penalty area; similarly for the Shots.
We observe in Fig. 8 the prevalent positions played in the three clusters. Figure 9 report a ternary plot with the membership degrees obtained with FCMd-MD-NC, for the three clusters. While there are several players that are fuzzy assigned, in particular to the first and the third cluster, over 50% of the units are allocated to a single cluster with a membership degree above 0.70.

Discussion and final remarks
In this paper, following the Partitioning Around Medoids (PAM) approach, we propose a robust fuzzy clustering with noise cluster and weighting system for mixed attributes. A simulation study is proposed. The model is used for analysing massive dataset in sport, i.e. for clustering football players based on their performance and positional attributes. The clustering model allows different types of variables, or attributes, to be taken into account. A weight is objectively assigned to the distance matrix associated to each set of attributes during the optimization process. The weights reflect the the relevance of each attribute type in the clustering results. A noise cluster neutralizes the effect of outliers. The simulation study has shown the ability of the model to weight properly the distance matrices of the attributes in the presence of outliers. The application to the clustering of football players on the basis of distance matrices of different attributes has shown the need to manage separately mixed type data. The obtained weights allow to understand which is the most The proposed clustering model has targeted some relevant issues in the research field. A mixed distance for the different attributes is considered; weights to distances related to different attribute types giving relevance to the variable types capable to increase the within cluster similarity are objectively provided by the model; a noise cluster represented by a noise prototype is introduced to achieve robustness with respect to outliers. For the practictioners, clustering of football players on the basis of playing characteristics, position, performance and success variables is relevant for clubs, either to drive team formation and selection of players, or for determining the value of a football player in the transfer window period (Behravan and Razavi 2021;Shelly et al. 2020;Narizuka and Yamazaki 2019). In football, performance assessment has typically been conducted to predict player's abilities, to rate player's performances, to drive their physical training or to explain a team's success (Mohr et al. 2003;Di Salvo et al. 2007;McHale et al. 2012). Future work will deal with other robust clustering solutions, the inclusion of financial variables of the player and of the team, the temporal aspect of the playing variables, if available, the interactions among team members and with the opposing players in the course of a game. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
Proof First, fixed w s , we determine the membership degrees u ic . We consider the Lagrangian function analyzing the objective function not splitted into C −1 clusters and the noise cluster: where u i = (u i1 , . . . , u ic , . . . , u iC ) and λ is the Lagrange multiplier. Therefore, we set the first derivatives of (11) with respect to u ic and λ equal to zero, yielding: From (12) we obtain: and, by considering (13): Finally, substituting (15) in (14) and taking into account the decomposition into the C − 1 clusters and the noise cluster we obtain u ic as in (4). Then, fixed u ic we derive w s . The Lagrangian function is: where w = (w 1 , . . . , w s , . . . , w S ) and ξ is the Lagrange multiplier. By setting the first derivatives of (16) with respect to w s and ξ equal to zero, we obtain respectively: From (17) we have: and using (18): Then, replacing (20) in (19), we obtain w s , as in (5).