Uvod

U ovom seminarskom radu istražujemo razvoj i karakteristike balističkih raketa kroz analizu podataka prikupljenih iz ruskih i kineskih sistema naoružanja.

Pri rešavanju zadatih problema, koristićemo podatke iz dosutpnog fajla rakete.ods. U zadatoj bazi podataka imamo sledeća obeležja raketa (čije je značenje jasno po samom nazivu):

NAPOMENA: rusko ime koristimo samo kao klasičan tekstualni podatak pa nam ovo nije pravo kategoričko obeležje (u teoriji jeste, ali u praksi ga ne koristimo praveći bilo kakve klase rakete), dok je NATO ime kodno te preko njega možemo podeliti rakete na ruske i kineske (prefiks SS- za ruske, a CSS- za kineske, biće detaljno objašnjeno u 6. zadatku) pa stoga možemo reći da je u pitanju kategoričko obeležje koje pored navedenog značaja koristimo kao i rusko tj. kao klasičan tekstualni podatak

Zadate podatke koristićemo u cilju analiziranja kako pomenuta obeležja utiču jedna na druga, da li to ima ikakve veze sa zemljom porekla raketa ili pak njihovom starošću, da li možemo da uočimo nekakve značajne pravilnosti u podacima, podelimo ih ili predvidimo. Naš pristup će biti matematički u prvih 5 zadataka, dok ćemo se u 6. zadatku orijentisati na istorijski i politički pristup zadatoj bazi.


Prvi zadatak (učitavanje podataka)

Naš prvi zadatak je zapravo priprema podataka, tačnije učitavanje istih u .ods formatu, te pretvaranje u data.frame.

Koristićemo paket “readODS” da ih učitamo, i funkciju “read_ods” iz tog paketa da ih pretvorimo u data.frame.

Dajemo i primer naših podataka u obliku tabele prvih 5 opservacija. 

#install.packages("readODS")
library(readODS)
library(knitr)

rakete_pocetna_baza = read_ods("C:\\Users\\sarac\\Downloads\\rakete.ods", 
                               range = NULL, col_names = TRUE)

attach(rakete_pocetna_baza)
kable(head(rakete_pocetna_baza, n = 5), 
      format = "html", 
      table.attr = "style='font-size: 12px; width: 70%'"
)
rusko_ime NATO_ime masa_kg duzina_m kalibar_m broj_podglava payload_kg tip_goriva max_domet_km nacin_lansiranja broj_faza_motora klasa_po_dometu
R-36M2 Voevoda SS-18 Satan 209600 32.2 3.05 10 8800 Tecno 16000 Silos 2 L
UR-100N SS-19 Stiletto 105600 27.0 2.50 6 3355 Tecno 10000 Silos 2 L
RT-2PM Topol SS-25 Sickle 45100 21.5 1.80 1 1000 Cvrsto 11000 RM TEL 3 L
RT-2PM2 Topol-M SS-27 Mod 1 47200 22.7 1.90 1 1200 Cvrsto 11000 RM TEL 3 L
RS-24 Yars SS-27 Mod 2 49600 23.0 2.00 3 1200 Cvrsto 12000 RM TEL 3 L



Priprema podataka za ostale zadatke

Za dalji rad (uključivaće pravljenje statističkih modela, razne podele podataka i funkcije primenjivane nad njima) biće nam potrebni unapred spremljeni podaci što podrazumeva:

rakete = rakete_pocetna_baza[, -c(1, 2)]

numerical_col = sapply(rakete, is.numeric)
rakete[numerical_col] = scale(rakete[numerical_col])

tip_goriva = factor(tip_goriva, levels = c("Cvrsto", "Tecno"), labels = c(1,2))
nacin_lansiranja = factor(nacin_lansiranja, levels = c("RM TEL", "Silos", "Mornaricki", "Avion"), labels = c(1,2,3,4))
klasa_po_dometu = factor (klasa_po_dometu, levels = c("S", "M", "L"), labels = c(1,2,3))

rakete$tip_goriva = as.numeric(tip_goriva)
rakete$nacin_lansiranja = as.numeric(nacin_lansiranja)
rakete$klasa_po_dometu = as.numeric(klasa_po_dometu)


Pogledajmo kako sada izgleda naša baza na primeru prvih 5 opservacija.

kable(head(rakete, n = 5), 
      format = "html", 
      table.attr = "style='font-size: 12px; width: 70%'"
)
masa_kg duzina_m kalibar_m broj_podglava payload_kg tip_goriva max_domet_km nacin_lansiranja broj_faza_motora klasa_po_dometu
3.1509705 2.2001557 2.0859515 1.5184599 3.4927289 2 2.0273365 2 0.0916496 3
1.2170642 1.5490973 1.3247589 0.5363071 0.8048411 2 0.8804534 2 0.0916496 3
0.0920514 0.8604778 0.3559684 -0.6913839 -0.3576889 1 1.0716006 1 1.2830942 3
0.1311014 1.0107221 0.4943671 -0.6913839 -0.2589602 1 1.0716006 1 1.2830942 3
0.1757300 1.0482831 0.6327657 -0.2003075 -0.2589602 1 1.2627478 1 1.2830942 3


Drugi zadatak (popunjavanje NA vrednosti)

U drugom zadatku cilj nam je popuniti NA vrednosti. Prvo pogledajmo koja obeležja ih imaju i koliko ih je, kako bi se odlučili za najsmilseniji način popunjavanja istih.

kable(colSums(is.na(rakete)), 
      format = "html",
      table.attr = "style='font-size: 15px; width: 70%'")
x
masa_kg 1
duzina_m 1
kalibar_m 0
broj_podglava 1
payload_kg 0
tip_goriva 0
max_domet_km 0
nacin_lansiranja 0
broj_faza_motora 0
klasa_po_dometu 0


Vidimo da NA vrednosti nema puno tj. postoji samo po jedna za obeležja masa_kg, duzina_m, broj_podglava.

Ako pogledamo sve varijable koje imamo, intuitivno nam deluje kao da su sve zavisne u manjoj ili većoj meri. Stoga ćemo ove vrednosti popuniti na sledeći način: koristićemo funkciju “mice” iz paketa “mice” koja popunjava ove vrednosti koristeći veze sa ostalim dostupnim poljima svih varijabli i na taj način nam omogućava realistično popunjavanje NA vrednosti.

NAPOMENA: sva obeležja koja sadrže NA vrednosti su numerička i već skalirana u pripremi podataka, ali nam ovo neće praviti problem pri imputaciji, samo je važno napomenuti da ni rezultati neće biti u stvarnom opsegu vrednosti. 

#install.packages("mice")
library(mice)

imputacija = mice(rakete, m = 5, method = 'pmm', seed = 500)
rakete = complete(imputacija, 1) # biramo prvu imputaciju


Treći zadatak (procena max dometa)

U trećem zadatku želimo da predvidimo maksimalan domet rakete (dakle to će biti naša ciljna promenljiva) na osnovu drugih prediktora, modelom po izboru. Nama su poznate vrednosti obeležja max_domet, ali praivćemo se da ih ne znamo, predvidićemo ih, a onda naše predviđene vrednosti uporediti sa stvarnim da bi videli da li smo uradili dobar posao.


Podela podataka

U procesu formiranja modela možemo koristiti sledeće skupove za podelu podataka:

1. Trening skup (na njemu se model “uči”)
2. Validacioni skup (međukorak za procenu performansi modela na neviđenim podacima i podešavanje hiperparametara modela, nije uvek neophodan)
3. Test skup (završna evaluacija tj. procena performansi modela na neviđenim podacima)

Vizuelizacija podele podataka


Za pouzdaniju procenu performansi modela, umesto posebnog validacionog skupa može se koristiti kros-validacija (unakrsna validacija). Unakrsna validacija sa k grupa podrazumeva sledeće korake:

  1. Podaci se dele na k foldova
  2. Za svaki od k koraka, jedan fold se koristi kao test skup, a preostalih k-1 kao trening skup
  3. Prosečna performansa modela na k test skupova daje konačnu ocenu perfromansi (dakle ovime smo evaluirali naš model)

U ovom zadatku podatke ćemo deliti samo na trening i test skup (bez validacionog). Razlog za ovo je jednostavan pristup problemu, nemamo puno podataka a modeli koje planiramo da koristimo (SVM i linearna regresija) su dovoljno jednostavni modeli da nam optimizacija hiperparametara (koje SVM ima) nije neophodna.

Kako je baza mala, nećemo rizikovati potpuno nasumičnom podelom, već ćemo koristiti funkciju “createDataPartition” iz paketa “caret”. Ona podelu vrši nasumično ali obezbeđuje da se zadrži približno ista distribucija ciljne promenljive u oba skupa.

library(caret)

set.seed(123)
train_indexes = createDataPartition(rakete$max_domet_km, p = 0.7, list = F) 
train_set = rakete[train_indexes, ]
test_set = rakete[-train_indexes, ]

Sada kada imamo podelu na trening i test skup možemo se odlučiti za prvi model koji ćemo isprobati.


Linearna regresija

Prvo koristimo tradicionalni model linearne regresije, kao najjednostavniji pristup. On je često, pa tako i sada, prvi intuitivni izbor, kada su očigledne direktne (i po našoj proceni linearne) veze među obeležjima (npr. u našem slučaju max_domet_km ima jasnu direktnu vezu sa masom, dužinom itd.).


Vizuelizacija lm modela


Sada nam je cilj da proverimo koliko je svaki od prediktora značajan za naš model i predviđanje ciljne promenljive max_domet. Proveru ćemo uraditi pregledom p-vrednosti za svaki od prediktora na trening skupu, jer p-vrednosti na test skupu nisu relevantne za treniranje modela niti za procenu značaja, a test skup zadržavamo kao potpuno neviđeni za validaciju performansi modela.

summary (lm(max_domet_km ~., data = train_set))
## 
## Call:
## lm(formula = max_domet_km ~ ., data = train_set)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49455 -0.09629  0.00787  0.16668  0.64291 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)      -0.27931    0.41787  -0.668   0.5123  
## masa_kg          -0.23053    0.50396  -0.457   0.6528  
## duzina_m          0.24045    0.24668   0.975   0.3426  
## kalibar_m         0.28976    0.28306   1.024   0.3196  
## broj_podglava    -0.02212    0.13764  -0.161   0.8741  
## payload_kg        0.42993    0.34230   1.256   0.2252  
## tip_goriva       -0.37295    0.19439  -1.919   0.0710 .
## nacin_lansiranja -0.06711    0.07262  -0.924   0.3676  
## broj_faza_motora  0.10651    0.11349   0.939   0.3604  
## klasa_po_dometu   0.43676    0.17097   2.555   0.0199 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3128 on 18 degrees of freedom
## Multiple R-squared:  0.9386, Adjusted R-squared:  0.9079 
## F-statistic: 30.59 on 9 and 18 DF,  p-value: 4.102e-09
library(car)
vif(lm(max_domet_km ~., data = train_set))
##          masa_kg         duzina_m        kalibar_m    broj_podglava 
##        68.790177        15.962495        20.337554         5.821232 
##       payload_kg       tip_goriva nacin_lansiranja broj_faza_motora 
##        41.464484         2.647774         1.501146         3.736724 
##  klasa_po_dometu 
##         6.187825

Vidimo da su p-vrednosti većine prediktora visoke, osim za tip_goriva i klasa_po_dometu, što bi moglo sugerisati da su ostali prediktori manje značajni.

Međutim, s obzirom na to da smo uvideli jaku korelisanost među prediktorima (multikolinearnost) koristeći funkciju “vif” i uočavajući njene pretežno visoke vrednosti, zaključujemo da visoke p-vrednosti mogu biti rezultat te međusobne povezanosti. Model tada teško razlikuje doprinos pojedinačnih prediktora, što može „prigušiti“ njihove statističke značajnosti.

Zbog toga ćemo se voditi intuicijom, smatrajući da svaki od prediktora može imati, manji ili veći, značaj za predviđanje max_domet, te ćemo ih sve zadržati u modelu. Ovakav pristup ne bi trebalo da mnogo naškodi našem modelu te nećemo rizikovati izbacivanjem prediktora.

lm_model = lm(max_domet_km ~., data = train_set)
lm_predictions = predict(lm_model, newdata = test_set)
lm_predictions
##          2          5         11         18         21         22         24 
##  0.9705408  1.0155622  0.5995547 -1.0883919  0.3406234 -0.9769980 -1.1350237 
##         30         31         33         35 
## -0.4461397 -0.1986037  0.1945715  1.2633674

Vidimo da smo zaista dobili predikcije maksimalnog dometa za 15 raketa koje se nalaze u test skupu. Pre nego što se upustimo u razmatranje koliko su ove predikcije zaista dobre, pokušajmo da uradimo isto drugim, nešto kompleksnijim modelom: Random Forest


Random Forest

Naš drugi pristup biće Random Forest. On je skup velikog broja stabala odlučivanja (decision trees) koje zajedno donose konacnu odluku.
Biramo ga jer je jak model koji može da uhvati složenije odnose, a ipak nije komplikovan i ne zahteva nužno korišćenje validacionog skupa.

Vizuelizacija RF modela


library(randomForest)

rf_model = randomForest(max_domet_km ~ ., data = train_set, ntree = 100)
rf_predictions = predict(rf_model, newdata = test_set)
rf_predictions
##           2           5          11          18          21          22 
##  0.94888566  0.97089657  0.47360932 -0.98102477  0.50966764 -0.92652807 
##          24          30          31          33          35 
## -0.85104518 -0.76670817 -0.42529189  0.02011261  1.71177480


Poređenje modela

Sada kada imamo dva različita modela koja imaju isti cilj možemo ih uporediti nekim od metrika koje se često koriste za procenu i upoređivanje modela (MSE, MAE, RMSE, \(R^2\)). Bolji je onaj model koji ima manju grešku i/ili veći koeficijent determinacije.

rmse_rf = sqrt(mean((test_set$max_domet_km - rf_predictions)^2))
rmse_lm = sqrt(mean((test_set$max_domet_km - lm_predictions)^2))

print(c(rmse_rf, rmse_lm))
## [1] 0.2692947 0.3245652


mae_rf = mean(abs(rf_predictions - test_set$max_domet_km))
mae_lm = mean(abs(lm_predictions - test_set$max_domet_km))

print(c(mae_rf, mae_lm))
## [1] 0.2094464 0.2576249


mae_rf = mean(abs(test_set$max_domet_km - rf_predictions))
mae_lm = mean(abs(test_set$max_domet_km - lm_predictions))

print(c(mae_rf, mae_lm))
## [1] 0.2094464 0.2576249


r2_rf = cor(test_set$max_domet_km, rf_predictions)^2
r2_lm = cor(test_set$max_domet_km, lm_predictions)^2

print(c(r2_rf, r2_lm))
## [1] 0.9372715 0.8946334

Nešto bolji je Random Forest model, mada su oba sličnih i ujedno dovoljno dobrih performansi.

Intuitivno ne uočavamo mnogo nelinearnih odnosa u podacima, pa je zato linearna regresija bila dobar izbor koji je imao i solidne performanse, jer je upravo za ovakve slučajeve i dizajnirana. Zato Random Forest nema izraženu prednost (mada ona postoji) nad linearnom regresijom. Tako je sasvim očekivano da su rezultati slični, a da je Random Forest eventualno malo bolji zbog svoje fleksibilnosti.

Ove modele možemo predstaviti i grafički. Očekujemo isti zaključak.

par(mfrow = c(1, 2))

plot(test_set$max_domet_km, rf_predictions, main = "RF: Stvarne vs Predviđene", 
     xlab = "Stvarne vrednosti", ylab = "Predviđene vrednosti", col = "blue")
abline(0, 1, col = "red")

plot(test_set$max_domet_km, lm_predictions, main = "LM: Stvarne vs Predviđene", 
     xlab = "Stvarne vrednosti", ylab = "Predviđene vrednosti", col = "darkgreen")
abline(0, 1, col = "red")



Četvrti zadatak (klasifikator pokretnosti lansera)

U četvrtom zadatku naš cilj je, za početak, napraviti klasifikator (model kojim se podaci svrstavaju u jednu od kategorija) koji pogađa da li raketa ima pokretan ili fiksan lanser.

Nakon toga treba da uočimo overfitting ukoliko je do njega došlo i u tom slučaju modifikujemo naš model.

Naše obeležje način_lansiranja uzima 4 vrednosti:

Mornarički tip lansera može biti i fiksan i pokretan u zavisnosti od toga da li se koristi u cilju lansiranja rakete sa broda/podmornice u pokretu, ili samo kako bi se došlo do udaljene, nekopnene tačke za stacionarno lansiranje. Kako bi olakšali dalji radi, opredelimo se za jedan, npr. fiksan tip mornaričkog lansera.

Pokretan Fiksan
Silos X
RM TEL X
Avion X
Mornarički X


Klasifikacija može biti višeklasna ili binarna. Stoga postoje dva pristupa našem problemu:

Radi lakše implementacije, opredelićemo se za drugi pristup, binarnu klasifikaciju.


Ciljna promenljiva

Na gore objašnjen način prvo formiramo novo obeležje koje ćemo nazvati tip_lansiranja i koje će ujedno biti i naša ciljna promenljiva. Pretvorićemo ga u faktor kako bi njegovo značenje bilo očigledno.
Nakon toga pogledajmo našu bazu na primeru prvih 5 opservacija nakon dodavanja kolone tip_lansiranja i naravno uklanjanja nacin_lansiranja s obzirom da nam on više nije potreban (sve što želimo da znamo o lanserima je trenutno to da li su fiksni ili pokretni, a znanje o tome nam je sačuvano u varijabli tip_lansiranja). Ipak, sačuvaćemo ga u jednoj promenljivoj kako bi kasnije mogli da vratimo bazu na staro.

rakete$tip_lansiranja = ifelse(rakete$nacin_lansiranja %in% c(2, 3), 0, 
                               ifelse(rakete$nacin_lansiranja %in% c(1, 4), 1, 0))
nacin_lansiranja_kopija = rakete$nacin_lansiranja
rakete = rakete[ ,-8]
kable(head(rakete, n = 5), 
      format = "html",
      table.attr = "style='font-size: 12px; width: 70%'")
masa_kg duzina_m kalibar_m broj_podglava payload_kg tip_goriva max_domet_km broj_faza_motora klasa_po_dometu tip_lansiranja
3.1509705 2.2001557 2.0859515 1.5184599 3.4927289 2 2.0273365 0.0916496 3 0
1.2170642 1.5490973 1.3247589 0.5363071 0.8048411 2 0.8804534 0.0916496 3 0
0.0920514 0.8604778 0.3559684 -0.6913839 -0.3576889 1 1.0716006 1.2830942 3 1
0.1311014 1.0107221 0.4943671 -0.6913839 -0.2589602 1 1.0716006 1.2830942 3 1
0.1757300 1.0482831 0.6327657 -0.2003075 -0.2589602 1 1.2627478 1.2830942 3 1

Sada kada smo odlučili kako ćemo pristupiti problemu i formirali našu novu ciljnu promenljivu, treba da podelimo naše podatke pre nego što počnemo pravljenje modela za koji se odlučimo.


Podela podataka

Za razliku od prošlog zadatka, ovde ćemo podatke podeliti u 3 skupa: trening(50%), validacioni(20%) i test(30%). Za uvođenje još i validacionog skupa opredeljujemo se jer nam on omogućava lakše prepoznavanje overfittinga (što nam je ovde potrebno) kao i modifikaciju modela u tom slučaju.

Takođe ćemo razdvojiti prediktore i ciljnu promenljivu radi lakše organizacije i fleksibilnosti, a i ukoliko dođe do overfittinga koristićemo funkciju “glmnet” za koju su nam potrebni ovako razdvojeni podaci.

validation_size = 0.2
test_size = 0.3
train_size = 0.5

set.seed(70)

train_val_idx = createDataPartition(rakete$tip_lansiranja, p = train_size + validation_size, 
                                    list = FALSE)
train_val_data = rakete[train_val_idx,]
test_data = rakete[-train_val_idx,]

x_train_val = train_val_data[, -10]
y_train_val = train_val_data[, 10]
x_test = test_data[, -10]
y_test = test_data[, 10]

set.seed(70)

train_idx = createDataPartition(train_val_data$tip_lansiranja, 
                                p = train_size/(validation_size + train_size), 
                                list = FALSE)
train_data = train_val_data[train_idx,]
val_data = train_val_data[-train_idx,]

x_train = train_data[,-10]
y_train = train_data[,10]
x_val = val_data[,-10]
y_val = val_data[,10]


Logistička regresija

Sada konačno možemo preći na pravljenje modela za koji smo se odlučili.

Logistička regresija je jedan od najjednostavnijih i najkorišćenijih klasifikacionih modela. Koristimo je kada je naša ciljna promenljiva kategoričkog tipa, kao u ovom slučaju.

logit_model = glm(tip_lansiranja~., data = train_data, family = "binomial")
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred


Da li je došlo do overfittinga?

Dobijamo upozorenja da je model predvideo verovatnoće jednake 0 ili 1. To znači da je model “previše siguran u svoje odluke” i ne prepoznaje nijanse u podacima tj. previše je dobro naučio o podacima trening skupa. Zato je ovo alarm da je verovatno došlo do overfittinga (preprilagođavanja) - pojava da se nauče obrasci sa trening skupa, ali se ne ispoljavaju dobro na novim podacima.

Proverimo na još neki način da li je ovo zaista slučaj.

Izračunaćemo tačnost modela na skupu podataka na kom je treniran (train_data) i na validacionom skupu koji nije još viđen (val_data). Koristićemo funkciju “mean”.

NAPOMENA: nećemo trenirati model na trening + validacija skupu i gledati njegove performanse na uniji treninga i validacije, pa na test skupu odvojeno da bi uočili overfitting (što bismo verovatno uspeli i tim putem). U tom slučaju bi koristili test skup da odlučimo koliko je naš model dobar ili nije, a želimo da ga sačuvamo kao konačni, neviđeni skup podataka za finalnu evaluaciju 

train_pred = predict(logit_model, newdata = train_data, type = "response")
train_pred = as.numeric(train_pred > 0.5)
train_acc = mean(train_pred == train_data$tip_lansiranja)
train_acc
## [1] 1
val_pred = predict(logit_model, newdata = val_data, type = "response")
val_pred = as.numeric(val_pred > 0.5)
val_acc = mean(val_pred == val_data$tip_lansiranja)
val_acc
## [1] 0.75

Vidimo da je tačnost modela na skupu na kom je treniran baš 1. Iako to na prvi pogled deluje logično (model savršeno predviđa podatke jer ih je već video) to nam zapravo ne odgovara. Cilj nije bio da model savršeno nauči podatke već da ih razume, da predvodi inteligencijom a ne memorijom. Stoga je nama zapravo pouzdaniji onaj model koji ima tačnost npr. 0.9 na skupu na kom je treniran a ne 1.

Dodatno, model je lošijih performansi na skupu neviđenih podataka tj. na validacionom skupu i ovaj zaključak ide u prilog tezi da naš model nije dubinski razumeo podatke i njihovu vezu već ih je samo naučio napamet.
Konačno, možemo zaključiti da je došlo do overfittinga.

Zbog svega navedenog odlučujemo se da model poboljšamo koliko možemo. Rešavamo overfitting regularizacijom što često ujedno dovodi i do smanjenja overconfidence-a.


Regularizacija

Overfitting ćemo prebroditi regularizacijom - tehnikom za smanjenje složenosti modela. Postoje 3 tipa regularizacije (u zavisnosti od argumenta alpha):

  • L1 (Lasso)
  • L2 (Ridge) - mi ćemo je koristiti
  • Elastic Net (kombinacija)

Važno je napomenuti da funkcija za pravljenje regularizovanog modela prima kao argument parametar regularizacije lambda. Da bi dobili što bolji model, ovaj parametar potrebno je optimizovati. Zato, umesto da pristupimo rešavanju problema overfittinga funkcijom “glmnet” (trenira model i primenjuje regularizaciju), primenićemo funkciju “cv.glmnet” (dodatno uključuje i unakrsnu validaciju na uniji trening skupa i validacionog skupa za pronalaženje najbolje vrednosti parametra lambda).

Postupak:

  1. Unakrsna validacija i optimizacija lambda (na validacionom skupu) korišćenjem funkcije “cv.glmnet” (Napomena: prosleđujemo funkciji eksplicitno opseg iz kog da traži optimalno lambda, nazvali smo ga lambdas i generisali ga kao jedan od najpoznatijih, dovoljno širokih opsega za ovaj postupak)
  2. Pravimo finalni, regularizovani model sa najboljim lambda iz 1. dela
  3. Taj model treniramo (na uniji trening i validacionog skupa)
  4. Evaluiramo model (na test skupu)
library(glmnet)
lambdas = exp(seq(-10,3, length.out = 100))

regularized_model_cv_L2 = cv.glmnet(as.matrix(x_train_val), as.matrix(y_train_val), 
                                 family="binomial", 
                                 type.measure = "class",
                                 standardize = FALSE,
                                 nfolds = 10, lambda = lambdas)
plot(regularized_model_cv_L2)

Na grafiku vidimo dve važne linije: prva nam predstavlja lambda za koje je greška (tj Binomial Deviance - metrika greške kod klasifikatora) najmanja, a druga nam predstavlja lambda sa 1 standardnom greškom.

Izabraćemo ono lambda koje minimizuje funkciju gubitka i napraviti i evaluirati naš finalni model na test skupu.

optimal_lambda = regularized_model_cv_L2$lambda.min

final_model_L2 = glmnet(x_train_val, y_train_val, 
                        lambda = optimal_lambda, family = "binomial", standardize = FALSE)
y_pred_reg_test = predict(final_model_L2, s = optimal_lambda, newx = as.matrix(x_test), type = "class")
test_acc_reg = mean(y_pred_reg_test == y_test)
test_acc_reg 
## [1] 0.6363636
y_pred_reg_val = predict(final_model_L2, s = optimal_lambda, newx = as.matrix(x_val), type = "class")
val_acc_reg = mean(y_pred_reg_val == y_val)
c(val_acc_reg, val_acc)
## [1] 0.875 0.750
y_pred_reg_train = predict(final_model_L2, s = optimal_lambda, newx = as.matrix(x_train), type = "class")
train_acc_reg = mean(y_pred_reg_train == y_train)
c(train_acc_reg, train_acc)
## [1] 0.95 1.00

Uspeli smo da smanjimo prilagođenost modela na skupu na kom je treniran i time smatramo da smo uspešno izvršili regularizaciju.

Sada kada konačno imamo naš finalni, regularizovan model, možemo pogledati šta “pogađa” koji tip lansera imaju naše rakete i to uporediti sa stvarnim vrednostima.

Takođe, možemo pogledati i kakve su performanse našeg modela preko nekih osnovnih metrika.

library("caret")

conf_matrix = confusionMatrix(factor(test_data$tip_lansiranja), factor(y_pred_reg_test))
conf_matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 0 1
##          0 3 3
##          1 1 4
##                                           
##                Accuracy : 0.6364          
##                  95% CI : (0.3079, 0.8907)
##     No Information Rate : 0.6364          
##     P-Value [Acc > NIR] : 0.6322          
##                                           
##                   Kappa : 0.2903          
##                                           
##  Mcnemar's Test P-Value : 0.6171          
##                                           
##             Sensitivity : 0.7500          
##             Specificity : 0.5714          
##          Pos Pred Value : 0.5000          
##          Neg Pred Value : 0.8000          
##              Prevalence : 0.3636          
##          Detection Rate : 0.2727          
##    Detection Prevalence : 0.5455          
##       Balanced Accuracy : 0.6607          
##                                           
##        'Positive' Class : 0               
## 
library(pROC)

roc_curve = roc(as.numeric(test_data$tip_lansiranja), as.numeric(y_pred_reg_test))
plot(roc_curve, main = "Roc kriva", col = "red", lwd = 3)

auc_value = auc(roc_curve)
text(0, 0.4, paste("AUC = ", round(auc_value, 2)), col = "navy", cex = 1.5)


Dodatak - unakrsna validacija

Kako tačnost na test skupu nakon evaluacije modela nije bila zadovoljavajuća (0.6), odlučili smo da primenom unakrsne validacije steknemo precizniji uvid u stvarne performanse naše regularizacije na nevidjenim podacima. Pretpostavljamo da je relativno niska početna tačnost mogla biti posledica malog broja opservacija u skupu podataka.

Koristićemo ugneždenu unakrsnu validaciju (nested cross-validation), pri čemu spoljašnja petlja svaki put izdvaja jedan fold kao test skup, dok se preostalih k - 1 folda koristi za trening i validaciju. Unutrašnja petlja zatim na tom trening + validacionom skupu izvodi unakrsnu validaciju pomoću funkcije “cv.glmnet”, kojom se bira optimalna vrednost regularizacionog parametra lambda. Performanse modela sa izabranim lambda na spoljašnjem test skupu beleže se u svakoj iteraciji, a najrealnija procena konačnih performansi dobija se kao prosečna tačnost na svim spoljašnjim test skupovima.

set.seed(123)

X = as.matrix(rakete[, -10])   # pretpostavljam da je 10. kolona ciljna
y = rakete[, 10]

k_outer = 10
outer_folds = createFolds(y, k = k_outer, list = TRUE)
lambdas = exp(seq(-10, 3, length.out = 100))

outer_accuracies = numeric(k_outer)

for (i in 1:k_outer) {
  
  test_idx = outer_folds[[i]]
  X_test = X[test_idx, ]
  y_test = y[test_idx]
  
  X_train_val = X[-test_idx, ]
  y_train_val = y[-test_idx]
  
  cv_model = cv.glmnet(X_train_val, y_train_val,
                        family = "binomial",
                        type.measure = "class",
                        standardize = FALSE,
                        nfolds =10,
                        lambda = lambdas,
                        alpha = 0)
  
  best_lambda = cv_model$lambda.min
  
  final_model = glmnet(X_train_val, y_train_val,
                        family = "binomial",
                        lambda = best_lambda,
                        alpha = 0,
                        standardize = FALSE)
  
  preds_test = predict(final_model, newx = X_test, type = "class", s = best_lambda)
  outer_accuracies[i] = mean(preds_test == y_test)
}

mean(outer_accuracies)
## [1] 0.775

Sada vidimo da je naša regularizacija bila uspešna - ne samo da smo prebrodili overfitting već model dobijen istom ima sasvim solidne performanse.

NAPOMENA: još jednom naglašavamo da je unakrsna validacija korišćena isključivo za dobijanje realnije procene performansi regularizovanog modela, uz maksimalno iskorišćenje dostupnih opservacija. Ona nije služila za dobijanje konačnog modela — u praktičnoj primeni koristili bismo inicijalni model sa nešto nižom tačnošću. Međutim, unakrsna validacija nam je pokazala da je zbog ograničenog broja opservacija prvobitno izmerena tačnost bila „lažno niska“

NAPOMENA: nećemo se sada detaljno baviti metrikama dobijenim unakrsnom validacijom, dovoljno nam je bilo da znamo da je naš model zadovoljavajući (ako bi hteli morali bi u svakoj iteraciji unutrašnje petlje zabeležiti i sve metrike dobijene na trenutnom test skupu, izračunati i njihovu srednju vrednost, pa prikazati takve)


Sada ćemo izbaciti prediktor koji smo sami pravili, tip_lansiranja i vratiti prediktor način_lansiranja. Vraćamo bazu na staro jer nam nije više potrebna modifikovana pri rešavanju novih problema.

rakete = subset(rakete, select = -c(tip_lansiranja))
rakete$nacin_lansiranja = nacin_lansiranja_kopija


Peti zadatak(klasterizacija)

U petom zdatku treba izbaciti sve rakete srednjeg dometa kao i prediktore max_domet i klasa_po_dometu, a zatim bez ikakvog znanja o daljem dometu raketa (sem što znamo da sada imamo samo one kratkog i dugog dometa) izvršiti klasterizaciju po istom.


Izbacivanje raketa srednjeg dometa i prediktora max_domet i klasa_po_dometu

Prediktor klasa_po_dometu uzima vrednosti S, M i L (koje smo mi kodirali redom u 1, 2 i 3). Pretpostavka je da ove vrednosti zapravo predstavljaju Small, Medium i Large domet, te da je dovoljno izbaciti iz baze rakete čije obeležje klasa_po_dometu uzima vrednost M (2 kod nas). Ipak, kako nemamo jasno utvrđen smisao ovog obeležja, možemo proveriti da li smo u pravu.

Proveravamo da li postoje jasne kategorije podataka tj. da li zaista sve rakete klase S imaju mali max_domet, one klase M osrednji i one klase L najveći. Koristićemo se “anova” testom za proveru.

  • H0: Nema razlike u srednjim vrednostima zavisne varijable (max_domet) između grupa.
  • H1: Postoji barem jedna grupa sa srednjom vrednošću koja se značajno razlikuje.
anova_model = aov(max_domet_km ~ as.factor(klasa_po_dometu), data = rakete)
summary(anova_model)
##                            Df Sum Sq Mean Sq F value   Pr(>F)    
## as.factor(klasa_po_dometu)  2 28.621  14.311   54.93 1.15e-11 ***
## Residuals                  36  9.379   0.261                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kako je p-vrednost izuzetno mala, to odbacujemo nultu hipotezu i mozemo zaključiti da je mala verovatnoća da su razlike u srednjim vrednostima zavisne varijable (max_domet) između grupa nezavisne varijable (L, S, M) nastale slučajno.

Pogledajmo kako ovo izgleda i grafički.

boxplot(max_domet_km ~ klasa_po_dometu, data = rakete,
        main = "Distribucija max_domet_km po klasa_po_dometu",
        xlab = "Klasa po dometu", ylab = "Max domet", col = "red")

Odavde se jasno vidi da SKORO sve rakete klase S imaju maksimalan domet veoma mali, one klase M srednji, a one klase L najveći. Vizuelno deluje kao da postoje neke opservacije koje mogu poremetiti naš zaključak o značenju klasa obeležja klasa_po_dometu (postoji preplitanje max_domet -a klase M i L iz klasa_po_dometu).

Proverimo ovo numerički.

opseg_S = c(min(subset(rakete, klasa_po_dometu == 1)$max_domet), max(subset(rakete, klasa_po_dometu == 1)$max_domet))
opseg_M = c(min(subset(rakete, klasa_po_dometu == 2)$max_domet), max(subset(rakete, klasa_po_dometu == 2)$max_domet))
opseg_L = c(min(subset(rakete, klasa_po_dometu == 3)$max_domet), max(subset(rakete, klasa_po_dometu == 3)$max_domet))

opseg_S
## [1] -1.0176382 -0.9354449
opseg_M
## [1] -0.91633023 -0.07528259
opseg_L
## [1] -0.1708562  2.4096309
sort(subset(rakete, klasa_po_dometu == 3)$max_domet)
##  [1] -0.17085618  0.07763517  0.21143820  0.21143820  0.49815899  0.55550314
##  [7]  0.55550314  0.55550314  0.68930618  0.88045337  1.07160056  1.07160056
## [13]  1.20540359  1.26274775  1.45389494  1.83618933  2.02733652  2.40963090


NAPOMENA: prisetimo se da su ove vrednosti skalirane, mada nam to ne pravi problem pri zaključivanju jer nam nisu važne realne brojke o kojima pričamo, već samo opseg svake klase i njihovo preplitanje.

Dakle, pravi nam problem što postoji jedna raketa klase L ciji je max_domet -0.17085618 (što ulazi u opseg za klasu M, jer -0.17085618 < -0.07528259 što je gornja granica opsega klase M).

Kako postoji samo jedno odstupanje od pravila koje smo uvideli, uvažićemo ga kao ispravno i zaključiti da su rakete srednjeg dometa one koje pripadaju klasi “M” obeležja klasa_po_dometu.

Konačno, upravo ovakve rakete izbacujemo, a nakon toga i prediktore klasa_po_dometu i max_domet.

indeksi_klasa_M = which(klasa_po_dometu == 2)
rakete_izbaceno = rakete[-indeksi_klasa_M, ]
rakete_izbaceno = subset(rakete_izbaceno, select = -c(max_domet_km, klasa_po_dometu))


Dalje je plan je da podacima pristupimo na 2 načina tj. da se koristimo sa 2 tipa klasterizacije. Nakon toga ćemo ih uporediti kako bi zaključili nešto više o pravom pristupu podacima i opredelili se za bolju klasterizaciju.


K - means klasterizacija

1. Odabir broja klastera K

Prvo, treba se oprediliti za broj klastera prema kom ćemo izvršiti podelu. Naš cilj je da uočimo pravilnost u dobijenim grupama, te da zaključimo koji klaster/klasteri predstavljaju rakete kratkog, a koji dugog dometa.

NAPOMENA: Zahtev za klasterizaciju prema dometu koji uzima 2 vrednosti (kratki i dugi) ne znači nužno da su nam potrebna samo 2 klastera. Moguće je da se desi da je najoptimalniji broj klastera veći i da je deo dobijenih klastera jedna “celina” koja predstavlja rakete kratkog, a drugi deo dugog dometa. Do ovoga može doći usled specifičnosti podataka ili nekih podgrupa koje postoje unutar raketa jednog određenog dometa.

Postoji nekoliko metoda za odabir broja K:

  • Elbow metoda
    Y - osa: suma kvadrata unutar klastera - “WCSS” (pokazuje koliko su podaci kompaktni unutar klastera za svako K)
    Odabir K: Mesto gde se kriva “lomi” (smanjenje WCSS se usporava) tj. dodavanje više klastera prestaje da značajno smanjuje WCSS
  • Silhouette analiza
    Y - osa: Silhouette koeficijent koji ocenjuje koliko su klasteri kompaktni i razdvojeni
    Odabir K: što veći koeficijent, ako je više K sličnih vrednosti koeficijenata, izabrati najmanje radi jednostavnosti
  • Gap Statistic
    Y - osa: Gap vrednost (razlika između WCSS za stvarne i nasumične podatke), pokazuje značajnost klasterizacije
    Odabir K: Prvo K za koje vrednost naglo poraste i postaje stabilna


2. K - means algoritam

Najbolji broj klastera K biraćemo jednom od navedenih metoda pa neka to bude i najpoznatija - Elbow metoda. Koristićemo se sa dva algoritma klasterizacije i odabrati nama intuitivniji. Za početak će to biti K - means algoritam.

K - means algoritam funkcioniše tako što se prvo (nasumično) odabere K centara klastera (“centroida”). Zatim se svaka opservacija dodeljuje klasteru čijoj je centroidi najbliža. Tada se računa centar novodobijenog klastera i u njega se premešta centroida.

library(factoextra)
library(cluster)
library(gridExtra)

elbow_plot = fviz_nbclust(rakete_izbaceno, kmeans, method = "wss") + 
  ggtitle("Elbow metoda - Odabir K")

elbow_plot

Za K = 3 kriva se “lomi” te je to lakat koji smo tražili. Za K = 8 dolazi do blagog skoka, ali ovo nije ništa značajno što može da poremeti naš zaključak već je najverovatnije samo posledica specifičnosti nekog podatka. 

set.seed(42)

final_model_cluster_kmeans = kmeans(rakete_izbaceno, centers = 3, nstart = 10)

# dodeljeni klasteri svakoj opservaciji
clusters_kmeans = final_model_cluster_kmeans$cluster
clusters_kmeans
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 17 18 20 22 23 24 27 28 29 34 35 
##  3  2  2  2  2  2  2  2  2  2  2  3  2  1  1  1  1  2  1  1  1  2  1  1  2  3 
## 36 38 39 
##  2  2  1
# prikaz koliko opservacija je svrstano u svaki od 3 klastera koje imamo
table(clusters_kmeans)
## clusters_kmeans
##  1  2  3 
## 10 16  3


Hierarchical Clustering

Hijerarhijsko klasterovanje je dobar alat za proveru i dodatno razumevanje strukture podataka. Gradi se “stablo” klastera (dendogram). Svaka iteracija počinje kao zaseban klaster, a zatim se iterativno spajaju najbliži klasteri dok ne ostane jedan najveći. Ovo možemo prikazati i vizuelno.

hc_result = hclust(dist(rakete_izbaceno), method = "ward.D2")
plot(hc_result, main = "Dendrogram", sub = "", xlab = "Opservacije", ylab = "Udaljenost")


Sada treba izabrati broj klastera K. Vizuelno uočavamo da postoje tri jasne grupe na dendogramu (tri jasne “grane”) pa pogledajmo kako izgleda dendogram ako bi broj klastera bio isti kao za K - means tj. K = 3. Na grafiku crvene granice pedstavljaju tri krajnja klastera koja smo zapravo i tražili. Nakon vizuelizacije, možemo se poslužiti i “elbow” metodom radi provere.

plot(hc_result, main = "Dendrogram", sub = "", xlab = "Opservacije", ylab = "Udaljenost")
rect.hclust(hc_result, k = 3, border = "darkred")

fviz_nbclust(rakete_izbaceno, FUN = hcut, method = "wss")


“Elbow” metoda nam potvrđuje da je naš vizuelni zaključak bio ispravan jer je lakat ponovo na K = 3.

Kao i kod K - means metode, vršimo ponovo podelu na 3 klastera i možemo pogledati koja opservacija pripada kom.

clusters_hc = cutree(hc_result, k = 3)

#dodeljeni klasteri svakoj opservaciji
clusters_hc
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 17 18 20 22 23 24 27 28 29 34 35 
##  1  1  2  2  2  2  2  2  2  2  2  1  2  2  3  3  3  2  3  3  3  2  3  3  2  1 
## 36 38 39 
##  2  2  3
# prikaz koliko opservacija je svrstano u svaki od 3 klastera koje imamo
table(clusters_hc)
## clusters_hc
##  1  2  3 
##  4 16  9


Poređenje klasterizacija

Sada kada znamo kom klasteru je dodeljena koja opservacija putem K - means i Hierarchical algoritma, možemo uporediti rezultate.

Pošto algoritmi klasterizacije nemaju “pravu” ciljnu varijablu ne možemo koristiti standardne metrike. Umesto toga koristićemo funkciju “silhouette” koja računa Silhouette score koji nam pokazuje koliko su objekti unutar klastera slični, a različiti od objekata u drugim klasterima. Ova mera uzima vrednosti od -1 do 1. Rezultat blizu 1 označava da je tačka dobro uklopljena u svoj klaster, rezultat blizu 0 da je na granici i rezultat blizu -1 da je tačka loše pozicionirana. Nakon što izračunamo srednju vrednost svih rezultata dobijamo procenu.

library(cluster)
sil_kmeans <- silhouette(clusters_kmeans, dist(rakete_izbaceno))
sil_hclust <- silhouette(clusters_hc, dist(rakete_izbaceno))

mean(sil_kmeans[, 3])
## [1] 0.4547102
mean(sil_hclust[, 3])
## [1] 0.4446867

Vidimo da su obe klasterizacije bile uspešne (rezultati oko 0.5 su solidno dobri), kao i da su gotovo jednako dobre. Prema ličnoj preferenci i boljem intuitivnom shvatanju izabraću Hijerarhijsko klasterovanje za dalji rad.


Da li možemo da zaključimo nešto o podeli raketa na kratak i dug domet?

Jasno, nama ništa ne garantuje da smo uspeli da dođemo do bilo kakve jasne podele na rakete prema njihovom dometu, koji je sve vreme nepoznat. Klasterizacija se koristila sličnostima i razlikama između raketa prema njihovim ostalim, poznatim obeležjima. Ono što mi intuitivno znamo jeste da ostala obeležja utiču i na domet raketa. Npr. ukoliko uvidimo jasno da su u nekim od 3 klastera koja imamo isključivo rakete slične po masi i dužini, recimo velike mase i dužine, to će nam omogućiti da neformalno zaključimo da upravo ovi klasteri predstavljaju ujedno i rakete dugog dometa (jer nam je jasno da baš takve rakete imaju dug domet). Na taj način nam homogenost klastera prema ostalim karakteristikama kao i njihova međusobna različitost omogućuju da zaključimo da svaki klaster predstavlja jedan “tip” raketa i da se najverovatnije ti tipovi odnose i na njihov domet (kratak ili dug, u zavisnost od vrednosti ostalih obeležja povezanih sa dometom).


Kako poznata obeležja utiču na domet (intuitivno shvatanje):

  • masa - veće (teže) rakete imaju veći domet jer mogu nositi jače motore i više goriva (ZNAČAJAN UTICAJ)
  • dužina - kao i za masu, duže rakete imaju veći domet (ZNAČAJAN UTICAJ)
  • kalibar - povezan je sa veličinom motora, pa veći kalibar znači i veći domet (ZNAČAJAN UTICAJ)
  • tip goriva - rakete sa tečnim gorivom imaju veći domet jer su efikasnije u korišćenju energije (SREDNJI UTICAJ)
  • način lansiranja - nije presudan (MALI UTICAJ)
  • payload (korisni teret) - veća težina tereta smanjuje kapacitet za gorivo pa smanjuje domet (SREDNJI UTICAJ)
  • broj podlgava - rakete sa više glava imaju veću stabilnost pa i domet (SREDNJI UTICAJ)
  • broj faza motora - višefazni motori povećavaju brzinu i efikasnost raketa pa i domet (SREDNJI UTICAJ)

Sada ćemo grafički prikazati rakete po 2 atributa, kako bi uvideli pravilnost u podeli na klastere i zaključili nešto više o dometu.

NAPOMENA: Prikazaćemo i legendu koja nam pokazuje koji od gore dodeljenih brojeva klastera (1,2,3) označava koju boju. Radi jednostavnosti koda uradićemo to samo prvi put, znajući da isto važi i za svaki naredni grafik.

clusters_hc_colors = as.factor(clusters_hc)

plot(rakete_izbaceno$masa_kg, rakete_izbaceno$duzina_m, 
     pch = 16, cex = 1.5, col = clusters_hc_colors)

colors_assigned = levels(clusters_hc_colors)

legend(x = 2.3, y = 0.3, 
       legend = paste("Klaster", 1:length(colors_assigned)), 
       fill = colors_assigned, 
       title = "Klasteri", 
       xpd = TRUE)

Svaka boja predstavlja jedan klaster (u legendi vidimo i koja koji), a svaka tačka jednu opservaciju. Vidimo da su opservacije prilično pravilno raspodeljene po klasterima u smislu njihove mase i dužine - rakete najmanje mase i dužine čine jedan klaster i tako redom.


plot(rakete_izbaceno$masa_kg, rakete_izbaceno$tip_goriva, 
     pch = 16, cex = 1.5, col = clusters_hc)

Setimo se da na prethodnom grafiku, tačke crne boje bile su rakete velike masei dužine. Sada vidimo da isti klaster (zovimo ga “Crni klaster”) karakteriše i čvrsto gorivo (omogućava veći domet). Možemo zaključiti da Crni klaster predstavlja rakete dugog dometa. Što se tiče druga dva klastera (“Zeleni klaster” i “Roze klaster”) ne oučavamo jasnu pravilnost prema tipu goriva (algoritam ih je “sekao” prema masi).

Zato pogledajmo još jedan grafik klasterizacije prema 2 atributa.


plot(rakete_izbaceno$duzina_m, rakete_izbaceno$kalibar_m, 
     pch = 16, cex = 1.5, col = clusters_hc)

Još jednom potvrđujemo da Crni klaster predstavlja rakete dugog dometa (rakete velike mase, dužine, kalibra i sa čvrstim gorivom). Vidimo da Zeleni klaster karakterišu rakete male mase, dužine i kalibra (setimo se da za tip goriva nismo uvideli jasnu pravilnost), te možemo zaključiti da on predstavlja rakete kratkog dometa.


Šta da zaključimo o Roze klasteru?
U ovom klasteru nalaze se rakete osrednje mase, dužine i kalibra (o tipu goriva nismo zaključili ništa) od raketa koje su ostale u bazi nakon izbacivanja raketa srednjeg dometa.

Za početak pogledajmo koje su to rakete u Roze klasteru (prisetimo se legende sa prvog grafika: Roze klaster je zapravo klaster br. 2).

#koristimo početnu bazu jer smo samo u njoj sačuvali nazive raketa
rakete_izbaceno_pocetna_baza = rakete_pocetna_baza[!rakete_pocetna_baza$klasa_po_dometu == "M", ]

rakete_izbaceno_pocetna_baza[clusters_hc == 2, "NATO_ime"]
## # A tibble: 16 × 1
##    NATO_ime              
##    <chr>                 
##  1 SS-25 Sickle          
##  2 SS-27 Mod 1           
##  3 SS-27 Mod 2           
##  4 SS-N-18 Stingray Mod 1
##  5 SS-N-18 Stingray Mod 2
##  6 SS-N-18 Stingray Mod 3
##  7 SS-N-23 Sineva Mod 2  
##  8 SS-N-23 Liner         
##  9 SS-N-32               
## 10 SS-N-27A Sizzler      
## 11 SS-N-27B Sizzler      
## 12 SS-X-31 Frontier      
## 13 SS-N-17 Snipe         
## 14 CSS-10                
## 15 CSS-N-14              
## 16 CSS-20

Vidimo da su nam u Roze klasteru pretežno sovjetske/ruske rakete (čak 11 njih, a svega 3 kineske). Stoga možemo zaključiti da nam i Roze klaster, kao i Crni, predstavlja rakete dugog dometa, a Zeleni rakete kratkog dometa.

NAPOMENA: prethodni zaključak biće detaljno objašnjen pred kraj 6. zadatka, kada se detaljnije uputimo u istoriju razvoja raketa i opravdamo zaključak da je klaster pretežno sovjetskih raketa najverovatnije klaster raketa dugog dometa

Takođe, u cilju donošenja (neformalnog) zaključka, možemo pogledati i masu raketa u Roze klasteru.

Dizajn raketa i njihov domet zavise od mnogo faktora, ali masa je definitivno veoma bitno obeležje koje utiče na domet raketa. Prema nezvaničnim podacima (dolazi često do preklapanja brojki, ali može se uhvatiti neka pravilnost) rakete dugog dometa obično teže između 30000 i 120000 kilograma, dok su one kratkog dometa između 1000 i 8000 kilograma.

rakete_izbaceno_pocetna_baza[clusters_hc == 2, "masa_kg"]
## # A tibble: 16 × 1
##    masa_kg
##      <dbl>
##  1   45100
##  2   47200
##  3   49600
##  4   35300
##  5   35300
##  6   38000
##  7   40300
##  8   40300
##  9   36800
## 10    1920
## 11    1570
## 12   36000
## 13   26900
## 14   42000
## 15   42000
## 16   80000

Vidimo da gotovo sve rakete iz Roze klastera teže preko 30000 kg. Ovaj podatak nam ide u prilog prethodnom zaključku o tome da se u ovom klasteru nalaze rakete dugog dometa. Naravno, dolazi do nekih odstupanja, ali to smo mogli i očekivati s obzirom na to smo znali da ovakva klasterizacija ne može savršeno podeliti rakete prema njihovom dometu kao nepoznatom obeležju.


Šesti zadatak

U poslednjem zadatku cilj nam je da razumemo naše podatke i pojmove koji se korite u njima na još neki način pored matematičkog. Primarni pristup biće istorijski i politički. Posebno ćemo, kako zadatak i nalaže, istražiti koliko se često javljaju rakete srednjeg dometa, koje su one starosti, kome pripadaju i zašto.

Prisetimo se kako je izgledala naša baza na samom početku (na primeru prvih 5 opservacija).

kable(head(rakete_pocetna_baza, n = 5), 
      format = "html", 
      table.attr = "style='font-size: 12px; width: 70%'"
)
rusko_ime NATO_ime masa_kg duzina_m kalibar_m broj_podglava payload_kg tip_goriva max_domet_km nacin_lansiranja broj_faza_motora klasa_po_dometu
R-36M2 Voevoda SS-18 Satan 209600 32.2 3.05 10 8800 Tecno 16000 Silos 2 L
UR-100N SS-19 Stiletto 105600 27.0 2.50 6 3355 Tecno 10000 Silos 2 L
RT-2PM Topol SS-25 Sickle 45100 21.5 1.80 1 1000 Cvrsto 11000 RM TEL 3 L
RT-2PM2 Topol-M SS-27 Mod 1 47200 22.7 1.90 1 1200 Cvrsto 11000 RM TEL 3 L
RS-24 Yars SS-27 Mod 2 49600 23.0 2.00 3 1200 Cvrsto 12000 RM TEL 3 L


Zemlje porekla raketa srednjeg dometa, godine proizvodnje

Dakle naš rad se bazirao ne raketama kojima znamo rusko ime i NATO ime. Međutim, nemamo eksplicitno naveden podatak u vlasništvu koje sile je svaka od njih. Ono što možemo da uradimo je da pogledamo NATO imena naših raketa.

kable(rakete_pocetna_baza[,2], 
      format = "html",
      table.attr = "style='font-size: 15px; width: 70%'")
NATO_ime
SS-18 Satan
SS-19 Stiletto
SS-25 Sickle
SS-27 Mod 1
SS-27 Mod 2
SS-N-18 Stingray Mod 1
SS-N-18 Stingray Mod 2
SS-N-18 Stingray Mod 3
SS-N-23 Sineva Mod 2
SS-N-23 Liner
SS-N-32
SS-X-30 Satan 2
SS-N-27A Sizzler
SS-N-27B Sizzler
SS-26 Stone
Dagger
SS-21 Scarab A
SS-21 Scarab B
SS-N-6 Serb
SS-X-31 Frontier
SS-20 Saber
SS-23 Spider
SS-1B Scud A
SS-1C Scud B
SS-1D Scud C
SS-1E Scud D
SS-N-17 Snipe
CSS-7
CSS-X-15
CSS-6
CSS-11
CSS-5
Guam Killer (unofficial)
CSS-10
CSS-4
CSS-N-14
CSS-5
CSS-20
CSS-?

Zašto gledamo baš NATO ime raketa?
Možemo uočiti da su NATO imena raketa kodna, tačnije sva imena počinju ili prefiksom SS- ili CSS- (izuzev dve rakete koje ćemo istražiti odvojeno a to su, prema NATO imenu, “Dagger” i “Guam Killer”).
NATO je sovjetske/ruske rakete imenovao prefiksom SS- (Surface-to-Surface), a kineske sa CSS- (Chinese Surface-to-Surface). Postojali su i drugi prefeksi NATO kodnog imenovanja, ali nam trenutno nisu od ključnog značaja. Izuzeci u NATO imenovanju su određene rakete čiji nazivi označavaju neke njihove posebne karakteristke ili su jednostavno nastala zbog medijskog interesa. Tako, imamo rusku raketu “Dagger” (jedna od najmodernijih oružja ruske vojske) i kinesku “Guam Killer” (poznata po dometu taman dovoljnom za gađanje američke baze na ostrvu Guam).


Prvo, dodajmo u našu bazu podataka još jedno obeležje koje će nam pomoći u daljem radu, a to je sila kojoj svaka raketa pripada (to ćemo uraditi prema gore objašnjenim NATO kodnim imenima). Ostavljamo prostora da postoje neke rakete pored “Dagger” i “Guam Killer” koje ne počinju sa SS- ili CSS-, a nismo ih uvideli (dodavanjem difoltne vrednosti “Nepoznato”).
Zatim, želimo da vidimo koliko ima ruskih, a koliko kineskih, raketa srednjeg dometa (one koje za klasa_po_dometu uzimaju vrednost “M”, što znamo iz prethodnog zadatka) i koje su to rakete, kako bi istražili kada su pravljene.

library(dplyr)

rakete_pocetna_baza = rakete_pocetna_baza %>%
  mutate(
    zemlja = case_when(
      grepl("^SS-", NATO_ime) ~ "Rusija",
      grepl("^CSS-", NATO_ime) ~ "Kina",
      NATO_ime == "Dagger" ~ "Rusija",
      NATO_ime == "Guam Killer (unofficial)" ~ "Kina",
      TRUE ~ "Nepoznato"
    )
  )

rakete_rusija_m = rakete_pocetna_baza %>%
  filter(zemlja == "Rusija", klasa_po_dometu == "M")

rakete_kina_m = rakete_pocetna_baza %>%
  filter(zemlja == "Kina", klasa_po_dometu == "M")

kable(rakete_rusija_m$NATO_ime, 
      format = "html", 
      table.attr = "style='font-size: 12px; width: 70%'"
)
x
Dagger
SS-N-6 Serb
SS-20 Saber
SS-1D Scud C
SS-1E Scud D 


kable(rakete_kina_m$NATO_ime, 
      format = "html", 
      table.attr = "style='font-size: 12px; width: 70%'"
)
x
CSS-6
CSS-11
CSS-5
Guam Killer (unofficial)
CSS-5

Pogledajmo sada kada su proizvedene sve ove rakete.

NAPOMENA: Imamo dve kineske rakete CSS-5 ali one nisu iste već je samo njihovo NATO ime isto (to su rakete Dong Feng 21 i Dong Feng 4, redom, prema ruskom imenu)

1. Dagger (R) - 2018
2. SS-N-6 Serb (R) - 1960-1961
3. SS-20 Saber (R) - 1976-1979
4. SS-1D Scud C (R) - 1970-1972
5. SS-1E Scud D (R) - 1976-1979


6. CSS-6 (K) - 1990
7. CSS-11 (K) - 1970 - 1971
8. CSS-5 (K) - 1985
9. Guam Killer (K) - 2010 - 2015 (nezvanično)
10. CSS-5 (K) - 1975

Vidimo da ne postoji nijedna ruska raketa srednjeg dometa proizvedena u periodu 1988-2018. Ovo nije slučaj kod kineskih raketa. Razmotrimo zašto.


Reagan - Gorbačev (INF sporazum)

1970-ih godina Sovjetski Savez je stavio u upotrebu raketu SS-20 Saber koju smo već pominjali. Ovo je bila raketa srednjeg dometa koja je lako mogla da pogodi Zapadnu Evropu sa sovjetskih teritorija, što je predstavljalo veliki izazov za NATO jer je povećalo nuklearnu pretnju.

Zato je NATO, pod uticajem Nemačke, doneo 1979. “Double Track Decision” što uključuje:

  • uklanjaje 1000 nuklearnih bojevih glava (ekplozivni uređaj na vrhu rakete zadužen za uništenje na cilju) iz Evrope
  • ako pregovori sa Sovjetskim Savezom ne uspeju, modernizacija i povećanje nuklearnih snaga NATO-a

Upravo je ova strategija i zabrinutost zbog velikog sovjetskog i kineskog nuklearnog rasta bila preduslov za INF sporazum potpisan 1987. između SAD-a (Ronald Reagan) i Sovjetskog Saveza (Mikhail Gorbachev).

Potpisivanje INF - a

Sporazum je strogo zabranio proivodnju i raspoređivanje raketa srednjeg dometa, čime su SAD i Sovjetski Savez ukinuli ovu pretnju i smanjili mogućnost nuklearnog sukoba u Evropi. Bio je značajan korak ka smanjenju napetosti u Hladnom ratu i stabilizaciji evropske sigurnosti, kao i smanjenju nuklearne opasnosti u svetu. Time je Sovjetski Savez bio u obavezi da sve postojeće rakete srednjeg dometa uništi i da ne proizvodi nove.

Za razliku od njih, Kina nije bila potpisnica i time nije bila deo ovog sporazuma. Nije bila ni u kakvoj obavezi. Kina je za to vreme usavršavala i povećala proivodnju raketa srednjeg dometa i tako stekla stratešku prednost, pogotovo u kontekstu njenog rastućeg uticaja u Aziji.


Bolton - Šoigu (sastanak 2018.)

Uprkos svim dogovorima, već od 2014. godine SAD optužuju Rusiju da krši INF sporazum. Naime, sporazum je zabranjivao isključivo kopnene balističke i krstareće rakete srednjeg dometa i nije se odnosio na vazdušne i pomorske. Prema našim podacima to bi značilo da su zabranjene sve rakete srednjeg dometa koje za tip_lansiranja uzimaju “Silos” ili “RM TEL”, ali ne i one koje uzimaju “Mornarički” i “Avion”. Na taj način Rusija je uspevala da doskoči INF sporazumu tako što je ipak proizvodila rakete srednjeg dometa koje su predstavljale pretnju, ali bez formalnog kršenja sporazuma.

Npr. raketa “Dagger” o kojoj je već bilo reči je prvi put javno predstavljena u martu 2018. ali je već krajem 2017. godine ušla u operativnu upotrebu. U pitanju je raketa srednjeg dometa koja bi bila zabranjena sporazumom da je lansirana sa kopna, ali je Rusija formalno kršenje istog preduhitrila tako što je “Dagger” raketa lanisirana sa aviona te nije pod INF zabranom.

Pored ovakvih doskočica Rusije postojao i konkretan problem, a to je bila SSC-8 (surface-to-Surface Cruise) raketa koju nemamo u bazi. U pitanju je kopnena krstareća raketa koja bi kao takva bila zabranjena INF-om ako bi bila srednjeg dometa. Problem je nastao jer su SAD tvrdile da je domet rakete SSC-8 2500km, a Rusija da ne prelazi 480km (ali nisu pružili dovoljno dokaza za ovu tvrdnju). Dodatno, Rusija je imala i altrenativu. To je bila raketa “Kalibr” koju takođe nemamo u bazi, a koja je preteča rakete SSC-8 samo što joj je lanser sa mornarice. Tako je izgledalo kao da Rsuija pronalazi “rupu” u INF sporazumu kako bi uspela da proizvede kopnenu krstareću raketu srednjeg dometa (SSC-8).

Zbog svega navedenog, u oktobru 2018. tadašnji američki savetnik za nacionalnu bezbednost J.Bolton sastao se sa ruskim ministrom odbrane S.Soigu i drugim zvaničnicima u Moskvi. Tema sastanka bila je INF sporazum, a Bolton je jasno preneo američku nameru da se povuče iz tog sporazuma, upravo zbog navodnog ruskog kršenja. Neposredno nakon toga, 20. oktobra 2018. je Donald Trump najavio zvanično povlačenje SAD iz INF sporazuma. Dodao je kao razlog i to što Kina (koja, podsetimo se, nije bila ograničena sporazumom) postaje ozbiljna pretnja američkim interisima na Pacifiku.

Konačno, 1. februara 2019. SAD su zvanično suspendovale svoje obaveze prema INF-u, a Rusija je učinila isto već dan kasnije. Do avgusta 2019. sporazum je u potpunosti prestao da važi. I SAD i Rusija ponovo su započe proizvodnju raketa srednjeg dometa svih tipova lansera.


Objašnjenje zaključka u prethodnom zadatku

Sada kada znamo više o istoriji razvoja raketa možemo bolje opravdati zaključak da je Roze klaster u 5. zadatku predstavljao rakete dugog, a ne kratkog dometa. Naime, još tada smo uvideli da se on sastoji pretežno od sovjetskih/ruskih raketa. Prisetimo se da smo u tom zadatku radili nad bazom u kojoj su bile izbačene sve rakete srednjeg dometa. Sada kada znamo da je Sovjetski Savez dugo imao zabranu proizvodnje raketa srednjeg dometa (tokom važenja INF sporazuma), možemo se osloniti na činjenicu da su u tom periodu proizvodili mnogo više raketa dugog dometa koje su imale bolje i značajnije performanse od onih kratkog dometa.

Za razliku od Kine koja je razvijala i rakete kratkog dometa (jer je bila fokusirana na regionalne sukobe u Pacifiku i jugoistočnoj Aziji), Rusija je više težila nečemu drugačijem. Cilj im je bio održavanje vojne moći sa SAD-om (to su postizali raketama dugog dometa jer su tako mogli da pokriju velike geografske udaljenosti) i nuklearno odvraćanje od napada (rakete dugog dometa bile su osnovni element ruske nuklearne strategije).

Takođe, Rusija je generalno kroz istoriju proizvodila više raketa dugog dometa (i pre INF sporazuma) u skladu sa njihovom strategijom i potrebama za nuklearnim odvraćanjem i globalnom dominacijom.

Stoga, veća brojnost ruskih raketa u Roze klasteru daje nam podlogu da (neformalno) zaključimo da se radi o raketama dugog dometa.