# loading all necessary libraries

library(readxl)
library(corrplot)
library(NbClust)
library(tidyverse)
library(cluster)    
library(factoextra)
library(dendextend)
library(stringr)
library(stats)
library(flexclust)
library(fpc)
library(clustertend)
library(ClusterR)
library(lubridate)
library(kableExtra)
library(gridExtra)
library(gtable)

Introduction

Global climate warming is one of the greatest challenges facing humanity today. It is caused primarily by human activity, and more specifically by the ever-increasing emissions of carbon dioxide (CO2) into the atmosphere. More and more governments of the world’s major economies are beginning to set goals to reduce global CO2 emissions. In order to reduce global CO2 emissions, it is also worth asking ourselves what characterizes the countries that emit the most CO2 and those that produce the least, and exactly what factors affect CO2 emissions the most.

In this study, a cluster analysis will be conducted to group countries with similar macroeconomic and non-macroeconomic characteristics. Its results should allow us to assess what common characteristics countries with similar CO2 emissions have, so that it will be possible to determine which factors most affect countries’ CO2 emissions. In addition, there will be an additional analysis examining the impact of non-macroeconomic factors (such as the Democracy Index) on CO2 emissions and on cluster formation. That is, whether the information on social factors allows us to draw more interesting conclusions from cluster analysis than from purely macroeconomic factors alone.

Review of the data

The data used in the paper comes from the World Bank’s World Development Indicators database and from The Economist Intelligence Unit’s website for the democracy index data. They cover 130 countries from around the world and are for the year 2014. The choice of 2014 is due to the lack of data on relevant indicators for later years in the World Bank database. All variables are described below and basic statistics for each variable are also provided.

# loading the dataset as data

data <- read_excel("CO2.xlsx")
options(scipen=999)

kable(head(data))
Country CO2 GDPpc GDPg Electric_power_cons Coal Gas Oil Resources_rents Urbanization Democracy_index Hydro_renewable_nucl Agriculture EU Education
Albania 1.668337 4578.633 1.774449 2309.3665 0.000000 0.00000 0.0000000 3.0758947 56.423 5.67 100.00000 19.990181 0 3.152284
Algeria 3.811645 5516.231 3.800000 1368.6215 0.000000 98.61150 0.9931198 24.1267274 70.221 3.83 0.39538 10.286397 0 4.467196
Angola 1.092216 5059.080 4.820000 310.0817 0.000000 0.00000 46.8248945 23.7586261 62.731 3.35 53.17511 7.547057 0 3.570000
Argentina 4.209096 12334.798 -2.512615 3074.7021 2.217284 48.19968 13.8401066 2.7452689 91.377 6.84 35.74293 5.156686 0 5.094119
Armenia 1.896240 4017.230 3.600000 1976.8645 0.000000 42.43871 0.0000000 0.3069489 63.112 4.13 57.56129 18.072036 0 2.100000
Australia 15.830421 62513.411 2.579017 10071.3990 61.164325 21.90974 2.0192215 4.7405285 85.602 9.01 14.90671 2.375955 0 4.864731 
# setting countries as row names and leaving only numeric variables in the dataset # - save as data_num 
data_num<-data[,2:15]
r_name<-data[,1]$Country
rownames(data_num) <- r_name
# Basic statistics of each variable in the dataset

summary(data_num)
##       CO2              GDPpc               GDPg         Electric_power_cons
##  Min.   : 0.0676   Min.   :   472.3   Min.   :-23.043   Min.   :   39.57   
##  1st Qu.: 1.6229   1st Qu.:  3442.1   1st Qu.:  1.712   1st Qu.:  918.85   
##  Median : 3.7139   Median :  7788.4   Median :  3.461   Median : 2598.33   
##  Mean   : 5.3529   Mean   : 18093.9   Mean   :  3.328   Mean   : 4427.75   
##  3rd Qu.: 7.2858   3rd Qu.: 23438.0   3rd Qu.:  5.243   3rd Qu.: 5654.10   
##  Max.   :36.8757   Max.   :123678.7   Max.   : 19.047   Max.   :53832.48   
##       Coal              Gas                 Oil           Resources_rents  
##  Min.   : 0.0000   Min.   :  0.00000   Min.   :  0.0000   Min.   : 0.0000  
##  1st Qu.: 0.0000   1st Qu.:  0.00576   1st Qu.:  0.2643   1st Qu.: 0.6079  
##  Median : 0.6451   Median : 13.40075   Median :  1.5157   Median : 2.2378  
##  Mean   :16.3496   Mean   : 26.88452   Mean   : 15.0654   Mean   : 7.1058  
##  3rd Qu.:27.9662   3rd Qu.: 47.15049   3rd Qu.: 15.2885   3rd Qu.: 8.4711  
##  Max.   :95.6766   Max.   :100.00000   Max.   :100.0000   Max.   :51.7565  
##   Urbanization    Democracy_index Hydro_renewable_nucl  Agriculture      
##  Min.   : 16.22   Min.   :1.740   Min.   :  0.000      Min.   : 0.03454  
##  1st Qu.: 51.78   1st Qu.:3.822   1st Qu.:  8.564      1st Qu.: 2.38154  
##  Median : 66.37   Median :6.255   Median : 35.802      Median : 5.23201  
##  Mean   : 63.96   Mean   :5.921   Mean   : 39.745      Mean   : 8.47905  
##  3rd Qu.: 79.39   3rd Qu.:7.713   3rd Qu.: 64.144      3rd Qu.:11.90835  
##  Max.   :100.00   Max.   :9.930   Max.   :100.000      Max.   :38.52045  
##        EU           Education     
##  Min.   :0.0000   Min.   : 0.850  
##  1st Qu.:0.0000   1st Qu.: 2.986  
##  Median :0.0000   Median : 4.182  
##  Mean   :0.2154   Mean   : 4.245  
##  3rd Qu.:0.0000   3rd Qu.: 5.107  
##  Max.   :1.0000   Max.   :13.010

As we can see from the statistics presented above, the variables in question have significantly different orders of magnitude by which, in order to perform cluster analysis, it will be necessary to standardize them. The first couple rows after the standardization of the variables are shown below.

data_z <- as.data.frame(lapply(data_num, scale))
rownames(data_z) <- r_name
kable(head(data_z))
CO2 GDPpc GDPg Electric_power_cons Coal Gas Oil Resources_rents Urbanization Democracy_index Hydro_renewable_nucl Agriculture EU Education
Albania -0.6412154 -0.5945006 -0.3805143 -0.3437715 -0.6676657 -0.8337592 -0.5795526 -0.3742667 -0.3703359 -0.1182093 1.8086308 1.3750123 -0.5219178 -0.6126878
Algeria -0.2682197 -0.5532582 0.1154673 -0.4964358 -0.6676657 2.2244409 -0.5413483 1.5807938 0.3078595 -0.9839676 -1.1811149 0.2158886 -0.5219178 0.1247714
Angola -0.7414767 -0.5733670 0.3652271 -0.6682159 -0.6676657 -0.8337592 1.2217569 1.5466070 -0.0602868 -1.2098176 0.4031284 -0.1113274 -0.5219178 -0.3784148
Argentina -0.1990521 -0.2533272 -1.4302558 -0.2195727 -0.5771188 0.6610387 -0.0471368 -0.4049730 1.3477132 0.4323000 -0.1201181 -0.3968588 -0.5219178 0.4763769
Armenia -0.6015540 -0.6191953 0.0664948 -0.3977300 -0.6676657 0.4823760 -0.5795526 -0.6314279 -0.0415600 -0.8428113 0.5347847 1.1458886 -0.5219178 -1.2028545
Australia 1.8233847 1.9538985 -0.1835057 0.9158528 1.8300947 -0.1542809 -0.5018751 -0.2196667 1.0638621 1.4533302 -0.7455408 -0.7290190 -0.5219178 0.3477265

A chart showing the correlations between each pair of variables (Pearson correlation coefficient was used) was used to examine the relationship between the variables. As we can read from the graph, some variables are more strongly correlated (positively or negatively) with the other variables. An example of such a variable is the CO2 variable, which is strongly positively correlated with the GDPpc, or Electric_power_cons variables and negatively with agriculture. The set also includes variables that are slightly correlated with the other variables, such as Oil or GDPg.

kor<-cor(data_z)
corrplot(kor, method="color", type="lower", tl.col = 'black', tl.cex = 0.75)

Prediagnostics

In order to examine whether and how well the data can be clustered the Hopkins ststistic was used. If Hopkins statistic is smaller than 0.5, then there are no visible clusters and data is uniformly distributed. If Hopkins statistics is bigger than 0.5, then dataset is significantly clusterable data. Graphical analysis was also used to determine if the data is clusterable, for this purpose the Ordered Dissimilarity Matrix was plotted.

get_clust_tendency(data_z, 2, graph=TRUE, gradient=list(low="red", mid="white", high="blue"), seed = 123)
## $hopkins_stat
## [1] 0.8516717
## 
## $plot

The Hopkins statistic for the dataset is 0.85, this means that we reject the null hypothesis, which states that the dataset is uniformly distributed. So, the dataset contains meaningful clusters. Similar conclusions can also be drawn from an analysis of the chart above. The blue color indicates a large distance (high dissimilarity) between objects of the dataset, and the red color indicates a small distance (low dissimilarity) between objects. In the above graph, both red and blue color blocks are visible, which allows us to conclude that the data tends to form into clusters.

Both the Hopkins statistic and the Ordered Dissimilarity Matrix analysis allow us to conclude that the dataset is clustered, so we can move on to cluster analysis

Cluster analysis

Selecting optimal number of clusters

The following methods were used to perform cluster analysis: k-means, PAM and hierarchical clustering. The CLARA method was omitted, as it is primarily applicable to large datasets

Before performing cluster analysis, it is necessary to first determine the optimal number of clusters into which the dataset should be divided. To do this, several different methods were used, such as the Silhouette index, the GAP statistic and the distance-based approach - a method that finds the optimal number of clusters depending on the distance between objects in dataset (the euclidean distance was chosen to preform this method). Rule of thumb and elbow method were not used because of they are too simplistic.

# Distance-based approach

opt1<-NbClust(data_z, distance="euclidean", min.nc=2, max.nc=10, method="complete", index="ch")
opt1$All.index
##       2       3       4       5       6       7       8       9      10 
##  6.3411 16.8244 14.4113 16.8271 17.5089 18.1827 16.5373 15.7827 16.8254
opt1$Best.nc
## Number_clusters     Value_Index 
##          7.0000         18.1827
# Silhouette index

# k-means
s1 <- fviz_nbclust(data_z, FUNcluster = kmeans, method = "silhouette") + theme_classic() 
# PAM
s2 <- fviz_nbclust(data_z, FUNcluster = cluster::pam, method = "silhouette") + theme_classic()
# hierarchical clustering
s3 <- fviz_nbclust(data_z, FUNcluster = hcut, method = "silhouette") + theme_classic()

grid.arrange(s1, s2, s3, ncol=2)

# GAP ststistic

# k-means
g1 <- fviz_nbclust(data_z, FUNcluster = kmeans, method = "gap_stat") + theme_classic() 
# PAM
g2 <- fviz_nbclust(data_z, FUNcluster = cluster::pam, method = "gap_stat") + theme_classic()
# hierarchical clustering
g3 <- fviz_nbclust(data_z, FUNcluster = hcut, method = "gap_stat") + theme_classic()

grid.arrange(g1, g2, g3, ncol=2)

As we can observe in the graphs showing the average Silhouette width depending on the number of clusters into which the dataset can be divided, for chosen methods the most optimal number of clusters is 3. The exception is the PAM algorithm, for which this number is 9, but choosing such a large number of clusters can cause risks of overfitting the data. A similar number of clusters is obtained from interpreting the graph representing the GAP statistic (k), where k is the number of clusters into which the set can be divided. The graphs show that the GAP statistic achieves its highest value for the highest number of clusters (9 or 10), but this is too large a number of clusters to divide a set of only 130 observations into. In the graphs for PAM and k-means, we can observe an “elbow point” where the Gap Statistic is maximized dal k = 3 or 4.

The last of the selected criteria for selecting the optimal number of clusters was a distance-based approach based on euclidean distance. According to this approach, the best choice would be a division into 7 clusters, but for 3 or 5 clusters a division would also be optimal. In summary, based on each of the criteria discussed, it can be decided that the best for each of the selected methods will be the division of the dataset into 3 clusters. We can then proceed to cluster analysis using each of the selected methods.

Clustering using K-means

data_clusters3 <- kmeans(data_z, 3)
data_clusters3$size
## [1] 49 62 19
fviz_cluster(data_clusters3, data = data_z, alpha=0.2, geom = c("point", "text"), labelsize = 6)

sil<-silhouette(data_clusters3$cluster, dist(data_z))
fviz_silhouette(sil)
##   cluster size ave.sil.width
## 1       1   49          0.22
## 2       2   62          0.23
## 3       3   19          0.18

Above, clustering of the dataset was carried out using the k-means method into 3 clusters. The graph above shows the dataset divided into 3 clusters, as you can see it is 2-dimensional although our dataset has 14 dimensions but it would be impossible to present it in this form. The dimensions shown on the set account for almost 50% of the variance in the entire set. As you can see from the graph, we managed to create three separate clusters that do not overlap. The silhouette chart, on the other hand, allows us to assess the quality of the clustering performed. It shows that the average silhouette width is 0.22. Each of the formed clusters also has an average silhouette width close to 0.2 so none of them is significantly “worse quality” than the others.

Clustering using PAM

pam_clusters3<-pam(data_z,3) 

#Number of observation in each cluster
as.data.frame(pam_clusters3$clustering) %>%
  group_by(pam_clusters3$clustering) %>%
  summarize(Count = n())
## # A tibble: 3 × 2
##   `pam_clusters3$clustering` Count
##                        <int> <int>
## 1                          1    40
## 2                          2    55
## 3                          3    35
fviz_cluster(pam_clusters3, data = data_z, alpha=0.2, geom = c("point", "text"), labelsize = 6)

sil_pam<-silhouette(pam_clusters3$clustering, dist(data_z))
fviz_silhouette(sil_pam)
##   cluster size ave.sil.width
## 1       1   40          0.22
## 2       2   55          0.05
## 3       3   35          0.29

Another of the clustering methods chosen in this study was PAM, the graph above shows the clusters obtained with this method on a two-dimensional graph. In it, we can observe that cluster 2 overlaps with the others, which may indicate an incorrect assignment of individual observations to this cluster. Analysis of the silhouette plot allows us to conclude that this method is slightly weaker than k-means, as the average width of the silhouette is 0.17. It can also be seen that one of the clusters is much weaker than the other two (the average width of the silhouette is equal to 0.05, while for the others it is above 0.2), this means that some of the observations from this cluster should be in other clusters.

Clustering using hierarchical clustering

hc_clusters3 <- eclust(data_z, "hclust", k = 3,
                 method = "ward.D2", graph = FALSE) 
hc_clusters3$size
## [1] 85 36  9
fviz_dend(hc_clusters3, rect = TRUE, show_labels = TRUE, cex = 0.5)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

fviz_cluster(hc_clusters3, data = data_z, alpha=0.2, geom = c("point", "text"), labelsize = 6)

sil_hc<-silhouette(hc_clusters3$cluster, dist(data_z))
fviz_silhouette(sil_hc)
##   cluster size ave.sil.width
## 1       1   85          0.15
## 2       2   36          0.32
## 3       3    9          0.33

The last method selected was hierarchical clustering. In order to visualize its results, a dendrogram was used. Two-dimensional graph were also presented, in order to more easily compare the results with the previous methods. The dendrogram shows that the first cluster is larger than the other two combined, and that it is less similar to them than they are to each other. The average silhouette width for this method is 0.21, which is the same as for k-means clustering. The silhouette plot also shows that the first (the largest of the clusters) has a much smaller silhouette width (0.15) than the other two clusters (above 0.3). (Note: the colors of the 2 and cluster do not match between the dendrogram and the 2-dimensional graph)

Short summary of cluster analisys

Based on the charts above, you can conduct a preliminary cluster analysis and try to find some common characteristics of countries in the same cluster. The most similar clusters to each other between the different methods are cluster 1 - k-means, cluster 3 - pam and cluster 2 - hc. The countries in each of these clusters can be described as mostly democratic highly developed countries. Cluster 3 for k-means and hc consists primarily of Middle Eastern countries with rich oil reserves. In contrast, the largest of the clusters for each of the methods comprises predominantly underdeveloped and developing countries.

Analysis of clustering results

In order to conduct a more detailed analysis of the clusters, the average values of each variable for each of the clusters formed by clustering method are presented above. When analyzing the above statistics, it should be borne in mind that for each clustering method, individual clusters have different numbers (e.g., for k-means, cluster 1 corresponds to cluster 3 for pam).

data$cluser_km <- data_clusters3$cluster
data$cluser_pam <- pam_clusters3$clustering
data$cluser_hc <- hc_clusters3$cluster
a <- aggregate(data, CO2 ~ cluser_km, mean)
b <- aggregate(data, CO2 ~ cluser_pam, mean)
c <- aggregate(data, CO2 ~ cluser_hc, mean)
total1 <- cbind(a, b)
total1 <- cbind(total1, c)

a <- aggregate(data, GDPpc ~ cluser_km, mean)
b <- aggregate(data, GDPpc ~ cluser_pam, mean)
c <- aggregate(data, GDPpc ~ cluser_hc, mean)
total2 <- cbind(a, b)
total2 <- cbind(total2, c)

a <- aggregate(data, GDPg ~ cluser_km, mean)
b <- aggregate(data, GDPg ~ cluser_pam, mean)
c <- aggregate(data, GDPg ~ cluser_hc, mean)
total3 <- cbind(a, b)
total3 <- cbind(total3, c)

a <- aggregate(data, Electric_power_cons ~ cluser_km, mean)
b <- aggregate(data, Electric_power_cons ~ cluser_pam, mean)
c <- aggregate(data, Electric_power_cons ~ cluser_hc, mean)
total4 <- cbind(a, b)
total4 <- cbind(total4, c)

a <- aggregate(data, Coal ~ cluser_km, mean)
b <- aggregate(data, Coal ~ cluser_pam, mean)
c <- aggregate(data, Coal ~ cluser_hc, mean)
total5 <- cbind(a, b)
total5 <- cbind(total5, c)

a <- aggregate(data, Gas ~ cluser_km, mean)
b <- aggregate(data, Gas ~ cluser_pam, mean)
c <- aggregate(data, Gas ~ cluser_hc, mean)
total6 <- cbind(a, b)
total6 <- cbind(total6, c)

a <- aggregate(data, Oil ~ cluser_km, mean)
b <- aggregate(data, Oil ~ cluser_pam, mean)
c <- aggregate(data, Oil ~ cluser_hc, mean)
total7 <- cbind(a, b)
total7 <- cbind(total7, c)

a <- aggregate(data, Resources_rents ~ cluser_km, mean)
b <- aggregate(data, Resources_rents ~ cluser_pam, mean)
c <- aggregate(data, Resources_rents ~ cluser_hc, mean)
total8 <- cbind(a, b)
total8 <- cbind(total8, c)

a <- aggregate(data, Urbanization ~ cluser_km, mean)
b <- aggregate(data, Urbanization ~ cluser_pam, mean)
c <- aggregate(data, Urbanization ~ cluser_hc, mean)
total9 <- cbind(a, b)
total9 <- cbind(total9, c)

a <- aggregate(data, Democracy_index ~ cluser_km, mean)
b <- aggregate(data, Democracy_index ~ cluser_pam, mean)
c <- aggregate(data, Democracy_index ~ cluser_hc, mean)
total10 <- cbind(a, b)
total10 <- cbind(total10, c)

a <- aggregate(data, Agriculture ~ cluser_km, mean)
b <- aggregate(data, Agriculture ~ cluser_pam, mean)
c <- aggregate(data, Agriculture ~ cluser_hc, mean)
total11 <- cbind(a, b)
total11 <- cbind(total11, c)

a <- aggregate(data, Hydro_renewable_nucl ~ cluser_km, mean)
b <- aggregate(data, Hydro_renewable_nucl ~ cluser_pam, mean)
c <- aggregate(data, Hydro_renewable_nucl ~ cluser_hc, mean)
total12 <- cbind(a, b)
total12 <- cbind(total12, c)

a <- aggregate(data, EU ~ cluser_km, mean)
b <- aggregate(data, EU ~ cluser_pam, mean)
c <- aggregate(data, EU ~ cluser_hc, mean)
total13 <- cbind(a, b)
total13 <- cbind(total13, c)

a <- aggregate(data, Education ~ cluser_km, mean)
b <- aggregate(data, Education ~ cluser_pam, mean)
c <- aggregate(data, Education ~ cluser_hc, mean)
total14 <- cbind(a, b)
total14 <- cbind(total14, c)
kable(total1)
cluser_km CO2 cluser_pam CO2 cluser_hc CO2
1 6.800509 1 1.223519 1 2.816187
2 1.809129 2 7.356201 2 7.755104
3 13.183401 3 6.924108 3 19.701772 
kable(total2)
cluser_km GDPpc cluser_pam GDPpc cluser_hc GDPpc
1 33662.856 1 3390.412 1 5774.517
2 3894.618 2 15830.843 2 41307.982
3 24276.794 3 38454.072 3 41587.117 
kable(total3)
cluser_km GDPg cluser_pam GDPg cluser_hc GDPg
1 1.908494 1 4.301232 1 3.845842
2 4.922061 2 3.417088 2 2.270958
3 1.790172 3 2.077376 3 2.671795 
kable(total4)
cluser_km Electric_power_cons cluser_pam Electric_power_cons cluser_hc Electric_power_cons
1 7438.079 1 922.0101 1 1729.447
2 1250.278 2 4381.8258 2 9005.819
3 7032.852 3 8506.4710 3 11599.425 
kable(total5)
cluser_km Coal cluser_pam Coal cluser_hc Coal
1 23.786668 1 3.633065 1 15.3754034
2 14.062451 2 24.714440 2 22.7066108
3 4.632959 3 17.737946 3 0.1219567 
kable(total6)
cluser_km Gas cluser_pam Gas cluser_hc Gas
1 17.09233 1 7.125256 1 23.97705
2 19.81800 2 49.069013 2 18.94869
3 75.19725 3 14.605196 3 86.08733 
kable(total7)
cluser_km Oil cluser_pam Oil cluser_hc Oil
1 7.225613 1 32.53185 1 21.237016
2 22.390151 2 7.46233 2 1.594995
3 11.382141 3 7.05151 3 10.659880 
kable(total8)
cluser_km Resources_rents cluser_pam Resources_rents cluser_hc Resources_rents
1 1.312829 1 7.225695 1 7.6489722
2 5.986187 2 10.968244 2 0.8851269
3 25.698752 3 0.899065 3 26.8579200 
kable(total9)
cluser_km Urbanization cluser_pam Urbanization cluser_hc Urbanization
1 76.28529 1 47.75620 1 56.68253
2 50.30518 2 68.09122 2 75.84606
3 76.71479 3 75.97760 3 85.11200 
kable(total10)
cluser_km Democracy_index cluser_pam Democracy_index cluser_hc Democracy_index
1 7.900408 1 4.845750 1 5.192588
2 5.111129 2 5.289818 2 8.198611
3 3.460526 3 8.142571 3 3.693333 
kable(total11)
cluser_km Agriculture cluser_pam Agriculture cluser_hc Agriculture
1 2.818767 1 16.261254 1 11.9238988
2 14.459487 2 6.618889 2 2.2758058
3 3.561522 3 2.508219 3 0.7573675 
kable(total12)
cluser_km Hydro_renewable_nucl cluser_pam Hydro_renewable_nucl cluser_hc Hydro_renewable_nucl
1 48.473406 1 54.37694 1 38.2754300
2 43.386336 2 17.97672 2 53.0954383
3 5.350697 3 57.22901 3 0.2184033 
kable(total13)
cluser_km EU cluser_pam EU cluser_hc EU
1 0.5714286 1 0.0 1 0.0235294
2 0.0000000 2 0.0 2 0.7222222
3 0.0000000 3 0.8 3 0.0000000 
kable(total14)
cluser_km Education cluser_pam Education cluser_hc Education
1 5.297073 1 3.922557 1 4.045955
2 3.645748 2 3.979542 2 4.882703
3 3.485328 3 5.029632 3 3.570075

The most similar clusters from each of the methods used are as follows: 1 for k-means, 3 for PAM and 2 for hierarchical clustering. After analyzing the statistics presented above, it can be determined that these clusters primarily include highly developed democratic countries (e.g. Norway, France, Canada). They are characterized, on average, by a high level of CO2 emissions per capita, a high level of electricity consumption and, on average, a higher share of renewable energy sources in total electricity generation than countries in the other clusters. These countries are also among the most developed, with on average the highest GDP per capita, low annual GDP growth, the highest level of urbanization and the lowest share of agriculture in GDP. They are also, on average, the most democratic countries and with the highest spending on education.

The second group of countries we can distinguish are the rich countries of the Middle East and Asia (e.g. Kuwait, Saudi Arabia, Singapore), which are represented primarily in the 3rd cluster for k-means and hierarchical clustering (these countries are also in the 2nd cluster obtained by the PAM method, but this cluster also contains many other countries and the average silhouette width is very low for it, so we will not interpret these results). They are characterized primarily by the highest CO2 consumption on average, high GDP per capita, and the highest share of natural gas in electricity generation while having a very low share of renewable energy sources. These countries also have the highest level of natural resource rents as a % of GDP. They are also the least democratic countries with low spending on education as a %.

The last group of countries we can distinguish includes primarily those we can describe as underdeveloped and developing countries (e.g. Ethiopia, Haiti, Paraguay). They are characterized by by far the lowest CO2 emissions per capita, the smallest GDP per capita and the lowest electricity consumption. They are also countries with, on average, the highest levels of economic growth, low levels of urbanization and a high share of agriculture in GDP (indicating that they are developing countries). Electricity generation in these countries is based primarily on renewable energy sources and oil. These countries are also on average less democratic than those in the first group of clusters discussed but more so than those in the second.

Macroeconomic indicators analysis

We now turn to the second part of the study’s analysis, which consists of comparing the clustering results with those of a set narrowed down only to purely macroeconomic variables (e.g., GDP per capita, GDP growth), with no “social” variables (e.g., Education level). This analysis will allow us to assess whether factors such as education level or democracy index have an impact on CO2 emissions, and whether countries with certain macroeconomic characteristics in a given group are also similar in other respects. We can therefore learn more about the characteristics of such countries and group them in a more interesting way.

data_ <- read_excel("CO2.xlsx")
data_<-data_[,c(2:8,12)]
data_m <- as.data.frame(lapply(data_, scale))
rownames(data_m) <- data[,1]$Country
kable(head(data_m))
CO2 GDPpc GDPg Electric_power_cons Coal Gas Oil Hydro_renewable_nucl
Albania -0.6412154 -0.5945006 -0.3805143 -0.3437715 -0.6676657 -0.8337592 -0.5795526 1.8086308
Algeria -0.2682197 -0.5532582 0.1154673 -0.4964358 -0.6676657 2.2244409 -0.5413483 -1.1811149
Angola -0.7414767 -0.5733670 0.3652271 -0.6682159 -0.6676657 -0.8337592 1.2217569 0.4031284
Argentina -0.1990521 -0.2533272 -1.4302558 -0.2195727 -0.5771188 0.6610387 -0.0471368 -0.1201181
Armenia -0.6015540 -0.6191953 0.0664948 -0.3977300 -0.6676657 0.4823760 -0.5795526 0.5347847
Australia 1.8233847 1.9538985 -0.1835057 0.9158528 1.8300947 -0.1542809 -0.5018751 -0.7455408

Diagnostic tests

hopkins(data_m, n=nrow(data_m)-1) 
## Warning in hopkins(data_m, n = nrow(data_m) - 1): Package `clustertend` is
## deprecated.  Use package `hopkins` instead.
## $H
## [1] 0.1669901
get_clust_tendency(data_m, 2, graph=TRUE, gradient=list(low="red", mid="white", high="blue"), seed = 123)
## $hopkins_stat
## [1] 0.8626502
## 
## $plot

Both the Hopkins statistic of 0.86 and the visible red and blue squares in the graph indicate that the dataset contains meaningful clusters. We can therefore proceed to further analysis

opt_m<-NbClust(data_m, distance="euclidean", min.nc=2, max.nc=10, method="complete", index="ch")
opt_m$All.index
##       2       3       4       5       6       7       8       9      10 
##  8.1624 19.2935 16.5169 17.1231 22.7409 20.6762 26.3413 24.4041 31.5724
opt_m$Best.nc
## Number_clusters     Value_Index 
##         10.0000         31.5724
s1_m <- fviz_nbclust(data_m, FUNcluster = kmeans, method = "silhouette") + theme_classic() 
g1_m <- fviz_nbclust(data_m, FUNcluster = kmeans, method = "gap_stat") + theme_classic() 
grid.arrange(s1_m, g1_m, ncol=2)

The graphs showing the average width of the silhouette and Gap stastistic indicate that we should choose 8 or 9 clusters. However, this is too large a number, which can lead to overfitting the data. We decide to choose 3 clusters (even if choosing 5 or 6 clusters would be more appropriate choice) , as this way it will be simpler for us to compare the results with those obtained from clustering the first dataset.

K-means

It was decided to choose the k-means algorithm because it performed best in clustering the underlying set. Only one algorithm was chosen, as it significantly simplifies the comparative analysis of the clustering of the two sets in question.

data_clusters3_m <- kmeans(data_m, 3)
data_clusters3_m$size
## [1] 22 94 14
fviz_cluster(data_clusters3_m, data = data_m, alpha=0.2, geom = c("point", "text"), labelsize = 6)

sil_m<-silhouette(data_clusters3_m$cluster, dist(data_m))
fviz_silhouette(sil_m)
##   cluster size ave.sil.width
## 1       1   22          0.12
## 2       2   94          0.26
## 3       3   14          0.30

The two-dimensional graph shows that of the 3 clusters formed, 1 and 3 overlap in part. The average width of the silhouette for this set is 0.2, with the highest value (0.28) for the largest of the clusters - cluster 1, and this is almost double that for the other clusters.

We can also observe, at a glance, that the perimeter of each cluster includes different countries from the clustering of the base set. For example, cluster 2 includes both highly developed EU countries and non-democratic Middle Eastern countries, which previously were always in separate clusters.

Comparative analysis

data_$cluser_km <- data_clusters3_m$cluster
kable(aggregate(data_, CO2 ~ cluser_km, mean))
cluser_km CO2
1 14.788400
2 2.853628
3 7.306406 
kable(aggregate(data_, GDPpc ~ cluser_km, mean))
cluser_km GDPpc
1 40450.939
2 7353.198
3 55077.431 
kable(aggregate(data_, GDPg ~ cluser_km, mean))
cluser_km GDPg
1 1.898081
2 3.861295
3 1.998413 
kable(aggregate(data_, Electric_power_cons ~ cluser_km, mean))
cluser_km Electric_power_cons
1 9098.962
2 2021.156
3 13245.823 
kable(aggregate(data_, Coal ~ cluser_km, mean))
cluser_km Coal
1 13.28083
2 17.62889
3 12.58219 
kable(aggregate(data_, Gas ~ cluser_km, mean))
cluser_km Gas
1 64.311859
2 20.814346
3 8.827027 
kable(aggregate(data_, Oil ~ cluser_km, mean))
cluser_km Oil
1 7.3947257
2 19.0380069
3 0.4463184 
kable(aggregate(data_, Hydro_renewable_nucl ~ cluser_km, mean))
cluser_km Hydro_renewable_nucl
1 9.149899
2 41.378572
3 76.852098

As we can see in the chart and statistics above, countries have been assigned to individual clusters in a different way from the baseline dataset. Thus, we can distinguish a set of rich countries using renewable energy sources (cluster 3), rich countries basing their energy mostly on natural gas (cluster 1) and developing countries (cluster 2). However, this way of grouping these countries causes us to lose a lot of important information about individual countries, so that countries that are similar on many levels end up in different clusters. For example, countries such as Luxembourg, the United States, or Australia are more similar to highly developed EU countries in terms of spending on education, or the level of democracy, which also affects the approach to climate change. These countries, therefore, should not be in the same group as the oil-rich countries of the Middle East. The situation is similar for countries such as Italy or Croatia , which are assigned to the cluster of developing countries instead of the democratic highly developed countries. Thus, it can be concluded that information on factors that are more ‘social’ and not directly related to CO2 emissions allows us to group countries into groups that better explain the approach to climate change by these countries.

Conclusion

In this study, a cluster analysis was conducted grouping countries with similar CO2 emissions characteristics. Using three selected clustering methods, it was possible to distinguish 3 groups of countries with different characteristics: highly developed democratic countries, rich countries of the Middle East and Asia, underdeveloped and developing countries. This analysis also made it possible to determine the impact of selected variables on CO2 emissions, for example, non-democratic countries with a high share of natural resource rents in GDP are among the largest producers of CO2.

An additional analysis examining the impact of non-macroeconomic factors on CO2 emissions was also conducted. It turned out that these factors allowed for a better assignment of countries to particular clusters, giving us a better picture of what similarities exist between countries within clusters.