Parciális rangsorok képzése hagyományos modell-, illetve adatredukciós módszerek alapján

1. Munkamenet (Tartalomjegyzék)

• Szükséges függvénykönyvtárak beolvasása
• Függvénykönyvtárak hivatkozásai
• Szükséges adatok beolvasása

2. Kétlépéses klaszterezhetőség vizsgálata

• Indikátorok (oszlopok) klaszterezhetősége
  • Hopkins mutató számítása
  • Hierarchikus klaszterezés az indikátorok esetében
  • Többdimenziós skálázás alkalmazása az indikátorokra
• Országok (sorok) klaszterezhetősége
  • Hopkins mutató számítása
  • Hierarchikus klaszterezés az indikátorok esetében
  • Többdimenziós skálázás alkalmazása az indikátorokra

3. Modellredukció

• Főkomponens-elemzés
  • R (erőforrások) főkomponens meghatározása
  • E (környezet) főkomponens meghatározása
  • C (kapcsolatok) főkomponenes meghatározása
  • O (kiment) főkomponenes meghatározása
  • Korrelációs vizsgálat a főkomponensekre vonatkozóan
• Faktor-elemzés

4. Adatredukció

• A faktorok 3 klaszterbe rendezése
• A faktorok 4 klaszterbe rendezése

5. Parciális rangsorok meghatározása

• Rangsorok közti korreláció 3 klaszter esetén
• Rangsorok közti korreláció 4 klaszter esetén

6. Összegzés

1. Szükséges függvénykönyvtárak és adatok betöltése

Szükséges függvénykönyvtárak beolvasása

#Korreláció-elemzés/megjelenítés
library(corrgram) #Korrelogram
library(corrplot) #Korreláció-ábra
source("http://www.sthda.com/upload/rquery_cormat.r") #Fejlett 2D-ábrák

#Faktor- és főkomponenselemzés, többdimenziós skálázás
library(psych) #Faktorelemzés, főkomponens-elemzés
library(factoextra) #További távolság/hasonlóság-mértékek
library(MASS) #isoMDS=többdimenziós skálázás

#Klaszterelemzés, klaszterezettség-elemzés
library(fclust) #VAT klaszterek
library(clustertend) #VAT, Hopkins mutatók
library(cluster) #Klaszterek kirajzolása
library(fpc) #Klaszterstatisztikák
library(mclust) #Többlépcsős klaszterelemzés

#Összetett modellek számítása, megjelenítése
library(lavaan) #CFA, SEM modellek számítása
library(semPlot) #Összetett modellek megjelenítése

#További statisztikai függvénykönyvtárak
library(Hmisc) #Fejlett statisztikai módszerek
library(seriation) #Sorbarendezés, szeriáció

#Adattranszformáció, táblázatmegjelenítés
library(reshape2) #Adattranszformáció
library("dplyr", character.only = TRUE) #Adattranszformáció
library(knitr) #Formázott táblázatok megjelenítése

#Térképes megjelenítés
library(maps) #Megjelenítés térképeken
library(mapdata) #Megjelenítés térképeken
library(maptools) #Megjelenítés térképeken

#Interaktív és 3D megjelenítések
library(shiny) #Interaktív megjelenítés
library(htmltools) #Interaktív megjelnítés
library(rgl) #3D megjelenítés
library(rglwidget) #3D interaktív megjelenítés
library(ggplot2) #Interaktív ábrák megjelenítése
library(ggrepel) #Interaktív ábrák megjelenítése
library(ggdendro) #Dendogram megjelenítése
library(d3heatmap) #Interaktív hőtérkép
library(igraph) #Gráf-megjelenítés
library(ggsci) #Simpson-paletta
library(plotly) #Plotly - interaktí diagram

Hivatkozások

#Korreláció-elemzés/megjelenítés
base::print(citation("corrgram"),style="text") #Korrelogram

## Wright K (2021). _corrgram: Plot a Correlogram_. R package version
## 1.14, <URL: https://CRAN.R-project.org/package=corrgram>.

base::print(citation("corrplot"), style ="text") #Korreláció-ábra

## Wei T, Simko V (2021). _R package "corrplot": Visualization of a
## Correlation Matrix_. (Version 0.88), <URL:
## https://github.com/taiyun/corrplot>.

#Faktor- és főkomponenselemzés, többdimenziós skálázás
base::print(citation("psych"),style="text") #Faktorelemzés, főkomponens-elemzés

## Revelle W (2021). _psych: Procedures for Psychological, Psychometric,
## and Personality Research_.  Northwestern University, Evanston,
## Illinois. R package version 2.1.3, <URL:
## https://CRAN.R-project.org/package=psych>.

base::print(citation("factoextra"),style="text") #További távolság/hasonlóság-mértékek

## Kassambara A, Mundt F (2020). _factoextra: Extract and Visualize the
## Results of Multivariate Data Analyses_. R package version 1.0.7, <URL:
## https://CRAN.R-project.org/package=factoextra>.

base::print(citation("MASS"),style="text") #isoMDS=többdimenziós skálázás

## Venables WN, Ripley BD (2002). _Modern Applied Statistics with S_,
## Fourth edition. Springer, New York. ISBN 0-387-95457-0, <URL:
## https://www.stats.ox.ac.uk/pub/MASS4/>.

#Klaszterelemzés, klaszterezettség-elemzés
base::print(citation("fclust"),style="text") #VAT klaszterek

## Ferraro M, Giordani P, Serafini A (2019). "fclust: An R Package for
## Fuzzy Clustering." _The R Journal_, *11*. <URL:
## https://journal.r-project.org/archive/2019/RJ-2019-017/RJ-2019-017.pdf>.

base::print(citation("clustertend"),style="text") #VAT, Hopkins mutatók

## Wright K, YiLan L, RuTong Z (2021). _clustertend: Check the Clustering
## Tendency_. R package version 1.5, <URL:
## https://CRAN.R-project.org/package=clustertend>.

base::print(citation("cluster"),style="text") #Klaszterek kirajzolása

## Maechler M, Rousseeuw P, Struyf A, Hubert M, Hornik K (2021). _cluster:
## Cluster Analysis Basics and Extensions_. R package version 2.1.2 - For
## new features, see the 'Changelog' file (in the package source), <URL:
## https://CRAN.R-project.org/package=cluster>.

base::print(citation("fpc"),style="text") #Klaszterstatisztikák

## Hennig C (2020). _fpc: Flexible Procedures for Clustering_. R package
## version 2.2-9, <URL: https://CRAN.R-project.org/package=fpc>.

base::print(citation("mclust"),style="text") #Többlépcsős klaszterelemzés

## Scrucca L, Fop M, Murphy TB, Raftery AE (2016). "mclust 5: clustering,
## classification and density estimation using Gaussian finite mixture
## models." _The R Journal_, *8*(1), 289-317. <URL:
## https://doi.org/10.32614/RJ-2016-021>.

#Összetett modellek számítása, megjelenítése
base::print(citation("lavaan"),style="text") #CFA, SEM modellek számítása

## Rosseel Y (2012). "lavaan: An R Package for Structural Equation
## Modeling." _Journal of Statistical Software_, *48*(2), 1-36. <URL:
## https://www.jstatsoft.org/v48/i02/>.

base::print(citation("semPlot"),style="text") #Összetett modellek megjelenítése

## Epskamp S (2019). _semPlot: Path Diagrams and Visual Analysis of
## Various SEM Packages' Output_. R package version 1.1.2, <URL:
## https://CRAN.R-project.org/package=semPlot>.

#További statisztikai függvénykönyvtárak
base::print(citation("Hmisc"),style="text") #Fejlett statisztikai módszerek

## Harrell Jr FE, Dupont wcfC, others. m (2021). _Hmisc: Harrell
## Miscellaneous_. R package version 4.5-0, <URL:
## https://CRAN.R-project.org/package=Hmisc>.

base::print(citation("seriation"),style="text") #Sorbarendezés, szeriáció

## Hahsler M, Buchta C, Hornik K (2020). _seriation: Infrastructure for
## Ordering Objects Using Seriation_. R package version 1.2-9, <URL:
## https://CRAN.R-project.org/package=seriation>.
## 
## Hahsler M, Hornik K, Buchta C (2008). "Getting things in order: An
## introduction to the R package seriation." _Journal of Statistical
## Software_, *25*(3), 1-34. ISSN 1548-7660, doi: 10.18637/jss.v025.i03
## (URL: https://doi.org/10.18637/jss.v025.i03), <URL:
## https://www.jstatsoft.org/v25/i03/>.

#Adattranszformáció, táblázatmegjelenítés
base::print(citation("reshape2"),style="text") #Adattranszformáció

## Wickham H (2007). "Reshaping Data with the reshape Package." _Journal
## of Statistical Software_, *21*(12), 1-20. <URL:
## http://www.jstatsoft.org/v21/i12/>.

base::print(citation("dplyr"),style="text") #Adattranszformáció

## Wickham H, François R, Henry L, Müller K (2021). _dplyr: A Grammar of
## Data Manipulation_. R package version 1.0.6, <URL:
## https://CRAN.R-project.org/package=dplyr>.

base::print(citation("knitr"),style="text") #Formázott táblázatok megjelenítése

## Xie Y (2021). _knitr: A General-Purpose Package for Dynamic Report
## Generation in R_. R package version 1.33, <URL:
## https://yihui.org/knitr/>.
## 
## Xie Y (2015). _Dynamic Documents with R and knitr_, 2nd edition.
## Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, <URL:
## https://yihui.org/knitr/>.
## 
## Xie Y (2014). "knitr: A Comprehensive Tool for Reproducible Research in
## R." In Stodden V, Leisch F, Peng RD (eds.), _Implementing Reproducible
## Computational Research_. Chapman and Hall/CRC. ISBN 978-1466561595,
## <URL: http://www.crcpress.com/product/isbn/9781466561595>.

#Térképes megjelenítés
base::print(citation("maps"),style="text") #Megjelenítés térképeken

## Becker OScbRA, Minka ARWRvbRBEbTP, Deckmyn. A (2018). _maps: Draw
## Geographical Maps_. R package version 3.3.0, <URL:
## https://CRAN.R-project.org/package=maps>.

base::print(citation("mapdata"),style="text") #Megjelenítés térképeken

## Becker OScbRA, Brownrigg. ARWRvbR (2018). _mapdata: Extra Map
## Databases_. R package version 2.3.0, <URL:
## https://CRAN.R-project.org/package=mapdata>.

base::print(citation("maptools"),style="text") #Megjelenítés térképeken

## Bivand R, Lewin-Koh N (2021). _maptools: Tools for Handling Spatial
## Objects_. R package version 1.1-1, <URL:
## https://CRAN.R-project.org/package=maptools>.

#Interaktív és 3D megjelenítések
base::print(citation("shiny"),style="text") #Interaktív megjelenítés

## Chang W, Cheng J, Allaire J, Sievert C, Schloerke B, Xie Y, Allen J,
## McPherson J, Dipert A, Borges B (2021). _shiny: Web Application
## Framework for R_. R package version 1.6.0, <URL:
## https://CRAN.R-project.org/package=shiny>.

base::print(citation("htmltools"),style="text") #Interaktív megjelnítés

## Cheng J, Sievert C, Chang W, Xie Y, Allen J (2021). _htmltools: Tools
## for HTML_. R package version 0.5.1.1, <URL:
## https://CRAN.R-project.org/package=htmltools>.

base::print(citation("rgl"),style="text") #3D megjelítés

## Murdoch D, Adler D (2021). _rgl: 3D Visualization Using OpenGL_. R
## package version 0.106.8, <URL: https://CRAN.R-project.org/package=rgl>.

base::print(citation("rglwidget"),style="text") #3D interaktív megjelenítés

## Murdoch D (2016). _rglwidget: 'rgl' in 'htmlwidgets' Framework_. R
## package version 0.2.1, <URL:
## https://CRAN.R-project.org/package=rglwidget>.

base::print(citation("ggplot2"),style="text") #Interaktív ábrák megjelnítése

## Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_.
## Springer-Verlag New York. ISBN 978-3-319-24277-4, <URL:
## https://ggplot2.tidyverse.org>.

base::print(citation("ggrepel"),style="text") #Interaktív ábrák megjelnítése

## Slowikowski K (2021). _ggrepel: Automatically Position Non-Overlapping
## Text Labels with 'ggplot2'_. R package version 0.9.1, <URL:
## https://CRAN.R-project.org/package=ggrepel>.

base::print(citation("ggdendro"),style="text") #Dendogram megjelenítése

## de Vries A, Ripley BD (2020). _ggdendro: Create Dendrograms and Tree
## Diagrams Using 'ggplot2'_. R package version 0.1.22, <URL:
## https://CRAN.R-project.org/package=ggdendro>.

base::print(citation("d3heatmap"),style="text") #Interaktív hőtérkép

## Cheng J, Galili T (2018). _d3heatmap: Interactive Heat Maps Using
## 'htmlwidgets' and 'D3.js'_. R package version 0.6.1.2, <URL:
## https://CRAN.R-project.org/package=d3heatmap>.

base::print(citation("igraph"),style="text") #Gráf-megjelenítés

## Csardi G, Nepusz T (2006). "The igraph software package for complex
## network research." _InterJournal_, *Complex Systems*, 1695. <URL:
## https://igraph.org>.

base::print(citation("ggsci"),style="text") #Simpson-paletta

## Xiao N (2018). _ggsci: Scientific Journal and Sci-Fi Themed Color
## Palettes for 'ggplot2'_. R package version 2.9, <URL:
## https://CRAN.R-project.org/package=ggsci>.

base::print(citation("plotly"),style="text") # Plotly - interaktív diagram

## Sievert C (2020). _Interactive Web-Based Data Visualization with R,
## plotly, and shiny_. Chapman and Hall/CRC. ISBN 9781138331457, <URL:
## https://plotly-r.com>.

Szükséges adatok beolvasása

Jelen html-ben az eredeti U21 adatokat használjuk.

load("U21_2014.RData")
rownames(U21_rank) <- rownames(U21_filtered)
orig_mtx <- as.matrix(U21_filtered)
#U21_filtered: 50 x 24 dataframe, amely az eredeti 2014-es U21-es adatokat tartalmazza
#U21_rank: 50 x 1 mátrix (vektor) U21 rangsor
#orig_mtx: 50 x 24 mátrix, eredeti U21 indikátorok az 50 országra

2. Kétlépéses klaszterezhetőség vizsgálata

Első lépésként megvizsgáljuk, hogy az adatok alkalmasak-e a klaszterezésre. A klaszterezhetőséget mind a változók, mind pedig az országok esetében ellenőrizni kell.

Indikátorok (oszlopok) klaszterezhetősége

Hopkins mutató számítása

A klaszterezhetőség ellenőrzésében a Hopkins mutató (H) nyújt segítséget. Ha a Hopkins mutató (H) nagyobb, mint 0,5, akkor nem érdemes klaszterezni. A változók klaszterezésénél a Spearman-korrelációt választjuk hasonlósági, az 1-Spearman-korrelációt távolsági mutatónak, mivel a rangkorreláció nem lineáris kapcsolatokat is kimutat.

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")
print(paste("H:= ",hopkins(dist.m, n=nrow(dist.m)-1, byrow = F, header = F)))

## [1] "H:=  0.311994877382919"

Mivel a Hopkins mutató kisebb, mint 0,5, ezért a változókat lehet klaszterezni.

Következő lépésként a klaszterek előzetes számának meghatározása érdekében kirajzoltatjuk a klaszterezettséget.

#Klaszterezettség meghatározása
clustend <- get_clust_tendency(scale(dist.m), nrow(dist.m)-1)
#Klaszterezettség kirajzolása
clustend$plot + 
  scale_fill_gradient(low = "yellow", high = "blue")

1. ábra A változók klaszterezettsége

Az 1. ábra alapján (ahol az adatmátrix cellái azok értékeinek nagysága szerint színezettek) 2-4 klaszter (ami itt indikátor csoport) meghatározása valószínűsíthető, de ezt további vizsgálatoknak kell megerősítenie.

Hierarchikus klaszterezés az indikátorok esetében

A következő lépésként hierarchikus klaszterezést alkalmazunk az indikátorokra vonatkozóan. A klaszterezés eredménye a 2. ábrán látható.

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="A hierarchikus klaszterezés eredménye az indikátorok esetében") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()) +
  coord_flip()

2. ábra A hierarchikus klaszterezés eredménye az indikátorok esetében

A hierarchikus klaszterezés is 2-4 klaszter (ebben az esetben indikátorcsoport) meglétét valószínűsíti. A távolságmátrix megjelenítése további segítséget nyújt a klaszterek számának meghatározásában. A 3. ábrán az indikátorokra vonatkozó távolságmátrix látható, ahol a kék cellák -1-hez közeli értékeket, míg a sötétvörös cellák 1-hez közeli értékeket jelentenek.

cormat <- rquery.cormat(dist.m, graphType="heatmap")

3. ábra Távolságmátrix az indikátorok esetében

Van Mechelen (2015) kritériumai alapján a klaszterek meghatározásánál az alábbi szempontok szerint kell eljárni.

1. Klasztereken belül az átlagos távolság (dissimilarity) alacsony legyen.
2. Klaszterek között az átlagos távolság (dissimilarity) nagy legyen.
3. A klaszterek stabilak (többszöri újrafuttás után sem kapunk lényegesen különböző eredményeket).
4. A klaszterek körülbelül hasonló elemszámúak legyenek.
5. A klaszterek száma alacsony.

A következő lépésben összehasonlítjuk azokat az eseteket, amikor 2, 3 illetve 4 klasztert alakítunk ki a változókra vonatkozóan.

cluster.stats(dist.m,cutree(hc, k=2))

## $n
## [1] 24
## 
## $cluster.number
## [1] 2
## 
## $cluster.size
## [1] 21  3
## 
## $min.cluster.size
## [1] 3
## 
## $noisen
## [1] 0
## 
## $diameter
## [1] 1.2498799 0.7726897
## 
## $average.distance
## [1] 0.5610211 0.6223547
## 
## $median.distance
## [1] 0.5672502 0.6203722
## 
## $separation
## [1] 0.68497 0.68497
## 
## $average.toother
## [1] 1.015888 1.015888
## 
## $separation.matrix
##         [,1]    [,2]
## [1,] 0.00000 0.68497
## [2,] 0.68497 0.00000
## 
## $ave.between.matrix
##          [,1]     [,2]
## [1,] 0.000000 1.015888
## [2,] 1.015888 0.000000
## 
## $average.between
## [1] 1.015888
## 
## $average.within
## [1] 0.5686878
## 
## $n.between
## [1] 63
## 
## $n.within
## [1] 213
## 
## $max.diameter
## [1] 1.24988
## 
## $min.separation
## [1] 0.68497
## 
## $within.cluster.ss
## [1] 4.111321
## 
## $clus.avg.silwidths
##         1         2 
## 0.4384823 0.3867175 
## 
## $avg.silwidth
## [1] 0.4320117
## 
## $g2
## NULL
## 
## $g3
## NULL
## 
## $pearsongamma
## [1] 0.6564944
## 
## $dunn
## [1] 0.5480287
## 
## $dunn2
## [1] 1.63233
## 
## $entropy
## [1] 0.3767702
## 
## $wb.ratio
## [1] 0.5597936
## 
## $ch
## [1] 10.44002
## 
## $cwidegap
## [1] 0.5204666 0.6203722
## 
## $widestgap
## [1] 0.6203722
## 
## $sindex
## [1] 0.68497
## 
## $corrected.rand
## NULL
## 
## $vi
## NULL

cluster.stats(dist.m,cutree(hc, k=3))

## $n
## [1] 24
## 
## $cluster.number
## [1] 3
## 
## $cluster.size
## [1] 15  3  6
## 
## $min.cluster.size
## [1] 3
## 
## $noisen
## [1] 0
## 
## $diameter
## [1] 1.0665562 0.7726897 0.9575082
## 
## $average.distance
## [1] 0.4433945 0.6223547 0.7098425
## 
## $median.distance
## [1] 0.4366579 0.6203722 0.6949694
## 
## $separation
## [1] 0.2946191 0.6849700 0.2946191
## 
## $average.toother
## [1] 0.7970161 1.0158882 0.7187454
## 
## $separation.matrix
##           [,1]      [,2]      [,3]
## [1,] 0.0000000 0.7845616 0.2946191
## [2,] 0.7845616 0.0000000 0.6849700
## [3,] 0.2946191 0.6849700 0.0000000
## 
## $ave.between.matrix
##           [,1]      [,2]      [,3]
## [1,] 0.0000000 1.0441516 0.6734484
## [2,] 1.0441516 0.0000000 0.9452299
## [3,] 0.6734484 0.9452299 0.0000000
## 
## $average.between
## [1] 0.8144531
## 
## $average.within
## [1] 0.5323766
## 
## $n.between
## [1] 153
## 
## $n.within
## [1] 123
## 
## $max.diameter
## [1] 1.066556
## 
## $min.separation
## [1] 0.2946191
## 
## $within.cluster.ss
## [1] 3.425579
## 
## $clus.avg.silwidths
##           1           2           3 
##  0.34829660  0.33917439 -0.05553832 
## 
## $avg.silwidth
## [1] 0.2461976
## 
## $g2
## NULL
## 
## $g3
## NULL
## 
## $pearsongamma
## [1] 0.5722906
## 
## $dunn
## [1] 0.276234
## 
## $dunn2
## [1] 0.9487293
## 
## $entropy
## [1] 0.9002561
## 
## $wb.ratio
## [1] 0.6536614
## 
## $ch
## [1] 8.082114
## 
## $cwidegap
## [1] 0.4325824 0.6203722 0.5920308
## 
## $widestgap
## [1] 0.6203722
## 
## $sindex
## [1] 0.2946191
## 
## $corrected.rand
## NULL
## 
## $vi
## NULL

cluster.stats(dist.m,cutree(hc, k=4))

## $n
## [1] 24
## 
## $cluster.number
## [1] 4
## 
## $cluster.size
## [1]  5 10  3  6
## 
## $min.cluster.size
## [1] 3
## 
## $noisen
## [1] 0
## 
## $diameter
## [1] 1.0665562 0.4296074 0.7726897 0.9575082
## 
## $average.distance
## [1] 0.6794132 0.2371534 0.6223547 0.7098425
## 
## $median.distance
## [1] 0.7161004 0.2366260 0.6203722 0.6949694
## 
## $separation
## [1] 0.2746400 0.2746400 0.6849700 0.2946191
## 
## $average.toother
## [1] 0.7358357 0.6848112 1.0158882 0.7187454
## 
## $separation.matrix
##           [,1]      [,2]      [,3]      [,4]
## [1,] 0.0000000 0.2746400 0.7845616 0.5761566
## [2,] 0.2746400 0.0000000 0.8014765 0.2946191
## [3,] 0.7845616 0.8014765 0.0000000 0.6849700
## [4,] 0.5761566 0.2946191 0.6849700 0.0000000
## 
## $ave.between.matrix
##           [,1]      [,2]      [,3]      [,4]
## [1,] 0.0000000 0.5818078 1.0129007 0.8540163
## [2,] 0.5818078 0.0000000 1.0597770 0.5831645
## [3,] 1.0129007 1.0597770 0.0000000 0.9452299
## [4,] 0.8540163 0.5831645 0.9452299 0.0000000
## 
## $average.between
## [1] 0.7571513
## 
## $average.within
## [1] 0.4956133
## 
## $n.between
## [1] 203
## 
## $n.within
## [1] 73
## 
## $max.diameter
## [1] 1.066556
## 
## $min.separation
## [1] 0.27464
## 
## $within.cluster.ss
## [1] 3.01833
## 
## $clus.avg.silwidths
##          1          2          3          4 
## -0.1273645  0.5697175  0.3391744 -0.1851225 
## 
## $avg.silwidth
## [1] 0.2069642
## 
## $g2
## NULL
## 
## $g3
## NULL
## 
## $pearsongamma
## [1] 0.5264653
## 
## $dunn
## [1] 0.2575017
## 
## $dunn2
## [1] 0.8196294
## 
## $entropy
## [1] 1.298077
## 
## $wb.ratio
## [1] 0.6545763
## 
## $ch
## [1] 6.723373
## 
## $cwidegap
## [1] 0.6839779 0.2432582 0.6203722 0.5920308
## 
## $widestgap
## [1] 0.6839779
## 
## $sindex
## [1] 0.27464
## 
## $corrected.rand
## NULL
## 
## $vi
## NULL

A fenti eredmények alapján azonban nem lehet eldönteni, hogy mennyi (2, 3 vagy 4) klaszter képezhető az indikátorokon. Tehát a klaszterek értelmezhetősége fog erről dönteni.

Többdimenziós skálázás alkalmazása az indikátorokra

A klaszterek értelmezhetőségének vizsgálatára a következő lépésben a többdimenziós skálázást alkalmazzuk az indikátorokra vonatkozóan. Elsőként két klasztert alakítunk ki, az erre vonatkozó eredmények a 4. ábrán láthatók.

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

## initial  value 0.234620 
## iter   5 value 0.102313
## iter  10 value 0.079063
## iter  15 value 0.058780
## iter  20 value 0.035771
## iter  25 value 0.022966
## iter  30 value 0.013233
## final  value 0.009420 
## converged

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 = "Klaszterek") +
  labs(x="", y="", title="Többdimenziós skálázás az indikátorokra vonatkozóan - 2 klaszter") + theme_bw()

4. ábra Többdimenziós skálázás az indikátorokra vonatkozóan - 2 klaszter esetén

Ahogy az a 4. ábrán is látszik, a két klaszter nehezen értelmezhető. Az első klaszterbe 19 indikátor került, míg a másodikba csak 3.

Az 5. ábrán a többdimenziós skálázás eredménye látható abban az esetben, ha 3 klasztert hozunk létre.

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

## initial  value 0.234620 
## iter   5 value 0.102313
## iter  10 value 0.079063
## iter  15 value 0.058780
## iter  20 value 0.035771
## iter  25 value 0.022966
## iter  30 value 0.013233
## final  value 0.009420 
## converged

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 = "Klaszterek") +
  labs(x="", y="", title="Többdimenziós skálázás az indikátorokra vonatkozóan - 3 klaszter") + theme_bw()

5. ábra Többdimenziós skálázás az indikátorokra vonatkozóan - 3 klaszter esetén

Az 5. ábra alapján a 3 klaszter (azaz 3 indikátorcsoport) is nehezen értelmezhető. Az első klaszterbe 15, a második klaszterbe 3, a harmadik klaszterbe pedig 6 indikátor került. A kapcsolatokat mérő indikátorok (C1, C2, C3 és C4) azonban már egy klaszterbe kerültek.

A következő lépésben megnézzük, hogy ha 4 klasztert hozunk létre, akkor hogyan alakul az indikátorok csoportosítása. Ennek eredménye a 6. ábrán látható.

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

## initial  value 0.234620 
## iter   5 value 0.102313
## iter  10 value 0.079063
## iter  15 value 0.058780
## iter  20 value 0.035771
## iter  25 value 0.022966
## iter  30 value 0.013233
## final  value 0.009420 
## converged

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 = "Klaszterek") +
  labs(x="", y="", title="Többdimenziós skálázás az indikátorokra vonatkozóan - 4 klaszter'") + theme_bw()

6. ábra Többdimenziós skálázás az indikátorokra vonatkozóan - 4 klaszter esetén

Ha négy klasztert hozunk létre, akkor az első klaszterbe 5, a másodikba 10, a harmadikba 3, a negyedikbe 6 indikátor kerül. Tisztán csak a C1-C4 változókat tartalmazó kapcsolatok változócsoportja jelenik meg a negyedik klaszterben, a többi klaszter esetében nagyon vegyes az indikátorok összetétele. Négy klaszter (változócsoport) is nehezen értelmezhető, ezért szükség van a változók számának redukálására, amit például főkomponens-elemzéssel tehetünk meg.

Országok (sorok) klaszterezhetősége

Hopkins mutató számítása

A klaszterezhetőség ellenőrzésében a Hopkins mutató (H) nyújt segítséget. Ha a Hopkins mutató (H) nagyobb, mint 0,5, akkor nem érdemes klaszterezni. Az országok klaszterezésénél a távolságon az Euklideszi távolságot értjük.

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))
print(paste("H:=",hopkins(dist.c, n=nrow(dist.c)-1, byrow = F, header = F)))

## [1] "H:= 0.220381143328058"

Mivel a Hopkins mutató kisebb, mint 0,5, ezért az országokat (azaz a sorokat) lehet klaszterezni.

Következő lépésként a klaszterek előzetes számának meghatározása érdekében kirajzoltatjuk a klaszterezettséget.

#Klaszterezettség
clustend <- get_clust_tendency(scale(dist.c), nrow(dist.c)-1)
#Klaszterezettség kirajzolása
clustend$plot + 
  scale_fill_gradient(low = "yellow", high = "blue")

7. ábra Az országok klaszterezettsége

A 7. ábra alapján 3 klaszter (azaz országcsoport) valószínűsíthető, de ezt további vizsgálatoknak kell megerősítenie.

###Hierarchikus klaszterezés az országok esetében

A következő lépésként hierarchikus klaszterezést alkalmazunk az országokra vonatkozóan. A klaszterezés eredménye a 8. ábrán látható.

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="A hierarchikus klaszterezés eredménye az országok esetében") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()) +
  coord_flip()

8. ábra A hierarchikus klaszterezés eredménye az országok esetében

A hierarchikus klaszterezés 3-4 klaszter (itt országcsoport) meglétét valószínűsíti.

A távolságmátrix megjelenítése további segítséget nyújt a klaszterek számának meghatározásában. A 9. ábrán az országokra vonatkozó távolságmátrix látható, ahol a kék cellák -1-hez közeli értékeket, míg a sötétvörös cellák 1-hez közeli értékeket jelentenek.

cormat <- rquery.cormat(dist.c, graphType="heatmap")

9. ábra Távolságmátrix az országok esetében

A 9. ábrán látható távolságmátrix alapján 3-4 klaszter valószínűsíthető. A következő lépésben összehasonlítjuk azokat az eseteket, amikor 3 vagy 4 klasztert alakítunk ki az országokra vonatkozóan.

cluster.stats(dist.c,cutree(hc, k=3))

## $n
## [1] 50
## 
## $cluster.number
## [1] 3
## 
## $cluster.size
## [1] 21 15 14
## 
## $min.cluster.size
## [1] 14
## 
## $noisen
## [1] 0
## 
## $diameter
## [1] 172.7304 187.5483 131.6234
## 
## $average.distance
## [1]  96.74970 114.40191  94.56826
## 
## $median.distance
## [1]  92.82482 111.94521  98.24285
## 
## $separation
## [1] 57.08949 64.51798 57.08949
## 
## $average.toother
## [1] 165.3137 174.3910 131.3474
## 
## $separation.matrix
##           [,1]      [,2]     [,3]
## [1,]   0.00000 136.35107 57.08949
## [2,] 136.35107   0.00000 64.51798
## [3,]  57.08949  64.51798  0.00000
## 
## $ave.between.matrix
##          [,1]     [,2]     [,3]
## [1,]   0.0000 200.0512 128.0951
## [2,] 200.0512   0.0000 135.9007
## [3,] 128.0951 135.9007   0.0000
## 
## $average.between
## [1] 157.7719
## 
## $average.within
## [1] 101.4346
## 
## $n.between
## [1] 819
## 
## $n.within
## [1] 406
## 
## $max.diameter
## [1] 187.5483
## 
## $min.separation
## [1] 57.08949
## 
## $within.cluster.ss
## [1] 258714.8
## 
## $clus.avg.silwidths
##         1         2         3 
## 0.2440727 0.1473629 0.1832766 
## 
## $avg.silwidth
## [1] 0.1980369
## 
## $g2
## NULL
## 
## $g3
## NULL
## 
## $pearsongamma
## [1] 0.566333
## 
## $dunn
## [1] 0.3043989
## 
## $dunn2
## [1] 1.119693
## 
## $entropy
## [1] 1.081972
## 
## $wb.ratio
## [1] 0.6429189
## 
## $ch
## [1] 24.42036
## 
## $cwidegap
## [1]  98.27018 125.24363 102.04286
## 
## $widestgap
## [1] 125.2436
## 
## $sindex
## [1] 61.91874
## 
## $corrected.rand
## NULL
## 
## $vi
## NULL

cluster.stats(dist.c,cutree(hc, k=4))

## $n
## [1] 50
## 
## $cluster.number
## [1] 4
## 
## $cluster.size
## [1] 21 11 14  4
## 
## $min.cluster.size
## [1] 4
## 
## $noisen
## [1] 0
## 
## $diameter
## [1] 172.7304 148.9883 131.6234 185.2957
## 
## $average.distance
## [1]  96.74970 106.13002  94.56826 128.13876
## 
## $median.distance
## [1]  92.82482 108.07670  98.24285 135.48427
## 
## $separation
## [1] 57.08949 60.81710 57.08949 60.81710
## 
## $average.toother
## [1] 165.3137 158.4569 131.3474 186.9005
## 
## $separation.matrix
##           [,1]      [,2]      [,3]     [,4]
## [1,]   0.00000 136.35107  57.08949 201.0139
## [2,] 136.35107   0.00000  64.51798  60.8171
## [3,]  57.08949  64.51798   0.00000 102.1068
## [4,] 201.01385  60.81710 102.10676   0.0000
## 
## $ave.between.matrix
##          [,1]     [,2]     [,3]     [,4]
## [1,]   0.0000 187.6245 128.0951 234.2247
## [2,] 187.6245   0.0000 124.8736 122.8686
## [3,] 128.0951 124.8736   0.0000 166.2251
## [4,] 234.2247 122.8686 166.2251   0.0000
## 
## $average.between
## [1] 155.9924
## 
## $average.within
## [1] 100.7137
## 
## $n.between
## [1] 863
## 
## $n.within
## [1] 362
## 
## $max.diameter
## [1] 185.2957
## 
## $min.separation
## [1] 57.08949
## 
## $within.cluster.ss
## [1] 247876.6
## 
## $clus.avg.silwidths
##           1           2           3           4 
##  0.24407268  0.07194809  0.15158446 -0.03399656 
## 
## $avg.silwidth
## [1] 0.158063
## 
## $g2
## NULL
## 
## $g3
## NULL
## 
## $pearsongamma
## [1] 0.5576146
## 
## $dunn
## [1] 0.3080993
## 
## $dunn2
## [1] 0.9588711
## 
## $entropy
## [1] 1.255947
## 
## $wb.ratio
## [1] 0.6456321
## 
## $ch
## [1] 17.30098
## 
## $cwidegap
## [1]  98.27018 106.99219 102.04286 174.91549
## 
## $widestgap
## [1] 174.9155
## 
## $sindex
## [1] 60.06623
## 
## $corrected.rand
## NULL
## 
## $vi
## NULL

A fenti eredmények alapján azonban nem lehet eldönteni, hogy hány (3 vagy 4) klaszter képezhető az országokon. Tehát az országok esetében is a klaszterek értelmezhetősége fog erről dönteni.

Többdimenziós skálázás alkalmazása az országokra

A klaszterek értelmezhetőségének vizsgálatára a következő lépésben a többdimenziós skálázást alkalmazzuk az országokra vonatkozóan. Elsőként három klasztert (azaz három országcsoportot) alakítunk ki, az erre vonatkozó eredmények a 10. ábrán láthatók.

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

## initial  value 0.728168 
## final  value 0.727938 
## converged

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 = "Klaszterek") +
  labs(x="", y="", title="Többdimenziós skálázás az országokra vonatkozóan - 3 klaszter") + theme_bw()

10. ábra Többdimenziós skálázás az országokra vonatkozóan - 3 klaszter esetén

Ahogy az a 10. ábrán is látható, az első klaszterbe 21, a másodikba 15, a harmadikba pedig 14 ország került.

A következő lépésben megnézzük, hogy ha 4 klasztert hozunk létre, akkor hogyan alakul az országok csoportosítása. Ennek eredménye a 11. ábrán látható.

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

## initial  value 0.728168 
## final  value 0.727938 
## converged

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 = "Klaszterek") +
  labs(x="", y="", title="Többdimenziós skálázás az országokra vonatkozóan - 4 klaszter") + theme_bw()

11. ábra Többdimenziós skálázás az országokra vonatkozóan - 4 klaszter esetén

Abban az esetben, ha négy klasztert hozunk létre, akkor az első klaszterbe 21, a másodikba 11, a harmadikba 14, míg a negyedikbe 4 ország került. Az előző esethez képest (amikor három klasztert hoztunk létre) az első és a harmadik klaszterben levő országok nem változtak. A negyedik klaszterbe azok az országok kerültek (Dánia, Svédország, Svájc, Egyesült Államok), amelyek az előző esetben a második klaszterhez tartoztak. A Van Mechelen (1993) által megfogalmazott elveket leginkább 4 klaszter teljesíti.

Klaszterezhetőség eredménye: Egyértelműen látszik, hogy mind a változók, mind az országok szerint lehet csoportosításokat végezni (lásd: Hopkins mutató). Ugyanakkor a kapott klaszterek értelmezhetősége meglehetősen nehézkes. A változócsoportok vegyesen tartalmaznak mutatókat, így szükség van az országok csoportosítása előtt modellredukcióra.

3. Modellredukció

Főkomponens-elemzés

Első lépésként a főkomponensek/faktorok (esetünkben indikátorcsoportok) számának meghatározására van szükség. Ehhez segítséget nyújt a 12. ábrán látható könyökdiagram.

scree(U21_filtered)

12. ábra Könyökdiagram

A 12. ábra alapján a főkomponens-elemzésnél maximum 5, faktorelemzésnél maximum 3 indikátorcsoportot érdemes megtartani.

Az U21 4 fő indikátorcsoportot tartalmaz: erőforrások (R1 - R5), környezet (E1 - E4), kapcsolatok (C1 - C6) és kimenetek (O1 - O9). A főkomponens elemzés esetén érdemes kiindulni ezekből az előre meghatározott indikátorcsoportokból. Ekkor azt várjuk, hogy 4 jól interpretálható főkomponenst kapunk, amely 4 indikátorcsoportot kódol. Az elemzés során megkapjuk, hogy az eredményül kapott főkomponensekre az eredeti változók mennyire illeszkednek jól.

R (erőforrások) főkomponens meghatározása

Elsőként az R főkomponenst határozzuk meg, amely az R1 - R5 indikátorokat tartalmazza, azaz azokat az indikátorokat, amelyek az egyes országok erőforrásait írják le.

Rs<-U21_filtered[,c("R1","R2","R3","R4","R5")]
KMO(cor(Rs))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(Rs))
## Overall MSA =  0.65
## MSA for each item = 
##   R1   R2   R3   R4   R5 
## 0.75 0.54 0.66 0.68 0.60

fitR <- principal(Rs,nfactors=1)
summary(fitR) # print variance accounted for 

## 
## Factor analysis with Call: principal(r = Rs, nfactors = 1)
## 
## Test of the hypothesis that 1 factor is sufficient.
## The degrees of freedom for the model is 5  and the objective function was  1.28 
## The number of observations was  50  with Chi Square =  58.63  with prob <  2.3e-11 
## 
## The root mean square of the residuals (RMSA) is  0.15

fitR$communality

##        R1        R2        R3        R4        R5 
## 0.4804492 0.3481106 0.7782275 0.8074245 0.8522280

pcaR <- prcomp(Rs, scale = TRUE)
eigs <- pcaR$sdev^2
varR <- eigs[1] / sum(eigs)
print(paste("Variance rate=",varR))

## [1] "Variance rate= 0.653287948596533"

Mivel a KMO értéke nagyobb, mint 0.5, ezért a főkomponens/faktorelemzés elvégezhető. Minden kommunalitás nagyobb, mint 0.25, amely azt jelenti, hogy az eredeti változók megfelelően illeszkednek az eredményváltozóra. Továbbá, a varianciahiányad 0.5-nél nagyobb, amely azt jelenti, hogy a modellredukció után is megmarad az információ legalább 50%-a. Ezek alapján a fentebb létrehozott R főkomponens megfelelő, felhasználható a további elemzések során.

E (környezet) főkomponens meghatározása

A következő lépésben az E főkomponenst hozzuk létre, amely az E1 - E4 indikátorokat tartalmazza, azaz azokat a mutatókat, amelyek a környezetet írják le.

Es<-U21_filtered[,c("E1","E2","E3","E4")]
KMO(cor(Es))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(Es))
## Overall MSA =  0.55
## MSA for each item = 
##   E1   E2   E3   E4 
## 0.61 0.56 0.53 0.54

fitE <- principal(Es,nfactors=1)
summary(fitE) # print variance accounted for 

## 
## Factor analysis with Call: principal(r = Es, nfactors = 1)
## 
## Test of the hypothesis that 1 factor is sufficient.
## The degrees of freedom for the model is 2  and the objective function was  0.28 
## The number of observations was  50  with Chi Square =  13.06  with prob <  0.0015 
## 
## The root mean square of the residuals (RMSA) is  0.21

fitE$communality

##        E1        E2        E3        E4 
## 0.3235923 0.1889422 0.5982086 0.6393162

pcaE <- prcomp(Es, scale = TRUE)
eigs <- pcaE$sdev^2
varE <- eigs[1] / sum(eigs)
print(paste("Variance rate=",varE))

## [1] "Variance rate= 0.43751482380739"

A fenti eredmények alapján látható, hogy az E2 mutató esetében a kommunalitás kisebb, mint 0.25 (0.189), tehát E2 nem illeszkedik megfelelően az eredményváltozóra. Ugyanakkor, ha a következő lépésben elhagynánk E2-t, E1 esetében lenne a kommunalitás kisebb, mint 0.25. A legtöbb változó akkor őrizhető meg, ha E3 változót hagyjuk el. Így a következő lépésben elhagyjuk E3-at.

Es<-U21_filtered[,c("E1","E2","E4")]
KMO(cor(Es))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(Es))
## Overall MSA =  0.55
## MSA for each item = 
##   E1   E2   E4 
## 0.53 0.54 0.61

fitE <- principal(Es,nfactors=1)
summary(fitE) 

## 
## Factor analysis with Call: principal(r = Es, nfactors = 1)
## 
## Test of the hypothesis that 1 factor is sufficient.
## The degrees of freedom for the model is 0  and the objective function was  0.22 
## The number of observations was  50  with Chi Square =  10.42  with prob <  NA 
## 
## The root mean square of the residuals (RMSA) is  0.25

fitE$communality

##        E1        E2        E4 
## 0.6167214 0.5175482 0.2970137

pcaE <- prcomp(Es, scale = TRUE)
eigs <- pcaE$sdev^2
varE <- eigs[1] / sum(eigs)
print(paste("Variance rate=",varE))

## [1] "Variance rate= 0.477094428557052"

E3 elhagyása után a KMO 0.5 feletti, minden változó esetében a kommunalitás nagyobb, mint 0.25. A varianciahányad értéke azonban 0.5 alatti (0.477), amely azt jelenti, hogy a modellredukció után az információ legalább 50%-a elveszik.

C (kapcsolatok) főkomponens meghatározása

A következő lépésben a C főkomponenst hozzuk létre, amely a C1 - C6 mutatókat tartalmazza, és az országok kapcsolatait írja majd le.

Cs<-U21_filtered[,c("C1","C2","C3","C4","C5","C6")]
KMO(cor(Cs))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(Cs))
## Overall MSA =  0.67
## MSA for each item = 
##   C1   C2   C3   C4   C5   C6 
## 0.57 0.66 0.37 0.70 0.82 0.78

fitC <- principal(Cs,nfactors=1)
summary(fitC) 

## 
## Factor analysis with Call: principal(r = Cs, nfactors = 1)
## 
## Test of the hypothesis that 1 factor is sufficient.
## The degrees of freedom for the model is 9  and the objective function was  0.53 
## The number of observations was  50  with Chi Square =  24.07  with prob <  0.0042 
## 
## The root mean square of the residuals (RMSA) is  0.14

fitC$communality

##        C1        C2        C3        C4        C5        C6 
## 0.4456108 0.5144444 0.1290487 0.6247154 0.6170022 0.5594279

pcaC <- prcomp(Cs, scale = TRUE)
eigs <- pcaC$sdev^2
varC <- eigs[1] / sum(eigs)
print(paste("Variance rate=",varC))

## [1] "Variance rate= 0.481708220917479"

C3 kommunalitása alacsonyabb, mint 0.25, tehát nem illeszkedik a főkomponensre, ezért ezt a mutatót a következő lépésben elhagyjuk.

Cs<-U21_filtered[,c("C1","C2","C4","C5","C6")]
KMO(cor(Cs))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(Cs))
## Overall MSA =  0.79
## MSA for each item = 
##   C1   C2   C4   C5   C6 
## 0.78 0.79 0.80 0.79 0.80

fitC <- principal(Cs,nfactors=1)
summary(fitC)

## 
## Factor analysis with Call: principal(r = Cs, nfactors = 1)
## 
## Test of the hypothesis that 1 factor is sufficient.
## The degrees of freedom for the model is 5  and the objective function was  0.23 
## The number of observations was  50  with Chi Square =  10.55  with prob <  0.061 
## 
## The root mean square of the residuals (RMSA) is  0.13

fitC$communality

##        C1        C2        C4        C5        C6 
## 0.5120676 0.5011992 0.5971696 0.6207606 0.5728886

pcaC <- prcomp(Cs, scale = TRUE)
eigs <- pcaC$sdev^2
varC <- eigs[1] / sum(eigs)
print(paste("Variance rate=",varC))

## [1] "Variance rate= 0.560817123312952"

A megmaradt változókon a KMO értéke magasabb, mint 0.5, minden változó esetében a kommunalitás nagyobb, mint 0.25, vagyis az eredeti változók megfelelően illeszkednek az eredményváltozóra. Továbbá, a megőrzött varianciahányad is 50% feletti (0.56), tehát a modellredukció után is megmarad az eredeti információ legalább 50%-a.

O (kimenet) főkomponens meghatározása

A legutolsó lépésben az O főkomponenst hozzuk létre, amely az O1 - O9 mutatókat tartalmazza, és az országok kimeneteit írja le.

Os<-U21_filtered[,c("O1","O2","O3","O4","O5","O6","O7","O8","O9")]
KMO(cor(Os))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(Os))
## Overall MSA =  0.73
## MSA for each item = 
##   O1   O2   O3   O4   O5   O6   O7   O8   O9 
## 0.40 0.73 0.77 0.79 0.63 0.83 0.74 0.85 0.64

fitO <- principal(Os,nfactors=1)
summary(fitO) 

## 
## Factor analysis with Call: principal(r = Os, nfactors = 1)
## 
## Test of the hypothesis that 1 factor is sufficient.
## The degrees of freedom for the model is 27  and the objective function was  3.28 
## The number of observations was  50  with Chi Square =  145.84  with prob <  2.9e-18 
## 
## The root mean square of the residuals (RMSA) is  0.14

fitO$communality

##           O1           O2           O3           O4           O5           O6 
## 0.1066000549 0.8603383885 0.7882600079 0.7698125877 0.6266131079 0.3970218624 
##           O7           O8           O9 
## 0.5495702936 0.7077575306 0.0005444639

pcaO <- prcomp(Os, scale = TRUE)
eigs <- pcaC$sdev^2
varO <- eigs[1] / sum(eigs)
print(paste("Variance rate=",varO))

## [1] "Variance rate= 0.560817123312952"

A fenti eredmények alapján az O1, és az O9 esetében is a kommunalitás alacsonyabb, mint 0.25, ugyanakkor a varianciahányad megfelelő (0.56). Két lehetőség kínálkozik: egy főkomponens, vagy két faktor használata.

Elsőként 1 főkomponenst hozunk létre, amelynek az eredményei lentebb láthatók.

Os<-U21_filtered[,c("O2","O3","O4","O5","O6","O7","O8")]
KMO(cor(Os))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(Os))
## Overall MSA =  0.79
## MSA for each item = 
##   O2   O3   O4   O5   O6   O7   O8 
## 0.74 0.79 0.81 0.75 0.82 0.78 0.85

fitO <- principal(Os,nfactors=1)
summary(fitO)

## 
## Factor analysis with Call: principal(r = Os, nfactors = 1)
## 
## Test of the hypothesis that 1 factor is sufficient.
## The degrees of freedom for the model is 14  and the objective function was  1.62 
## The number of observations was  50  with Chi Square =  73.02  with prob <  5.4e-10 
## 
## The root mean square of the residuals (RMSA) is  0.1

fitO$communality

##        O2        O3        O4        O5        O6        O7        O8 
## 0.8860009 0.7906704 0.7907994 0.5663786 0.4074405 0.5578915 0.7229660

pcaO <- prcomp(Os, scale = TRUE)
eigs <- pcaC$sdev^2
varO <- eigs[1] / sum(eigs)
print(paste("Variance rate=",varO))

## [1] "Variance rate= 0.560817123312952"

A KMO értéke nagyobb, mint 0.5, minden változó esetében a kommunalitás nagyobb, mint 0.25, a varinciahányad értéke 0.56, vagyis a modellredukció után is megmarad az eredeti információ legalább 50%-a. Két fontos mutató (O1 és O9) azonban kiesett a modellből. Az O1 mutató a felsőoktatási intézmények által publikált cikkek számát jelenti, míg az O9 a diplomás munkanélküliek aránya. Ezeket a mutatókat is érdemes és szükséges lenne megőrizni a modellredukció után.

A második lehetőség, hogy több mutatót őrizzünk meg, 2 faktor használata. 2 faktorra vonatkozó eredmények lentebb olvashatók.

Os<-U21_filtered[,c("O1","O2","O3","O4","O5","O6","O7","O8","O9")]
KMO(cor(Os))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(Os))
## Overall MSA =  0.73
## MSA for each item = 
##   O1   O2   O3   O4   O5   O6   O7   O8   O9 
## 0.40 0.73 0.77 0.79 0.63 0.83 0.74 0.85 0.64

fitO2 <- principal(Os,nfactors=2)
summary(fitO2)

## 
## Factor analysis with Call: principal(r = Os, nfactors = 2)
## 
## Test of the hypothesis that 2 factors are sufficient.
## The degrees of freedom for the model is 19  and the objective function was  1.87 
## The number of observations was  50  with Chi Square =  81.98  with prob <  8.4e-10 
## 
## The root mean square of the residuals (RMSA) is  0.11

fitO2$communality

##        O1        O2        O3        O4        O5        O6        O7        O8 
## 0.8438830 0.9223750 0.7955037 0.8210897 0.8665949 0.4008533 0.5506741 0.7255365 
##        O9 
## 0.3083455

fitO2$loadings

## 
## Loadings:
##    RC1    RC2   
## O1  0.138  0.908
## O2  0.959       
## O3  0.886  0.104
## O4  0.905       
## O5  0.670  0.646
## O6  0.629       
## O7  0.732  0.124
## O8  0.850       
## O9         0.547
## 
##                  RC1   RC2
## SS loadings    4.656 1.579
## Proportion Var 0.517 0.175
## Cumulative Var 0.517 0.693

Ha két faktort hozunk létre, akkor valamennyi kommunalitás magasabb, mint 0.25, azonban az O5 mutatóról nem lehet megállapítani, hogy melyik faktorhoz tartozik. Ezért a következő lépésben O5-t elhagyjuk, amely az adott ország 3 legjobb egyetemének Shanghai rangsorban szereplő pontértékét mutatja.

Os<-U21_filtered[,c("O1","O2","O3","O4","O6","O7","O8","O9")]
KMO(cor(Os))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(Os))
## Overall MSA =  0.78
## MSA for each item = 
##   O1   O2   O3   O4   O6   O7   O8   O9 
## 0.39 0.75 0.83 0.78 0.79 0.80 0.84 0.55

fitO2 <- principal(Os,nfactors=2)
summary(fitO2) 

## 
## Factor analysis with Call: principal(r = Os, nfactors = 2)
## 
## Test of the hypothesis that 2 factors are sufficient.
## The degrees of freedom for the model is 13  and the objective function was  1.21 
## The number of observations was  50  with Chi Square =  53.62  with prob <  7.1e-07 
## 
## The root mean square of the residuals (RMSA) is  0.1

fitO2$communality

##        O1        O2        O3        O4        O6        O7        O8        O9 
## 0.3429981 0.9340280 0.7946173 0.8450090 0.5310990 0.5954249 0.7515102 0.7198578

fitO2$loadings

## 
## Loadings:
##    RC1    RC2   
## O1  0.131  0.571
## O2  0.966       
## O3  0.891       
## O4  0.911 -0.123
## O6  0.615  0.390
## O7  0.716  0.288
## O8  0.848  0.180
## O9 -0.111  0.841
## 
##                  RC1   RC2
## SS loadings    4.195 1.319
## Proportion Var 0.524 0.165
## Cumulative Var 0.524 0.689

Ezáltal a KMO értéke nagyobb, mint 0.5, és minden mutató esetében a kommunalitás nagyobb, mint 0.25, azaz az eredeti változók megfelelően illeszkednek az eredményváltozóra. A fenti eredmények alapján a két faktor nem tartalmaz közös indikátorokat.

Az első faktorhoz tartoznak az alábbi indikátorok:

• O2: Népességre vetített publikációk száma
• O3: Hivatkozások átlagos száma a Karolinska Intézet normalizált impakt faktora alapján
• O4: Shanghai TOP 500-ban szereplő intézmények pontjaival súlyozott átlag/népesség
• O6: Felsőfokú továbbtanulási ráta
• O7: A 24-64 éves korosztályban a diplomások aránya
• O8: Kutatók száma/népesség

A második faktorhoz tartoznak az alábbi indikátorok:

• O1: A felsőoktatási intézmények által publikált cikkek száma
• O9: Diplomás munkanélküliek aránya a 25-64 éves korosztályban/diplomával nem, csak középfokú végzettséggel rendelkezők munkanélküliségi rátája

Így az O főkomponens létrehozásához képest a második faktorban mind az O1 mind pedig az O9 mutató megőrzésre került. Mivel főkomponens elemzést végeztünk R, C, E, O mutatókra, így az esetek akár korrelálhatnak is. Ezért a következő lépésben ellenőrizzük a közöttük levő korrelációt.

Korrelációs vizsgálat a főkomponensekre vonatkozóan

A korrelációs vizsgálat eredményei a 13. ábrán olvashatók. A 13. ábra alapján az R, O és C komponensek szoros szignifikáns kapcsolatban vannak.

factors<-cbind(fitR$scores,fitE$scores,fitC$scores,fitO$scores)
colnames(factors) <- c("R","E","C","O")
cor(factors)

##           R         E         C         O
## R 1.0000000 0.1659944 0.7354010 0.8207528
## E 0.1659944 1.0000000 0.2930795 0.2203930
## C 0.7354010 0.2930795 1.0000000 0.8482095
## O 0.8207528 0.2203930 0.8482095 1.0000000

corrgram(factors,order = TRUE,main="Komponensek közötti korrelációs kapcsolatok",lower.panel=panel.conf)

13. ábra A főkomponensek közötti korrelációs kapcsolatok

A korrelációs vizsgálatot elvégezzük az R, E, C főkomponensekre és az O első és második faktorára (O_fact1 és O_fact2) vonatkozóan. A korrelációs vizsgálat eredménye a 14. ábrán látható.

factors<-cbind(fitR$scores,fitE$scores,fitC$scores,fitO2$scores)
colnames(factors) <- c("R","E","C","O_fact1","O_fact2")
cor(factors)

##                  R         E          C      O_fact1       O_fact2
## R        1.0000000 0.1659944  0.7354010 8.350790e-01 -5.855570e-02
## E        0.1659944 1.0000000  0.2930795 2.054947e-01  1.568645e-01
## C        0.7354010 0.2930795  1.0000000 8.631607e-01 -1.107934e-01
## O_fact1  0.8350790 0.2054947  0.8631607 1.000000e+00  1.713747e-16
## O_fact2 -0.0585557 0.1568645 -0.1107934 1.713747e-16  1.000000e+00

corrgram(factors,order = TRUE,main="Komponensek és a faktorok közötti korrelációs kapcsolatok",lower.panel=panel.conf)

14. ábra A főkomponensek és a faktorok közötti korrelációs kapcsolatok

A 14. ábra alapján az R, O_fact1 és C komponensek szoros szignifikáns kapcsolatban vannak. Ugyanakkor az R, C, O_fact1 összevonása után 3 változócsoport azonosítható: az egyik, amelyet az erősen korreláló R, C, O_fact1 összevonásával kapunk, a másik az O_fact2, a harmadik az E.

Faktorelemzés

A főkomponens-elemzés után faktorelemzést végzünk a mutatókon. Első lépésként a faktorelemzés alkalmazhatóságát vizsgáljuk, amiben a könyökdiagram kirajzoltatása nyújt segítséget.

scree(U21_filtered)

15. ábra Könyökdiagram

KMO(cor(U21_filtered))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(U21_filtered))
## Overall MSA =  0.73
## MSA for each item = 
##   R1   R2   R3   R4   R5   E1   E2   E3   E4   C1   C2   C3   C4   C5   C6   O1 
## 0.70 0.51 0.82 0.75 0.79 0.56 0.35 0.69 0.82 0.76 0.63 0.37 0.83 0.83 0.76 0.36 
##   O2   O3   O4   O5   O6   O7   O8   O9 
## 0.80 0.80 0.79 0.75 0.73 0.66 0.74 0.30

Ahogy az a 15. ábrán látszódik, faktorelemzés esetén 3 független változócsoport használata célravezető. (Ugyanezt láttuk már a 12. ábrán is.)

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 

## 
## Factor analysis with Call: principal(r = U21s, nfactors = 3, rotate = "varimax")
## 
## Test of the hypothesis that 3 factors are sufficient.
## The degrees of freedom for the model is 207  and the objective function was  11.54 
## The number of observations was  50  with Chi Square =  440.6  with prob <  5.8e-19 
## 
## The root mean square of the residuals (RMSA) is  0.09

fit$communality

##        R1        R2        R3        R4        R5        E1        E2        E3 
## 0.5654944 0.5743389 0.8261170 0.7944043 0.8819861 0.3833381 0.1206436 0.5076963 
##        E4        C1        C2        C3        C4        C5        C6        O1 
## 0.5389281 0.3994398 0.7234222 0.2496209 0.5761145 0.6911836 0.6506483 0.5851133 
##        O2        O3        O4        O5        O6        O7        O8        O9 
## 0.9083897 0.8677446 0.8564676 0.7466957 0.3476795 0.5736386 0.7121123 0.3973431

fit$loadings

## 
## Loadings:
##    RC1    RC2    RC3   
## R1  0.397 -0.639       
## R2  0.439 -0.239 -0.570
## R3  0.906              
## R4  0.849 -0.265       
## R5  0.914 -0.211       
## E1  0.158  0.191  0.567
## E2         0.147  0.302
## E3  0.337  0.593  0.204
## E4  0.516  0.512  0.106
## C1  0.576         0.255
## C2  0.541 -0.159  0.637
## C3  0.154 -0.111  0.462
## C4  0.669  0.311  0.178
## C5  0.815 -0.153       
## C6  0.787  0.163       
## O1  0.227  0.495 -0.537
## O2  0.942         0.138
## O3  0.920         0.140
## O4  0.908         0.175
## O5  0.731  0.305 -0.345
## O6  0.526  0.208 -0.165
## O7  0.644  0.101 -0.386
## O8  0.836              
## O9         0.626       
## 
##                   RC1   RC2   RC3
## SS loadings    10.059 2.231 2.188
## Proportion Var  0.419 0.093 0.091
## Cumulative Var  0.419 0.512 0.603

biplot(fit)

16. ábra A három faktor jellemzői

A 16. ábra a három faktor jellemzőit mutatja. Ahogy az a fenti eredmények alapján is látszik, alacsony kommunalitású (pl. C3) és közös indikátorok (pl. R2) is vannak a kialakult faktoroknál. Ahhoz, hogy a faktorok egymáshoz viszonyított helyzetéről, illetve az egyes faktorokhoz tartozó mutatókról átfogó képet kapjunk, célszerű a faktorok (RC1, RC2, RC3) kezdeti állapotát 3D-ben kirajzoltatni. Ez a 3D-s ábra látható a 17. ábrán. Az ábra forgatásával láthatóvá válik, hogy az egyes mutatók hogyan csoportosulnak, mely faktorokhoz tartoznak. A fenti eredményekkel összhangban itt is látható, hogy az első faktorhoz tartozik például az R3, R4 és R5 mutató.

save <- options(rgl.useNULL = TRUE)
open3d()

## null 
##    1

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()

17. ábra A három létrehozott faktor kezdeti állapota 3D-ben

Az alacsony kommunalitású és a közös indikátorok eltávolítása után a 18. ábrán láthatóak a három faktor jellemzői.

columns<-c("R3","R4","R5","E1","E2","C1","C5","C6","O1","O4","O8","O9")
U21s<-U21_filtered[,columns]
KMO(cor(U21s))

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(U21s))
## Overall MSA =  0.81
## MSA for each item = 
##   R3   R4   R5   E1   E2   C1   C5   C6   O1   O4   O8   O9 
## 0.79 0.83 0.80 0.40 0.47 0.93 0.85 0.83 0.41 0.85 0.87 0.49

fit <- principal(U21s,nfactors=3,rotate = "varimax")
summary(fit) 

## 
## Factor analysis with Call: principal(r = U21s, nfactors = 3, rotate = "varimax")
## 
## Test of the hypothesis that 3 factors are sufficient.
## The degrees of freedom for the model is 33  and the objective function was  2.21 
## The number of observations was  50  with Chi Square =  93.17  with prob <  1.2e-07 
## 
## The root mean square of the residuals (RMSA) is  0.08

fit$communality

##        R3        R4        R5        E1        E2        C1        C5        C6 
## 0.8262154 0.8290226 0.9132030 0.6894734 0.5995707 0.4921809 0.6855794 0.6979789 
##        O1        O4        O8        O9 
## 0.6054390 0.8427117 0.7250600 0.6423689

fit$loadings

## 
## Loadings:
##    RC1    RC2    RC3   
## R3  0.888         0.190
## R4  0.905              
## R5  0.954              
## E1  0.134  0.815       
## E2 -0.126  0.764       
## C1  0.616  0.202 -0.269
## C5  0.824              
## C6  0.796 -0.161  0.197
## O1  0.148 -0.114  0.755
## O4  0.906  0.136       
## O8  0.829         0.169
## O9 -0.108  0.296  0.737
## 
##                  RC1   RC2   RC3
## SS loadings    5.787 1.450 1.312
## Proportion Var 0.482 0.121 0.109
## Cumulative Var 0.482 0.603 0.712

biplot(fit)

18. ábra A végső faktorstruktúra jellemzői

A 19. ábrán pedig a végső három faktor egymáshoz viszonyított állapota, illetve az egyes faktorokhoz tartozó mutatók láthatók 3D-ben.

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()

19. ábra A végső három faktor állapota 3D-ben

A fenti 3D-s ábra szemlélteti a három egymástól jól elkülöníthető mutatócsoportot.

A 20. ábrán látható a struktúramodell, amely grafikusan mutatja, hogy mely mutatók mely faktorhoz tartoznak. Az ábrán a nyilak vastagsága és erőssége azt szemlélteti, hogy az egyes mutatók mennyiben járulnak hozzá az adott faktorhoz.

Az első faktorhoz (f1) tartoznak az alábbi mutatók:

• R3: Éves összes kiadás dollárban/hallgató (nappali tagozatos egyenértékben számolva)
• R4: Felsőoktatás K+F kiadásai a GDP százalékában
• R5: Felsőoktatás K+F kiadásai/népesség
• C1: Nemzetközi hallgatók aránya
• C5: Tudástranszfer az egyetemek és az üzleti vállalkozások között
• C6: Egyetemi és ipari szereplők közösen publikált cikkeinek aránya
• O4: Shanghaj TOP 500-ban szereplő intézmények pontjaival súlyozott átlag/népesség
• O8: Kutatók száma/népesség

A második faktorhoz (f2) tartoznak az alábbi indikátorok:

• E1: Női hallgatók aránya
• E2: Női oktatók aránya a felsőoktatásban

A harmadik faktorhoz (f3) tartoznak az alábbi mutatók:

• O1: A felsőoktatási intézmények által publikált cikkek száma
• O9: Diplomás munkanélküliek aránya a 25-64 éves korosztályban / diplomával nem, csak középfokú végzettséggel rendelkezők munkanélküliségi rátája

Az első faktor (f1) egyaránt tartalmaz kiadás-jellegű, kapcsolatokra vonatkozó, illetve output-jellegű mutatókat. A második faktor (f2) a nők szerepét mutatja a felsőoktatásban, a harmadik faktor (f3) pedig outputra vonatkozó indikátorokat tartalmaz.

model='f1=~R3+R4+R5+C1+C5+C6+O4+O8
       f2=~E1+E2
       f3=~O1+O9
       f1~~0*f2 + 0*f3
       f2~~0*f3'
fitsem <- cfa(model,U21_filtered,std.lv = TRUE)
semPaths(fitsem,residuals=FALSE,"std",edge.color = "black",color="yellow")

20. ábra A végső struktúramodell

Dimenzió-redukció összegzése: Amennyiben cél, hogy független változócsoportokat alakítsunk ki, akkor mind a főkomponens, mind a faktorelemzéssel 3 változócsoportot definiálhatunk. Ezek közül a két “tiszta” változócsoport csak 2-2 mutatót tartalmaz, míg a harmadik vegyes csoport R, C, O mutatókból alakít ki egy csoportot. Komoly probléma, hogy a 24 indikátorból csak 12 őrződik meg.

#4. Adatredukció

##A faktorok 3 klaszterbe rendezése A következőkben a fentebb kapott három faktort klaszterekbe rendezzük. Az előzetes klaszterezettség vizsgálat alapján 3-4 klaszter volt valószínűsíthető. Így első lépésként három klasztert hozunk létre.

set.seed(2)
km3 <- kmeans(fit$scores,3)$cluster
cluster.stats(dist(fit$scores),km3)

## $n
## [1] 50
## 
## $cluster.number
## [1] 3
## 
## $cluster.size
## [1] 18 20 12
## 
## $min.cluster.size
## [1] 12
## 
## $noisen
## [1] 0
## 
## $diameter
## [1] 2.180009 5.704993 4.494681
## 
## $average.distance
## [1] 1.146773 1.610924 1.901035
## 
## $median.distance
## [1] 1.108597 1.356117 1.479310
## 
## $separation
## [1] 0.7804456 0.7390499 0.7390499
## 
## $average.toother
## [1] 2.479894 2.470991 2.636661
## 
## $separation.matrix
##           [,1]      [,2]      [,3]
## [1,] 0.0000000 0.7804456 0.8288624
## [2,] 0.7804456 0.0000000 0.7390499
## [3,] 0.8288624 0.7390499 0.0000000
## 
## $ave.between.matrix
##          [,1]     [,2]     [,3]
## [1,] 0.000000 2.373189 2.657735
## [2,] 2.373189 0.000000 2.617694
## [3,] 2.657735 2.617694 0.000000
## 
## $average.between
## [1] 2.520423
## 
## $average.within
## [1] 1.513456
## 
## $n.between
## [1] 816
## 
## $n.within
## [1] 409
## 
## $max.diameter
## [1] 5.704993
## 
## $min.separation
## [1] 0.7390499
## 
## $within.cluster.ss
## [1] 77.04899
## 
## $clus.avg.silwidths
##         1         2         3 
## 0.4922942 0.3117005 0.2117167 
## 
## $avg.silwidth
## [1] 0.3527181
## 
## $g2
## NULL
## 
## $g3
## NULL
## 
## $pearsongamma
## [1] 0.4333483
## 
## $dunn
## [1] 0.1295444
## 
## $dunn2
## [1] 1.248367
## 
## $entropy
## [1] 1.076819
## 
## $wb.ratio
## [1] 0.6004771
## 
## $ch
## [1] 21.33511
## 
## $cwidegap
## [1] 0.9296514 2.7278667 2.2474481
## 
## $widestgap
## [1] 2.727867
## 
## $sindex
## [1] 0.764972
## 
## $corrected.rand
## NULL
## 
## $vi
## NULL

Az alábbi forgatható ábrán (21. ábra) három, egymástól jól elkülönülő klasztert láthatunk. Az egyes klaszterekbe tartozó országok a következők:

Első klaszter	Második klaszter	Harmadik klaszter
Ausztrália	Argentína	Chile
Ausztria	Brazília	Görögország
Belgium	Bulgária	India
Kanada	Kína	Indonézia
Dánia	Horvátország	Irán
Finnország	Csehország	Olaszország
Franciaország	Németország	Japán
Hong Kong	Magyarország	Dél-Korea
Írország	Malajzia	Mexikó
Izrael	Lengyelország	Szaúd Arábia
Hollandia	Románia	Tajvan
Új-Zéland	Oroszország	Törökország
Norvégia	Szerbia
Portugália	Szlovákia
Szingapúr	Szlovénia
Svédország	Dél-Afrika
Svájc	Spanyolország
Egyesült Királyság	Thaiföld
	Ukrajna
	Egyesült Államok

Az első klaszterbe 18, a másodikba 20, a harmadikba pedig 12 ország került.

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()

21. ábra A létrehozott három klaszter

A 22. ábrán térképen ábrázolva láthatók, hogy mely országok kerültek egy klaszterbe. A térképen pirossal vannak jelölve azok az országok, amelyek az 1. klaszterbe, kékkel, amelyek a 2. klaszterbe, zölddel pedig azok, amelyek a 3. klaszterbe tartoznak.

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])

22. ábra A három klaszter térképen ábrázolva

##A faktorok 4 klaszterbe rendezése

A következő lépésként a faktorokat 4 klaszterbe rendezzük.

set.seed(2)
km4 <- kmeans(fit$scores,4)$cluster
cluster.stats(dist(fit$scores),km4)

## $n
## [1] 50
## 
## $cluster.number
## [1] 4
## 
## $cluster.size
## [1] 22  4 17  7
## 
## $min.cluster.size
## [1] 4
## 
## $noisen
## [1] 0
## 
## $diameter
## [1] 2.581561 3.503032 2.180009 3.714822
## 
## $average.distance
## [1] 1.228931 2.257800 1.134808 1.939406
## 
## $median.distance
## [1] 1.237384 2.276079 1.096520 1.818473
## 
## $separation
## [1] 0.3953279 0.6144058 0.3953279 0.7669507
## 
## $average.toother
## [1] 2.389488 3.094787 2.484270 3.014312
## 
## $separation.matrix
##           [,1]      [,2]      [,3]      [,4]
## [1,] 0.0000000 0.6144058 0.3953279 0.7669507
## [2,] 0.6144058 0.0000000 1.1363751 2.1025880
## [3,] 0.3953279 1.1363751 0.0000000 1.3234252
## [4,] 0.7669507 2.1025880 1.3234252 0.0000000
## 
## $ave.between.matrix
##          [,1]     [,2]     [,3]     [,4]
## [1,] 0.000000 2.807635 2.142717 2.749849
## [2,] 2.807635 0.000000 3.167920 3.819656
## [3,] 2.142717 3.167920 0.000000 3.167064
## [4,] 2.749849 3.819656 3.167064 0.000000
## 
## $average.between
## [1] 2.612725
## 
## $average.within
## [1] 1.378705
## 
## $n.between
## [1] 831
## 
## $n.within
## [1] 394
## 
## $max.diameter
## [1] 3.714822
## 
## $min.separation
## [1] 0.3953279
## 
## $within.cluster.ss
## [1] 53.58231
## 
## $clus.avg.silwidths
##         1         2         3         4 
## 0.3877404 0.1005759 0.4489660 0.2472150 
## 
## $avg.silwidth
## [1] 0.3659104
## 
## $g2
## NULL
## 
## $g3
## NULL
## 
## $pearsongamma
## [1] 0.5644199
## 
## $dunn
## [1] 0.1064191
## 
## $dunn2
## [1] 0.9490284
## 
## $entropy
## [1] 1.205341
## 
## $wb.ratio
## [1] 0.5276886
## 
## $ch
## [1] 26.73279
## 
## $cwidegap
## [1] 0.8288624 2.7278667 0.9296514 2.2474481
## 
## $widestgap
## [1] 2.727867
## 
## $sindex
## [1] 0.5028333
## 
## $corrected.rand
## NULL
## 
## $vi
## NULL

Ahogy az a 23. ábrán látható, egymástól jól elkülönülő 4 klasztert kaptunk. Az első klaszterbe 22, a másodikba 4, a harmadikba 17, a negyedikbe pedig 7 ország került, amelyek az alábbiak:

Első klaszter	Második klaszter	Harmadik klaszter	Negyedik klaszter
Argentína	Kína	Ausztrália	India
Brazília	Németország	Ausztria	Indonézia
Bulgária	Magyarország	Belgium	Irán
Chile	Egyesült Államok	Kanada	Japán
Horvátország		Dánia	Dél-Korea
Csehország		Finnország	Tajvan
Görögország		Franciaország	Törökország
Olaszország		Hong Kong
Malajzia		Írország
Mexikó		Izrael
Lengyelország		Hollandia
Portugália		Új-Zéland
Románia		Norvégia
Oroszország		Szingapúr
Szaúd Arábia		Svédország
Szerbia		Svájc
Szlovákia		Egyesült Királyság
Szlovénia
Dél-Afrika
Spanyolország
Thaiföld
Ukrajna

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()

23. ábra A létrehozott négy klaszter

A 24. ábrán térképen ábrázolva látható a 4 klaszter.

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])

24. ábra A négy klaszter térképen ábrázolva

5. Parciális rangsorok meghatározása

A parciális (országklasztereken belüli) rangsorok meghatározásánál egyrészt felhasználjuk a kiszámított faktorsúlyokat, másrészt pedig az U21 által használt eredeti súlyokat. A következő lépésben meghatározzuk az U21 által használt eredeti súlyokat.

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")

A következő lépésben kiszámításra kerülnek a parciális rangsorok, amelyek az 1. táblázatban láthatók. A táblázat utolsó két oszlopában jelzett km3, illetve km4 jelzi azt, hogy az egyes országok 3, illetve 4 klaszter esetén melyik klaszterbe kerültek. Az U21_overall_Rank az eredeti U21 rangsor rangsorszámait tartalmazza, az FA_Rank a faktorscore-ok alapján képzett rangsort, az U21_Rank_selected pedig azokat a helyezéseket, amelyek már csak a faktoroknak megfelelő indikátorokat tartalmazzák.

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="**1. táblázat** Parciális rangsor")

1. táblázat Parciális rangsor

	ID	U21_overall_Rank	FA_weights	FA_Rank	Score_selected	U21_Rank_selected	km3	km4
Argentina	1	41	3.716323e-02	39	53.20702	16	2	1
Australia	2	9	1.226081e-01	20	46.85638	23	1	3
Austria	3	12	1.205554e-01	22	77.76420	6	1	3
Belgium	4	13	2.847575e-01	6	41.03496	27	1	3
Brazil	5	38	1.626992e-01	14	41.07081	26	2	1
Bulgaria	6	40	5.989465e-02	30	22.48808	47	2	1
Canada	7	4	2.566109e-01	7	77.78374	5	1	3
Chile	8	33	1.120821e-02	43	20.57509	48	3	1
China	9	34	6.702021e-02	28	39.96442	29	2	2
Croatia	10	44	6.752780e-02	27	26.79678	42	2	1
Czech Republic	11	26	1.535241e-01	15	61.93087	9	2	1
Denmark	12	3	1.276323e-01	17	56.33062	12	1	3
Finland	13	5	4.592789e-01	2	95.87306	2	1	3
France	14	18	2.193017e-01	9	40.09535	28	1	3
Germany	15	14	4.278289e-01	3	73.57370	7	2	2
Greece	16	32	3.890460e-02	38	28.52945	41	3	1
Hong Kong	17	15	1.276280e-01	18	50.23970	21	1	3
Hungary	18	29	1.664197e-01	12	31.49285	38	2	2
India	19	50	1.537641e-03	47	37.16362	31	3	4
Indonesia	20	48	1.284827e-03	48	17.93802	50	3	4
Iran	21	49	4.994809e-03	45	35.34617	33	3	4
Ireland	22	17	2.919132e-01	5	41.66170	25	1	3
Israel	23	19	1.652540e-01	13	54.70782	13	1	3
Italy	24	27	4.908183e-02	34	32.20203	36	3	1
Japan	25	20	1.944635e-04	49	61.59416	10	3	4
Korea, Rep. (South)	26	21	1.681845e-05	50	31.96663	37	3	4
Malaysia	27	28	5.352892e-02	31	50.55125	20	2	1
Mexico	28	46	9.292853e-03	44	20.50180	49	3	1
Netherlands	29	7	1.669516e-01	11	68.99285	8	1	3
New Zealand	30	16	7.233421e-02	26	39.72090	30	1	3
Norway	31	11	2.270396e-01	8	78.31016	4	1	3
Poland	32	31	8.456104e-02	24	25.39432	43	2	1
Portugal	33	24	1.182719e-01	23	52.54867	18	1	1
Romania	34	39	4.490066e-02	35	24.88180	44	2	1
Russia	35	36	7.997149e-02	25	44.89472	24	2	1
Saudi Arabia	36	30	2.028512e-02	41	32.57981	35	3	1
Serbia	37	34	5.163543e-02	33	53.96639	14	2	1
Singapore	38	10	5.246536e-02	32	48.01649	22	1	3
Slovakia	39	37	1.219388e-01	21	50.87108	19	2	1
Slovenia	40	25	1.228653e-01	19	30.92769	39	2	1
South Africa	41	45	4.384428e-02	36	29.91428	40	2	1
Spain	42	23	1.420620e-01	16	32.97279	34	2	1
Sweden	43	2	2.060980e-01	10	100.00000	1	1	3
Switzerland	44	6	6.600255e-02	29	52.70262	17	1	3
Taiwan	45	22	3.518628e-02	40	58.63325	11	3	4
Thailand	46	42	4.127924e-02	37	24.36351	45	2	1
Turkey	47	47	3.550801e-03	46	35.56118	32	3	4
Ukraine	48	43	1.884380e-02	42	23.44289	46	2	1
United Kingdom	49	8	3.468806e-01	4	80.62078	3	1	3
United States	50	1	4.693511e-01	1	53.27870	15	2	2

A következő lépésben a rangsorok közötti rangkorrelációk kerülnek kiszámításra.

Rangsorok közötti korreláció 3 klaszter esetén

A 25-28. ábrák szemléltetik a rangsorok között fennálló rangkorrelációk mértékét. A 25. ábra az összes országra, a 26-28. ábra külön-külön az 1-3. klaszter országaira vonatkozóan.

Ha az összes országot vizsgáljuk (25. ábra), az eredeti U21 rangsor rangsorszámai (U21_overall_Rank) és az általunk kiszámított azon rangsor között, amely már csak az egyes faktoroknak megfelelő indikátorokat tartalmazzák (U21_Rank_selected) szoros pozitív kapcsolat figyelhető meg. Ez azt jelenti, hogy a két rangsorban a rangszámok hasonlóak egymáshoz. Ha csak az 1. és a 2. klaszterbe tartozó országok rangszámait vizsgáljuk, hasonló eredményt kapunk (lásd 26-27. ábra).

par(mfrow=c(2,2))
RANKS<-c("U21_overall_Rank","FA_Rank","U21_Rank_selected")
corrgram(base[,RANKS],cor.method = "spearman",order = TRUE,main="Összes országra",lower.panel=panel.cor)

25. ábra A parciális rangsorok közötti rangkorrelációk három klaszter esetében - az összes országra

corrgram(base[km3==1,RANKS],cor.method = "spearman",order = TRUE,main=paste("1. klaszter országaira, ",sum(km3==1)," ország"),lower.panel=panel.cor)

26. ábra A parciális rangsorok közötti rangkorrelációk három klaszter esetében - az 1. klaszter országaira

corrgram(base[km3==2,RANKS],cor.method = "spearman",order = TRUE,main=paste("2. klaszter országaira, ",sum(km3==2)," ország"),lower.panel=panel.cor)

27. ábra A parciális rangsorok közötti rangkorrelációk három klaszter esetében - a 2. klaszter országaira

corrgram(base[km3==3,RANKS],cor.method = "spearman",order = TRUE,main=paste("3. klaszter országaira, ",sum(km3==3)," ország"),lower.panel=panel.cor)

28. ábra A parciális rangsorok közötti rangkorrelációk három klaszter esetében - a 3. klaszter országaira

Rangsorok közötti korrelálció 4 klaszter esetén

A rangkorrelációk kiszámításra kerültek abban az esetben is, ha négy klasztert hozunk létre. Az 1.,2. és 3. klaszter esetében közepesen szoros kapcsolat van az eredeti U21-es rangszámok, és az általunk kiszámításra került rangszámok között.

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("1. klaszter országaira, ",sum(km4==1)," ország"),lower.panel=panel.cor)

29. ábra A parciális rangsorok közötti rangkorrelációk négy klaszter esetében - az 1. klaszter országaira

corrgram(base[km4==2,RANKS],cor.method = "spearman",order = TRUE,main=paste("2. klaszter országaira, ",sum(km4==2)," ország"),lower.panel=panel.cor)

30. ábra A parciális rangsorok közötti rangkorrelációk négy klaszter esetében - a 2. klaszter országaira

corrgram(base[km4==3,RANKS],cor.method = "spearman",order = TRUE,main=paste("3. klaszter országaira, ",sum(km4==3)," ország"),lower.panel=panel.cor)

31. ábra A parciális rangsorok közötti rangkorrelációk négy klaszter esetében - a 3. klaszter országaira

corrgram(base[km4==4,RANKS],cor.method = "spearman",order = TRUE,main=paste("4. klaszter országaira, ",sum(km4==4)," ország"),lower.panel=panel.cor)

32. ábra A parciális rangsorok közötti rangkorrelációk négy klaszter esetében - a 4. klaszter országaira

6. Összegzés

A következőkben összefoglaljuk a kétlépéses klaszterezés eredményeit.

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))

Az alábbi, 33. ábrán az eredeti U21-es adatok hőtérképe látható. A cellák színezése az egyes mutatókhoz tartozó értékektől függ. Minél sötétebb vörös az adott cella, annál kisebb az adott mutató értéke, minél pedig sötétebb kék a cella, annál magasabb. Az ábrán az országok a 3 klaszter szerint vannak rendezve.

p <- plot_ly(x=colnames(ORDORIG3), y=rownames(ORDORIG3), z = as.matrix(ORDORIG3), type = "heatmap",colorscale="RdBu",reversescale=TRUE)
p

33. ábra Az eredeti U21 adatok hőtérképe - a 3 klaszter szerint rendezve

A 34. ábrán a létrehozott három klaszterbe tartozó országok láthatók, illetve azok az indikátorok, amelyek a három főkomponenshez tartoznak. Azok az indikátorok amelyek egyetlen főkomponenshez sem tartoznak, vörös színű cellával jelennek meg. Látható, hogy az indikátorok fele (12) kiesett az elemzés során.

p <- plot_ly(x=colnames(ORDERED3), y=rownames(ORDERED3), z = as.matrix(ORDERED3), type = "heatmap",colorscale="RdBu",reversescale=TRUE)
p

34. ábra A létrehozott három klaszter

p <- plot_ly(x=colnames(ORDORIG4), y=rownames(ORDORIG4), z = as.matrix(ORDORIG4), type = "heatmap",colorscale="RdBu",reversescale=TRUE)
p

35. ábra Az eredeti U21 adatok hőtérképe - a 4 klaszter szerint rendezve

A 35. ábra - a 33. ábrához hasonlóan - mutatja az országok eredeti U21 adatainak hőtérképet, azzal a különbséggel, hogy a függőleges tengelyen felsorolt országok most nem a 3, hanem a 4 klaszternek megfelelően csoportosítottak.

A 36. ábrán a négy létrehozott klaszter látható, illetve azok az indikátorok, amelyek a három főkomponenshez tartoznak.

p <- plot_ly(x=colnames(ORDERED4), y=rownames(ORDERED4), z = as.matrix(ORDERED4), type = "heatmap",colorscale="RdBu",reversescale=TRUE)
p

36. ábra A létrehozott négy klaszter

Összegzés: A kétlépéses klaszterezés ilyen kevés adat esetén kétséges. Maga a módszer nem engedi meg a klaszterek átlapolódását, ugyanakkor lehet, hogy adott ország valamely tényezőiben jól, míg másokban rosszul teljesít. A klaszterezés eredményeképpen nagyon kis klasztereket is kaptunk, amelyeket nehéz elemezni, értelmezni.

A tanulmányhoz csatolt kiegészítő fájlok (jelen html is) tudományos célra szabadon felhasználhatók. Kérjük az ezekre való hivatkozáskor magát a Közgazdasági Szemlében megjelent tanulmányt hivatkozni.

Az eredeti cikk

Kosztyán, Z. T., Banász, Z., Csányi, V. V., & Telcs, A. (2019). Felsőoktatási ligák, parciális rangsorok képzése biklaszterezési eljárásokkal. Közgazdasági Szemle, 66(9), 905-931.

További hivatkozások

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