Add files via upload

notebook.Rmd

@@ -0,0 +1,131 @@
+---
+title: "Global Terrorism Database - analiza"
+output: pdf_document
+---
+
+```{r echo=FALSE, message=FALSE}
+## Czynności przygotowawcze - załadowanie potrzebnych pakietów oraz utworzenie i wypełnienie bazy danych
+library(RSQLite)
+library(tidyverse)
+library(ggpubr)
+df <- read.csv("globalterrorismdb_0919dist.csv")
+conn <- dbConnect(RSQLite::SQLite(), "gtd.db")
+#dbWriteTable(conn, "Attacks", df)
+```
+
+# 1. Wstęp
+
+Raport zawiera analizę danych pochodzących z Global Terrorism Database (GTD) - bazy utrzymywanej przez amerykańskie National Consortium for the Study of Terrorism and Responses to Terrorism (START) na uniwersytecie Marylandu w College Park, zawierającej listę ponad 190 tys. ataków terrorystycznych od 1970 r. do końca 2018 r (pomijając rok 1993). Każdy opisany jest przez najwyżej 135 zmiennych, określających m.in datę ataku, jego lokalizację, cel, rodzaj użytej broni lub narzędzi czy ilość poszkodowanych. Pełen opis bazy znajduje się w dołączonym dokumencie 'Codebook: Inclusion Criteria and Variables'. 
+
+# 2. Analiza wstępna/eksploracyjna
+
+## 2.1 Ilość ataków rocznie wg regionu
+
+Najpierw sprawdźmy, jak ilość ataków jest uzależniona od regionu świata, pozwoli nam to stwierdzić, które miejsca są obecnie najbardziej zagrożone atakiem terrorystycznym.
+
+```{r echo=FALSE}
+tmp <- df %>% group_by(region_txt, iyear) %>% summarize(attack_count = n())
+g <- ggplot(tmp, aes(iyear, attack_count)) + geom_line() + facet_wrap(region_txt ~., nrow = 4, ncol = 3, scales = "free_y") + labs(x = "Rok") + labs(y = "Ilość ataków")
+print(g)
+```
+Na wykresie zauważyć można nagły wzrost ilości ataków na Bliskim Wschodzie i Afryce północnej, Afryce subsaharyjskiej oraz Azji południowej ok. roku 2014 - jest to najprawdopodobniej efekt wzrostu aktywności ugrupowań islamistycznych, zwłaszcza tzw. Państwa Islamskiego; widać także chwilowy wzrost ilości ataków w Europie wschodniej, co prawdopodobnie związane jest z pojawieniem się prorosyjskich ugrupowań separatystycznych na terenach wschodniej Ukrainy. W pozostałych regionach albo obserwujemy spadek, albo też ilość ataków terrorystycznych nie jest obecnie duża.
+
+## 2.2 Najaktywniejsze ugrupowania terrorystyczne w Europie Wschodniej oraz na Bliskim Wschodzie, Afryce północnej i subsaharyjskiej, i Azji południowej i połudnowo-wschodniej od roku 2010
+
+```{r echo = FALSE, fig.width=7.5,fig.height=8}
+tmp <- dbGetQuery(conn, "SELECT iyear, region_txt, gname FROM Attacks WHERE iyear >= 2010 AND (region = 9 OR
+                  region = 10 OR region = 11 OR region = 6 OR region = 5) AND NOT gname = 'Unknown'")
+tmp2 <- tmp %>% group_by(region_txt, gname) %>% summarize(attack_count = n())
+EE <- tmp2 %>% filter(region_txt == "Eastern Europe") %>% arrange(desc(attack_count)) %>% head(3) 
+ME <- tmp2 %>% filter(region_txt == "Middle East & North Africa") %>% arrange(desc(attack_count)) %>% head(3)
+SSA <- tmp2 %>% filter(region_txt == "Sub-Saharan Africa") %>% arrange(desc(attack_count)) %>% head(3)
+SEA <- tmp2 %>% filter(region_txt == "Southeast Asia") %>% arrange(desc(attack_count)) %>% head(3)
+SA <- tmp2 %>% filter(region_txt == "South Asia") %>% arrange(desc(attack_count)) %>% head(3)
+tmp3 <- rbind(EE,ME,SSA,SEA,SA) %>% mutate(Index = 4 - with_order(order_by = attack_count, fun = row_number, x = attack_count ))
+
+g <- ggplot() + geom_col(data = tmp3, aes(x = reorder(gname, Index), y = attack_count)) + facet_wrap(region_txt ~., nrow = 5, ncol = 1, scales = "free") +  labs(x = "Organizacja") + labs(y = "Ilość ataków")
+print(g)
+```
+
+Wykres potwierdza przypusczenia co do źródła ataków, pierwsze miejsca zajmują organizacje islamistyczne, w Europie wschodniej separatyści, natomiast w Azji południowo-wschodniej najwięcej ataków zorganizowała komunistyczna New People's Army z Filipin, jednak drugie i trzecie miejsce przypada organizacjom islamistycznym.
+
+## 2.3 Sumaryczna ilość ofiar (zabitych lub rannych) w zależności od zastosowanej broni
+
+```{r echo=FALSE, fig.height=4}
+tmp <- dbGetQuery(conn,"SELECT weaptype1_txt, nwound, nkill FROM Attacks WHERE success = 1 AND NOT (nwound IS NULL OR nkill IS NULL OR weaptype1_txt = 'Unknown')")
+
+tmp2 <- tmp %>% group_by(weaptype1_txt) %>% summarize(victim_count = sum(nwound, nkill))
+tmp2[11,1] <- "Vehicle"
+
+g <- ggplot(tmp2, aes(x=reorder(weaptype1_txt,-victim_count), y=victim_count)) + geom_col(aes(fill = weaptype1_txt)) + theme(axis.text.x=element_text(angle=90, hjust=1), legend.position="none") + labs(x = "Rodzaj broni") + labs(y = "Ilość ofiar") + geom_text(aes(label = victim_count), vjust = -0.1) 
+print(g)
+```
+Wzięto pod uwagę jedynie udane ataki (parametr *success* równy 1), których liczba ofiar i rodzaj użytego uzbrojenia są znane. Największą liczbę ofiar spowodowały ładnuki wybuchowe - ponad 2 razy więcej niż reszta rodzajów broni razem wzięta, co może świadczyć zarówno o skuteczności tego narzędzia jak i o popularności. Żeby to sprawdzić, należy policzyć średnią ilość ofiar na atak, a także odchylenie standardowe, co pozwoli stwierdzić jak bardzo różne efekty daje zastosowanie danego rodzaju broni.
+
+## 2.4 Średnia i odchylenie standardowe ilości ofiar w zależności od zastosowanej broni
+
+```{r echo=FALSE, fig.width=9,fig.height=4}
+tmp2 <- tmp %>% group_by(weaptype1_txt) %>% summarize(victim_count = sum(nwound, nkill), avg = mean(nwound+nkill), sd=sd(nwound+nkill)) %>% filter(victim_count > 0)
+
+tmp2[9,1] <- "Vehicle"
+
+g1 <- ggplot(tmp2, aes(x=reorder(weaptype1_txt,-avg), y=avg)) + geom_col(aes(fill = weaptype1_txt)) + theme(axis.text.x=element_text(angle=90, hjust=1), legend.position="none") + labs(x = "Rodzaj broni") + labs(y = "Średnia ilość ofiar") + geom_text(aes(label = sprintf("%0.2f", round(avg, digits = 2))), vjust = -0.1)
+
+g2 <- ggplot(tmp2, aes(x=reorder(weaptype1_txt,-avg), y=sd)) + geom_col(aes(fill = weaptype1_txt)) + theme(axis.text.x=element_text(angle=90, hjust=1), legend.position="none") + labs(x = "Rodzaj broni") + labs(y = "Odchylenie standardowe") + geom_text(aes(label = sprintf("%0.2f", round(sd, digits = 2))), vjust = -0.1, hjust= 0.4)
+
+print(ggarrange(g1,g2,ncol=2,nrow = 1))
+x <- dbGetQuery(conn,"SELECT * FROM Attacks WHERE nwound = 10878")
+```
+Wartość średniej i odchylenia standardowego dla kategorii *Vehicle* sugeruje występowanie wartości odstających w danych, zawyżających wynik na tyle, że uzyskane wartości nie pozwalają na jednoznaczne określenie relacji pomiędzy poziomami skuteczności różnych rodzajów uzbrojenia. Po przejrzeniu danych udało się znaleźć przyczynę takiego stanu rzeczy - jest to uwzględnienie informacji nt. ataku na WTC z 11.09.2001 r. Po odrzuceniu odpowiadających mu wpisów otrzymujemy:
+
+```{r echo=FALSE, fig.width=9,fig.height=4}
+tmp <- dbGetQuery(conn,"SELECT weaptype1_txt, nwound, nkill FROM Attacks WHERE success = 1 AND NOT (nwound IS NULL OR nkill IS NULL OR weaptype1_txt = 'Unknown' OR nwound = 10878)") #Ilość rannych w ataku na WTC określono na 10878 na samolot
+tmp2 <- tmp %>% group_by(weaptype1_txt) %>% summarize(victim_count = sum(nwound, nkill), avg = mean(nwound+nkill), sd=sd(nwound+nkill)) %>% filter(victim_count > 0)
+
+tmp2[9,1] <- "Vehicle"
+
+g1 <- ggplot(tmp2, aes(x=reorder(weaptype1_txt,-avg), y=avg)) + geom_col(aes(fill = weaptype1_txt)) + theme(axis.text.x=element_text(angle=90, hjust=1), legend.position="none") + labs(x = "Rodzaj broni") + labs(y = "Średnia ilość ofiar") + geom_text(aes(label = sprintf("%0.2f", round(avg, digits = 2))), vjust = -0.1)
+
+g2 <- ggplot(tmp2, aes(x=reorder(weaptype1_txt,-avg), y=sd)) + geom_col(aes(fill = weaptype1_txt)) + theme(axis.text.x=element_text(angle=90, hjust=1), legend.position="none") + labs(x = "Rodzaj broni") + labs(y = "Odchylenie standardowe") + geom_text(aes(label = sprintf("%0.2f", round(sd, digits = 2))), vjust = -0.1, hjust= 0.4)
+
+print(ggarrange(g1,g2,ncol=2,nrow = 1))
+x <- dbGetQuery(conn,"SELECT * FROM Attacks WHERE nwound = 10878")
+```
+Teraz widać, że najwyższą średnią ofiar przynoszą ataki z zastosowaniem broni chemicznej lub biologicznej, ataki przy użyciu materiałów wybuchowych są dopiero na czwartym miejscu, co oznacza, że wysoka łączna ilość ofiar ataków z zastosowaniem materiałów wybuchowych wynika jedynie z popularności tej metody ataku.
+
+# 3. Modelowanie statystyczne
+
+## 3.1 Zależność między ilością zabitych w ataku a ilością rannych
+
+```{r echo=FALSE, fig.height=3}
+tmp <- tmp %>% filter(nkill < 600) %>% filter(nwound < 2000) #odcięcie wartości odstających
+g <- ggplot(tmp, aes(nwound, nkill)) + geom_point() + geom_smooth(method = "lm") + labs(x ="Liczba rannych", y = "Liczba zabitych")
+linearModel <- lm(nkill ~ nwound, data=tmp)
+print(g)
+```
+Zauważyć można wzrost liczby zabitych przy rosnącej liczbie rannych - wytłumaczyć to można wzrostem skali ataku. Parametry prostej regresji:
+
+```{r echo=FALSE}
+summary(linearModel)
+```
+
+## 3.2 Analiza rozkładu ilości zabitych w ataku terrorystycznym
+```{r echo=FALSE, fig.height=4}
+testData <- dbGetQuery(conn,"SELECT region_txt, iyear, imonth, iday, nwound, nkill FROM Attacks WHERE success = 1 AND NOT (nwound IS NULL OR nkill IS NULL OR weaptype1_txt = 'Unknown') and nkill > 0")
+testData2 <- testData %>% group_by( nkill) %>% summarize(attack_count = n())
+g <- ggplot(testData2, aes(nkill , attack_count)) + geom_line() + labs(y = "Ilość ataków", x = "Ilość zabitych" )
+g2 <- ggplot(testData2, aes(log(nkill) , log(attack_count))) + geom_point() + geom_smooth(method = "lm", se=FALSE) + labs(y = "Logarytm ilości ataków", x = "Logarytm ilości zabitych" )
+print(ggarrange(g,g2, ncol = 2, nrow = 1))
+```
+Wykresy przedstawiają ilość ataków w zależności od liczby zabitych, przy czym wzięto pod uwagę jedynie te o niezerowej liczbie ofiar. Na drugim wykresie widać, że zaleźność między logarytmem ilości ataków a logarytmem ilości zabitych jest bliska liniowej, co sugeruje, że rozkład ilości zabitych w ataku terrorystycznym jest zgodny z prawem potęgowym (power law), czyli, że funkcja prawdopodobieństwa rozkładu ma postać $p(x) = P(X=x) = Cx^{-\alpha}$. W celu sprawdzenia dopasowania rozkładu wykorzystano pakiet *poweRlaw* i wchodzącej w jego skład funkcji *bootstrap_p*. Uzyskana w ramach testu wartość *p*:
+```{r echo=FALSE}
+library(poweRlaw)
+m <- displ$new(testData$nkill)
+m$pars <- estimate_pars(m)
+bs_p <- bootstrap_p(m, no_of_sims=100, threads=2)
+print(bs_p$p)
+```
+pozwala na przyjęcie hipotezy zerowej. Estymowana wartość parametru $\alpha$ wynosi:
+```{r echo=FALSE}
+print(m$pars)
+```