clustering_project

Schools in Poland - Clustering project

The aim of this project is to perform kmeans clustering on dataset that contains information on average score of 8th grade exam in Poland together with some variables like school financing, number of teachers in each county (powiat), etc. The result is a cartogam with counties that are similar to each other. The data set was prepared by me using data from OKE (Okręgowe Komisje Egzaminacyjne) and GUS.

Data manipulation and overview

Setting seed and downloading necessary packages.

set.seed(1210)

#install.packages("factoextra")
#install.packages("flexclust")
#install.packages("fpc")
#install.packages("clustertend")
#install.packages("cluster")
#install.packages("ClusterR")
library(factoextra)
Warning: package 'factoextra' was built under R version 4.3.2
Loading required package: ggplot2
Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(flexclust)
Warning: package 'flexclust' was built under R version 4.3.2
Loading required package: grid
Loading required package: lattice
Loading required package: modeltools
Loading required package: stats4
library(fpc)
Warning: package 'fpc' was built under R version 4.3.2
library(clustertend)
Package `clustertend` is deprecated.  Use package `hopkins` instead.
library(cluster)
Warning: package 'cluster' was built under R version 4.3.2
library(ClusterR)
Warning: package 'ClusterR' was built under R version 4.3.2

Now, we need to load the data set and check the structure. X1 variable is the amount of books per school in a county, X2 is the expenditure per school in a county, X3 is the number of full-time teacher positions in a school and X4 is the amount of students in a class. O1 is the average exam score of 8th graders in Poland. The variables x1, x2, x3 and x4 are for the year 2021 and the exam score is from 2022 as, in literature it is assumed, that inputs from year t have an influence on output of year t+1.

primary2022E8 <- read.csv("C:/Users/mateu/Downloads/primary2022.csv")
primary2022<-primary2022E8[1:380, 1:6]

variables <- primary2022[2:6]

str(variables)
'data.frame':   380 obs. of  5 variables:
 $ x1_ksiegozbior_na_szkole: num  9700 12060 7756 9484 10469 ...
 $ x2_wydatki_na_szkole    : num  2394259 2542486 4114696 2428588 2404206 ...
 $ x3_etaty_na_szkole      : num  15.8 17.1 23.6 15.1 14.7 ...
 $ x4_uczniowie_na_odział  : num  17 17 17 16 15 16 15 15 16 17 ...
 $ o1_wyniki_e8            : num  53.3 56.7 56 53.7 55.3 ...

Kmeans clustering

Dataset has been checked for NAs beforehand. Applying standarization for clustering.

variables_z <- as.data.frame(lapply(variables, scale))

head(variables_z)
  x1_ksiegozbior_na_szkole x2_wydatki_na_szkole x3_etaty_na_szkole
1                0.1937364           -0.3698476        -0.37244964
2                0.8296968           -0.1497750        -0.09514117
3               -0.3300446            2.1844746         1.20954649
4                0.1356062           -0.3188791        -0.51248762
5                0.4011379           -0.3550783        -0.59528296
6                0.3503851           -0.6827499        -0.88925925
  x4_uczniowie_na_odział o1_wyniki_e8
1             0.45837234   -0.6810717
2             0.45837234   -0.1066095
3             0.45837234   -0.2188132
4            -0.08765114   -0.5998561
5            -0.63367461   -0.3302558
6            -0.08765114   -0.4586431

The dataset is not big, nor very complicated. Kmeans algorithim is sufficient.

primary_clusters <-kmeans(variables_z, 2)

primary_clusters$size
[1] 275 105

Creating the clusters.

fviz_cluster(list(data=variables_z, cluster=
              primary_clusters$cluster), 
             ellipse.type = "convex", geom="point",
             stand=FALSE, pallette='jco', ggtheme=theme_classic())

Variance is explained in 80% which is sufficient.

Checking the silhouette plot.

primary_cluster_test <- kmeans(variables_z, 2)
sil<-silhouette(primary_cluster_test$cluster, dist(variables_z))
fviz_silhouette(sil)
  cluster size ave.sil.width
1       1  105          0.28
2       2  275          0.53

The silhouette plot looks good. There are many observations close to 0.7 and little below 0.

opt2<-Optimal_Clusters_KMeans(variables_z, max_clusters=10, plot_clusters = TRUE)

I have decided to stick to 2 clusters.

Cartogram creation

Some data manipulation has been already done in excel. But we still need to make e correction to data. Initial territorial values did not have 0 in front if the length was less equal to 6.

primary2022$TERYT <- ifelse(nchar(primary2022$TERYT)==6,  paste0("0", primary2022$TERYT), primary2022$TERYT)

Now, I am making a data set containing only territorial values and clusters. It will be easier to plot it on cartogram this way.

mergedterytcluster<- cbind(primary2022$TERYT, primary_clusters$cluster)
colnames(mergedterytcluster)<-c("JPT_KOD_JE", "cluster")
mergedterytcluster <- as.data.frame(mergedterytcluster)
library(sf)
Warning: package 'sf' was built under R version 4.3.2
Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
shp<-st_read("C:/Users/mateu/Downloads/powiaty/powiaty.shp")
Reading layer `powiaty' from data source 
  `C:\Users\mateu\Downloads\powiaty\powiaty.shp' using driver `ESRI Shapefile'
Simple feature collection with 380 features and 36 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 171677.5 ymin: 133223.4 xmax: 861895.8 ymax: 775019.1
Projected CRS: ETRF2000-PL / CS92
shpdf<-as.data.frame(shp)

Now we need to truncate last two zeros from TERYT column to match it with JPT_KOD_JE

mergedterytcluster$JPT_KOD_JE<-substr(mergedterytcluster$JPT_KOD_JE, 1, 4)
shp_merged <- merge(shpdf, mergedterytcluster, by="JPT_KOD_JE")
str(shp_merged)
'data.frame':   380 obs. of  38 variables:
 $ JPT_KOD_JE: chr  "0201" "0202" "0203" "0204" ...
 $ gml_id    : chr  NA NA NA NA ...
 $ JPT_SJR_KO: chr  "POW" "POW" "POW" "POW" ...
 $ JPT_POWIER: chr  "130386" "47851" "44326" "73816" ...
 $ JPT_NAZWA_: chr  "powiat bolesławiecki" "powiat dzierżoniowski" "powiat głogowski" "powiat górowski" ...
 $ JPT_ORGAN_: chr  NA NA NA NA ...
 $ JPT_JOR_ID: chr  "13410" "13410" "13410" "13410" ...
 $ WERSJA_OD : chr  "20211210" "20120926" "20120926" "20120926" ...
 $ WERSJA_DO : chr  "0" "0" "0" "0" ...
 $ WAZNY_OD  : chr  "20120926" "20120926" "20120926" "20120926" ...
 $ WAZNY_DO  : chr  "0" "0" "0" "0" ...
 $ JPT_KOD__1: chr  NA NA NA NA ...
 $ JPT_NAZWA1: chr  NA NA NA NA ...
 $ JPT_ORGAN1: chr  "NZN" "NZN" "NZN" "NZN" ...
 $ JPT_WAZNA_: chr  "NZN" "NZN" "NZN" "NZN" ...
 $ ID_BUFORA_: chr  "250003" "13888" "13888" "13888" ...
 $ ID_BUFORA1: chr  "0" "0" "0" "0" ...
 $ ID_TECHNIC: chr  "829341" "829318" "829319" "829320" ...
 $ IIP_PRZEST: chr  "PL.PZGIK.200" "PL.PZGIK.200" "PL.PZGIK.200" "PL.PZGIK.200" ...
 $ IIP_IDENTY: chr  "275bb110-6b8f-4f24-873d-abbb862c02a6" "6e4b2d8b-4559-4329-a4c9-eefdb5c5eb89" "f6753ea2-1749-4460-9a52-ae8524c5b67b" "65b5a197-3d86-401c-8cea-ea80c3b4bf7c" ...
 $ IIP_WERSJA: chr  "2021-12-10T14:23:10+01:00" "2012-09-27T07:36:28+02:00" "2012-09-27T07:36:28+02:00" "2012-09-27T07:36:28+02:00" ...
 $ JPT_KJ_IIP: chr  "EGIB" "EGIB" "EGIB" "EGIB" ...
 $ JPT_KJ_I_1: chr  "0201" "0202" "0203" "0204" ...
 $ JPT_KJ_I_2: chr  NA NA NA NA ...
 $ JPT_OPIS  : chr  NA NA NA NA ...
 $ JPT_SPS_KO: chr  "UZG" "UZG" "UZG" "UZG" ...
 $ ID_BUFOR_1: chr  "0" "0" "0" "0" ...
 $ JPT_ID    : chr  "1361679" "829318" "829319" "829320" ...
 $ JPT_POWI_1: chr  "0" "0" "0" "0" ...
 $ JPT_KJ_I_3: chr  NA NA NA NA ...
 $ JPT_GEOMET: chr  "0" "0" "0" "0" ...
 $ JPT_GEOM_1: chr  "0" "0" "0" "0" ...
 $ SHAPE_LENG: chr  "2.85878211456" "1.5716746878" "2.08474100874" "2.02724245867" ...
 $ SHAPE_AREA: chr  ".16812415183" ".0609316661" ".05759524697" ".0958298518" ...
 $ REGON     : chr  "23082162900000" "89071781700000" "39064721600000" "41105046800000" ...
 $ RODZAJ    : chr  "powiat" "powiat" "powiat" "powiat" ...
 $ geometry  :sfc_MULTIPOLYGON of length 380; first list element: List of 1
  ..$ :List of 1
  .. ..$ : num [1:4015, 1:2] 250582 250592 250636 250669 250695 ...
  ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
 $ cluster   : chr  "1" "1" "2" "1" ...
shp_merged$cluster <- as.numeric(shp_merged$cluster)

Final result

library(RColorBrewer)
library(ggplot2)
ggplot(shp_merged$geometry)+
  geom_sf(aes(fill=as.factor(shp_merged$cluster))) + scale_fill_brewer(name="clusters", type = "seq", palette = "Dark2")

As we can see the city counties and counties located nearby have been classified as 1 cluster. The second one are more rural areas. We can see there is a significant difference when it comes to education system in Poland.

mergeddfcluster<- cbind(primary2022, primary_clusters$cluster)
colnames(mergeddfcluster)<-c("TERYT", "x1_ksiegozbior_na_szkole", "x2_wydatki_na_szkole", "x3_etaty_na_szkole","x4_uczniowie_na_odział", "o1_wyniki_e8", "cluster")

tapply(mergeddfcluster[2:6], mergeddfcluster$cluster, summary)
$`1`
 x1_ksiegozbior_na_szkole x2_wydatki_na_szkole x3_etaty_na_szkole
 Min.   : 3238            Min.   :1488178      Min.   :10.17     
 1st Qu.: 6244            1st Qu.:2069560      1st Qu.:13.57     
 Median : 7469            Median :2284001      Median :14.72     
 Mean   : 7795            Mean   :2321420      Mean   :15.05     
 3rd Qu.: 8859            3rd Qu.:2546151      3rd Qu.:16.47     
 Max.   :14770            Max.   :3537927      Max.   :21.98     
 x4_uczniowie_na_odział  o1_wyniki_e8  
 Min.   :11.00          Min.   :43.78  
 1st Qu.:14.00          1st Qu.:52.09  
 Median :15.00          Median :54.97  
 Mean   :15.36          Mean   :55.40  
 3rd Qu.:16.00          3rd Qu.:58.30  
 Max.   :18.00          Max.   :74.78  

$`2`
 x1_ksiegozbior_na_szkole x2_wydatki_na_szkole x3_etaty_na_szkole
 Min.   : 5320            Min.   :2418552      Min.   :17.95     
 1st Qu.: 8026            1st Qu.:3067690      1st Qu.:20.97     
 Median :11252            Median :3403088      Median :23.65     
 Mean   :12086            Mean   :3486557      Mean   :24.34     
 3rd Qu.:15237            3rd Qu.:3898023      3rd Qu.:26.64     
 Max.   :34157            Max.   :5622205      Max.   :36.55     
 x4_uczniowie_na_odział  o1_wyniki_e8  
 Min.   :15.00          Min.   :49.87  
 1st Qu.:18.00          1st Qu.:58.07  
 Median :18.00          Median :62.21  
 Mean   :18.26          Mean   :62.37  
 3rd Qu.:19.00          3rd Qu.:66.59  
 Max.   :20.00          Max.   :77.97  

As we can see from descriptive statistics counties classified as 2, which are the city counties or located nearby cities have, on average, all variables higher. It means that those counties have higher expenditures per school, more books in schools, have more teachers, classes are more populated and they have higher average exam scores.