---
title: "Partial ranking with two-step clustering"
output:
  html_document:
    highlight: haddock
    theme: lumen
    toc: yes
    toc_float:
      collapsed: no
  html_notebook:
    fig_caption: yes
    highlight: haddock
    theme: lumen
    toc: yes
    toc_float:
      collapsed: no
  pdf_document:
    keep_tex: yes
    latex_engine: lualatex
    toc: yes
  word_document:
    toc: yes
---

Workflow

1. Loading necessary libraries and data
  + Loading libraries 
  + Print citations
  + Loading necessary data
2. Investigation of two-step clustering
  + Cluster tendency (columns)
    + Calculating Hopkins index
    + Hierarchical clustering for indicators
    + Apply multidimensional scaling to indicators
  + Cluster tendency (rows)
    + Calculating Hopkins index
    + Hierarchical clustering for indicators
    + Apply multidimensional scaling to indicators
3. Model reduction
  + Principal component analysis
    + Determination of PC **R** (resources)
    + Determination of PC **E** (environment)
    + Determination of PC **C** (connectivity)
    + Determination of PC **O** (output)
    + Correlation test for principal components
  + Factor analysis
4. Data reduction
  + Organize the factors into 3 clusters
  + Organize the factors into 4 clusters
5. Determining partial rankings
  + Correlation between rankings for 3 clusters
  + Correlation between rankings for 4 clusters
6. Summary and conclusion


# 1. Loading necessary libraries and data

## Loading libraries 
```{r, warning=FALSE, message=FALSE}
# Correlation analysis
library(corrgram) # Correlogram
library(corrplot) # Correlation plots
source("http://www.sthda.com/upload/rquery_cormat.r") # Correlation matrix representation

# Dimension reduction methods
library(psych) # Factor Analysis, Principal Component Analysis
library(factoextra) # Distance metrics
library(MASS) # Multi-dimensional Scaling

# Cluster analysis
library(fclust) # Cluster analysis
library(clustertend) # Hopkins coefficients
library(cluster) # Cluster plots
library(fpc) # Cluster statistics
library(mclust) # Multi-step clustering

# Complex models
library(lavaan) # CFA, SEM models
library(semPlot) # Plot models

# Further statistical libraries
library(Hmisc) # Misceleaus statistical functions
library(seriation) # Seriation

# Data manipulation
library(reshape2) # Data transformation
library("dplyr", character.only = TRUE) # Data manipulation
library(knitr) # Formatted tables

# Maps
library(maps) # Mapping fuctions
library(mapdata) # Mapping data
library(maptools) # Mapping tools

# Interactive 3D figures

library(shiny) # Interactive figures 
library(htmltools) # HTML widgets
library(rgl) #3D shaping
library(rglwidget) #3D shaping
library(ggplot2) # Complex figures
library(ggrepel) # Complex figures
library(ggdendro) # Dendogram
library(d3heatmap) # Heat map
library(igraph) # Graph representation
library(ggsci) # Simpson palette
library(plotly) # Plotly - interactive plot
```

## Citations

```{r, warning=FALSE, message=FALSE}
# Correlation analysis
base::print(citation("corrgram"),style="text") # Correlogram
base::print(citation("corrplot"),style="text") # Correlation plots

# Dimension reduction methods
base::print(citation("psych"),style="text") # Factor Analysis, Principal Component Analysis
base::print(citation("factoextra"),style="text") # Distance metrics
base::print(citation("MASS"),style="text") # Multi-dimensional Scaling

# Cluster analysis
base::print(citation("fclust"),style="text") # Cluster analysis
base::print(citation("clustertend"),style="text") # Hopkins coefficients
base::print(citation("cluster"),style="text") # Cluster plots
base::print(citation("fpc"),style="text") # Cluster statistics
base::print(citation("mclust"),style="text") # Multi-step clustering

# Complex models
base::print(citation("lavaan"),style="text") # CFA, SEM models
base::print(citation("semPlot"),style="text") # Plot models

# Further statistical libraries
base::print(citation("Hmisc"),style="text") # Misceleaus statistical functions
base::print(citation("seriation"),style="text") # Seriation

# Data manipulation
base::print(citation("reshape2"),style="text") # Data transformation
base::print(citation("dplyr"),style="text") # Data manipulation
base::print(citation("knitr"),style="text") # Formatted tables

# Maps
base::print(citation("maps"),style="text") # Mapping fuctions
base::print(citation("mapdata"),style="text") # Mapping data
base::print(citation("maptools"),style="text") # Mapping tools

# Interactive 3D figures

base::print(citation("shiny"),style="text") # Interactive figures 
base::print(citation("htmltools"),style="text") # HTML widgets
base::print(citation("rgl"),style="text") #3D shaping
base::print(citation("rglwidget"),style="text") #3D shaping
base::print(citation("ggplot2"),style="text") # Complex figures
base::print(citation("ggrepel"),style="text") # Complex figures
base::print(citation("ggdendro"),style="text") # Dendogram
base::print(citation("d3heatmap"),style="text") # Heat map
base::print(citation("igraph"),style="text") # Graph representation
base::print(citation("ggsci"),style="text") # Simpson palette
```

## Loading necessary data

The original U21 data is used.

```{r, warning=FALSE, message=FALSE}
load("U21_2014.RData")
rownames(U21_rank) <- rownames(U21_filtered)
orig_mtx <- as.matrix(U21_filtered)
#U21_filtered: 50 by 24 dataframe of original, U21 data
#U21_rank: 50 by 1 matrix (vector) U21 ranking
#orig_mtx: 50 by 24 matrix, original U21 indicators of 50 countries
```

# 2. Two-step clustering

As a first step, we examine whether the data are suitable for clustering. Clustering should be checked for both variables and countries.

## Cluster tendency of indicators

### Calculating Hopkins indicator

The *Hopkins index (H)* helps to control clustering. If the *Hopkins index (H)* is greater than 0.5, it is not worth clustering.

When clustering the indicators, we choose *Spearman correlation* as a similarity indicator and *1-Spearman correlation* as a distance indicator, because the rank correlation also shows nonlinear relationships.

```{r,warning=FALSE, message=FALSE}
vars <- U21_filtered[,c("R1","R2","R3","R4","R5","E1","E2","E3","E4","C1","C2","C3","C4","C5","C6","O1","O2","O3","O4","O5","O6","O7","O8","O9")]
rownames(vars)<-rownames(U21_filtered)
countries<-t(vars)
colnames(countries) <- rownames(U21_filtered)

 dist.m<-1-cor(vars,method="spearman")
paste("H:= ",hopkins(dist.m, n=nrow(dist.m)-1, byrow = F, header = F))
```

Since the Hopkins index is less than 0.5, the indicators can be clustered.

The next step is to plot cluster tendency in order to determine the preliminary number of clusters.

```{r, warning=FALSE, message=FALSE, fig.width = 7, fig.height = 5, fig.cap="**Figure 1. ** Cluster tendency of indicators"}
# Determine cluster tendency
clustend <- get_clust_tendency(scale(dist.m), nrow(dist.m)-1)
# Plot cluster tendency
clustend$plot + 
  scale_fill_gradient(low = "yellow", high = "blue")
```

The **Figure 1** shows, that 2-4 clusters (which are the groups of indicators here) can be specified, but this should be confirmed by further studies.

### Hierarchical clustering for indicators

The next step is to use hierarchical clustering for the indicators. The result of the clustering is shown in **Figure 2**.

```{r,  warning=F, message=FALSE,fig.width = 7, fig.height = 5, fig.cap="**Figure 2.** Hierarchical clustering for the indicators"}
hc <- hclust(dist(dist.m), method = "complete")
hcdata <- dendro_data(hc, type="rectangle")
o.cols <- hcdata$labels$label
ggplot(hcdata$segments) + 
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend))+
  geom_text(data = hcdata$labels, aes(x, y, label = label), 
            hjust = 1.3, vjust=0.5, size = 4) +
  labs(x="", y="", title="Hierarchical clustering for the indicators") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()) +
  coord_flip()
```

Hierarchical clustering also makes 2-4 clusters (in this case a group of indicators) probable. Displaying the distance matrix provides additional help in determining the number of clusters. **Figure 3** shows the distance matrix for the indicators, where the blue cells represent values close to -1 and the dark red cells represent values close to 1.

```{r, message=FALSE,warning=FALSE, fig.width = 5, fig.height = 5, fig.cap="**Figure 3** Distance matrix for the indicators"}
cormat <- rquery.cormat(dist.m, graphType="heatmap")
```

Van Mechelen et al. (1993) seggested the following criteria for clustering:

1. Within-cluster dissimilarities should be small.

2. Between-cluster dissimilarities should be large.

3. Clusters should be fitted well by certain homogeneous probability models such as the Gaussian or a uniform distribution on a convex set, or, if appropriate, by linear, time series or spatial process models.

4. Members of a cluster should be well represented by its centroid.

5. The dissimilarity matrix of the data should be well represented by the clustering (i.e., by the ultrametric induced by a dendrogram, or by defining a binary metric "in same cluster/in different clusters").

6. Clusters should be stable.

7. Clusters should correspond to connected areas in data space with high density.

8. The areas in data space corresponding to clusters should have certain characteristics (such as being convex or linear).

9. It should be possible to characterize the clusters using a small number of variables.

10. Clusters should correspond well to an externally given partition or values of one or more variables that were not used for computing the clustering.

11. Features should be approximately independent within clusters.

12. All clusters should have roughly the same size.

13. The number of clusters should be low.

Cluster statistics are calculated for 2-4 clusters.
```{r, warning=FALSE, message=FALSE}
cluster.stats(dist.m,cutree(hc, k=2))
cluster.stats(dist.m,cutree(hc, k=3))
cluster.stats(dist.m,cutree(hc, k=4))
```

Unfortunately, based on the above results, it is not possible to decide how many (2, 3 or 4) clusters can be formed on the indicators. So the interpretability of the clusters will decide this problem.

### Multidimensional scaling for indicators

To examine the interpretability of the clusters, the next step is to use multidimensional scaling (MDS) for the indicators. First, we form two clusters, the results of which are shown in **Figure  4**.


```{r, warning=F, message=F, fig.cap="**Figure 4.** Multidimensional scaling for indicators considering 2 clusters"}
mds.isomds <- isoMDS(dist(dist.m), tol=1e-3, k=10, maxit = 200)

plotdata <- mds.isomds$points
plotdata <- as.data.frame(plotdata)
plotdata$names <- rownames(mds.isomds$points)
plotdata$cut <- cutree(hc, k=2)
ggplot(plotdata, aes(V1, V2, label=names)) + 
  geom_point(aes(colour=factor(cut)), size=2.3) +
  geom_text_repel(aes(colour=factor(cut)), size=4) + 
  scale_colour_discrete(name = "Clusters") +
  labs(x="", y="", title="Multidimensional scaling for indicators considering 2 clusters") + theme_bw()
```

As shown in **Figure 4**, the two clusters are difficult to interpret. The first cluster included 19 indicators, while the second included only 3 indicators. 


**Figure 5** shows the result of multidimensional scaling in the case of specifying 3 clusters. 

```{r, warning=F, message=F, fig.cap="**Figure 5.** Multidimensional scaling for indicators considering 3 clusters"}
mds.isomds <- isoMDS(dist(dist.m), tol=1e-3, k=10, maxit = 200)

plotdata <- mds.isomds$points
plotdata <- as.data.frame(plotdata)
plotdata$names <- rownames(mds.isomds$points)
plotdata$cut <- cutree(hc, k=3)
ggplot(plotdata, aes(V1, V2, label=names)) + 
  geom_point(aes(colour=factor(cut)), size=2.3) +
  geom_text_repel(aes(colour=factor(cut)), size=4) + 
  scale_colour_discrete(name = "Clusters") +
  labs(x="", y="", title="Multidimensional scaling for indicators considering 3 clusters") + theme_bw()
```

Based on **Figure 5**, the 3 clusters (i.e. 3 groups of indicators) are also difficult to interpret. There were 15 indicators in the first cluster, 3 in the second cluster and 6 in the third cluster. However, the indicators measuring the relationships (C1, C2, C3 and C4) have already been clustered.

In the next step, we will look at how to group the indicators if we create 4 clusters. The result is shown in **Figure 6**.

```{r, warning=F, message=F, fig.cap="**Figure 6.** Multidimensional scaling for indicators considering 3 clusters"}
mds.isomds <- isoMDS(dist(dist.m), tol=1e-3, k=10, maxit = 200)

plotdata <- mds.isomds$points
plotdata <- as.data.frame(plotdata)
plotdata$names <- rownames(mds.isomds$points)
plotdata$cut <- cutree(hc, k=4)
ggplot(plotdata, aes(V1, V2, label=names)) + 
  geom_point(aes(colour=factor(cut)), size=2.3) +
  geom_text_repel(aes(colour=factor(cut)), size=4) + 
  scale_colour_discrete(name = "Clusters") +
  labs(x="", y="", title="Multidimensional scaling for indicators considering 4 clusters") + theme_bw()
```

If four clusters are specified, 5 indicators will be added to the first cluster, 10 to the second, 3 to the third, and 6 to the fourth. Purely only the variable group of the relationships containing the variables C1-C4 appears in the fourth cluster, in the case of the other clusters the composition of the indicators is very mixed. Four clusters (groups of variables) are also difficult to interpret, so it is necessary to reduce the number of variables, which can be done, for example, by principal component analysis.

## Cluster tendency of countries

### Calculating Hopkins indicator

The *Hopkins index (H)* helps to control clustering. If the *Hopkins index (H)* is greater than 0.5, it is not worth clustering. In the case of clustering rows (countries), distance is the Euclidean distance.

```{r,warning=FALSE, message=FALSE}
vars <- U21_filtered[,c("R1","R2","R3","R4","R5","E1","E2","E3","E4","C1","C2","C3","C4","C5","C6","O1","O2","O3","O4","O5","O6","O7","O8","O9")]
rownames(vars)<-rownames(U21_filtered)
countries<-t(vars)
colnames(countries) <- rownames(U21_filtered)
dist.c<-as.matrix(dist(vars))
paste("H:=",hopkins(dist.c, n=nrow(dist.c)-1, byrow = F, header = F))
```

Since the Hopkins index is less than 0.5, countries (i.e., rows) can be clustered.

The next step is to plot the clustering in order to determine the preliminary number of clusters.

```{r,warning=FALSE,message=FALSE, fig.width = 7, fig.height = 5, fig.cap="**Figure 7.** cluster tendency of countries"}
# Cluster tendency
clustend <- get_clust_tendency(scale(dist.c), nrow(dist.c)-1)
# Plot the cluster tendency
clustend$plot + 
  scale_fill_gradient(low = "yellow", high = "blue")
```

Based on **Figure 7**, 3-4 clusters (ie a group of countries) are likely, but this needs to be confirmed by further studies.

### Hierarchical clustering of countries

The next step is to apply hierarchical clustering for countries. The result of the clustering is shown in **Figure 8**.


```{r, warning=F, message=FALSE, fig.width = 8, fig.height = 7,fig.cap="**Figure 8.** Result of hierarchical clustering for countries"}
hc <- hclust(dist(dist.c), method = "complete")
hcdata <- dendro_data(hc, type="rectangle")
o.rows <- hcdata$labels$label

ggplot(hcdata$segments) + 
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend))+
  geom_text(data = hcdata$labels, aes(x, y, label = label), 
            hjust = 1, vjust=0.7, size = 3) +
  labs(x="", y="", title="Hierarchical clustering for countries") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()) +
  coord_flip()
```

Hierarchical clustering makes the existence of 3-4 clusters (here a group of countries) probable.

Displaying the distance matrix provides additional help in determining the number of clusters. **Figure  9** shows the distance matrix for countries, with blue cells representing values close to -1 and dark red cells representing values close to 1.

```{r, message=FALSE,warning=FALSE, fig.width = 5, fig.height = 5, fig.cap="**Figure 9.** Distance matrix for countries"}
cormat <- rquery.cormat(dist.c, graphType="heatmap")
```

Based on the distance matrix shown in **Figure 9**, 3-4 clusters are probable. In the next step, we compare the cases where we create 3 or 4 clusters for the countries.

```{r, warning=FALSE, message=FALSE}
cluster.stats(dist.c,cutree(hc, k=3))
cluster.stats(dist.c,cutree(hc, k=4))
```
Nevertheless, based on the results, it is not possible to decide how many (3 or 4) clusters can be specified for the countries. Then the number of clusters should be decided by their interpretability. 

### Multidimensional scalig for countries

To examine the interpretability of clusters, the next step is to use multidimensional scaling for countries. First, we specify three clusters (i.e. three groups of countries), the results of which are shown in **Figure 10**.

```{r, warning=F, message=FALSE, fig.cap="**Figure 10.** Multidimensional scaling for countries considering 3 clusters"}
mds.isomds <- isoMDS(dist(dist.c), tol=1e-3, k=10, maxit = 200)

plotdata <- mds.isomds$points
plotdata <- as.data.frame(plotdata)
plotdata$names <- rownames(mds.isomds$points)
plotdata$cut <- cutree(hc, k=3)

ggplot(plotdata, aes(V1, V2, label=names)) + 
  geom_point(aes(colour=factor(cut)), size=2.3) +
  geom_text_repel(aes(colour=factor(cut)), size=4) + 
  scale_colour_discrete(name = "Clusters") +
  labs(x="", y="", title="Multidimensional scaling for countries considering 3 clusters") + theme_bw()
```

**Figure 10** shows that the first cluster included 21 countries, the second 15 countries and 14 countries.

In the next step, we will look at how to specify 4 clusters. The result is shown in **Figure 11**.

```{r, warning=F, message=FALSE, fig.cap="**Figure 11.** Multidimensional scaling for countries considering 4 clusters"}

mds.isomds <- isoMDS(dist(dist.c), tol=1e-3, k=10, maxit = 200)

plotdata <- mds.isomds$points
plotdata <- as.data.frame(plotdata)
plotdata$names <- rownames(mds.isomds$points)
plotdata$cut <- cutree(hc, k=4)

ggplot(plotdata, aes(V1, V2, label=names)) + 
  geom_point(aes(colour=factor(cut)), size=2.3) +
  geom_text_repel(aes(colour=factor(cut)), size=4) + 
  scale_colour_discrete(name = "Clusters") +
  labs(x="", y="", title="Multidimensional scaling for countries considering 4 clusters") + theme_bw()
```

In the case of specifying 4 clusters, 21 countries were included in the first cluster, 11 in the second, 14 in the third, and 4 in the fourth. Compared to the previous case (when we specified three clusters), the countries in the first and in the third clusters did not change. The fourth cluster included the countries (Denmark, Sweden, Switzerland, United States) that belonged to the second cluster in the previous case. The principles formulated by Van Mechelen (1993) are mostly fulfilled by 4 clusters.

**The results of cluster tendency:** It is clear that clustering can be made by both variables and countries (see Hopkins index). However, the resulting clusters are quite difficult to interpret. The groups of variables contain mixed indicators, so model reduction is needed before grouping countries.

# 3. Model reduction

## Principal Component Analysis

The first step is to determine the number of principal components / factors. The **Figure 12** shows the scree plot, which can be used to specify the number of latent variables (i.e. principal components).

```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 12 scree plot**"}
scree(U21_filtered)
```
Based on **Figure 12**, maximum 5 principal component (PC), while maximum 3 factor should be kept.

U21 contains 4 main groups of indicators: resources (R1 - R5), environment (E1 - E4), connectivity (C1 - C6) and outputs (O1 - O9). In the case of principal component analysis, it is worth starting from these predefined groups of indicators. 

### *R (resources)* PC

First, we define the principal component *R*, which contains the indicators R1-R5, that describe the resources of each country.

```{r, warning=FALSE, message=FALSE}
Rs<-U21_filtered[,c("R1","R2","R3","R4","R5")]
KMO(cor(Rs))
fitR <- principal(Rs,nfactors=1)
summary(fitR) # print variance accounted for 
fitR$communality
pcaR <- prcomp(Rs, scale = TRUE)
eigs <- pcaR$sdev^2
varR <- eigs[1] / sum(eigs)
paste("Variance rate=",varR)
```

Since the KMO value is greater than 0.5, principal component / factor analysis can be performed. All communality is greater than 0.25, which means that the original variables fit the principal component correctly. Furthermore, the variance ratio is greater than 0.5, which means that at least 50% of the information remains even after model reduction. Based on these, the *R* principal component created above is appropriate and can be used in further analyzes.

### *E (environment)* PC

Next, we specify the principal component *E*, which contains the indicators E1-E4, that describe the environmental indicators of each country.

```{r, warning=FALSE, message=FALSE}
Es<-U21_filtered[,c("E1","E2","E3","E4")]
KMO(cor(Es))
fitE <- principal(Es,nfactors=1)
summary(fitE) # print variance accounted for 
fitE$communality
pcaE <- prcomp(Es, scale = TRUE)
eigs <- pcaE$sdev^2
varE <- eigs[1] / sum(eigs)
paste("Variance rate=",varE)
```

Based on the above results, it can be seen that the communality of E2 is less than 0.25 (0.189), so E2 does not fit the PC properly. However, if we exclude E2 in the next step, the communality value of E1 will be decreased less to 0.25. Most indicators can be saved, if E3 will be excluded, instead of E2. Thus, in the next step, we omit E3.

```{r, warning=FALSE, message=FALSE}
Es<-U21_filtered[,c("E1","E2","E4")]
KMO(cor(Es))
fitE <- principal(Es,nfactors=1)
summary(fitE) 
fitE$communality
pcaE <- prcomp(Es, scale = TRUE)
eigs <- pcaE$sdev^2
varE <- eigs[1] / sum(eigs)
paste("Variance rate=",varE)
```


After excluding E3, the communality value is greater than 0.25 for all variables. The KMO value is greater than 0.5. However, the value of the variance ratio is less than 0.5 (0.477), which means that at least 50% of the information is lost after the model reduction.

### *C (connenctions)* PC

Next, we specify the principal component *C*, which contains the indicators C1-C6, that describe the connectivity and communications indicators of each country.

```{r, warning=FALSE, message=FALSE}
Cs<-U21_filtered[,c("C1","C2","C3","C4","C5","C6")]
KMO(cor(Cs))
fitC <- principal(Cs,nfactors=1)
summary(fitC) 
fitC$communality
pcaC <- prcomp(Cs, scale = TRUE)
eigs <- pcaC$sdev^2
varC <- eigs[1] / sum(eigs)
paste("Variance rate=",varC)
```

The communality value of C3 is lower than 0.25, therefore, this variable should be excluded in the next step.

```{r, warning=FALSE, message=FALSE}
Cs<-U21_filtered[,c("C1","C2","C4","C5","C6")]
KMO(cor(Cs))
fitC <- principal(Cs,nfactors=1)
summary(fitC)
fitC$communality
pcaC <- prcomp(Cs, scale = TRUE)
eigs <- pcaC$sdev^2
varC <- eigs[1] / sum(eigs)
paste("Variance rate=",varC)
```

Both the KMO and the variance rate is greater than 0.5.

### *O (outputs)* PC

At the last step, we specify the principal component *O*, which contains the indicators O1-O9, that describe the research output indicators of each country.

```{r, warning=FALSE, message=FALSE}
Os<-U21_filtered[,c("O1","O2","O3","O4","O5","O6","O7","O8","O9")]
KMO(cor(Os))
fitO <- principal(Os,nfactors=1)
summary(fitO) 
fitO$communality
pcaO <- prcomp(Os, scale = TRUE)
eigs <- pcaC$sdev^2
varO <- eigs[1] / sum(eigs)
paste("Variance rate=",varO)
```

Based on the above results, for both O1 and O9, the communality is lower than 0.25, however, the variance ratio is adequate (0.56). There are two options: use one principal component or two factors.

#### OPT1 - Find one principal component

```{r, warning=FALSE, message=FALSE}
Os<-U21_filtered[,c("O2","O3","O4","O5","O6","O7","O8")]
KMO(cor(Os))
fitO <- principal(Os,nfactors=1)
summary(fitO)
fitO$communality
pcaO <- prcomp(Os, scale = TRUE)
eigs <- pcaC$sdev^2
varO <- eigs[1] / sum(eigs)
paste("Variance rate=",varO)
```

The value of KMO is greater than 0.5, for all indicators the communality is greater than 0.25, the variance rate is 0.56, i.e. at least 50% of the original information remains even after model reduction. However, two important indicators (O1 and O9) were dropped from the model. The O1 indicator is the number of articles published by higher education institutions, while the O9 is the proportion of unemployed graduates. These indicators are also worthwhile and would need to be preserved after model reduction.

#### OPT2 - Find two (or more) factors

Second opportunity to save more indicators with using two factors.

```{r, warning=FALSE, message=FALSE}
Os<-U21_filtered[,c("O1","O2","O3","O4","O5","O6","O7","O8","O9")]
KMO(cor(Os))
fitO2 <- principal(Os,nfactors=2)
summary(fitO2)
fitO2$communality
fitO2$loadings
```

In the case of two (principal) factors, akk communality values are greater than 0.25; however, the correlation between factor 1 and O5, and the factor 2 and O5 is very similar (|cor(RC1,O5)|<2|cor(RC2,O5)| or |cor(RC1,O5)|<|cor(RC2,O5)|+0.25). Therefore O5 cannot be matched to any factor. In this case we say that O5 is a common indicator, which has to be excluded from the model.

```{r, warning=FALSE, message=FALSE}
Os<-U21_filtered[,c("O1","O2","O3","O4","O6","O7","O8","O9")]
KMO(cor(Os))
fitO2 <- principal(Os,nfactors=2)
summary(fitO2) 
fitO2$communality
fitO2$loadings
```

Both the KMO (0.78) and the variance rate (0.689) is greater than 0.5. All communality values are greater than 0.25. And there is no common indicator.

The first factor contains the folowing indicators:

  + $O_2$: Total number of articles that are produced by higher education institutions per capita
  + $O_3$: Average impact of articles, as measured by citations in 2014 of articles that were published in previous years using the Karolinska Institute normalized impact factor
  + $O_4$: Depth of world-class universities in a country. This is calculated as an average of the institutions' score of a country that is listed in the top 500 of the Shanghai ranking, divided by the country's population
  + $O_6$: Enrollment in tertiary education as a % of the eligible population, which is defined as the 5-year age group after secondary education
  + $O_7$: % of the population aged 25-64 with a tertiary qualification
  + $O_8$: Number of researchers (full-time equivalent) in the nation per population
  
The second factor contains the following indicators:

  + $O_1$: Total number of journal articles  that are produced by higher education institutions
  + $O_9$: Unemployment rates among tertiary-educated aged 25-64 years compared with unemployment rates for those with only upper-secondary or post-secondary non-tertiary education


Thus, compared to the specification of the principal component $O$, both the $O_1$ and $O_9$ indicators in the second factor were saved. 

Since we specified principal components all for $R$, $C$, $E$, $O$ indicators, the components can even be correlated. Therefore, in the next step, the correlation between them have to be checked.


### Correlation analysis for PCs

**Figure 13** shows the results of the correlation analysis. The **Figure 13** shows that there are strong correlations among principal components $R$, $O$, and $C$.

```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 13.** Correlation between principal components"}
factors<-cbind(fitR$scores,fitE$scores,fitC$scores,fitO$scores)
colnames(factors) <- c("R","E","C","O")
cor(factors)
corrgram(factors,order = TRUE,main="Correlation between principal components",lower.panel=panel.conf)
```

Correlation analysis is also performed among the PCs $R$, $E$, and $C$ and two specified factors, such as $O_{fact_{1}}$, and $O_{fact_{2}}$. The result of the correlation analysis is shown in **Figure 14**.

```{r, warning=FALSE, message=FALSE,fig.cap="**Figure 14.** Correlations between PCs and factors"}
factors<-cbind(fitR$scores,fitE$scores,fitC$scores,fitO2$scores)
colnames(factors) <- c("R","E","C","O_fact1","O_fact2")
cor(factors)
corrgram(factors,order = TRUE,main="Correlations between PCs and factors",lower.panel=panel.conf)
```

Based on **Figure 14** latent variables R, O_fact1, and C are strongly correlated each other. Correlation analysis specifies three latent variable groups: (1) $R$, $C$, $O_{fact_{1}}$, (2) $O_{fact_{2}}$, (3) $E$.

## Factor analysis

After principal component analysis, factor analysis (FA) is performed on the indicators. As a first step, we examine the applicability of factor analysis, which is aided by scree-plots and calculating the KMO value.

```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 15.** scree plot"}
scree(U21_filtered)
KMO(cor(U21_filtered))
```

**Figure 15** suggests us to specify 3 independent groups of variables (see the same results in **Figure 12**.)


```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 16.** Biplot of 3 factors"}
U21s<-U21_filtered[,c("R1","R2","R3","R4","R5","E1","E2", "E3", "E4","C1","C2","C3","C4","C5","C6","O1","O2","O3","O4","O5","O6","O7","O8","O9")]
fit <- principal(U21s,nfactors=3,rotate = "varimax")
summary(fit) # print variance accounted for 
fit$communality
fit$loadings
biplot(fit)
```

**Figure 16** shows the properties of the three factors. As can be seen from the above results, there are also low communality (e.g. $C_3$) and common indicators (e.g. $R_2$). In order to get a comprehensive picture of the relative position of the factors and the indicators belonging to each factor, it is expedient to plot the initial state of the factors ($RC_1$, $RC_2$, $RC_3$) in 3D. This 3D figure is shown in **Figure 17**. By rotating the figure, you can see how each indicator is grouped, to which factors it belongs. Consistent with the above results, it can be seen here that the first factor includes, for example, $R_3$, $R_4$ and $R_5$.


```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 17.** Factors and indicators in 3D"}
save <- options(rgl.useNULL = TRUE)
open3d()
options(save)
plot3d(fit$loadings)
text3d(fit$loadings,texts = as.vector(rownames(fit$loadings)))
coords <- NULL
for (i in 1:nrow(fit$loadings)) {
  coords <- rbind(coords, rbind(c(0,0,0),fit$loadings[i,]))
}
lines3d(coords, col="red", lwd=4)
rglwidget()
```


Excluding common and low communality value indicators the result is shown in the **Figure 18.**.

```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 18.** shows the final factor structure"}
columns<-c("R3","R4","R5","E1","E2","C1","C5","C6","O1","O4","O8","O9")
U21s<-U21_filtered[,columns]
KMO(cor(U21s))
fit <- principal(U21s,nfactors=3,rotate = "varimax")
summary(fit) 
fit$communality
fit$loadings

biplot(fit)
```

**Figure 19** shows the final factor structure in 3D.

```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 19** final factor structure in 3D"}
plot3d(fit$loadings)
text3d(fit$loadings,texts = as.vector(rownames(fit$loadings)))
coords <- NULL
for (i in 1:nrow(fit$loadings)) {
  coords <- rbind(coords, rbind(c(0,0,0),fit$loadings[i,]))
}
lines3d(coords, col="red", lwd=4)
rglwidget()
```

The 3D figure above illustrates the three distinct groups of indicators.

**Figure 20** shows the structure equation model, which graphically shows which indicators belong to which factor. In the figure, the thickness and strength of the arrows illustrate the extent to which each indicator contributes to a given factor.

The list of indicators, which belong to the factor 1:

  + $R_3$: Annual exp. per student (full-time equivalent) by tertiary education institutions in USD, PPP
  + $R_4$: Exp. in tertiary education institutions for R&D as a % of GDP
  + $R_5$: Exp. in tertiary education institutions for R&D per head of population at USD, PPP
  + $C_1$: % of international students in tertiary education
  + $C_5$: Responses to question "Knowledge transfer is highly developed between companies and universities", which was asked of business executives in the annual survey by IMD World Development Centre, Switzerland
  + $C_6$: % of university research publications that are co-authored with industry researchers
  + $O_4$: Depth of world-class universities in a country. This is calculated as an average of the institutions' score of a country that is listed in the top 500 of the Shanghai ranking, divided by the country's population
  + $O_8$: Number of researchers (full-time equivalent) in the nation per population
  
The list of indicators, which belong to the factor 2:

  + $E_1$: % of female students in tertiary education
  + $E_2$: % of female academic staff in tertiary institutions
  
The list of indicators, which belong to the factor 3:

  + $O_1$: Total number of journal articles  that are produced by higher education institutions
  + $O_9$: Unemployment rates among tertiary-educated aged 25-64 years compared with unemployment rates for those with only upper-secondary or post-secondary non-tertiary education

The first factor ($RC_1$) contains both expenditure-type, relationship-type and output-type indicators. The second factor ($RC_2$) shows the role of women in higher education, and the third factor ($RC_3$) contains output indicators.

```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 20.** The final structure ecuation model"}
model='RC1=~R3+R4+R5+C1+C5+C6+O4+O8
       RC2=~E1+E2
       RC3=~O1+O9
       RC1~~0*RC2 + 0*RC3
       RC2~~0*RC3'
fitsem <- cfa(model,U21_filtered,std.lv = TRUE)
semPaths(fitsem,residuals=FALSE,"std",edge.color = "black",color="yellow")
```

**The summary of the evaluatiion of dimension reduciton:** If the goal is to specify independent groups of variables, then 3 groups of variables can be defined with both principal component and factor analysis. Of these, two "pure" groups of variables contain only 2-2 indicators, while the third mixed group forms a group of R, C, O indicators. A serious problem is that only 12 of the 24 indicators are preserved.

# 4. Data reduction

## Specifiying 3 clusters on 3 factors

Based on the preliminary clustering study, 3-4 clusters were probable for countries (rows). Thus, as a first step, we specify three clusters.

```{r, warning=FALSE, message=FALSE}
set.seed(2)
km3 <- kmeans(fit$scores,3)$cluster
cluster.stats(dist(fit$scores),km3)
```

**Figure 21*) shows the results of clustering based on 3 factors

18 countries belong to the first,20 countries belong to the second, and 12 countries belong to the third cluster.
 
```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 21.** the results of three clusters"}
  plot3d(fit$loadings)
text3d(fit$loadings,texts = rownames(fit$loadings))
coords <- NULL
for (i in 1:nrow(fit$loadings)) {
  coords <- rbind(coords, rbind(c(0,0,0),fit$loadings[i,]))
}
lines3d(coords, col="red", lwd=4)

text3d(fit$scores,texts = rownames(U21_filtered),col=km3)
rglwidget()
```

**Figure 22.** shows which countries belong to which clusters. These countries are separated align with the three factors.

```{r, warning=FALSE, message=FALSE,fig.cap="**Figure 22.** The specified three clusters on the map"}
data(wrld_simpl)
countries_km1 = wrld_simpl@data$NAME %in% c("Australia","Austria","Belgium","Canada","Denmark","Finland","France","Hon Kong","Ireland","Israel","Netherlands","New Zealand","Norway","Portugal","Singapore","Sweden","Switzerland","United Kingdom")
countries_km2 = wrld_simpl@data$NAME %in% c("Argentina","Brazil","Bulgaria","China", "Croatia", "Czech Republic", "Germany", "Hungary","Malaysia","Poland","Romania","Russia","Serbia","Slovakia","Slovenia","South Africa","Spain","Thailand","Ukraine","United States")
countries_km3 = wrld_simpl@data$NAME %in% c("Chile", "Greece", "India", "Indonesia", "Iran", "Italy", "Japan", "Mexico", "Saudi Arabia", "Taiwan", "Turkey", "Korea, Republic of")
countries_km2 = countries_km2 *2
countries_km3 = countries_km3 *3
countries_map = countries_km1+countries_km2+countries_km3+1
plot(wrld_simpl, col = c(gray(.80), "red", "blue", "green")[countries_map])

```

## Specifiying 4 clusters on 3 factors

Next step is to specify four clusters on 3 factors.

```{r, warning=FALSE, message=FALSE}
set.seed(2)
km4 <- kmeans(fit$scores,4)$cluster
cluster.stats(dist(fit$scores),km4)
```
**Figure 23** shows, that the specified four clusters are well separated. There are 22 countries in the first cluster, 4 in the second, 17 in the third and 7 in the fourth.

```{r, warning=FALSE, message=FALSE, fig.cap="**Figure 23.** specified four clusters on three factors"}
  plot3d(fit$loadings)
text3d(fit$loadings,texts = rownames(fit$loadings))
coords <- NULL
for (i in 1:nrow(fit$loadings)) {
  coords <- rbind(coords, rbind(c(0,0,0),fit$loadings[i,]))
}
lines3d(coords, col="red", lwd=4)

text3d(fit$scores,texts = rownames(U21_filtered),col=km4)
rglwidget()
```

**Figure 24** shows the map of clustered countries.

```{r, warning=FALSE, message=FALSE,fig.cap="**Figure 24.** The specified four clusters on the map"}
library(maptools)
data(wrld_simpl)
countries_km1 = wrld_simpl@data$NAME %in% c("Argentina","Brazil","Bulgaria","Chile","Croatia","Czech Republic","Greece", "Italy","Malaysia","Mexico","Poland","Portugal","Romania","Russia","Saudi Arabia","Serbia","Slovakia","Slovenia","South Africa","Spain","Thailand","Ukraine")
countries_km2 = wrld_simpl@data$NAME %in% c("China","Germany","Hungary","United States")
countries_km3 = wrld_simpl@data$NAME %in% c("Australia","Austria","Belgium","Canada","Denmark","Finland","France","Hong Kong", "Ireland","Israel","Netherlands","New Zealand","Norway","Singapore","Sweden","Switzerland","United Kingdom")
countries_km4 = wrld_simpl@data$NAME %in% c("India","Indonesia","Iran","Japan","Taiwan","Turkey","Korea, Republic of")
countries_km2 = countries_km2 *2
countries_km3 = countries_km3 *3
countries_km4 = countries_km4 *4
countries_map = countries_km1+countries_km2+countries_km3+countries_km4+1
plot(wrld_simpl, col = c(gray(.80), "red", "blue", "green", "orange")[countries_map])
```



# 5. Specifying partial rankings within specified clusters.

When specifying the partial rankings (within clusters of countries), we should either use the calculated factor weights or the original weights used by U21.

Firt, we calculate the original weights used by U21.

```{r,message=FALSE,warning=FALSE}
weights=matrix(0,1,24)
colnames(weights)<-c("R1","R2","R3","R4","R5","E1","E2","E3","E4","C1","C2","C3","C4","C5","C6","O1","O2","O3","O4","O5","O6","O7","O8","O9")
weights[1,] <- c(0.05,0.05,0.05,0.025,0.025,0.02,0.02,0.02,0.14,0.04,0.04,0.02,0.02,0.04,0.04,0.13,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03)
rownames(weights)<-c("weights")
```

In the next step, the partial rankings are calculated, which (see **Table 1**). In the last two columns of the **Table 1**, variables *km3* and *km4* indicate cluster memberships for clusters 3 and 4, respectively. *U21_overall_Rank* contains the rankings of the original U21 ranking, *FA_Rank* the ranking based on factor scores, and *U21_Rank_selected* is the rankings that only contain indicators, which corresponds to the factors.

```{r,warning=FALSE, message=FALSE}
base<-as.data.frame(1:nrow(U21_rank))
colnames(base)<-paste("ID")
#base$countries<-rownames(U21_rank)
base$U21_overall_Rank<-U21_rank
base$FA_weights <- as.matrix(pnorm(fit$scores[,1])*pnorm(fit$scores[,2])*pnorm(fit$scores[,3]))
base$FA_Rank<-rank(-base$FA_weights)
base$Score_selected <- rowSums(U21_filtered[,columns]*weights[,columns])*100/max(rowSums(U21_filtered[,columns]*weights[,columns]))
base$U21_Rank_selected<-rank(-base$Score_selected)
base$km3 <-km3
base$km4 <-km4
rownames(base)<-rownames(U21_rank)
kable(base, caption="**Table 1** Partial rankings")
```

In the next step, the rank correlations between the rankings are calculated. 

## Rank correlations, in the case of 3 clusters

**Figures 25-28** show rank correlations between rankings. **Figure 25** shows rank correlations for all countries, while **Figures 26-28** show rank correlations separately for the clusters 1-3.

If we examine all countries (see **Figure 25**), there is a close positive relationship between the rankings of the original U21 ranking (*U21_overall_Rank*) and the ranking we calculated, which contains only the indicators corresponding to each factor (*U21_Rank_selected*). If we examine only the rankings of the countries belonging to clusters 1 and 2, a similar result is obtained (see **Figures 26-27**).

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 25.** Partial rankings fo all countries"}
par(mfrow=c(2,2))
RANKS<-c("U21_overall_Rank","FA_Rank","U21_Rank_selected")
corrgram(base[,RANKS],cor.method = "spearman",order = TRUE,main="For all countries",lower.panel=panel.cor)

```

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 26.** Partial correlactions between calculated and the original rankings on cluster 1"}
corrgram(base[km3==1,RANKS],cor.method = "spearman",order = TRUE,main=paste("Partial rank correlations of",sum(km3==1)," countries from cluster 1"),lower.panel=panel.cor)
```

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 27.** Partial correlactions between calculated and the original rankings on cluster 2"}
corrgram(base[km3==2,RANKS],cor.method = "spearman",order = TRUE,main=paste("Partial rank correlations of",sum(km3==2)," countries from cluster 2"),lower.panel=panel.cor)
```

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 28.** Partial correlactions between calculated and the original rankings on cluster 3"}
corrgram(base[km3==3,RANKS],cor.method = "spearman",order = TRUE,main=paste("Partial rank correlations of",sum(km3==3)," countries from cluster 3"),lower.panel=panel.cor)
```

## Rank correlations, in the case of 4 clusters

Rank correlations are also calculated when four clusters are created. For 1, 2, and 3 clusters, there are moderately strong relationship between the original U21 rankings and calculated rankings.

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 29.** Partial correlactions between calculated and the original rankings on cluster 1"}
par(mfrow=c(2,2))
RANKS<-c("U21_overall_Rank","FA_Rank","U21_Rank_selected")
corrgram(base[km4==1,RANKS],cor.method = "spearman",order = TRUE,main=paste("Partial rank correlations of",sum(km4==1)," countries from cluster 1"),lower.panel=panel.cor)
```

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 30.** Partial correlactions between calculated and the original rankings on cluster 2"}
corrgram(base[km4==2,RANKS],cor.method = "spearman",order = TRUE,main=paste("Partial rank correlations of",sum(km4==2)," countries from cluster 2"),lower.panel=panel.cor)
```

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 31.** Partial correlactions between calculated and the original rankings on cluster 3"}
corrgram(base[km4==3,RANKS],cor.method = "spearman",order = TRUE,main=paste("Partial rank correlations of",sum(km4==3)," countries from cluster 3"),lower.panel=panel.cor)
```

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 32.** Partial correlactions between calculated and the original rankings on cluster 4"}
corrgram(base[km4==4,RANKS],cor.method = "spearman",order = TRUE,main=paste("Partial rank correlations of",sum(km4==4)," countries from cluster 4"),lower.panel=panel.cor)
```

# 6. Summary and conclusions

The results of the two-step clustering are summarized below.

```{r,warning=FALSE, message=FALSE}
MTX3<-as.data.frame(matrix(0,nrow(U21_filtered),ncol(U21_filtered)))
MTX4<-as.data.frame(matrix(0,nrow(U21_filtered),ncol(U21_filtered)))
rownames(MTX3)<-rownames(U21_filtered)
colnames(MTX3)<-colnames(U21_filtered)
rownames(MTX4)<-rownames(U21_filtered)
colnames(MTX4)<-colnames(U21_filtered)
MTX3[,columns]<-km3
MTX4[,columns]<-km4
cols<-c('0'="#FFFFFF",'1'="#99FF66",'2'="#66FF33",'3'="#33CC00",'4'="#009900")
SER3<-c(seriate(dist(as.matrix(MTX3))),seriate(dist(t(as.matrix(MTX3)))))
SER4<-c(seriate(dist(as.matrix(MTX4))),seriate(dist(t(as.matrix(MTX4)))))
ORDERED3<-MTX3[get_order(SER3,dim=1),get_order(SER3,dim=2)]
ORDORIG3<-U21_filtered[get_order(SER3,dim=1),get_order(SER3,dim=2)]
ORDERED4<-MTX4[get_order(SER4,dim=1),get_order(SER4,dim=2)]
ORDORIG4<-U21_filtered[get_order(SER4,dim=1),get_order(SER4,dim=2)]
rownames(ORDERED3)<-rownames(MTX3[get_order(SER3,dim=1),get_order(SER3,dim=2)])
rownames(ORDERED4)<-rownames(MTX4[get_order(SER4,dim=1),get_order(SER4,dim=2)])
rownames(ORDORIG3)<-rownames(MTX3[get_order(SER3,dim=1),get_order(SER3,dim=2)])
rownames(ORDORIG4)<-rownames(MTX4[get_order(SER4,dim=1),get_order(SER4,dim=2)])
par(mfrow=c(2,2))
```

**Figure 33** shows the heat map of the original U21 data. The coloring of the cells depends on the values for each indicator. The darker red cells indicate smaller values, and the darker blue cells indicate higher values. In the **Figure 33**, the countries are groupped by the set of 3 clusters. Cluster 1 contains countries from Belgium to Canada, cluster 2 contains countries from Hungary to Brazil, cluster 3 containes countries from Greece to Indonesia.

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 33.** Heat map of U21 indicators grouped by the set of tree clusters"}
p <- plot_ly(x=colnames(ORDORIG3), y=rownames(ORDORIG3), z = as.matrix(ORDORIG3), type = "heatmap",colorscale="RdBu",reversescale=TRUE)
p
```

**Figure 34** shows the countries belonging to the three clusters, as well as the indicators that belong to the three factors. Indicators that do not belong to any factor are displayed with a dark red cell. It can be seen that half (12) of the indicators were dropped during the analysis. Countries in the first cluster appear in peach color, those in the second in light blue, and those in the third in dark blue.

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 34.** selected countries and indicators"}
p <- plot_ly(x=colnames(ORDERED3), y=rownames(ORDERED3), z = as.matrix(ORDERED3), type = "heatmap",colorscale="RdBu",reversescale=TRUE)
p
```


```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 35.** Heat map of U21 indicators grouped by the set of four clusters"}
p <- plot_ly(x=colnames(ORDORIG4), y=rownames(ORDORIG4), z = as.matrix(ORDORIG4), type = "heatmap",colorscale="RdBu",reversescale=TRUE)
p
```

**Figure 35** shows the heat map of the original U21 data. The coloring of the cells depends on the values for each indicator. The darker red cells indicate smaller values, and the darker blue cells indicate higher values. In the **Figure 35**, the countries are groupped by the set of 4 clusters. Cluster 1 contains countries from Brasil to Argentina, cluster 2 contains countries from Hungary to Germany, cluster 3 contains countries from Belgium to Austria, cluster 4 containes countries from Iran to Taiwan.

**Figure 36.** shows the four clusters specified and the indicators that belong to the three factors.

```{r,warning=FALSE, message=FALSE, fig.cap="**Figure 36.** selected countries and indicators"}
p <- plot_ly(x=colnames(ORDERED4), y=rownames(ORDERED4), z = as.matrix(ORDERED4), type = "heatmap",colorscale="RdBu",reversescale=TRUE)
p
```

**Summary:** *Applying two-step clustering is questionable, first for the low sample. The method itself does not allow clusters to overlap, however, it may perform well in some factors in one country and poorly in others. As a result of clustering, we also obtained very small clusters, which are difficult to analyze and interpret.*

**The additional files attached to the study (including this html) are free to use for scientific purposes. Please cite the original paper**

# Original paper

[Kosztyán, Z. T., Banász, Z., Csányi, V. V., & Telcs, A. (2019). Rankings or leagues or rankings on leagues?-Ranking in fair reference groups. Tertiary Education and Management, 25(4), 289-310.](https://link.springer.com/content/pdf/10.1007/s11233-019-09028-x.pdf)

# Further references

Van Mechelen, I., J. Hampton, R. S. Michalski, and P. Theuns (1993). Categories and Concepts - Theoretical Views and Inductive Data Analysis. Academic Press, London.