---
title: Applying bi-clustering algorithms for specifying partial rankings based on
  U21 and related data
author: "Zsolt T. Kosztyan, Zsuzsanna Banasz, Vivien Valeria Csanyi, Andras Telcs"
output:
  html_document:
    fig_height: 8
    fig_width: 11
    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 libraries and data
  + Loading necessary libraries
  + Loading necessary data sources
2. Estimation of the number of bi-clusters
  + Seriation of normalized data
  + Estimation of the number of bi-clusters
3. Applying iBBiG on normalized data
  + Calculation of top league bi-cluster
  + Calculation of F-tests of the top league bi-cluster
  + Heat maps of the top league
  + Checking stability of the top league bi-cluster
  + Calculation of partial rankings
  + Comparison of partial rankings to the adequate part of the original ranking
  + Venn diagrams of U21's League A and the League A of the extended database
4. Applying iBBiG on reversed normalized data
  + Calculation of lower league bi-cluster
  + Calculation of F-tests of the lower league bi-cluster
  + Heat maps of the lower league
  + Checking stability of the lower league bi-cluster
  + Calculation of partial rankings
  + Comparison of partial rankings to the adequate part of the original ranking
  + Venn diagrams of U21's League C and the League C of the extended database
5. Applying BicARE
  + Calculation of middle league bi-cluster
  + Calculation of F-tests of the middle league bi-cluster
  + Heat maps of the middle league
  + Checking stability of the middle league bi-cluster
  + Calculation of partial rankings
  + Comparison of partial rankings to the adequate part of the original ranking
  + Venn diagrams of U21's League B and the League B of the extended database
6. Specifying overlaps
  + Venn diagrams of indicators
  + Venn diagrams of countries
7. Summing up the conclusions


#Loading libraries and data

##Loading necessary libraries
```{r,warning=FALSE,message=FALSE}

#Seriation and biclustering
library(seriation) #For seriate a matrix
library(biclust) #General biclustering functions
library(iBBiG) #iBBiG biclustering
library(BicARE) #BicARE biclustering
library(BcDiag) #Diagnostics of biclusters

#Additional statistic packages
library(factoextra) #Additional distance functions
library(VennDiagram) #VennDiagram
library(venneuler) #VennDiagram

#For formatted and interactive tables and figures
library(knitr) #Formatted tables, figures
```

Citations

```{r,warning=FALSE,message=FALSE}
#Seriation and biclustering
print(citation("seriation"),style="text") #For seriate a matrix
print(citation("biclust"),style="text") #General biclustering functions
print(citation("iBBiG"),style="text") #iBBiG biclustering
print(citation("BicARE"),style="text") #BicARE biclustering
print(citation("BcDiag"),style="text") #Diagnostics of biclusters

#Additional statistic packages
print(citation("factoextra"),style="text") #Additional distance functions
print(citation("VennDiagram"),style="text") #VennDiagram
print(citation("venneuler"),style="text") #VennDiagram

#For formatted and interactive tables and figures
print(citation("knitr"),style="text") #Formatted tables, figures
```

##Loading necessary data

In the case of this analysis an extended database is used. The extended database contains data from the U21 ranking (R1-5, E1-4, C1-6, O1-9), the Global Competitiveness Index (GCI) and the indicators of its two pillars related to higher education and R&D and GDP from World Bank. New indicators are the following:

5th pillar of the GCI (Higher education and training):

  + GCI_64: Secondary education enrollment, gross %
  + GCI_65: Tertiary education enrollment, gross %
  + GCI_66: Quantity of education
  + GCI_67: Quality of education system (from 1 to 7, where 7 is the best)
  + GCI_68: Quality of Math and Science education (from 1 to 7, where 7 is the best)
  + GCI_69: Quality of Management schools (from 1 to 7, where 7 is the best)
  + GCI_70: Internet access in schools (from 1 to 7, where 7 is the best)
  + GCI_71: Quality of education
  + GCI_72: Availability of research and training services (from 1 to 7, where 7 is the best)
  + GCI_73: Extent to staff training (from 1 to 7, where 7 is the best)
  + GCI_74: On-the-job training
  + GCI_75: Higher education and training
  
12th pillar of the GCI (Innovation):
  
  + GCI_149: Capacity for innovation (from 1 to 7, where 7 is the best)
  + GCI_150: Quality of scientific research institutions (from 1 to 7, where 7 is the best)
  + GCI_151: Company spending on R&D (from 1 to 7, where 7 is the best)
  + GCI_152: University-industry collaboration in R&D (from 1 to 7, where 7 is the best)
  + GCI_153: Governmental procurement of advanced tech products (from 1 to 7, where 7 is the best)
  + GCI_154: Availability of scientists and engineers (from 1 to 7, where 7 is the best)
  + GCI_155: PCT patents, application/million population
  + GCI_156: Innovation
  
Overall GCI:  
  
  + GCI_158: Global Competitiveness Index
  
World Bank/World Development Indicators/GDP:
  
  + WB_WDI_223: GDP per capita, PPP (constant 2011 international $)

```{r,warning=FALSE,message=FALSE}
load("U21_AND_MORE_BIC.RData")
#ALL: 50 by 24 dataframe of normalized data, with two letters indicator names
#U21: 50 by 24 dataframe of normalized data, with original indicator names
#LESS_NORM: 50 by 46 dataframe of normalized U21 and Global Competitive Index indicators for 50 countries
#LESS_ORIG: 50 by 46 dataframe of original U21 and Global Competitive Index indicators for 50 countries
#orig_mtx: 50 by 46 matrix of original U21 and Global Competitive Index indicators for 50 countries
#orig_mtx_u21: 50 by 24 matrix of original U21 indicators for 50 countries
#mtx: 50 by 46 matrix of normalized U21 and Global Competitive Index indicators for 50 countries
#mtx_u21: 50 by 34 matrix of normalized U21 indicators for 50 countries
#U21_rank: 50 by 1 matrix (vector) of U21'S country ranks
#weights: 1 by 46 matrix (vector) of equally weights (in this case all weights are 1/46)
```

#Estimation of the number of bi-clusters

##Seriation of normalized data

Firstly, to obtain a preliminary picture of the possible bi-clusters, seriation of the normalized data is done.

```{r,warning=FALSE,message=FALSE}
ALL_SER <- c(seriate(dist(mtx,"euclidean"),method="MDS"),seriate(get_dist(t(mtx),"spearman"),method="MDS"))
ALL_ORDERED<-mtx[get_order(ALL_SER,dim=1),get_order(ALL_SER,dim=2)]
rownames(ALL_ORDERED) <- rownames(mtx[get_order(ALL_SER,dim=1),get_order(ALL_SER,dim=2)])
```

##Estimating the number of bi-clusters

After seriation, plotting the heat map of the normalized data (see **Fig. 1**) can help to detect the possible bi-clusters. On **Fig. 1** the blue cells indicate lower values of the indicators, while red cells mean higher values. 

```{r,warning=FALSE,message=FALSE,fig.width = 9, fig.height = 7, fig.cap = "**Fig. 1** Heat map of the normalized and seriated matrix"}

hmap(as.matrix(ALL_ORDERED),col=bluered(100),showdist="both")
```

After seriation, two-two bigger homogeneous blocks can be identified based on **Fig. 1**. Two blocks of the red cells on the top left corner on **Fig. 1** indicates the top leagues, while the two bigger blue blocks, which indicates the remaining (lower) leagues, can be discovered in the bottom of the figure. The dendrogram of two-way clustering also shows that regarding rows and columns two-two main blocks can be specified.


#Applying iBBiG on normalized data

Applying the iBBiG method on the normalized data can specify the top league bi-clusters. First of all, the number of bi-clusters needed to be determined.

```{r,warning=FALSE,message=FALSE}
res <- iBBiG(binaryMatrix=binarize(mtx,threshold = NA),nModules = 2,alpha=0.3,pop_size = 100,mutation = 0.08,stagnation = 50,selection_pressure = 1.2,max_sp = 15,success_ratio = 0.6)
summary(res)
```

The results above show that the iBBiG algorithm found 2 bi-clusters. In the first bi-cluster there are 24 rows (countries) and 36 columns (indicators), in the second bi-cluster there are 25 rows (countries) and 7 columns (indicators). The next step is to calculate F-tests to determine which bi-cluster is significant.

```{r,message=FALSE,warning=FALSE}
Obs.FStat <- NULL
for (i in c(1:res@Number)){
  Obs.FStat[[i]] <- computeObservedFstat(x=mtx,bicResult=res,number=i)
}
kable(Obs.FStat,caption = "**Table 1** The results of row and column effects for two bi-clusters")
```
In **Table 1** the results of the F-tests can be seen. In the case of the first bi-cluster, both the row and column effect are significant. On the other hand, in the case of the second bi-cluster the row effect is not significant at *p = 0.01* level; thus the second bi-cluster is not significant. Only the first bi-cluster is significant which is the top league (League A).

##Calculating top league bi-cluster

The next step is to identify the top league bi-cluster.

```{r,warning=FALSE,message=FALSE}
res <- iBBiG(binaryMatrix=binarize(mtx,threshold = 0.50),nModules = 1,alpha=0.3,pop_size = 100,mutation = 0.08,stagnation = 50,selection_pressure = 1.2,max_sp = 15,success_ratio = 0.8)
summary(res)
```

The results of the iBBiG algorithm can be seen above. The algorithm found one bi-cluster in which there are 24 rows (countries) and 36 columns (indicators). The following step is to calculate the F-tests to determine whether the bi-cluster is significant or not.

##Calculating F-tests of the top league bi-cluster

In **Table 2** the results of the row and column effects (F-tests) can be read. Both the row and column effects are significant, so the top league bi-cluster is significant.

```{r,message=FALSE,warning=FALSE}
Obs.FStat <- NULL
for (i in c(1:res@Number)){
  Obs.FStat[[i]] <- computeObservedFstat(x=mtx,bicResult=res,number=i)
}
kable(Obs.FStat,caption = "**Table 2** The results of row and column effects for the top league")
```

To see the differences between those countries that are included in the top league and those that are not in the top league, one can plot the means, medians, variances, and median absolute deviations (MAD) for rows (countries) (see **Fig. 2**). (The horizontal axis of **Fig. 2** shows the countries, and the vertical axis depicts the mentioned parameters on normalized (0-1) data of indicators.)  On **Fig. 2** the red line indicates those countries that are included in the top league, while the black line shows the remaining countries. The countries in the top league have higher means and medians.

```{r,message=FALSE,warning=FALSE,fig.width = 7, fig.height = 7, fig.cap ="**Fig. 2** Mean, median, variance and MAD for the top league countries"}
exploreBic(dset=mtx,bres=res,mname='biclust',pfor='all',gby='genes',bnum=1)
```

One can do the same for the columns (indicators). **Fig. 3** shows the differences between the means, medians, variances, and median absolute deviations (MAD) for the indicators. (The horizontal axis of **Fig. 3** shows the indicators, and the vertical axis depicts the mentioned parameters on normalized (0-1) data of countries.) The red line indicates those indicators that are in the top league, while the black line shows the remaining indicators.

```{r,message=FALSE,warning=FALSE,fig.width = 7, fig.height = 7, fig.cap="**Fig. 3** Mean, median, variance and MAD for the top league indicators"}
exploreBic(dset=mtx,bres=res,mname='biclust',pfor='all',gby='conditions',bnum=1)
```

##Heat maps of the top league

**Fig. 4** depicts the heat map of the top league. The top league contains 24 countries and 36 indicators. The red cells mean higher values in the case of indicators, while blue cells symbol lower values.

```{r,warning=FALSE,message=FALSE, fig.cap= "**Fig. 4** Heat map of the top league"}
drawHeatmap(x=mtx,bicResult=res,number=1,beamercolor = TRUE,paleta=bluered(100))
```

**Fig. 5** shows the heat map of the results of iBBiG on normalized data. The top league can be seen in the upper left corner of **Fig.5**, which is enlarged on **Fig.4**.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap= "**Fig. 5** Heat map of the output after applying iBBiG on normalized data"}
drawHeatmap(x=mtx,bicResult=res,local=FALSE,number=1,beamercolor = TRUE,paleta=bluered(100))
```


##Checking stability of the top league bi-cluster

The next step is to check the stability of the top league bi-cluster. The results of the bootstrap can be seen below.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9}
Bootstrap <- diagnoseColRow(x=mtx,bicResult=res,number=1,
  nResamplings=100,replace=TRUE)
Bootstrap
```

Since the *bootstrapPvalueRow* and the *bootstrapPvalueCol* are lower than 0.01, the bi-cluster is stable.

##Calculating partial rankings

The next step is to calculate partial rankings based on the indicators and countries in the top league. The calculation of each score of each indicator of the countries are done according to the method used in the original U21 ranking. The original weights of U21's indicators are used, and the total scores for the countries are calculated using the selected indicators in the given bi-cluster. **Table 3** shows the result of partial ranking in the case of the top league.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9}
C<-biclust::bicluster(orig_mtx,res,number=1)
B<-as.data.frame(C[[1]])
selectedweight<-weights[,colnames(B)]
B$Overall_Score<-rowSums(B*selectedweight)*100/max(rowSums(B*selectedweight))
B$Rank<-rank(-B$Overall_Score)
kable(B,caption = "**Table 3** League A by iBBiG",digits = 2)
```

##Comparing partial rankings to the adequate part of the original ranking

The final step is to compare the results of the partial rankings to the corresponding part of the original U21 ranking. The results of the Spearman's rank correlation show a moderate positive relationship between the two rankings.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9}
cor(B$Rank,U21_rank[rownames(B),],method="spearman")
cor.test(B$Rank,U21_rank[rownames(B),],method="spearman")
```


##Venn diagram of U21's *League A* and the *League A* of the extended database

An interesting question is how the *League A* of the U21 and the *League A* of the extended database (U21 + GCI + GDP) are overlap. In order to answer this question, one can make a Venn diagram of the two leagues to detect overlaps. The first step before plotting the Venn diagrams (one for the indicators and one for the countries), is to calculate the top league bi-cluster of the U21 database.

```{r,warning=FALSE,message=FALSE}
res_U21 <- iBBiG(binaryMatrix=binarize(mtx_u21,threshold = 0.50),nModules = 1,alpha=0.05,pop_size = 100,mutation = 0.08,stagnation = 50,selection_pressure = 1.2,max_sp = 15,success_ratio = 0.8)
T<-biclust::bicluster(orig_mtx,res,number = 1)
T_U21<-biclust::bicluster(orig_mtx_u21,res_U21,number = 1)
League_A<-T[[1]]
League_A_U21<-T_U21[[1]]
```

**Fig.6** shows the Venn diagram of the indicators in the two different top leagues. The green circle contains those indicators that are in the U21's top league *(League A)* while in the red circle the extended database's top league's indicators can be found. 15 indicators are the same in the two top leagues.

```{r,warning=FALSE,message=FALSE,fig.width = 5, fig.height = 9, fig.cap="**Fig. 6** Venn diagram of the indicators in the top leagues (original vs. extended database)"}
draw.pairwise.venn(area1=nrow(as.matrix(colnames(League_A))),area2=nrow(as.matrix(colnames(League_A_U21))),cross.area=nrow(as.matrix(intersect(colnames(League_A),colnames(League_A_U21)))),fill=c("red","green"),category=c("Extended","U21"),main="Venn Diagram of indicators on League As")
```

One can draw a Venn diagram for the countries in the two top leagues. On **Fig. 7** the green circle contains those countries that are in the U21's top league while in the red circle the extended database's top league's countries can be found. 21 countries are the same in the two top leagues.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap="**Fig. 7** Venn diagram of the countries in the top leagues (original vs. extended database)"}
draw.pairwise.venn(area1=nrow(as.matrix(rownames(League_A))),area2=nrow(as.matrix(rownames(League_A_U21))),cross.area=nrow(as.matrix(intersect(rownames(League_A),rownames(League_A_U21)))),fill=c("red","green"),category=c("Extended","U21"),main="Venn Diagram of countries on League As")
```

#Applying iBBiG on reversed normalized data

To identify the lower league(s), the reverse data matrix is needed to be calculated firstly. Applying the iBBiG method on the reversed normalized data can specify the lower league bi-clusters.

```{r,warning=FALSE,message=FALSE}
rmtx <- 1-mtx
rres <- iBBiG(binaryMatrix=binarize(rmtx,threshold = NA),nModules = 2,alpha=0.3,pop_size = 100,mutation = 0.08,stagnation = 50,selection_pressure = 1.2,max_sp = 15,success_ratio = 0.6)
summary(rres)
```

The results above show that the iBBiG algorithm found 2 bi-clusters on the reversed data. In the first bi-cluster there are 29 rows (countries) and 39 columns (indicators), in the second bi-cluster there are 22 rows (countries) and 9 columns (indicators). The next step is to calculate F-tests to determine which bi-cluster is significant.


```{r,message=FALSE,warning=FALSE} 
Obs.FStat <- NULL 
for (i in c(1:rres@Number)){
  Obs.FStat[[i]] <- computeObservedFstat(x=rmtx,bicResult=rres,number=i) 
  } 
kable(Obs.FStat,caption = "**Table 4** The results of row and column effects of the reversed data for two bi-clusters")
```

In **Table 4** the results of the F-tests can be seen. In the case of the first bi-cluster, both the row and column effect are significant. On the other hand, in the case of the second bi-cluster the row effect is not significant; thus the second bi-cluster is not significant. Only the first bi-cluster is significant which is the lower league (League C).

##Calculating lower league bi-cluster

The next step is to identify the lower league bi-cluster.

```{r,warning=FALSE,message=FALSE,echo=FALSE} 
  rres <- iBBiG(binaryMatrix=binarize(rmtx,threshold = NA),nModules = 1,alpha=0.3,pop_size = 100,mutation = 0.08,stagnation = 50,selection_pressure = 1.2,max_sp = 15,success_ratio = 0.8)
  score<-rres@Clusterscores
summary(rres) 
```

The results of the iBBiG algorithm on the reversed data can be seen above. The algorithm found one bi-cluster in which there are 29 rows (countries) and 39 columns (indicators). The following step is to calculate the F-tests to determine whether the bi-cluster is significant or not.

##Calculating F-tests of the lower league bi-cluster

In **Table 5** the results of the row and column effects (F-tests) can be read. Both the row and column effects are significant, so the lower league bi-cluster is significant.

```{r,message=FALSE,warning=FALSE} 
Obs.FStat <- NULL 
for (i in c(1:rres@Number)){
  Obs.FStat[[i]] <- computeObservedFstat(x=rmtx,bicResult=rres,number=i) 
  } 
kable(Obs.FStat,caption = "**Table 5** The results of row and column effects of the lower league")
```

To see the differences between those countries that are included in the lower league and those that are not in the lower league, one can plot the means, medians, variances, and median absolute deviations (MAD) for rows (countries) (see **Fig. 8**). On **Fig. 8** the red line indicates those countries that are included in the lower league, while the black line shows the remaining countries. The countries in the lower league have lower means and medians.

```{r,message=FALSE,warning=FALSE,fig.width = 7, fig.height = 7, fig.cap="**Fig. 8** Mean, median, variance and MAD for lower league countries "}
exploreBic(dset=mtx,bres=rres,mname='biclust',pfor='all',gby='genes',bnum=1)
```
One can do the same for the columns (indicators). **Fig. 9** shows the differences between the means, medians, variances, and median absolute deviations (MAD) for the indicators. The red line indicates those indicators that are in the lower league, while the black line shows the remaining indicators.

```{r,message=FALSE,warning=FALSE,fig.width = 7, fig.height = 7, fig.cap="**Fig. 9** Mean, median, variance and MAD for lower league indicators"}
exploreBic(dset=mtx,bres=rres,mname='biclust',pfor='all',gby='conditions',bnum=1)
```


##Heat maps of the lower league

**Fig. 10** depicts the heat map of the lower league. The lower league contains 29 countries and 39 indicators. The red cells mean higher values in the case of indicators, while blue cells symbol lower values.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 7, fig.cap="**Fig. 10** Heat map of the lower league"}
drawHeatmap(x=mtx,bicResult=rres,number=1,paleta=bluered(100),beamercolor = TRUE)
``` 

**Fig. 11** shows the heat map of the results of iBBiG on reversed normalized data. The lower league can be seen in the upper left corner of **Fig. 11**, which is enlarged on **Fig.10**.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap="**Fig. 11** Heat map of the output of after applying iBBiG on reversed normalized data"}
drawHeatmap(x=mtx,bicResult=rres,local=FALSE,number=1,paleta=bluered(100),beamercolor = TRUE) 
```

##Checking stability of the lower league bi-cluster

The next step is to check the stability of the lower league bi-cluster. The results of the bootstrap can be seen below.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9}
Bootstrap <- diagnoseColRow(x=mtx,bicResult=rres,number=1,
  nResamplings=100,replace=TRUE)
Bootstrap
```

Since the *bootstrapPvalueRow* and the *bootstrapPvalueCol* are lower than 0.01, the bi-cluster is stable.

##Calculating partial rankings

The next step is to calculate partial rankings based on the indicators and countries in the lower league. The calculation of each score of each indicator of the countries are done according to the method used in the original U21 ranking. The original weights of U21's indicators are used, and the total scores for the countries are calculated using the selected indicators in the given bi-cluster. **Table 6** shows the result of partial ranking in the case of the lower league.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9}
C<-biclust::bicluster(orig_mtx,rres,number=1)
B<-as.data.frame(C[[1]])
selectedweight<-weights[,colnames(B)]
B$Overall_Score<-rowSums(B*selectedweight)*100/max(rowSums(B*selectedweight))
B$Rank<-rank(-B$Overall_Score)
kable(B,caption = "**Table 6** League C by reverse iBBiG",digits = 2)
```

##Comparing partial rankings to the adequate part of the original ranking

The final step is to compare the results of the partial rankings to the corresponding part of the original U21 ranking. The results of the Spearman's rank correlation show a strong positive relationship between the two ranking.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9}
cor(B$Rank,U21_rank[rownames(B),],method="spearman")
cor.test(B$Rank,U21_rank[rownames(B),],method="spearman")
```


##Venn diagram of U21's *League C* and the *League C* of the extended database

An interesting question is how the *League C* of the U21 and the *League C* of the extended database (U21 + GCI + GDP) are overlap. In order to answer this question, one can make a Venn diagram of the two leagues to detect overlaps. The first step before plotting the Venn diagrams (one for indicators and one for countries), is to calculate the lower league bi-cluster of the U21 database.

```{r,warning=FALSE,message=FALSE}
rres_U21 <- iBBiG(binaryMatrix=binarize(1-mtx_u21,threshold = 0.50),nModules = 1,alpha=0.08,pop_size = 100,mutation = 0.08,stagnation = 50,selection_pressure = 1.2,max_sp = 15,success_ratio = 0.8)
T<-biclust::bicluster(orig_mtx,rres,number = 1)
T_U21<-biclust::bicluster(orig_mtx_u21,rres_U21,number = 1)
League_C<-T[[1]]
League_C_U21<-T_U21[[1]]
```

**Fig.12** shows the Venn diagram of the indicators in the two different lower leagues. The green circle contains those indicators that are in the U21's lower league *(League C)* while in the red circle the extended database's lower league's indicators can be found. 19 indicators are the same in the two lower leagues.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap="**Fig. 12** Venn diagram of the indicators in the lower leagues"}
draw.pairwise.venn(area1=nrow(as.matrix(colnames(League_C))),area2=nrow(as.matrix(colnames(League_C_U21))),cross.area=nrow(as.matrix(intersect(colnames(League_C),colnames(League_C_U21)))),fill=c("red","green"),category=c("Extended","U21"),main="Venn Diagram of indicators on League Cs")
```

One can draw a Venn diagram for the countries in the two lower leagues. On **Fig. 13** the green circle contains those countries that are in the U21's lower league while in the brown circle the extended database's lower league's countries can be found. An interesting finding is that all of the countries in the extended database (29) are part of the U21 countries. 

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap="**Fig. 13** Venn diagram of the countries in the lower leagues"}
draw.pairwise.venn(area1=nrow(as.matrix(rownames(League_C))),area2=nrow(as.matrix(rownames(League_C_U21))),cross.area=nrow(as.matrix(intersect(rownames(League_C),rownames(League_C_U21)))),fill=c("red","green"),category=c("Extended","U21"),main="Venn Diagram of countries on League Cs")
```

#Bi-clustering by correlations - Applying BiCARE algorithm

With the BicARE method the middle league(s) can be specified. In order to use the BiCARE method, the first step is to convert the BicARE object to a Biclust object.

```{r,message=FALSE,warning=FALSE}
bicare2biclust <- function(x){
	if(class(x)=="Biclust"){
		return(x)
	} else if(class(x)=="biclustering"){
		Parameters <- list(numberofbicluster=x$param[1,2],residuthreshold=x$param[2,2],genesinitialprobability=x$param[3,2],samplesinitialprobability=x$param[4,2],numberofiterations=x$param[5,2],date=x$param[6,2])
		RowxNumber <- t(x$bicRow==1)
		NumberxCol <- x$bicCol==1
		Number <- as.numeric(dim(RowxNumber)[2])
		info <- list()
		return(new("Biclust",Parameters=Parameters,RowxNumber=RowxNumber,NumberxCol=NumberxCol,Number=Number,info=info))
	}
}
```

In the next step the BicARE algorithm determines the number of the bi-clusters.

```{r,warning=F,message=FALSE}
BICARE_res <- FLOC(Data=mtx,k=2,pGene=1.0,pSample=1.0,r=NULL,N=8,M=6,t=500,blocGene=NULL,blocSample=NULL)
bres <- bicare2biclust(BICARE_res)
summary(bres)
bres@Parameters$residuthreshold
```

The results above show that the BicARE algorithm found 2 bi-clusters. In the first bi-cluster there are 35 rows (countries) and 6 columns (indicators), in the second bi-cluster there are 27 rows (countries) and 6 columns (indicators). The next step is to calculate F-tests to determine which bi-cluster is significant.

```{r,message=FALSE,message=FALSE}
RObs.FStat <- NULL
for (i in c(1:bres@Number)){
  RObs.FStat[[i]] <- computeObservedFstat(x=mtx,bicResult=bres,number=i)
}
kable(RObs.FStat,caption = "**Table 7** The results of row and column effects of the middle league for two bi-clusters")
```
In **Table 7** the results of the F-tests can be seen. Although both bi-clusters are significant, in the following only the first bi-cluster is considered because based on the U21 data only one bi-cluster was found and can be compared with the extended database's results. 
The first bi-cluster will be considered as the middle league (League B).

##Calculating middle league bi-cluster

The next step is to identify the middle league bi-cluster.

```{r,warning=F,message=FALSE}
BICARE_res <- FLOC(Data=mtx,k=1,pGene=1.0,pSample=1.0,r=NULL,N=8,M=6,t=500,blocGene=NULL,blocSample=NULL)
bres <- bicare2biclust(BICARE_res)
summary(bres)
bres@Parameters$residuthreshold
```

The results of the BicARE algorithm can be seen above. The algorithm found one bi-cluster in which there are 27 rows (countries) and 6 columns (indicators). The following step is to calculate the F-tests to determine whether the bi-cluster is significant or not.

##Calculating F-tests of the middle league bi-cluster

In **Table 8** the results of the row and column effects (F-tests) can be read. Both the row and column effects are significant, so the middle league bi-cluster is significant.

```{r,message=FALSE,message=FALSE}
RObs.FStat <- NULL
for (i in c(1:bres@Number)){
  RObs.FStat[[i]] <- computeObservedFstat(x=mtx,bicResult=bres,number=i)
}
kable(RObs.FStat,caption = "**Table 8** The results of row and column effects of the middle league")
```

To see the differences between those countries that are included in the middle league and those that are not in the middle league, one can plot the means, medians, variances, and median absolute deviations (MAD) for rows (countries) (see **Fig. 14**). On **Fig. 14** the red line indicates those countries that are included in the middle league, while the black line shows the remaining countries. **Fig. 14** shows that in the case of the countries in the middle league the variances differ the most.


```{r,message=FALSE,warning=FALSE,fig.width = 7, fig.height = 7, fig.cap="**Fig. 14** Mean, median, variance and MAD for middle league countries "}
exploreBic(dset=mtx,bres=bres,mname='biclust',pfor='all',gby='genes',bnum=1)
```

One can do the same for the columns (indicators). **Fig. 15** shows the differences between the means, medians, variances, and median absolute deviations (MAD) for the indicators. The red line indicates those indicators that are in the middle league, while the black line shows the remaining indicators. 

```{r,message=FALSE,warning=FALSE,fig.width = 7, fig.height = 7, fig.cap="**Fig. 15** Mean, median, variance and MAD for middle league indicators"}
exploreBic(dset=mtx,bres=bres,mname='biclust',pfor='all',gby='conditions',bnum=1)
```

##Heat maps of the middle league

**Fig. 16** depicts the heat map of the middle league. The middle league contains 27 countries and 6 indicators. The red cells mean higher values in the case of indicators, while blue cells symbol lower values.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 7, fig.cap="**Fig. 16** Heat map of the middle league"}
drawHeatmap(x=mtx,bicResult=bres,number=1,paleta=bluered(100),beamercolor = TRUE) 
``` 

**Fig. 17** shows the heat map of the results of BicARE on normalized data.The middle league can be seen in the upper left corner of **Fig. 17**, which is enlarged on **Fig.16**.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap="**Fig. 17** Heat map of the output after applying BicARE on normalized data"}
drawHeatmap(x=mtx,bicResult=bres,local=FALSE,number=1,paleta=bluered(100),beamercolor = TRUE) 
```


##Checking stability of the middle league bi-cluster

The next step is to check the stability of the middle league bi-cluster. The results of the bootstrap can be seen below.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9}
Bootstrap <- diagnoseColRow(x=mtx,bicResult=bres,number=1,
  nResamplings=100,replace=TRUE)
Bootstrap
```

Since the *bootstrapPvalueRow* and the *bootstrapPvalueCol* are lower than 0.01, the bi-cluster is stable.


##Calculating partial rankings

The next step is to calculate partial rankings based on the indicators and countries in the middle league. The calculation of each score of each indicator of the countries are done according to the method used in the original U21 ranking. The original weights of U21?s indicators are used, and the total scores for the countries are calculated using the selected indicators in the given bi-cluster. **Table 9** shows the result of partial ranking in the case of the middle league.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9}
C<-biclust::bicluster(orig_mtx,bres,number=1)
B<-as.data.frame(C[[1]])
selectedweight<-weights[,colnames(B)]
B$Overall_Score<-rowSums(B*selectedweight)*100/max(rowSums(B*selectedweight))
B$Rank<-rank(-B$Overall_Score)
kable(B,caption = "**Table 9** League B by BicARE",digits = 2)
```

##Comparing partial rankings to the adequate part of the original ranking

The final step is to compare the results of the partial rankings to the corresponding part of the original U21 ranking. The results of the Spearman?s rank correlation show a moderate positive relationship between the two ranking.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9}
cor(B$Rank,U21_rank[rownames(B),],method="spearman")
cor.test(B$Rank,U21_rank[rownames(B),],method="spearman")
```


##Venn diagram of U21's *League B* and the *League B* of the extended database

An interesting question is how the *League B* of the U21 and the *League B* of the extended database (U21 + GCI + GDP) are overlap. In order to answer this question, one can make a Venn diagram of the two leagues to detect overlaps. The first step before plotting the Venn diagrams(one for indicators and one for countries), is to calculate the middle league bi-cluster of the U21 database.

```{r,warning=FALSE,message=FALSE}
BICARE_res_U21 <- FLOC(Data=as.matrix(ALL),k=1,pGene=1,pSample=1,r=1e-16,N=17,M=6,t=10000,blocGene=NULL,blocSample=NULL)
bres_U21 <- bicare2biclust(BICARE_res_U21)

T<-biclust::bicluster(orig_mtx,bres,number = 1)
T_U21<-biclust::bicluster(orig_mtx_u21,bres_U21,number = 1)
League_B<-T[[1]]
League_B_U21<-T_U21[[1]]
```

**Fig.18** shows the Venn diagram of the indicators in the two different top leagues. The green circle contains those indicators that are in the U21's middle league *(League B)* while in the red circle the extended database's middle league's indicators can be found. There is no common indicator in the two middle leagues.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap="**Fig. 18** Venn diagram of the indicators in the middle leagues"}
draw.pairwise.venn(area1=nrow(as.matrix(colnames(League_B))),area2=nrow(as.matrix(colnames(League_B_U21))),cross.area=nrow(as.matrix(intersect(colnames(League_B),colnames(League_B_U21)))),fill=c("red","green"),category=c("Extended","U21"),main="Venn Diagram of indicators on League Bs")
```

One can draw a Venn diagram for the countries in the two middle leagues. On **Fig. 19** the green circle contains those countries that are in the U21's middle league while in the red circle the extended database's middle league's countries can be found. 8 countries are the same in the two middle leagues.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap="**Fig. 19** Venn diagram of the countries in the middle leagues"}
draw.pairwise.venn(area1=nrow(as.matrix(rownames(League_B))),area2=nrow(as.matrix(rownames(League_B_U21))),cross.area=nrow(as.matrix(intersect(rownames(League_B),rownames(League_B_U21)))),fill=c("red","green"),category=c("Extended","U21"),main="Venn Diagram of countries on League Bs")
```


#Specifying overlaps

Bi-clusters might (or might not) overlap. Analyzing the overlaps can give a subtler picture of countries and indicators.

##Venn diagram of indicators

**Fig. 20** shows the Venn diagram of the extended database's (U21 + GCI + GDP) indicators. The red circle contains those indicators that are in the upper league (League A), in the blue circle those indicators can be found that are in the lower league (League C). The green circle indicates the middle league's (League B) indicators.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap="**Fig. 20** Venn diagram of the indicators"}
T<-biclust::bicluster(orig_mtx,res,number = 1)
League_A<-T[[1]]
R<-biclust::bicluster(orig_mtx,rres,number = 1)
League_C<-R[[1]]
B<-biclust::bicluster(orig_mtx,bres,number = 1)
League_B<-B[[1]]
draw.triple.venn(area1=nrow(as.matrix(colnames(League_A))),area2=nrow(as.matrix(colnames(League_B))),area3=nrow(as.matrix(colnames(League_C))),n12=nrow(as.matrix(intersect(colnames(League_A),colnames(League_B)))),n13=nrow(as.matrix(intersect(colnames(League_A),colnames(League_C)))),n23=nrow(as.matrix(intersect(colnames(League_B),colnames(League_C)))),n123=nrow(as.matrix(intersect(intersect(colnames(League_A),colnames(League_B)),colnames(League_C)))),fill=c("red","green","blue"),category=c("League A","League B","League C"),main="Venn Diagram of indicators of the extended database")
```

One interesting finding is, that all of the indicators in the middle league are part of the lower league's indicators.

##Venn diagram of countries

**Fig. 21** shows the Venn diagram of the extended database's (U21 + GCI +GDP) countries. The red circle contains those countries that are in the upper league (League A), in the the blue circle those countries can be found that are in the lower league (League C). The green circle indicates the middle league's (League B) countries.

```{r,warning=FALSE,message=FALSE,fig.width = 7, fig.height = 9, fig.cap="**Fig. 21** Venn diagram of the countries"}
draw.triple.venn(area1=nrow(as.matrix(rownames(League_A))),area2=nrow(as.matrix(rownames(League_B))),area3=nrow(as.matrix(rownames(League_C))),n12=nrow(as.matrix(intersect(rownames(League_A),rownames(League_B)))),n13=nrow(as.matrix(intersect(rownames(League_A),rownames(League_C)))),n23=nrow(as.matrix(intersect(rownames(League_B),rownames(League_C)))),n123=nrow(as.matrix(intersect(intersect(rownames(League_A),rownames(League_B)),rownames(League_C)))),fill=c("red","green","blue"),category=c("League A","League B","League C"),main="Venn Diagram of countries of the extended database")
```

If the extended database (U21 + GCI + GDP) is used there are no countries which belong to all of the three leagues at the same time.

#Summing up the conclusions

The above proposed bi-clustering methods can help to identify which countries can be compared based on which indicators in higher education system rankings. With the iBBiG method the upper and lower league(s) *(League(s) A and League(s) C)* can be determined and the BicARE method can specify the middle league(s) *(League(s) B)*. In this case an extended database was used during the bi-clustering. The exentended database contains not just indicators of the U21 ranking, but indicators of the Global Competitive Index (GCI) and an indicator from the World Bank (WB/GDP).
The results show that using this extended database do not change significantly the countries and the indicators in the upper and lower league compared to the results of the leagues of the U21 database. On the other hand, in the case of the middle league there are differences between the results of the extended database and the U21 database.