New skills in symbolic data analysis for official statistics

2.1Motivation

In our program Clamix [6] for clustering symbolic data, they are represented by discrete distributions (n,𝐩) where 𝐩 is an empirical probability distribution and n is the number of original data units summarized by 𝐩. This representation has two important properties

In this paper, we will elaborate on the second observation.

For example, let us consider the population pyramids of the world’s countries. How to join the population pyramids of China and Vanuatu (see Fig. 1)?

When comparing two countries A and B we compare the shapes of their probability distributions 𝐩A and 𝐩B. But to determine the correct probability distribution of their union A∪B we need to know also the sizes nA and nB of countries A and B

2.2Aggregation

In an analysis of large data sets the aggregation is a standard way for reducing the size (complexity) of the data. Recently some books dealing with the theoretical and algorithmic background of the traditional aggregation (replacing values of a variable over a group by a single value) were published [7, 8, 9, 10, 11].

Data analysis programs provide aggregation functions such as means (arit, geom, harm, median, modus), min, max, product, bounded sum, counting, etc. [12]. Special care has to be given to variables measured in different measurement scales.

In theoretical discussion the traditional aggregation functions are usually “normalized” to the interval [0,1] – they take real arguments in [0,1]k and produce a value in [0,1], and satisfy the conditions: f⁢(𝟎)=0, f⁢(𝟏)=1, and monotonicity 𝐱⩽𝐲⇒f⁢(𝐱)⩽f⁢(𝐲). Often, in applications, also idempotency and symmetry are required.

The applications of traditional aggregation functions are used, besides determining a representative value for a group of measurements, mainly to combine partial criteria into a single criterion (multicriteria optimization and decision-making) or to express the membership degree in combined fuzzy sets.

A problem with traditional aggregation is that often too much information is discarded, thus reducing the precision of the obtained results.

A much better, preserving more information, summarization of original data can be achieved by representing aggregated data using selected types of complex data such as symbolic objects [2, 13], compositions [14], functional data [15], etc. In the Symbolic Data Analysis (SDA) framework, much work is devoted to the summarization process, for example, the function classic.to.sym in RSDA [16], and SODAS or SYR software.

2.3Mergeable summaries

In complex data analysis the measured values over a selected subset of units A are aggregated into a complex object Σ⁢(A) and not into a single value. Most of the aggregation theory does not apply directly. In our contribution, we present an attempt to start building a theoretical background of complex aggregation.

An interesting question is, which complex data types are compatible with the merging of disjoint sets of units

Selecting a name for this kind of summary we were inclined towards the term hierarchical or mergeable summary. Searching on Google we learned that the term mergeable summary was already proposed and elaborated by [17]. They enable parallelization in big data algorithms and stream processing. The summarization in big data is not deterministic and allows some errors. A summary is mergeable if the error and space (size of the summary) do not increase after the merge.

In this paper, we will discuss exactly mergeable summaries “without errors”.

2.4Exactly mergeable summaries

A summary Σ⁢(A) is an exactly mergeable summary if and only if it requires a fixed space of small size and satisfies the relation Eq. (1).

We can consider merging as a partially defined binary operation Σ⁢(A)*Σ⁢(B)=F⁢(Σ⁢(A),Σ⁢(B)). For mutually disjoint subsets A, B, and C we have

2.5Conclusions

In measurement theory [18, 19] measurement scales are divided into absolute, ratio, interval, ordinal, and nominal. The corresponding “best representatives” are count, geometric mean, average (arithmetic mean), median, and mode. The count is an exactly mergeable summary (2.4.1.1). So are the geometric mean (2.4.1.7) and the average (2.4.1.5), provided that we keep also the size of the corresponding set of units.

Median and mode are not exactly mergeable. A good exactly mergeable alternative is to use the corresponding frequency distribution (histogram or bar chart). In the case of a large number of categories, less frequent categories can be combined into a common category. By using the frequency distribution also for aggregation of numerical (ratio and interval) variables, we get a uniform representation for all types of variables.

3.1Context

This study concerns the Portuguese Labour Force Survey (LFS), analysing data from the 1^st trimester of 2008 and the 4^th trimester of 2010. We only consider people who were unemployed at the time of the survey (had no job and were looking for one), and focus on the Activity Time (in years)(AT) and Unemployment Time (in months) (UT). Disregarding records with missing values, and keeping only those from mainland Portugal (i.e. excluding Madeira and Azores), we end up with 1150 observations in 2008 and 1569 in 2010.

These micro-data were then gathered, in each case, on the basis of Gender (Mas, Fem), Region (North, Centre, Lisbon and Tagus Valley (LTV), South), Age-Group (Young: 15–24, Prime: 25–44, Mature: 45 and above) and Education (Basic or less, Secondary, Higher), leading to 58 sociological groups in 2008 (T1) and 68 in 2010 (T4) (as some of the 72 possible combinations do not occur) and which constitute the statistical units to be analysed.

We note that although the individuals at micro data level are not the same in 2008 and 2010, the aggregate units formed correspond to the same sub-populations – e.g., Young Women, from the North, with Secondary Education – and data are hence comparable at aggregate level.

The objective of this study is to cluster the aggregate units in each year, and compare the obtained partitions, trying to get insights about the dynamics between 2008 and 2010. For this purpose, we rely on the parametric model for interval-valued variables proposed in [20] and the model-based clustering methodology developed in [21]. Data aggregation as well as all analysis are done with R package MAINT.Data [22, 23].

3.2LFS Interval data

Let S={s1,…,sn}, be the set of n units under analysis. An interval-valued variable is defined by an application

where B is the set of all intervals of an underlying set O⊆I⁢R.

Let I be an n×p data array representing the values of p interval-valued variables on S. Each row of I corresponds to an element si∈S, represented by a p-dimensional vector of intervals, Ii=(Ii⁢1,…,Ii⁢p)=([li⁢1,ui⁢1],…,[li⁢p,ui⁢p]),i=1,…,n, as in Table 1.

In our case, for each aggregate unit, in each year, the minimum and maximum values of each of the Activity Time (AT) and Unemployment Time (UT) were recorded. As a result, each group is described by two intervals, that represent the within range of variation of the Activity Time and Unemployment Time in the corresponding year. Table 2 displays some rows of the 2008 and 2010 data arrays.

Comparing the interval data arrays for 2008 and 2010, we could observe that in several cases the UT interval became much wider within the two years, with the maximum value for UT showing a large increase. This is specially the case for groups with higher education levels, being not so frequent and clear in groups with basic education. Examples of such cases are the Mature Women from the South with Secondary education (Fem-South-Mat-Sec), for whom the UT interval went from [3,26] to [1,131], or Prime Man from the Centre with Superior education (Mas-Cen-Pri-Sup), for whom it changed from [7,13] to [1,171].

Figure 2.

Interval data for 2008 (left) and 2010 (right), groups with superior education (top), secondary education (center), and basic education (bottom).

Figure 2 displays the interval-valued data separately for the two years under analysis and the three education levels. We note that in 2008 the UT intervals (along the horizontal axis) differ considerably across education levels, getting larger as education level decreases. However, these intervals are much wider in 2010 than in 2008 for groups with superior education (upper figures), somehow for groups with secondary education (middle figures), but not so much for groups with basic education (lower figures). This effect reduces the differences in the UT intervals across education levels observed in 2008. In both years, the Activity Time variability increases as education level decreases (intervals along the vertical axis get wider). However, we do not observe remarkable changes in AT from 2008 to 2010, for any of the three education levels.

3.3Model

The value of an interval-valued variable Yj for each si∈S is usually defined by the lower and upper bounds li⁢j and ui⁢j of Ii⁢j=Yj⁢(si). For modelling purposes, however, we consider an alternative parameterisation, representing each interval Yj⁢(si) by the MidPoint ci⁢j=li⁢j+ui⁢j2 and Range ri⁢j=ui⁢j-li⁢j of Ii⁢j.

The Gaussian model (see [20]) assumes a multivariate Normal distribution for the MidPoints C and the logs of the Ranges, R*=ln⁡(R),(C,R*)∼𝒩2⁢p(μ,𝚺), with μ=[μCt,μR*t]t and 𝚺=(𝚺C⁢C𝚺C⁢R*𝚺R*⁢C𝚺R*⁢R*) where μC and μR* are p-dimensional column vectors of the mean values of, respectively, the MidPoints and the Log-Ranges, and 𝚺C⁢C,𝚺C⁢R*,𝚺R*⁢C and 𝚺R*⁢R* are p×p matrices with their variances and covariances.

We note that the model does not allow considering observations with degenerate intervals, where the range is null.

This model allows for the application of classical inference methods; however one should keep in mind that the MidPoint ci⁢j and the Range ri⁢j of the value of an interval-valued variable Ii⁢j=Yj⁢(si) pertain to one same variable, and must therefore be considered together. Specific configurations of the global covariance matrix allow taking into account the link that may exist between MidPoints and Ranges of the same or different variables. We consider the following configurations:

From the Normality assumption it obviously follows that imposing non-correlations with Log-Ranges is equivalent to imposing non-correlations with Ranges. In cases C⁢2, C⁢3 and C⁢4, 𝚺 may be written as a block diagonal matrix, after a possible rearrangement of rows and columns.

The mean vector and the variance-covariance matrix may be estimated by maximum likelihood. In the restricted configurations C⁢2, C⁢3, and C⁢4, estimation may be done block-wise (see [20]).

3.4Model-based clustering

Model-based Clustering considers the data as coming from a distribution that is a mixture of several components [24, 25, 26]. Each component is then associated with a cluster, characterized by a conditional density/mass function, and has a probability or “weight”. When the conditional probability is specified as the multivariate Gaussian, the model will be a finite mixture of multivariate Normals, known as the Gaussian mixture model.

The model parameters for each component, and the membership (posterior) probabilities of each unit, must be estimated, this is commonly accomplished by the Expectation-Maximisation (EM) algorithm [27]. This algorithm alternates an expectation (E) step, where the expectation of the log-likelihood at the current parameter estimates is computed, and a maximisation (M) step, where parameters are estimated by maximising the expected log-likelihood found in the E step.

Model-based Clustering of interval-valued data has been developed in [21], considering the Gaussian model described above (see Section 3.3), where the EM algorithm has been suitably adapted to the likelihood maximisation for the different covariance configurations.

Figure 3.

BIC values for 2008 (left) and 2010 (right).

The finite mixture model with k components for a data vector 𝐱 is defined as

where all weights πℓ are positive and π1+…+πk=1; θℓ are the parameters of the conditional distribution of component ℓ. In the Gaussian case, we consider two alternatives: a homoscedastic setup, where the covariance matrix is constant across components, and the conditional distribution is given by N⁢(μℓ,𝚺), and a heteroscedastic setup with one covariance matrix per component, with conditional distribution N⁢(μℓ,𝚺ℓ). Maximum likelihood parameter estimation involves the maximisation of the log-likelihood function

with 𝚪=(π1,…,πk,μ𝟏,…,μ𝐤,𝚺) in the homoscedastic case, and 𝚪=(π1,…,πk,μ𝟏,…,μ𝐤,𝚺𝟏,…,𝚺𝐤) in the heterocedastic case.

In Model-based Clustering of interval data, 𝐱𝐢=[𝐜𝐢t,𝐫𝐢*t]t is defined as the 2⁢p-dimensional vector comprising all the MidPoints and Log-Ranges for si. In the unrestricted case, the M-step formulas for 𝚺^ or 𝚺^ℓ are the classical ones; for the restricted configurations, 𝚺^ or 𝚺^ℓ are obtained maximising the likelihood for each block separately (see [20]).

The covariance configuration and the number of components k are selected as those that minimise the Bayesian Information Criterion (BIC) [28].

3.5Clustering of the LFS data

The method presented above was applied to the LFS data described in Section 3.1, separately for 2008 and 2010. However, we had to disregard six units in 2008 and two in 2010, since they presented a degenerate interval in at least one of the two variables. The units removed in 2008 were: F-Center-Young-Sup, F-LTV-Young-Sup, F-North-Young-Sup, F-South-Mature-Sup, M-Center-Prime-Sec, M-Center-Young-Sup, and in 2010: F-Center-Young-Sec and M-North-Young-Sup.

Figure 3 shows the BIC values for k=2,…,8, all four covariance configurations, and both homoscedastic and heteroscedastic setups.

We notice that the lowest BIC value is attained for configuration C⁢2 and a heteroscedastic setup in both cases, with k=5 in 2008 and k=4 in 2010.

The best partition obtained for 2008 is displayed in Table 3, and Table 4 gathers the mean values per component for the four indicators. Figure 4 displays the parallel coordinate plot for this partition, providing some insights for its characterisation.

Figure 4.

Parallel coordinate plot for 2008, best model.

Component 5 comprises Mature units with Basic education, and is characterised by high Activity and Unemployment Times, both with high variability (as conveyed by the log-ranges).

Component 1 is also composed by Mature units, now mostly with some education, with high Activity Times, lower Unemployment Times, both with not so high variability.

Component 3 is formed by Young and a few Prime units, it is characterised by low Activity and Unemployment Times, both with low variability.

Component 4 is mostly composed by Prime units with some education, and a few Young units. Its characterisation is similar to that of Component 3, but less pronounced; however, Unemployment time shows higher variability.

Finally, Component 2 comprises Prime units mostly with Basic education; it shows intermediate Activity and Unemployment Times, both with relatively high variability.

Table 5 describes the best partition obtained for 2010, Table 6 provides the corresponding mean values, and Fig. 5 displays the parallel coordinate plot.

Figure 5.

Parallel coordinate plot for 2010, best model.

Component 4 gathers Young units and is characterised by low Activity and Unemployment Times, both with low variability. It somehow corresponds to Component 3 of 2008, although now with no Prime units. In comparison, and as expected, Activity Time is on average lower and with lower variability. However, Unemployment Time is now on average higher, and with higher variability, which may be a consequence of the European sovereign debt crisis.

Component 1 is essentially formed by Mature units, irrespective of education level, coming from CP1 and CP5 of the 2008 partition. Activity and Unemployment Times are high, as expected, both with relatively high variability.

Component 2 is mostly composed by Prime units, with a few Mature. It is characterised by the high variability of both Activity and Unemployment Times, and high Unemployment Time.

Finally, Component 3, similarly to Component 4 of 2008, gathers mainly Prime and some Young units, with intermediate Activity and Unemployment Times, and respective variabilities.

We note that, unlike what happened in 2008, in 2010 education seems not to play a role in the formation of clusters. In particular, in 2008 there is a clear separation of Mature groups in two clusters, mainly distinguished in terms of education, while in the 2010 partition most Mature are grouped in one single cluster.

As a final remark, we observe that, although the individuals at micro data level are not the same in 2008 and 2010, the aggregate units formed correspond to the same sub-populations, allowing for a comparative analysis of the resulting partitions.

In both years, given the two variables used, age seems to be a driving force in the formation of clusters; in 2008 education also appears to play a role, but not so much in 2010.

It is noteworthy that, although some cluster correspondences may be identified, the partition is essentially not stable from 2008 to 2010, the sociological units gather in a different way. This may be the result of the change in the Portuguese labour market created by the economic crisis.

4.1Application on the data of international measurement of reading achievement among young students

In the first example, we are interested in the concordances and discordances of countries based on the proportion of students who scored at or above a high (high or advanced) level on the traditional paper reading assessment in the Progress in International Reading Literacy Study (PIRLS [31]). The PIRLS is an international assessment and research project designed to measure the reading achievement of fourth graders and the instructional practices of schools and teachers. The reading achievement scale is derived from several variables that measure the quality of reading. More detailed information can be found on the website of the cited reference (PIRLS 2016 User Guide, Chapter 4, p. 63). We focus our study on the 2016 survey data from fifty participating countries.

In this case, individuals are students and classes are countries. The collection P is the set of all 50 countries participating in the study. The x category includes high and advanced levels of traditional reading achievement. For each country c, fc⁢(x) indicates the proportion of students who achieved a high or advanced level of traditional reading achievement. To obtain the function gx⁢(c,P), we divided the interval [0,1] into five equidistant subintervals and, based on the values fc′⁢(x),c′∈P, obtained the distribution of countries shown in Table 7. Figure 6 presents the positioning of the countries in the plane with axes fc⁢(x) and gx⁢(c,P).

As it can be seen from Table 7, almost half of the countries have more than 40% and up to 60% of students achieving at least a high level of reading literacy on paper, so we can say that each of these countries is concordant with the collection of all countries. This interpretation is captured by the measure Sconc1. In Fig. 6, these countries are positioned in the upper right place.

Figure 6.

Positioning of the countries based on their values of the functions f and g for the traditional (paper) reading.

The second chosen s-concordant measure Sconc2 can be interpreted as the area of the rectangles in the plane in Fig. 6, for each country c determined by the values fc⁢(x) and gx⁢(c,P). The larger area of the rectangle corresponding to the country expresses a higher concordance of this country with the collection of all countries. Among the 23 countries with a proportion of well-skilled students between 0.4 and 0.6, Chinese Taipei had the highest proportion, and the value of its Sconc2 is the largest for this country. Thus, taking into account both a good fit between the country and observed literacy achievement (in our case high or advanced literacy) which results in high fc⁢(x), and the frequent occurrence of such fits among countries, we can say that Chinese Taipei is the country for which the concordance with a collection of all countries is the highest.

There are seven countries where the proportion of well-skilled students is above 0.6 (positioned at the bottom right of Fig. 6). Based on their good fit with the category (high fc⁢(x)) and the rare frequency of such good fits, these most discordant countries can be identified by the highest s-discordance value Sdisc3. In PIRLS 2016 data, the Russian Federation has the highest proportion of well-skilled students in paper reading. Since there are only 7 out of 50 countries in our data that have more han 60% of well-skilled students, the Russian Federation is the most discordant country in the collection of all participating countries.

4.2Application on textual data with the comparison with Tf-Idf

In the second example, we explore the use of the discordance measure as an alternative to the well-known Tf-Idf (term frequency - inverse document frequency) measure in Text Mining. The vector space model for representing text collections of documents is represented by a term-document matrix in which the documents are represented by columns and the index terms used to index the document collection are represented by rows. The basic idea of Tf-Idf (see [32]) is to characterize a category (here presence of the term x) of a class of documents by its “relevance" compared to the other classes. The relevance of the term x to the class is high if

There are several variants of the Tf-Idf measure. We will use its basic form, which is for the term x and class of documents c defined as

where

n⁢(x,c) is the number of documents in the class c in which the term x occurs

|c| is the number of documents in the class c

K is the total number of classes in the collection of documents

k is the number of classes in the collection of documents for which there is at least one document in which the term x occurs

Note that the definition of Tf-Idf does not take into account the differences between classes due to the number of documents in which the term occurs. In the definition of the s-discordance, however, these differences are included in the function gx⁢(c,P)

We illustrate the difference between these measures in the collection of 15 documents (book titles) from the field of data mining (DM documents), linear algebra (LA documents), and one document combining these two fields (Table 8) [33].

The list of index terms consists of terms that appear in at least two documents, with so called stop words or words commonly used in a language sorted out, and word variants mapped in their base form. In this way, a list of 16 index terms was created.

The values of the function gx⁢(c,P) are determined using five equidistant subintervals [0,0.2], (0.2,0.4], (0.4, 0.6], (0.6, 0.8], (0.8, 1.0]. In our case of two classes gx⁢(c,P) can take only two values: 1, if fc⁢(x) falls in the same interval for both classes; and 0.5, otherwise. In order to compare the values of Sdisc3⁢(c,P;x) and Tf-Idf(x,c), these values are normalized to the interval [0,1], i.e. divided by the maximum possible values of these functions in the case of two classes. Calculated values for both classes and for each of 16 terms are presented in Table 9. Identifying a term as characteristic for a class if it has a higher value of observed measures, we can see that both measures recognize relevant terms for classes of DM and LA documents. Values of measures Sdisc3⁢(𝐷𝑀,P;x) and Tf-Idf(x,𝐷𝑀) are greater for terms x related to the field of data mining (such as classification, clustering, data, document, information, mining, retrieval, and text) while measures Sdisc3⁢(𝐿𝐴,P;x) and Tf-Idf(x,𝐿𝐴) are greater for terms x related to the field of linear algebra (algebra, linear, matrix, space, and vector).

To distinguish relevant terms for classes based on defined measures of s-discordance and Tf-Idf we define α-relevance for these measures. For a given threshold α>0, Tf-Idf is considered α-relevant if it is at least α, and analogously, the Sdics3 value is called α-relevant if it is at least α. By choosing α=0.3, we show some differences between the measures in the recognition of relevant terms shown in Table 9. In this case, the Sdisc3 measure recognized relevant terms more successfully than Tf-Idf, because unlike Tf-Idf, Sdisc3 recognized the terms information and retrieval as relevant for the class DM and also recognized the terms space and vector as relevant for the class LA. Reason for that is that the Tf-Idf measure equalizes classes in which a specific term appears at least once, while the s-discordance measure makes a difference according to the frequency of the term appearance in classes.

4.3Discussion and further work

We have shown two examples of possible applications of the new measures s-concordance and s-discordance from different contexts.

In the first example, these measures were used to measure the concordance and discordance of a single country with the collection of all countries. We used data from a PIRLS 2016 survey where we focused on at least a high level of traditional paper reading assessment. Measures of s-concordance and s-discordance are used to compare the country to the collection of all countries. Since in the definitions of s-concordance and s-discordance we use the distribution of classes included with the function g, which in application is usually based on the empirical distribution, a more detailed study of its influence is needed in the future.

In the second example, we used the measure of s-discordance to identify relevant terms that are characteristic of a particular class of documents. We compared our results with those obtained using the Tf-Idf measure. In our case, the s-discordance measure identified more relevant terms for classes.

The Tf-Idf measure detects relevant terms for a class if that term occurs only in that class, while the s-discordance measure detects relevant terms for a class if the frequency of their occurrence in that class is greater than in other classes. Because of this property, the s-discordance measure could be used to extend a weighting of terms when representing a document in a vector space model to improve classification performance. It could also be used in sentiment analysis to automatically capture a lexicon by calculation of s-discordance measure of index terms for classes of documents with positive and negative sentiment.