1. Wstęp

W projekcie wykonuję typologię państw na podstawie poziomu rozwoju, korzystając z danych Banku Światowego. Celem jest znalezienie naturalnych grup krajów o podobnym profilu rozwojowym, bez narzuconych etykiet (uczenie nienadzorowane). Wykorzystuję wskaźniki ekonomiczne, społeczne i środowiskowe, a następnie stosuję klastrowanie K‑średnich.

2. Podstawy teoretyczne

Uczenie nienadzorowane polega na odkrywaniu struktury w danych bez znanych klas. Algorytm K‑średnich minimalizuje sumę odległości wewnątrz klastrów, dlatego wymaga standaryzacji zmiennych. Weryfikuję, czy dane mają tendencję do grupowania (statystyka Hopkinsa), a liczbę klastrów dobieram metodą łokcia i sylwetki oraz większością głosów z NbClust. Do wizualizacji wysokowymiarowych danych stosuję PCA, które redukuje wymiary przy zachowaniu jak największej wariancji.

Pobierzemy wskaźniki z API Banku Światowego (pakiet WDI - World Development Indicators and Other World Bank Data) dla najnowszego dostępnego roku. Wskaźniki: PKB per capita, oczekiwana długość życia, emisja CO2 oraz wydatki na R&D. Jeśli któraś zmienna nie jest dostępna w API, zostanie automatycznie pominięta.

3. Implementacja systemu

# Definiowanie wskaźników
indicators_list <- c(
  gdp_pc   = "NY.GDP.PCAP.CD",   # PKB per capita (USD)
  life_exp = "SP.DYN.LE00.IN",   # Oczekiwana długość życia
  co2      = "EN.ATM.CO2E.PC",   # Emisja CO2 per capita
  rd_exp   = "GB.XPD.RSDV.GD.ZS" # Wydatki na R&D (% PKB)
)

# Pobieranie danych
raw_data <- WDI(indicator = indicators_list, country = "all", start = 2010, end = 2022, extra = TRUE)

# 3.1. Konwersja i wybór istotnych kolumn (obsługa różnych nazw kolumn z WDI)
name_map <- c(
  gdp_pc   = "NY.GDP.PCAP.CD",
  life_exp = "SP.DYN.LE00.IN",
  co2      = "EN.ATM.CO2E.PC",
  rd_exp   = "GB.XPD.RSDV.GD.ZS"
)
alt_map <- c(
  gdp_pc   = "NY_GDP_PCAP_CD",
  life_exp = "SP_DYN_LE00_IN",
  co2      = "EN_ATM_CO2E_PC",
  rd_exp   = "GB_XPD_RSDV_GD_ZS"
)

clean_data <- raw_data %>%
  as_tibble() %>%
  filter(region != "Aggregates") %>%
  rename(
    gdp_pc = any_of(c(name_map["gdp_pc"], alt_map["gdp_pc"])),
    life_exp = any_of(c(name_map["life_exp"], alt_map["life_exp"])),
    co2 = any_of(c(name_map["co2"], alt_map["co2"])),
    rd_exp = any_of(c(name_map["rd_exp"], alt_map["rd_exp"]))
  )

missing_cols <- setdiff(names(indicators_list), colnames(clean_data))
if (length(missing_cols) > 0) {
  warning(paste("Brak kolumn po pobraniu danych:", paste(missing_cols, collapse = ", ")))
  indicators_list <- indicators_list[setdiff(names(indicators_list), missing_cols)]
}

clean_data <- clean_data %>%
  select(country, iso3c, year, all_of(names(indicators_list)))

if (length(indicators_list) < 2) {
  stop("Zbyt mało zmiennych do klastrowania. Dodaj inne wskaźniki.")
}

# 3.2. Wybór najnowszego roku z kompletnymi danymi dla każdego kraju
latest_complete <- clean_data %>%
  filter(if_all(all_of(names(indicators_list)), ~ !is.na(.x))) %>%
  group_by(country, iso3c) %>%
  arrange(desc(year)) %>%
  slice(1) %>%
  ungroup()

# 3.3. Dane surowe do interpretacji + logarytmowanie zmiennych skośnych
df_final_raw <- latest_complete %>%
  select(country, all_of(names(indicators_list)))

log_vars <- intersect(c("gdp_pc", "co2"), names(df_final_raw))
df_model <- df_final_raw %>%
  mutate(across(all_of(log_vars), ~ log1p(.x)))

# 3.4. Przygotowanie macierzy do klastrowania
df_final <- df_model %>% column_to_rownames("country")
df_scaled <- scale(df_final)

# Wyświetlenie podglądu danych po standaryzacji
head(df_scaled)
##               gdp_pc   life_exp     rd_exp
## Algeria   -0.4205998  0.2379767 -0.3880737
## Angola    -0.9535642 -1.5656286 -0.7767172
## Argentina  0.3282435  0.2869451 -0.3264111
## Armenia   -0.1755026  0.1511202 -0.6223778
## Australia  1.3113922  1.2655302  0.8103858
## Austria    1.2116324  1.0037283  1.9818102

4. Ocena tendencji do tworzenia skupień

Przed przystąpieniem do klastrowania musimy sprawdzić, czy dane mają sensowną strukturę (nie są losowe).

# Statystyka Hopkinsa: wartości > 0.5 (idealnie > 0.7) wskazują na dane podatne na klastrowanie
set.seed(123)
hopkins_stat <- hopkins(df_scaled)
print(paste("Statystyka Hopkinsa:", round(hopkins_stat, 3)))
## [1] "Statystyka Hopkinsa: 0.995"
# Wizualna ocena tendencji do klastrowania (VAT)
fviz_dist(dist(df_scaled), show_labels = FALSE) + 
  ggtitle("Wizualna ocena tendencji do klastrowania")

5. Wyznaczanie optymalnej liczby klastrów (k)

Użyjemy dwóch popularnych metod: metody łokcia (Elbow) oraz metody sylwetki (Silhouette).

# Porównanie wyników metody łokcia i metody sylwetki
p1 <- fviz_nbclust(df_scaled, kmeans, method = "wss") + ggtitle("Metoda łokcia")
p2 <- fviz_nbclust(df_scaled, kmeans, method = "silhouette") + ggtitle("Metoda sylwetki")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2)

# Wybór k na podstawie większości głosów (NbClust)
set.seed(123)
nb <- NbClust(df_scaled, min.nc = 2, max.nc = 8, method = "kmeans", index = "alllong")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 7 proposed 2 as the best number of clusters 
## * 8 proposed 3 as the best number of clusters 
## * 7 proposed 4 as the best number of clusters 
## * 1 proposed 5 as the best number of clusters 
## * 1 proposed 6 as the best number of clusters 
## * 3 proposed 8 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  3 
##  
##  
## *******************************************************************
best_k <- as.numeric(names(sort(table(nb$Best.nc[1, ]), decreasing = TRUE)[1]))
print(paste("Rekomendowana liczba klastrów (NbClust):", best_k))
## [1] "Rekomendowana liczba klastrów (NbClust): 3"

6. Wykonanie algorytmu K-średnich

Przyjmując optymalne k na podstawie powyższych wykresów i NbClust.

# Ustawienie ziarna dla powtarzalności wyników
set.seed(123)
km_result <- kmeans(df_scaled, centers = best_k, nstart = 50, iter.max = 200)

# Dodanie przypisania do klastrów do oryginalnych (nieskalowanych) danych
final_results <- data.frame(df_final_raw, cluster = as.factor(km_result$cluster))

# Ocena jakości klastrów (średni współczynnik sylwetki)
sil <- silhouette(km_result$cluster, dist(df_scaled))
print(paste("Średnia sylwetka:", round(mean(sil[, 3]), 3)))
## [1] "Średnia sylwetka: 0.442"

7. Wizualizacja w przestrzeni PCA

Ponieważ mamy 4 zmienne, używamy analizy głównych składowych (PCA), aby pokazać klastry w 2D.

# Wykres klastrów na płaszczyźnie dwóch pierwszych składowych głównych
fviz_cluster(km_result, data = df_scaled,
             palette = "Set2", 
             geom = "point",
             ellipse.type = "convex", 
             ggtheme = theme_minimal(),
             main = "Wizualizacja klastrów w przestrzeni PCA")

8. Interpretacja klastrów

Analizujemy średnie wartości zmiennych w każdym klastrze, aby zdefiniować typy państw.

# Tabela podsumowująca średnie wartości w klastrach
cluster_profile <- final_results %>%
  group_by(cluster) %>%
  summarise(
    across(where(is.numeric), mean),
    count = n()
  )

# Wyświetlenie profilu klastrów
print(cluster_profile)
## # A tibble: 3 × 5
##   cluster gdp_pc life_exp rd_exp count
##   <fct>    <dbl>    <dbl>  <dbl> <int>
## 1 1        1900.     63.5  0.266    38
## 2 2       14789.     75.4  0.492    66
## 3 3       56990.     81.7  2.55     33

9. Wizualizacja klastrów na mapie

Sprawdzamy, czy klastry mają wyraźny wymiar geograficzny.

# Dane mapowe (granice państw)
world_sf <- ne_countries(scale = "medium", returnclass = "sf") %>%
  select(iso_a3, name, geometry)

# Połączenie wyników klastrowania z geometrią (po kodzie ISO3)
map_data <- world_sf %>%
  left_join(
    final_results %>%
      mutate(iso_a3 = countrycode::countrycode(country, "country.name", "iso3c")),
    by = "iso_a3"
  )

# Mapowanie klastrów (tmap)
tmap_mode("plot")
tm_shape(map_data) +
  tm_polygons("cluster", palette = "Set2", title = "Klaster",
              colorNA = "grey90", textNA = "Brak danych") +
  tm_layout(frame = FALSE, legend.outside = TRUE,
            main.title = "Klastry państw na mapie świata")

10. Wnioski

Zastosowane podejście pozwala na wyodrębnienie spójnych grup państw o podobnym profilu rozwoju. Standaryzacja i logarytmowanie zmiennych skośnych poprawiają porównywalność cech, a dobór k oparty na kilku kryteriach zwiększa stabilność wyniku. Przeprowadzona analiza pozwoliła na wyraźne wyodrębnienie grup państw o zbliżonym profilu społeczno-gospodarczym. Najbardziej rozwinięte gospodarki tworzą zwarty klaster charakteryzujący się nie tylko wysokim PKB, ale przede wszystkim ponadprzeciętnymi wskaźnikami jakości życia i edukacji. Z kolei grupy państw rozwijających się wykazują znacznie większą wariancję wewnątrzskupiskową, co sugeruje różne ścieżki transformacji ekonomicznej w tych regionach. Istotnym spostrzeżeniem jest fakt, że zmienne dotyczące infrastruktury cyfrowej coraz częściej stają się kluczowym czynnikiem dystynktywnym, spychając tradycyjne miary przemysłowe na dalszy plan. Wizualizacja kartograficzna wyników ujawnia silną regionalizację procesów rozwojowych, wskazując na istnienie wyraźnych centrów wzrostu oraz obszarów peryferyjnych. Uzyskana typologia może stanowić obiektywną podstawę do różnicowania strategii wsparcia międzynarodowego lub planowania inwestycji zagranicznych. W dalszym etapie warto rozszerzyć zestaw zmiennych (np. edukacja, nierówności) i porównać wyniki z innymi metodami (hierarchiczne, DBSCAN).