Visualizing Proﬁles of Large Datasets of Weighted and Mixed Data

: This work provides a procedure with which to construct and visualize proﬁles, i.e., groups of individuals with similar characteristics, for weighted and mixed data by combining two classical multivariate techniques, multidimensional scaling (MDS) and the k -prototypes clustering algorithm. The well-known drawback of classical MDS in large datasets is circumvented by selecting a small random sample of the dataset, whose individuals are clustered by means of an adapted version of the k -prototypes algorithm and mapped via classical MDS. Gower’s interpolation formula is used to project remaining individuals onto the previous conﬁguration. In all the process, Gower’s distance is used to measure the proximity between individuals. The methodology is illustrated on a real dataset, obtained from the Survey of Health, Ageing and Retirement in Europe (SHARE), which was carried out in 19 countries and represents over 124 million aged individuals in Europe. The performance of the method was evaluated through a simulation study, whose results point out that the new proposal solves the high computational cost of the classical MDS with low error.


Introduction
One of the most important goals in visualizing data is to get a sense of how near or far objects are from each other. Often, this is done with a scatter plot, because the Euclidean distance is the only one that our brain can easily interpret. However, scatter plots cannot always be obtained from raw data, nor is the Euclidean distance always the appropriate one to be computed on raw data. This may be the case when comparing a high number of variables, where a dimension reduction is usually necessary to better see the proximities between objects, or when working with more complex datasets, such as weighted mixed data or functional data, where other distances are preferred to the Euclidean one. For instance, survey data coming from macro-surveys at national and cross-national levels are rather complex datasets of weighted and mixed data. They are composed of variables of different natures, such as binary, multi-state categorical and numerical variables; and as result of a multi-stage sampling methodology, they each include a weighting variable, so that each individual represents a group of different size for the target population. Another added complexity may be their large or very large sample size (10 4

or larger).
Multidimensional scaling (MDS) is one of the most extended methodologies to analyze and visualize the profile structure of data, and can address some of those problems. This dimensionality reduction technique takes a dissimilarity or distance matrix as the input and produces a pictorial representation of the data in a Euclidean space, similar to a scatter plot. An important limitation when working with large and very large datasets is that it relies on the eigendecomposition of the full distance matrix between objects or individuals, thereby requiring large quantities of memory and very long computing times. Paradis (2018) [1] proposed an approach to avoid the limitations of the standard MDS procedure, which is based on a random selection of a small number of observations and the application of standard MDS with one or two dimensions. In a second step, the remaining Note that our proposal starts by clustering the individuals in the original space, and next, we use MDS to visualize the clustering in the Euclidean space. For that reason we use a clustering algorithm able to cope with mixed data, k-prototypes (although other methods can be used [3,4]). Another possibility would be to start with the MDS representation and next do the clustering on the Euclidean space using k-means clustering. In any case, when working with large datasets, a "fast" MDS is required in order to reduce both computational time and memory use. Thus, in any case, the idea of applying the methodology to a small portion of the data and to project the remaining individuals on the MDS-map remains. In this work, we explore the former; that is, we cluster the individuals in the original space instead of doing the clustering on the MDS-map, since an additional motivation was to incorporate Gower's distance in the k-prototypes clustering algorithm.
The performance of our proposal was evaluated through a simulation study based on a real dataset of weighted and mixed data that came from the Survey of Health, Ageing and Retirement in Europe (SHARE), carried out in 19 countries. The analysis was applied to 60,020 individuals, who represent a target population of more than 124 million Europeans aged 55 years or over.
The work proceeds as follows: In Section 2 we review weighted MDS; in Section 3 we present the proposed methodology; in Section 4.1 we apply the method to data coming from SHARE database; Section 4.2 contains the simulation studies; and we conclude in Section 5.

Materials and Methods
In this section we give a general overview of classical MDS for weighted mixed data, introduce some useful notation and present Gower's interpolation formula.
The purpose of MDS is to construct a set of points in a Euclidean space whose interdistances are either equal (classical MDS) or approximately equal (ratio, interval or ordinal MDS) to those in a given matrix of dissimilarities, so that the interpoint distances approximate the interobject dissimilarities. That is, given an n × n matrix D (2) ij is the squared dissimilarity between objects i, j, for i, j = 1, . . . , n, the objective of MDS is to search for a configuration of n points on a set of orthogonal axes, so that the l 2 -distances between the coordinates of these n points coincide with the corresponding entries in D (2) . These coordinates are called a Euclidean configuration/map/representation or MDS configuration/map/representation of D (2) . Several possible measures of approximation between interpoint distances and interobject dissimilarities can be used, each yielding a different MDS configuration. In this work, these coordinates are obtained via spectral decomposition. General contextual references are [5][6][7][8].
Sample surveys often employ multistage sampling schemes which involve unequal selection probabilities at some or all stages of the sampling process. We refer to this situation as a weighted context, where each individual can represent a population group of different size. In this framework, classical methods of data analysis which assume simple random sampling may no longer be valid, and weighting may appear as the only or best alternative. Albarrán et al. (2015) [9] reviewed the extension of classical MDS concepts to the weighted context.

Weighted MDS
Let {x i , i = 1, . . . , n} be n p-dimensional vectors which contain the observations or measurements of p variables for n different individuals and D (2) be the matrix of squared distances between n individuals, with entries δ 2 (x i , x j ), 1 ≤ i, j ≤ n. Remember that the information contained in the p variables can be either of a quantitative or qualitative nature, or both, hence it is crucial to select an appropriate dissimilarity function in the computation of D (2) in order to incorporate all the statistical information contained in the data.
Additionally, since each individual in the dataset can represent a group of a different size of the target population (weighted context), we have a vector of weights w = (w 1 , . . . , w n ) , such that w i > 0, for i = 1, . . . , n, and 1 w = 1, where 1 is the n × 1 vector of 1 s. Note that in the case of simple random sampling-that is, when each individual in the dataset represents a group of the same size of the target population-w i = 1/n, for i = 1, . . . , n, and classical MDS formulae are recovered.
Suppose that we are interested in obtaining an MDS representation of D (2) , provided that D (2) satisfies the Euclidean requirement.
Given w, we define an n × n diagonal matrix D w = diag(w), whose diagonal is the vector of weights, and J w = I − 1w is the w-centering matrix, where I is the n × n identity matrix.
The w-centering matrix J w is an orthogonal projector with respect to D w , idempotent (that is, J 2 w = J w ) and self-adjoint with respect to D w (that is, J w D w = D w J w ). In the weighted context, the doubly w-centered inner product matrix is given by whose standardized version is given by which is called standardized inner product matrix. The condition that D (2) satisfies the Euclidean requirement is equivalent to imposing that G w is positive semidefinite, which means that there exists a matrix Y w such that G w = Y w Y w . In the weighted context, Y w is called a w-centered Euclidean representation of D (2) and satisfies the following two properties: The squared l 2 -distances between the rows of Y w coincide with the corresponding entries in D (2) ; that is, for each pair of individuals i, j, we have that Matrix Y w is the w-weighted MDS representation of D (2) and is computed by means of the spectral decomposition of matrix F w defined in Equation (1). That is, given F w = UΛU , where Λ is a diagonal matrix with the eigenvalues of F w , ordered in descending order, and U is the corresponding matrix of eigenvectors (in column), whose rows are the principal coordinates of n individuals, and its columns are the principal axes of this representation.

Gower's Distance
Gower's similarity coefficient [10] is one of the most popular similarity measures and perhaps the easiest way to obtain a distance measure when working with mixed data. It is the Pitagorean sum of three similarity coefficients, one for each type of variable. In particular, it uses Jaccard's coefficient for binary variables, the simple matching coefficient for multi-state categorical variables and range-normalized city block distance for quantitative variables. Given two p-dimensional vectors x i and x j , Gower's similarity coefficient is defined as: where p 1 is the number of quantitative variables, R h is the range of the h-th quantitative variable, a is the number of positive matches, d is the number of negative matches for the p 2 binary variables, α is the number of matches for the p 3 multi-state categorical variables and p = p 1 + p 2 + p 3 . The entries of matrix D (2) are computed as: and Gower (1971) [10] proved that Equation (3) satisfies the Euclidean requirement.
Other metrics for mixed data can be considered, although the k-prototypes algorithm should be modified accordingly. For instance, a more robust metric that can overcome some of the shortcomings of Gower's is related metric scaling (RelMS) by Cuadras (1998) [11], which was used in [9] to obtain robust profiles in weighted and mixed datasets. However, in this work, we prefer to illustrate our methodology by using Gower's metric due to the computational complexity of RelMS.

Gower's Interpolation Formula
A very useful tool to project new data points onto a given MDS configuration is Gower's interpolation formula [8], which was extended to the weighted context by [12]. The following proposition can be derived from the Theorem 1 in [12]. This formula is a key tool in the visualization algorithm that we propose. Proposition 1. Let E be a set of n individuals; w = (w 1 , . . . , w n ) a vector of weights, such that w i > 0, for i = 1, . . . , n, and 1 w = 1; and D (2) be a matrix of squared distances between the n individuals satisfying the Euclidean requirement and Y w the w-weighted metric scaling representation of D (2) .
Given a new individual n + 1 for whom the squared distances to n individuals of E are known, δ = (δ 2 n+1,1 , . . . , δ 2 n+1,n ), its principal coordinates can be computed as: where g w = diag(G w ) is a row vector containing the diagonal elements of G w , D w = diag(w), G w = Y w Y w and Λ is the diagonal matrix containing the eigenvalues of matrix F w defined in (1).
Proof. The squared distance between the individual n + 1 and any individual i ∈ E is given by δ 2 n+1,i = (y n+1 − y i )(y n+1 − y i ) = y n+1 y n+1 − 2y n+1 y i + y i y i . In matrix notation, we have that Post multiplying expression (5) by D w Y w and after operating, we have that: Note that y n+1 2 1 D w Y w = 0 since 1 D w = w and w Y w = 0. Therefore, the principal coordinates of individual n + 1 are given by:

Methodology
In this section we discuss a methodology for visualizing profiles for large datasets of weighted and mixed data.
Among all the approaches proposed for visualizing data, MDS is one of the most common techniques. However, we find the classical MDS algorithm a limited tool when visualizing large datasets, since it requires very large CPU time or large computing memory when dealing with large distance matrices. Delicado and Pachón-García (2020) [13] showed both: that the time needed to compute MDS as a function of the sample size increases notably when using cmdscale() R function (in stats package by R Development Core Team) and that at least 400 MB of RAM memory is required to store the distance matrix when there are 10,000 observations.
There have been several attempts to solve the scalability problem, such as steerable multidimensional scaling [14], incremental MDS [15], relative MDS [16], FastMap [17], MetricMap [18], landmark MDS [19], the diagonal majorization algorithm [20] and uniform manifold approximation and projection [21]. SteerMDS proposed by Williams et al. [14] is based on a spring-mass model, introduced by Chalmers [22] (see also [23] for a samplingbased variant of the algorithm). These methods calculate lower-dimensional coordinates by iteratively minimizing a cost or stress function that is proportional to the distance between the current coordinates and the given dissimilarities. Incremental MDS is similar to the previous methods, although it focuses on overall shape instead of local details. Relative MDS combines MDS with the learning vector quantization clustering method. FastMap, MetricMap and Landmark MDS approximate classical MDS by solving MDS for a subset of the data and fit the remainder to the solution. Platt [24] studied these methods and concluded that Landmark MDS was the fastest and most accurate of the them. The diagonal majorization algorithm is a modification of the Guttman majorization algorithm [25] that is used to minimize the stress function, and is able to save computing time taking into account several factors (see [26]). Finally, the uniform manifold approximation and projection is a learning technique for dimension reduction that can be used for visualization, which is indirectly related to MDS and more closely to Isomap.
The visualization method that we propose is based on classical MDS applied to a random portion of the dataset plus the projection of the remaining individuals via Gower's interpolation formula. Thus, our proposal shares the idea behind several of the existing methods of applying MDS to a portion of the dataset, but differs from them in the projection tool used to obtain the final MDS representation.

The Visualization Algorithm
We start by summarizing the proposed method, and later we remark on some important aspects concerning the clustering algorithm and the feasible implementation of Gower's interpolation formula.
The starting point of the algorithm is a large dataset of weighted and mixed data; that is, our dataset was composed of several variables of different natures (binary, multi-state categorical and numerical variables) plus a weighting variable containing the individual weights. Remember that weights are given exogenously and are related to the survey sampling technique. No pre-processing of the data is needed, except for the normalization of weights to sum to 1.

1.
Select a small random sample, using the weights to produce a more informative sample, that is, trying to follow as much as possible the sampling scheme. Depending on the size of dataset, this selection can be 2.5, 5 or 10% of total observations. Let us denote this small sample by X n×p , where n is the number of individuals and p the number of variables.

2.
Compute the distance matrix between the rows of X n×p using Gower's distance Formula (3).

3.
Carry out the k-prototypes clustering algorithm in order to find the different clusters and label the individuals accordingly. Determine the number of clusters in the dataset by the "elbow" rule.

4.
Obtain the principal coordinates of the labeled individuals through weighted MDS.

5.
Compute the representatives (or centroids) of the clusters. This can be done by calculating the weighted mean or weighted median of those point-coordinates belonging to the same cluster in the MDS configuration. 6.
Project the rest of the individuals (the remaining 97.5, 95 or 90%) onto the MDS configuration using Gower's interpolation formula. 7.
Finally, from the MDS configuration, assign the new points to an existing cluster based on the closest centroid (according to Euclidean distance) and label/color them accordingly.
Once all points have been assigned to a cluster, it is possible to visualize the clusters on the MDS configuration, and thus, to see the proximities between them. Finally, a profile is defined as the "average" member of each cluster. To do so, for categorical variables the mode can be computed and the mean and the median are good options for quantitative variables.

Some Important Remarks
Next, we introduce a few observations on the steps presented above.

On the Clustering Algorithm
Classical hierarchical clustering can handle also mixed data by using Gower's coefficient. However, it is well known that it struggle as sample size increases. Reference [2] proposed a clustering algorithm called k-prototypes which is quite efficient, O(T + 1)kn, where n is the number of observations, k the number of clusters and T the number of iterations; and it has been previously used by [27] for profile construction. The k-prototypes algorithm is based on a dissimilarity measure that takes into account both quantitative and categorical variables. Let X r 1 , . . . , X r p 1 , X c (p 1 +1) , . . . , X c p be the variables available in the dataset, where X r h stands for a quantitative variable and X c h for a categorical one; then the dissimilarity between two individuals x i , x j ∈ R p can be measured by where δ(p, q) = 0 for p = q and δ(p, q) = 1 for p = q, and γ ≥ 0 is a coefficient that measures the influences of numeric and categorical variables. Note that when γ = 0, clustering only depends on numeric variables. As it happens in any non-hierarchical clustering method, the number of clusters, k, must be determined in advance. To do so, a variety of techniques exist, and sometimes determining the optimal number of clusters is an inherently subjective measure that depends on the goal of the analysis. Due to the large size of the dataset, we decided to use the "elbow" method, instead of the average silhouette width or other time-consuming criteria. To apply the "elbow" method, we ran the algorithm for different values of k and calculated the cost function for each run. Then, we plotted the cost function in a line graph; and the point where a turning point (or "elbow") was observed, that is, the point at which the cost function levels off, was selected as the optimal k value. Recently, Aschenbruck and Szepannek (2020) [28] examined the transferability of cluster validation indices to mixed data and evaluated them through simulation studies. They concluded that the average silhouette width was the most suitable with respect to both runtime and determination of the correct number of clusters. However, these conclusion rely on rather small datasets (≤400 individuals).
In this work, we introduce two particularities to the k-prototypes algorithm so that it can cope with Gower's distance and weighted datasets.
First, instead of the d 2 measure described in (6), we used Gower's distance Formula (3), which is a very popular dissimilarity measure for mixed data and satisfies the Euclidean requirement [10]. The second particularity introduced refers to weighted datasets. Since the ultimate difference from the standard algorithm is in centroid calculation, weighted averages of quantitative variables and weighted modes of categorical variables are used, instead of standard means and modes.
In what follows, we summarize the adapted version of the k-prototypes algorithm: • Initial prototypes selection. Select k distinct individuals from the dataset as the initial centroids. • Initial allocation. Each individual of the dataset is assigned to the closest prototype's cluster, according to distance (3). • Reallocation. The prototypes for the previous and current clusters of the individuals must be updated, taking into account individual weights. This repeats until there is no reallocation of individuals.
Some authors pointed out possible inaccurate clustering results when using Hamming distance and proposed other alternatives ( [29,30]). In our simulations, we did not experiment with such situations (see Section 4.2 for graphical representations of the cost function).
However, we propose to adapt the k-prototypes algorithm, although its convergence properties in combination with using Gower's distance were not further investigated.

On Gower's Interpolation Formula
As mentioned before, MDS's final configuration is obtained by projecting the remaining observations via Gower's interpolation formula on the initial configuration. But how can this be implemented in a feasible way? The idea is the following: We call M m×p the remaining observations after selecting X n×p .
• Split M m×p row-wise into partitions M 1 , . . . , M , equally sized, with perhaps the exception of M , which can be smaller. The number of partitions is set to be (n + m)/ , where × is the size of the largest distance matrix that a computer can calculate efficiently [13]. • Apply Gower's interpolation formula to each matrix M j (j = 1, . . . , ) and store the coordinates. The application of Gower's interpolation formula to a matrix M j whose rows are m "new" individuals is rather straightforward from Formula (4).
With the same notation as in Section 2.3, let ∆ be the m × n matrix whose rows contain the squared distances of the "new" m individuals to the n individuals of E . Then, the principal coordinates of these "new" m individuals can be computed by: where 1 m is a m × 1 vector of ones. This is simple, but strongly advantageous, since one of the main problems in MDS is memory consumption when computing distance matrices. Gower's interpolation formula allows us to iteratively get the MDS configurations without facing memory problems, because we are reducing the size of the corresponding matrices that contain the squared distances to the n individuals of the existing MDS configuration. This aspect is studied in Section 4.2.
All the functions listed above require a distance matrix as the main argument to work with. In case data are not in the distance/dissimilarity matrix format, R-functions dist, daisy and gower.dist may be of help. Moreover, some of the previous packages provide their own functions for calculating distances.
With respect to MDS configurations, the wcmdscale() function is used in this work. It is based on function cmdscale (base package of R-stats) and it can use point weights. Points with high weights will have stronger influences on the result than those with low weights. Setting equal weights will give ordinary multidimensional scaling.
Concerning the k-prototypes algorithm, the R package clustMixType by [37] contains the function kproto() needed to perform this technique. As mentioned before, we modified this function to achieve our desired goal, that is, to cope with Gower's distance and weighted datasets.

Results
In this section we present a real data application and a simulation study to evaluate the performance of the visualization algorithm.

Application
Prior to analyzing the performance of the proposed methodology, we applied the algorithm described in Section 3 to a dataset from the Survey of Health, Ageing and Retirement in Europe (SHARE). SHARE is a research infrastructure, formed by a panel database of micro data, for studying the effects of health, social, economic and environmental policies over the lifetimes of European citizens and beyond. It was founded in 2002 and is coordinated centrally at the Munich Center for the Economics of Aging (MEA), the Max Planck Institute for Social Law and Social Policy (Munich, Germany). The SHARE interview is ex-ante harmonized and all aspects of the data generation process, from sampling to translation, and from fieldwork to data processing, have been conducted according to strict quality standards. As a result, SHARE has the advantage of encompassing cross-national variations in public health and socioeconomic living conditions of European individuals, and becoming a major pillar of the European research area, with over 9000 researchers registered as SHARE users. See their website http://www.share-project.org/ for further details. This work uses Wave 6 of SHARE, which was conducted in 2015 in 18 European countries and Israel. It asked questions ranging from an individual's financial situation to his/her self-perception of health.

Description of the Dataset
The dataset to be analyzed consists of 60,020 observations and 13 variables, and includes a weighting variable that scales to represent over 124 million elderly individuals in Europe. Descriptions of variables are in Table 1. It is important to remark that the last four correspond to health and wellbeing indexes and were not in the original dataset, but created by [27] from the aggregations of 30 variables. Higher values of the indices correspond to situations of greater vulnerability. Although the process of creating those indices and other details about the dataset are described in their work, here we give a brief summary of the indices for better interpretation of the profiles.
The dependency index summarizes the loss of personal autonomy. It includes variables that reflect difficulties in performing activities of daily living, instrumental activities of daily living, mobility limitations and so on.
The self-perception of health index is quite a subjective measure that captures some aspects of subjective wellbeing. This index includes variables related to what an individual thinks about their health rather than their physical reality, including how satisfied the interviewee is with their life, whether they are feeling depression, etc.
The physical health and nutrition index captures the risk of an individual to suffering or developing serious health problems. It contains information related to body mass index, grip strength, nutrition and chronic conditions of the individual.
The mental agility index captures the mental acuteness of the respondents and is related to cognitive functions. It contains the results of tests around numeracy, orientation and linguistic fluency.
Next, we proceed with the visualization and construction of the profiles for this dataset of weighted and mixed data. Gender "Male", "Female" CT Ages "55-60", "61-65", "66-75", "76+" B Employment status "Employed", "Not working" B Marital status "Has no spouse", "Has a spouse" CT Education "No education", "Primary", "Secondary", "University" B Household in financial distress "Yes", "No" CT Household receives "Payments and no benefits", "No benefits benefits or has payments? and no payments", "Benefits and payments", "Payments and no benefits" C Dependency index Form 0 to 10 C Physical health and nutrition index From 0 to 10 C Self-perception of health index From 0 to 10 C Mental agility index From 0 to 10 B = binary, CT = categorical, C = continuous.

Visualization of Profiles and Findings
The result of applying the visualization algorithm is shown in Figure 2. In order to select the number of clusters, the algorithm was run for k = 2, . . . , 10. The cost function showed an "elbow" at around k = 3 or k = 4. We investigated having 3 or 4 profiles and selected k = 4, since k = 3 led to overly broad results. Panel (a) contains the MDS configuration based on Gower's interpolation formula computed from a portion of 2.5% of the data. Red triangles correspond to the mapped points of the random selection, and gray crosses stand for the projected ones. In panel (b) we show the pictorial representation of the clusters, where blue triangles stand for cluster centroids and circles represent their confidence regions, whose radii were computed as the 90th percentile of the (Euclidean) distance between each point and the corresponding centroid. We can observe beforehand two distinct groups of individuals, that is, a set of points grouped on the left (cluster 3) and another bunch of crowded points on the right, which splits into three small clusters. As will be seen later, cluster 3 corresponds to the least disadvantaged profile, whereas clusters 2 and 4 contain the most vulnerable individuals, according to the descriptive variables.
When original variables (and not just a matrix of distances or similarities) are available, it may be of interest to determine the influences of these original variables on the MDS dimensions, i.e., to determine which variables explain the most the homogeneity within groups. It was therefore necessary to calculate some correlation coefficient (or association measure) between the principal coordinates and the variables. We used Pearson's correlation coefficient for continuous variables, Spearman's correlation coefficient for ordinal variables and Cramer's V measure of association for nominal variables. Results are shown in Table 2, where it can be seen that categorical variables, such as gender, job (Employment status), fdistress (household in financial distress) and paybene (household receives benefits or has payments) have great influences on the axes. For instance, the first principal coordinate is mostly determined by the variables job (employment status) and paybene (household receives benefits or has payments), whereas fdistress (household in financial distress) and gender are influential to the second and third principal coordinates. Quantitative variables do not seem to have many impacts on the axes, although the variable mental (mental agility index) influences the second coordinate somewhat. This was expected somehow, since Gower's metric tends to give more weight to qualitative variables.   A graphical way to see the influences of the original variables on the construction of the profiles is to assign colors to the categories of the original variables (or groups of values, in the case of the quantitative ones) and represent the individuals colored accordingly in the MDS configuration. In Figure 3 we show the most representative projections of the MDS configuration regarding those variables more correlated/associated with the principal axes. In the following we summarize the main characteristics of the four estimated profiles shown in Figure 2 (right panel). The results were acquired by taking the weighted mode for qualitative variables and the weighted mean and median for quantitative ones. See also Figures 4 and 5 and  From the previous description we can sort the profiles from the most to the least disadvantaged, in terms of levels of health and socioeconomic wellbeing. In particular, P4 defines the group with the lowest levels of health and wellbeing, closely followed by P2. On the other hand, P1 defines a medium-low social vulnerability profile and P3 a profile of low risk of social vulnerability.
From Figure 4 and Table 3 we observe that profiles P4 and P2 more resemble each other compared to both other profiles across "dependency," "mental agility" and "physical health," whereas "self-perception of health" is more balanced. Remember that higher values of the indices correspond to situations of greater vulnerability.  P2 and P4 are the most vulnerable, although P2 is still better than P4. These profiles are especially vulnerable, since they have high percentages of individuals belonging to the 76 years or older group and they rank mostly bad in the indices, mainly due to the loss of physical or mental/intellectual autonomy produced by gradual aging. However, P4 is the one that ranked worst in all wellbeing indices, with mean differences of 0.31-1.39 points with respect to P2. Moreover, people in P4 were more likely to live alone than people in P2 (49.95% in front of 39.81%). Therefore, people in P4 were of the type that require more assistance and or extensive help in order to carry out common everyday actions. Additionally, we see that 96.69% of households in P4 were in financial distress, in front of 63.88% of the households in P2. Still comparing P2 and P4, another interesting finding is that the P4 individuals were more pessimistic with respect to their health. At least, this is what the distribution of the "self-perception of health" index seems to indicate (with a median difference of three points). This variable expresses what individuals think about their health: whether they are satisfied with their life, feeling depression, etc. Clearly, the least disadvantaged group was P3, which heavily skewed towards younger individuals, since almost 60% of the individuals were aged between 55 and 60, and those over 76 made up less than 5% of this profile. In addition, 57.27% of them were women, and they tend to work until an advanced age. In addition, the indices indicate good health for this group, which may be related to the fact that continuing working in later life has been proved to be correlated with positive health outcomes [38]. Low values in wellbeing indices are consistent with the fact that people in this profile have no health-related benefits.
Finally, we see that education is one of the factors that most influences the profiles' differences. Notice that most disadvantaged profiles have in common lower education levels (secondary, primary or none), whereas the most advantaged always have a high percentage of members who attended university. This finding supports the idea that society's disparities are explained by differences in education among individuals. Education is often used as a proxy for socioeconomic status and its impact on later-life health outcomes is well researched [39,40].

Profiles across Europe
The variable "country," which represents the country of residence of the individuals, was not included in the process described so far. However, it seemed interesting to see how the profiles are distributed across Europe, particularly for the most disadvantaged ones. Thus, for each country we computed the percentage of people belonging to each profile. Results are shown in Figure 6, where it can be seen that a combination of the least disadvantaged profiles is predominant in most European countries, except for Portugal, Estonia and Poland.
A more interesting question is to find out whether SDG-3 of United Nations 2030 Agenda, that is, to ensure healthy lives and promote wellbeing for all at all ages, is fulfilled in these developed countries. To do so, we analyzed the percentages of the most disadvantaged profiles in these EU countries and found that they are concentrated in the Southern, Central and Eastern European countries. In particular, in countries such as Portugal, Spain, Italy and Poland it is estimated that over 20% of their population of 55 years old or older belong to P4. Regarding P2, it is estimated to be over 20% in 66.6% of the EU countries, geographically distributed in Central-Eastern European countries.

Simulation Study
Recall from Section 3 that we introduced Gower's interpolation approach as a way to solve the scalability problem. However, is this true? Have we dealt successfully with the problem? We consider those questions next. The aim of this section is to evaluate the discrepancies between two MDS configurations, the classical one computed from the complete dataset and the one obtained with our algorithm. Besides, we analyze the time required to compute both configurations.
In this simulation study, the starting point is the dataset; hence, the computing times reported in this work are not comparable with those reported in [1], where the starting points were the first two principal coordinates.

Design of the Simulation Study
The analysis was carried out on samples of the dataset used in Section 4.1, and it was structured as follows: • Sample sizes. To evaluate the elapsed time, a total of nine different sample sizes were used: n = 500, 1000, 5000, 10,000, 20,000, 30,000, 40,000, 50,000 and 60,000. Discrepancies between two MDS configurations were evaluated through a total of six different sample sizes: n = 500, 1000, 2000, 3000, 4000 and 5000. • Portion of data. Recall that the first step of the algorithm was to select a small sample from the data. We wished to see whether there exists a significant difference in using 2.5, 5 or 10% as the initial portion. • Each scenario was the combination of a sample size and a portion of data and was repeated 100 times.
For each repetition, we computed elapsed time; the errors in configuration eigenvalues, measured in terms of mean squared error (MSE); and the cophenetic correlation between Euclidean configurations, as a measure of distortion between two distance matrices [41].
We used a personal computer to run the simulation analysis, whose technical specifications were: computer processing unit: Intel ; core™ i5-4200U CPU @ 1.60 GHz 2.30 GHz; RAM: 4 GB (Santa Clara, CA, USA).
Tables 4 and 5 contain the mean values per scenario.

Time to Compute MDS
As stated before, in this section we analyze the divergence of the time required to compute MDS configurations when interpolating data points and when using all points (complete MDS from now on) to construct the coordinates. Table 4 contains the simulation study results. Clearly, we see that MDS based on Gower's interpolation was, by far, the fastest one. Note that the table does not display any result for the complete MDS approach from 10 4 observations onward. This was because of memory limitation problems. It could not even store the distance matrix. Nevertheless, we see that for small sample sizes we obtained a greater elapsed time for the complete MDS case, although we were not considering the time required to get the distance matrix (if we considered this, the difference would become even larger). As sample size increases, so do the time and memory needed. For example, for sample sizes of 500 to 1000, the elapsed time was almost six times higher for the complete MDS than for Gower's approach. Another observation is that as sample size increases, it is more interesting to consider a portion of 2.5 or 5% rather than 10%. We will see in the next section whether there is a significant difference in choosing 2.5% instead of 10%.

Discrepancies in MDS Configurations
Here we compare the distortion of both MDS configurations by calculating the cophenetic correlation between them. That is, we want to know how different the configuration obtained through Gower's interpolation formula is from the classical MDS (complete MDS). We also analyze how much the eigenvalues of both approaches differ by means of the mean square error (MSE), computed on the normalized positive eigenvalues of the two MDS configurations. As stated before, due to memory limitations, we could not do the comparisons for sample sizes greater than 10 4 observations. Results are shown in Table 5, where we can see, first, that the normalized eigenvalues do not differ significantly (average MSE of 0.045 and median of 0.034), and second, that inter-distances between both configurations are preserved (average cophenetic correlation of 0.789 and median of 0.824). Overall, this simulation study reveals that discrepancies between configurations decrease as the sample size increases, and that, for a given sample size, they tend to decrease when we consider a higher sample portion. Although one can think of selecting a greater sample portion to reduce them, we must keep in mind that this requires more time to get the final MDS configuration.

Cost Function
Finally, we illustrate the convergence of the cost function of the k-prototypes algorithm, for k = 4 and several of the scenarios described above. In Figure 7 we depict the mean and median values of the corresponding cost functions for which convergence was achieved in a rather small number of iterations.

Conclusions
In this paper we presented a methodology for visualizing profiles of large datasets of weighted and mixed data that combines two classical multivariate techniques, MDS and a clustering algorithm. The clustering algortihm is used to partition the individuals in the original space, and once individuals are labeled accordingly, MDS is used to visualize the clusters in a Euclidean space, where it is easier to explore the proximities among clusters. The scalability problem is solved by means of Gower's interpolation formula. Finally, profiles are obtained as the "average" member of each cluster, where the mode is considered for categorical variables and the mean or the median for numerical ones.
In this work we have illustrated this methodology by means of the k-prototypes clustering algorithm, although other clustering procedures able to cope with mixed data can be used [4]. Additional motivations for using k-prototypes were to adapt it to the weighted context and to use Gower's metric as the dissimilarity measure. However, some authors pointed out possible inaccurate clustering results when using Hamming distance and proposed other alternatives ( [29,30]). However, we did not experiment with such situations in our case study.
We tested the procedure through a simulation study, where we evaluated computational costs (elapsed time, errors in configuration eigenvalues) and discrepancies between classical and Gower's interpolated MDS configurations (cophenetic correlation). The results show that MDS based on Gower's interpolation formula solves the main issue of classical MDS (high computational cost) with few errors.
We applied the proposed methodology to find several profiles of healthy life and wellbeing in the European Union, using the respondents to the Survey of Health, Ageing and Retirement in Europe survey data. We found that the most vulnerable group (people with the poorest health and wellbeing) contained elderly people who were less educated, living alone, had financial problems, had low personal autonomy and had impairments in cognitive abilities. This profile is more likely to be found in Southern and Eastern European countries.
Although Gower's metric is widely used when dealing with mixed datasets, it presents several shortcomings. For example, it gives more weight to categorical variables than to quantitative ones; it does not take into account the possible associations or correlations between variables; and it is not robust against atypical data (see [9,42,43]).
An interesting direction for future research would be to consider other metrics, such as the hybrid dissimilarity coefficient by Jian and Song [30], and their modification of the k-prototypes algorithm, although it is more time consuming than the classical k-prototypes. Nevertheless, this hybrid dissimilarity still ignores the association between variables. Another possibility would be to use related metric scaling (RelMS) to tailor the metric and incorporate it into the clustering algorithm. RelMS was introduced in [11,44] and provides more robust and stable configurations than Gower's, but it has a high computational cost, which should be reduced so that it can be workable for large datasets. Another interesting direction for future research would be to explore the possibility of combining our methodology with clustering algorithms based on archetype and archetypoid analysis, which focuses on extreme individuals instead of centroids ( [45,46]).  Acknowledgments: The authors are grateful to Pedro Delicado for his useful comments and suggestions.

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