report.Rmd

---
title: "**Analiza kriminala i socio-ekonomskih faktora**"
author:
- Jan Grgić
- Bernard Spiegl
- Korina Šimičević
- Dunja Šmigovec
output:
  pdf_document: default
header-includes:
  \usepackage{float}
  \floatplacement{figure}{H}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Uvod

Cilj ovog projekta je analizirati i vizualno prikazati podatke o kriminalu na području Chicaga bilježene u periodu od godinu dana.

## Učitavanje paketa

```{r echo=T, error=FALSE, warning=FALSE, message=FALSE}
library(dplyr)
library(corrplot)
library(car)
library(MASS)
library(lmtest)
library(het.test)
```

## Učitavanje podataka

```{r}
chicago_year = read.csv("data/Crimes_-_One_year_prior_to_present.csv")
chicago_pc = read.csv("data/Chicago_poverty_and_crime.csv")
dim(chicago_pc)
dim(chicago_year)
```

## Kategorije kriminala u datasetu
```{r}
categories <- unique(chicago_year$PRIMARY.DESCRIPTION)
categories
length(categories)
```

## Prikaz frekvencija različitih kategorija zločina unazad godinu dana

U nastavku (Figure \ref{fig:category_frequencies}) možemo vidjeti da je su najčešći lakši oblici zločina što je bilo i za očekivati.
```{r category_frequencies, fig.cap = "Frekventnost kategorija unazad godinu dana", out.width='0.9\\linewidth'}
table <- table(chicago_year$PRIMARY.DESCRIPTION)
data <- as.data.frame.table(table)


df <- data[order(data$Freq, decreasing = TRUE),]
op <- par(mar = c(12, 4, 1, 2))
barplot(
  df$Freq,
  names.arg = df$Var1,
  las = 2,
  cex.names = 0.6,
  col = "skyblue",
  ylim = range(pretty(c(0, df$Freq))),
  main = "Frequency of different crime activities"
)

```

## Prikaz kategorija zločina koji su se pojavili manje od 2000 puta unazad godinu dana

Promotrimo li kategorije zločina koje s frekvencijama manjim od 2000 (Figure \ref{fig:sub_2000}) uočavamo možemo vidjeti međusoban odnos frekventnosti težih kriminalnih radnji.
```{r sub_2000, fig.cap = "Kategorije kriminala s frekvencijom <2000", out.width='0.9\\linewidth'}
dfs <- df[df$Freq < 2000, ]
op <- par(mar = c(13, 4, 1, 2))
barplot(
  dfs$Freq,
  names.arg = dfs$Var1,
  las = 2,
  cex.names = 0.7,
  col = "skyblue",
  ylim = range(pretty(c(0, 2000))),
  main = "Crimes with a frequency <2000"
)


```
## Učitavanje podataka o vremenu

```{r}
#dohvaća listu datuma
dates <-
  strptime(c(chicago_year$DATE..OF.OCCURRENCE),
           format = "%m/%d/%Y %H:%M:%S %p",
           tz = "America/Chicago")

#dohvaća listu sati
hours <-
  as.numeric(format(
    strptime(
      chicago_year$DATE..OF.OCCURRENCE,
      format = "%m/%d/%Y %I:%M:%S %p",
      tz = "America/Chicago"
    ),
    format = "%H"
  ))
table(hours)
```

## Podaci o učestalosti kriminala ovisno o dobu dana

Promotrimo najprije ovisnost frekventnosti kriminalnih aktivnosti o trenutku u danu (Figure \ref{fig:frequency_hist1}).\
Inicijalna pretpostavka većine je vjerojatno da su kriminalne aktivnosti najfrekventnije u noćnim satima, ali histogram nam pokazuje drukčije.

```{r frequency_hist1, fig.cap = "Frekventnost kriminalnih radnji u periodu od 24 sata", out.width='0.9\\linewidth'}
hist(
  hours + 1,
  breaks = 0:24,
  main = "Frequency of criminal activities in 24 hours",
  xlab = "Hour",
  ylab = "Frequency",
  col = "light blue",
  xaxt = "n"
)
axis(side = 1,
     at = seq(0, 24, 2),
     tick = T)
rug(seq(1, 23, 2), ticksize = -0.03, side = 1)
```

Izdvojimo sad broj kriminalnih djela u 3 različita doba dana: noć koja traje od 22h do 6h, jutro koje traje od 6h do 14h i popodne koje smo uzeli da traje od 14h do 22h.
```{r}
day_parts<-ifelse(hours>=22|hours<6,rr2<-3,
                     ifelse(hours>=6&hours<14,rr2<-1,
                            rr2<-2))
mytable<-table(day_parts)
print(mytable)
```

Nakon što smo izdvojili podatke, možemo nacrtati histogram po dijelovima dana.

```{r frequency_hist_day, fig.cap = "Frekventnost kriminalnih radnji po dijelovima dana", out.width='0.9\\linewidth'}
hist(
  day_parts,
  breaks = 0:3,
  main = "Frequency of criminal activities for different parts of the day",
  xlab = "Day part",
  ylab = "Frequency",
  col = "light blue",
  xaxt = "n"
)
axis(side = 1,
     at = seq(0.5, 3, 1),
     labels = c("morning", "afternoon", "night")
     )
```
Iz histograma je očito da se najviše kriminalnih djela događa između 14 i 22h. Sada HI kvadrat testom možemo provjeriti uniformnost ove distribucije, odnosno vidjeti je li ta razlika statistički značajna ili nije. Pretpostavit ćemo da je distribucija uniformna.
Testirajmo hipotezu H0 jesu li kriminalna djela jednako distribuirana po svim djelovima dana, uz hipotezu H1 da nisu jednako distribuirana uz razinu značajnosti alpha = 0.05.

```{r chi squared test}
res<-chisq.test(as.integer(table(day_parts)), p = c(1/3, 1/3, 1/3))
res
```
S obzirom da je p-value < 2.2e-16 izrazito manja od alpha = 0.05, odbacujemo H0 u korist H1 te bismo mogli zaključiti da ova distribucija nije uniformna, odnosno da broj kriminalnih djela nije isti u svakom dobu dana. Naravno, treba pripaziti i na veličinu podataka. S obzirom na velik broj podataka, može se dogoditi da se odbaci H0 u korist H1 kad to nije trebalo napraviti, ali s obzirom na izgled dijagrama, slažemo se s odbacivanjem H0 u korist H1.

## Podaci o učestalosti kriminala ovisno o mjesecu u godini

Razmišljajući o prethodnom problemu, zaključujemo da bi bilo zanimljivo proučiti i ovisnost učestalosti kriminala o mjesecu u godini. U sljedećem dijelu, provjerit ćemo radi li se u ovom slučaju o uniformnoj razdiobi.
Izdvojimo za početak broj kriminala vezan uz svaki mjesec.

```{r}
#dohvaća listu mjeseci
months <-
  as.numeric(format(
    strptime(
      chicago_year$DATE..OF.OCCURRENCE,
      format = "%m/%d/%Y %I:%M:%S %p",
      tz = "America/Chicago"
    ),
    format = "%m"
  ))
table(months)
```
Podaci izgledaju kao da dolaze iz uniformne razdiobe, prikažimo ih histogramom bolje vizualizacije radi.

```{r criminal_action_per_month, fig.cap = "Frekventnost kriminalnih radnji svaki mjesec u godini", out.width='0.9\\linewidth'}
hist(
  months,
  breaks = 0:12,
  main = "Frequency of criminal activities every month in a year",
  xlab = "Month",
  ylab = "Frequency",
  col = "light blue",
  xaxt = "n"
)
axis(side = 1,
     at = seq(0.5, 12, 1),
     labels = c(1:12)
     )
```
Iz histograma distribucija izgleda donekle uniformno. Vidimo ponešto manje brojke u ožujku i travnju. Opet ćemo iskoristiti HI kvadrat test da provjerimo uniformnost ove distribucije.
Testirajmo hipotezu H0 jesu li kriminalna djela jednako distribuirana po svim mjesecima, uz hipotezu H1 da nisu jednako distribuirana uz razinu značajnosti alpha = 0.05.


```{r chi squared test 2}

```{r}

res<-chisq.test(as.integer(table(months)), p = c(rep(1/12, each=12)))
res
```
U ovom slučaju dobivamo nešto više od 3 puta manju vrijednost X-squared uz 11 stupnjeva slobode. p-value je poprilično mala te bismo mogli odbaciti H0 u korist H1. S obzirom da 3. i 4. mjesec odstupaju zbog pandemije koronavirusa, pokušajmo ih izbaciti iz podataka, te ponoviti test.

```{r}
months_filtered<-ifelse(months==1,rr2<-1,
                     ifelse(months==2,rr2<-2,
                            ifelse(months==5,rr2<-5,
                                   ifelse(months==6,rr2<-6,
                                          ifelse(months==7,rr2<-7,
                                                 ifelse(months==8,rr2<-8,
                                                        ifelse(months==9,rr2<-9,
                                                               ifelse(months==10,rr2<-10,
                                                                      ifelse(months==11,rr2<-11,
                                                                             ifelse(months==12,rr2<-12,
                            rr2<-NA))))))))))
table_months <- table(months_filtered)
table_months

```
```{r chi squared test 3}

res<-chisq.test(as.integer(table_months), p = c(rep(1/10, each=10)))
res
```
X-squared je sad pao za još 4 puta u odnosu na prošli test. Iako je p vrijednost i dalje mala, odučili smo ne odbaciti H0 zbog velikog skupa podataka. Smatramo da je distribucija po mjesecima uniformna te da je u 2020. godini bilo drugačije zbog lockdowna u ožujku i travnju, ali ostali mjeseci imaju približno uniformnu distribuciju.


## Učestalost krađa i kriminala vezanih uz narkotike

Dijagram (Figure \ref{fig:category_frequencies}) prikazuje učestalost različitih kategorija zločina unazad godinu dana. Najveću učestalost prema grafu ima kategorija krađe "THEFT", dok je kategorija "NARCOTICS" tek na 9. mjestu. Uz kategoriju krađe osim "THEFT" dodali smo "BURGLARY" i "ROBBERY" jer im je cilj krađa. Možemo reći da su specijalizacija krađe. Testom o jednoj proporciji provjeravamo je li razlika učestalosti krađa i kriminala vezanih uz narkotike statistički signifikantna, odnosno je li učestalost krađa doista veća od učestalosti kriminala vezanih za narkotike. Pretpostavit ćemo da se radi o normalnoj distribuciji.

Testiramo hipotezu H0: udio krađa je 0.5, nasuprot alternativi H1: udio krađa je veći od 0.5. Za razinu značajnosti alpha uzeli smo 0.05. 

```{r theft_narcotics, warning=FALSE, fig.cap = "Učestalost krađa i kriminala vezanih uz narkotike", out.width='0.9\\linewidth'}

theft_crimes <- c("THEFT", "ROBBERY", "BURGLARY")
theft <- chicago_year[chicago_year$PRIMARY.DESCRIPTION %in% theft_crimes,]
narcotic_crimes <- c("NARCOTICS", "OTHER NARCOTIC VIOLATION")
all_narocitcs <- chicago_year[chicago_year$PRIMARY.DESCRIPTION %in% narcotic_crimes,]
x <- nrow(theft) # broj krađa
n <- nrow(theft) + nrow(all_narocitcs) # veličina uzorka

res <- prop.test(x, n, p = 0.5, correct = TRUE,
              alternative = "greater")

res

```

Dobili smo p-vrijednost 2.2e-16, koja je puno manja od alpha = 0.05, iz tog razloga odbacujemo H0 u korist H1 te zaključujemo da su krađe doista učestalije. U prilog ovom zaključku ide i procjena vjerojatnosti krađe koja je 0.8836476 i nalazi se u 95%-tnom intervalu povjerenja 0.8816474 < 0.8836476 < 1.0000000.

## Podaci o ukradenoj vrijednosti

Nakon što smo vidjeli da je krađa najčešći oblik kriminala, zanimala nas je vrijednost ukradenih stvari. Dijagram (Figure \ref{fig:theft_crimes_values}) pokazuje da je učestalost većih krađa (preko $500) manja. Ponovo ćemo testom o jednoj proporciji provjeriti je li to stvarno istina. 

```{r theft_crimes_values, warning=FALSE, fig.cap = "Podaci o vrijednosti krađe", out.width='0.9\\linewidth'}


secondary_theft_types <- c("$500 AND UNDER", "OVER $500")
theft_secondary <- chicago_year[chicago_year$SECONDARY.DESCRIPTION %in% secondary_theft_types,]

table <- table(theft_secondary$SECONDARY.DESCRIPTION)
data <- as.data.frame.table(table)

barplot(
data$Freq,
names.arg = data$Var1,
cex.names = 0.6,
col = "skyblue",
ylim = range(pretty(c(0, data$Freq))),
main = "Theft value"
)
```

Testiramo hipotezu H0: udio krađa iznad \$500 je 0.5, nasuprot alternativi H1: udio krađa iznad \$500 je manji od 0.5. Za razinu značajnosti alpha uzeli smo 0.05.


```{r theft_value2, warning=FALSE, fig.cap = "Učestalost krađa", out.width='0.9\\linewidth'}

over500 <- c("OVER $500")

theft_more <- chicago_year[chicago_year$SECONDARY.DESCRIPTION %in% over500,]

x <- nrow(theft_more) # broj krađa većih od $500
n <- nrow(theft_secondary) # veličina uzorka

res <- prop.test(x, n, p = 0.5, correct = TRUE,
              alternative = "less")

res

```
Dobili smo p-vrijednost 2.2e-16, koja je znatno manja od alpha = 0.05, iz tog razloga odbacujemo H0 u korist H1 te zaključujemo da su krađe iznad \$500 manje učestale od krađa ispod \$500. Vjerojatnost krađe iznad $500 je 0.3745399. 


## Pad kriminala vezan uz koronavirus
Sljedeći graf prikazuje fascinantu činjenicu o smanjenju broja kriminalnih djela nakon početka pandemije koronavirusa i lockdowna. Prikaz je informativan te nećemo ulaziti u detaljniju analizu.

```{r theft_vs_narcotics, fig.cap = "Pad kriminala nakon pandemije koronavirusa - informativno", out.width='0.9\\linewidth'}

theft <- chicago_year[chicago_year$PRIMARY.DESCRIPTION == "THEFT",]
theft_dates <-
  table(as.Date(
    strptime(theft$DATE..OF.OCCURRENCE, format = "%m/%d/%Y %I:%M:%S %p", tz =
               "America/Chicago"),
  ))

narcotics <-
  chicago_year[chicago_year$PRIMARY.DESCRIPTION == "NARCOTICS",]
narcotics_dates <-
  table(as.Date(
    strptime(narcotics$DATE..OF.OCCURRENCE, format = "%m/%d/%Y %I:%M:%S %p", tz =
               "America/Chicago"),
  ))

all_crimes <-
  table(as.Date(
    strptime(chicago_year$DATE..OF.OCCURRENCE, format = "%m/%d/%Y %I:%M:%S %p", tz =
               "America/Chicago"),
  ))

plot(
  theft_dates,
  type = "l",
  main = "Criminal acts decrease during lockdown",
  xlab = "Date",
  ylab = "Number of crimes",
  ylim=c(0,750)
)
lines(all_crimes, col = "green")
lines(theft_dates, col = "red")
lines(narcotics_dates, col = "blue")
legend("topright", c("All crimes", "Theft", "Narcotics"), fill = c(rgb(0, 1, 0, 0.5), rgb(1, 0, 0, 0.5), rgb(0, 0, 1, 0.5)))


```


## Zavisnosti različitih socio-ekonomskih faktora
Prije nego što krenemo na veze između socio-ekonomskih faktora i kategorije kriminala, prvo nas zanima kakve su zavisnosti između različitih socio-ekonomskih faktora. U tu svrhu, napravili smo korelacijsku matricu.

```{r}
matrixData <- chicago_pc[, c(-1, -2, -3, -4)]
matrix <- cor(matrixData)
matrix
```
Uzmimo na primjer varijablu Below.Poverty.Level i pogledajmo kako se ona ponaša.
Najveću pozitivnu korelaciju smo dobili za par varijabli Below.Poverty.Level i Unemployment (0.76). Dosta dobru negativnu korelaciju dobili smo za par varijabli Below.Poverty.Level i Per.Capita.Income (-0.53). Možemo zaključiti da je nezaposlenost veliki razlog zašto netko živi ispod granice siromaštva. Iz negativne korelacije između Below.Poverty.Level i Per.Capita.Income mogli bismo pretpostaviti da su kvartovi eventualno grupirani po socijalnom statusu te bogatiji kvartovi imaju manje siromašnih ljudi.

Za varijable Unemployment i Crowded.Housing dobivamo slabu pozitivinu korelaciju (0.14). Za te smo varijable očekivali pozitivnu korelaciju. Naime, ako netko nema posao onda vjerojatno nema ni dovoljno za stanarinu koja bi odgovarala veličini prostora koji je potreban za zdrav život. Međutim broj koji smo dobili je puno manji od očekivanog.

Kako smo dobili broj koji je blizu 0, ovaj primjer je dobar kandidant da napravimo korelacijski test.

```{r}
cor.test(chicago_pc$Unemployment, chicago_pc$Crowded.Housing)
```
P-vrijednost testa je 0.2105. To znači da u slučaju da varijable Unemployment i Crowded.Housing nisu u nikakvoj korelaciji (tj. korelacija je 0), vjerojatnost da dobijemo pozitivnu korelaciju od 0.144 i više je 21%. Možemo posumnjati u pozitivnu korelaciju ove dvije varijable. 

Od drugih očekivanih vrijednosti imamo Per.Capita.Income i No.High.School.Diploma (-0.71). Ima smisla da su kvartovi siromašniji što je manji stupanj edukacije ljudi koji žive tamo.
Također ako manji Dependecy znači zdraviju i bogatiju ekonomiju, onda također očekujemo da će varijable Dependecy i Per.Capita.Income biti jako negativno korelirane što i jesu (-0.76).

Najveća vrijednost koja nije korelacija varijable sa samom sobom vidimo u slučaju No.High.School.Diploma i Crowded.Housing (0.91). Ovo je dosta zanimljiv i na prvi pogled jako neočekivani rezultat koji smo dobili. Sljedće što smo napravili je test.

```{r}
cor.test(chicago_pc$No.High.School.Diploma, chicago_pc$Crowded.Housing)
```
Vidimo da je 95%-tni interval između 0.85 i 0.93. Zaključujemo da su ove dvije varijable zaista jako zavisne. Iz tog razloga, možemo zaključiti da manjak obrazovanja najviše pridonosi prenatrpanim domovima. Međutim još uvijek ne znamo objasniti zašto je to tako. 

Ovo je dobar primjer da probamo napraviti linearnu regresiju i pogledamo kako to sve izgleda (Figure \ref{fig:regression}).

```{r}
no.high.school.diploma <-
  as.matrix(chicago_pc[8]) #No.High.School.Diploma
crowded.housing <- as.matrix(chicago_pc[6]) #Crowded.Housing

reg1 <- lm(crowded.housing  ~ no.high.school.diploma, chicago_pc)

summary(reg1)

```
```{r regression1, fig.cap = "Prikaz linearne regresije (Crowded Housing)", out.width='0.9\\linewidth'}
plot(
  no.high.school.diploma,
  crowded.housing ,
  pch = 16,
  cex = 1.3,
  col = "lightgreen",
  xlab = "Percentage of people without high school diploma",
  ylab = "Percentage of people living in crowded houses",
)
abline(lm(crowded.housing  ~ no.high.school.diploma), col = "pink", lwd = 2)
```

## Veze između socio-ekonomskih faktora i kategorije kriminala

## Kategorija Firearm related

Kategorija kriminala koju najprije promatramo je Firearm.related. Prvo ćemo napraviti korelacijsku matricu da dobijemo dojam kako se brojke kreću i od kojih varijabli bismo mogli krenuti raditi regresijski model.

```{r}
matrixData2 <- chicago_pc[, c(-1, -2, -3)]
matrix2 <- cor(matrixData2)
matrix2[,1]
```

Firearm.related i Unemployment imaju najveći korelacijski koeficijent pa ćemo započeti graditi regresijski model s tom varijablom.

```{r}
unemployment <- as.matrix(chicago_pc[10]) #Unemployment
firearm.related <- as.matrix(chicago_pc[4]) #Firearm.related
reg <- lm(firearm.related ~ unemployment, chicago_pc)

summary(reg)
```

```{r regression2, fig.cap = "Prikaz linearne regresije (Firearm)", out.width='0.9\\linewidth'}
plot(
  unemployment,
  firearm.related,
  pch = 16,
  cex = 1.3,
  col = "lightgreen",
  xlab = "Percentage of unemployed people",
  ylab = "Firearm related crimes per 100 000",
)
abline(lm(firearm.related ~ unemployment), col = "pink", lwd = 2)
```
Slika i regresija koju smo dobili nije loša, međutim odmah primjećujemo da imamo stršećih vrijednosti. Također vidimo da nam koeficijent b0 nije statistički značajan. Primjenjujući transformaciju nad varijablom Unemployment, pokušat ćemo dobiti nešto bolju regresiju.


```{r}
unemployment <- as.matrix(chicago_pc[10]) #Unemployment
firearm.related <- as.matrix(chicago_pc[4]) #Firearm.related
reg <- lm(firearm.related ~ log(unemployment), chicago_pc)

summary(reg)
```
Ako primijenimo logaritamsku transformaciju, koeficijent b0 postaje jako značajan, a i koeficijent determinacije se povećao, te F-test daje puno bolje rezultate. 


```{r regression, fig.cap = "Prikaz linearne regresije (Firearm)", out.width='0.9\\linewidth'}
plot(
  log(unemployment),
  firearm.related,
  pch = 16,
  cex = 1.3,
  col = "lightgreen",
  xlab = "Percentage of unemployed people",
  ylab = "Firearm related crimes per 100 000",
)
abline(lm(firearm.related ~ log(unemployment)), col = "pink", lwd = 2)
```
Na slici se još uvijek primjećju stršeće vrijednosti. 

Sljedeće dvije varijable koje su u velikoj korelaciji su Below.Poverty.Level i Dependency. Obje varijable su visoko korelirane s varijablom Unemployment što smo promatrali prije. 

```{r}
unemployment <- as.matrix(chicago_pc[10]) #Unemployment
dependecy <- as.matrix(chicago_pc[7]) #Dependency
firearm.related <- as.matrix(chicago_pc[4]) #Firearm.related
reg <-
  lm(firearm.related ~ log(unemployment) + dependecy, chicago_pc)

summary(reg)
```
 Ako dodamo varijablu Dependency u model, on se nije puno poboljšao, štoviše varijabla je statistički neznačajna.
 
```{r}

unemployment <- as.matrix(chicago_pc[10]) #Unemployment
dependecy <- as.matrix(chicago_pc[7]) #Dependency
firearm.related <- as.matrix(chicago_pc[4]) #Firearm.related
below.poverty.level <- as.matrix(chicago_pc[5]) #Below.Poverty.Level
reg <-
  lm(firearm.related ~ log(unemployment) + dependecy + below.poverty.level,
     chicago_pc)

summary(reg)
```
Te ako dodamo i varijablu Below.Poverty.Level dobivamo sličnu stvar što je i očekivano. Sve tri varijable su dosta zavisne te ne dobivamo puno novih informacija.
 
Iz tog razloga ćemo izbaciti te varijable iz modela.

Tražimo neku drugu varijablu koja nije jako kolerirana s našom prvom odabranom varijablom i tražimo onu koja daje trenutačno najbolji model. Dolazimo do rezultata da je to varijabla No.High.School.Diploma.
```{r}
unemployment <- as.matrix(chicago_pc[10]) #Unemployment
no.high.school.diploma <-
  as.matrix(chicago_pc[8]) #No.High.School.Diploma
firearm.related <- as.matrix(chicago_pc[4]) #Firearm.related
reg <-
  lm(firearm.related ~ log(unemployment) + no.high.school.diploma,
     chicago_pc)

summary(reg)
```
Tražimo postoji li još kakva varijabla koja bi mogla pridonjeti modelu.

Ako vizualiziramo parove Per.Capita.Income, Firearm.related možemo uočiti pravilnost. U siromašnijim kvartovima, kriminal koji uključuje vatreno oružje je puno češće nego u bogatijim kvartovima. Štoviše, uočavamo eksponencijalni pad pa ćemo primijeniti transformaciju podataka nad tom varijablom.
 
```{r fig.cap = "Odnos oružanih kriminala i GDP per capita", out.width='0.9\\linewidth'}
plot(
  chicago_pc$Per.Capita.Income,
  firearm.related,
  pch = 16,
  cex = 1.3,
  col = "lightgreen",
  xlab = "Per capita income",
  ylab = "Firearm related crimes per 100 000",
)
```

```{r}
unemployment <- as.matrix(chicago_pc[10]) #Unemployment
no.high.school.diploma <-
  as.matrix(chicago_pc[8]) #No.High.School.Diploma
firearm.related <- as.matrix(chicago_pc[4]) #Firearm.related
per.capita.income <- as.matrix(chicago_pc[9]) #Per.Capita.Income
reg <-
  lm(
    log(firearm.related) ~ log(unemployment) + no.high.school.diploma + log(per.capita.income),
    chicago_pc
  )

summary(reg)

```
Dodavanjem ostalih dviju varijabli Crowded.Housing i Below.Poverty.Level nismo dobili bolji model pa s ovime završavamo gradnju modela.

Sada je vrijeme da vidimo zadovoljala li model sve pretpostavke.

Provjeravamo homogenost varijanci. Kako bi vrijedila ta pretpostavka, trebamo dobiti što bolju moguću horizontalnu liniju. 

```{r fig.cap = "Prikaz standardiziranih reziduala", out.width='0.9\\linewidth'}
plot(reg, 3, col="lightgreen", pch=16)
```
Horizontalna linija je dobar znak da homogenost varijance vrijedi. \
Provjerimo sada to i testom.

```{r}
bptest(reg)
```
Test prolazi na razini značajnosti od 0.05.\
Dalje nas zanima normalnost reziduala.
```{r fig.cap = "Q-Q Plot", out.width='0.9\\linewidth'}
plot(reg, 2, col ="lightgreen", pch=16)
```
Reziduali ne prate najbolje liniju koju trebaju pa možemo pretpostaviti da reziduali nisu normalno distribuirani. To znači da naš model ne bi davao pouzdane rezultate kada bismo ga koristili za predikciju vrijednosti.\
Našu pretpostavku možemo još provjeriti Kolmogorov-Smirnovljevim testom.

```{r}
ks.test(rstandard(reg), 'pnorm')
```
Vidimo da model pada na KS testu te zaključujemo da reziduali nisu normalno distribuirani.\

Zadnje što provjeravamo je nezavisnost pomoću Durbin-Watsonovog testa.

```{r}
dwtest(reg)
```
Vidimo da reziduali nisu autokorelirani što znači da su nezavisni.

### Zaključak

Na kraju zaključujemo da na kriminal koji uključuje vatreno oružje utječe nezaposlenost, manjak više edukacije i GDP doprinos kvarta. Rezultati koji smo dobili nisu ništa neobičajni već konzistentni s onime što bi povezali s kriminalom. Loše socio-ekonomske prilike stvaraju atmosferu za kriminal.

## Kategorija Assault Homicide

Promotrimo sada kategoriju Assault Homicide i njezinu korelaciju s drugim faktorima. 

```{r}
matrixHomicide <- chicago_pc[,c(-1, -2, -4)]
corMatrixHomicide <- cor(matrixHomicide)
corMatrixHomicide[,1]

```
Kao i kod kriminala vezanih za vatrena oružja vidimo jaku korelaciju s nezaposlenosti, što sugerira da regresijski model ponovo gradimo od te dvije varijable.

```{r}
unemployment <- as.matrix(chicago_pc[10])
homicide <- as.matrix(chicago_pc[3])
regression <- lm(homicide ~ unemployment, chicago_pc)

summary(regression)
```
Kako bismo dobili bolje rezultate na naše podatke možemo primijeniti neku od transformacija.

```{r}
regression_sqrt <- lm(homicide ~ sqrt(unemployment), chicago_pc)
regression_log <- lm(homicide ~ log(unemployment), chicago_pc)

summary(regression_sqrt)
summary(regression_log)
```

Usporedbom rezultata uočavamo da je sqrt transformacija prikladnija u ovoj situaciji, također primijećujemo da smo dobili bolju signifikantnost varijable.  
Vizualizirajmo dobivenu regresiju.

```{r fig.cap = "Prikaz linearne regresije (Homicide)", out.width='0.9\\linewidth'}
plot(
  sqrt(unemployment),
  homicide,
  pch = 16,
  cex = 1.3,
  col = "cornflowerblue",
  xlab = "Unemployment percentage",
  ylab = "Homicide crimes per 100 000"
)
abline(regression_sqrt, col = "lightgreen", lwd = 2)
```
Sljedeća varijabla po razini korelacije Below poverty level. \
Ako promotrimo ovisnot stope ubojstava sa stopom siromaštva uočavamo dosta jasnu poveznicu.

```{r homicide_gdp, fig.cap = "Odnos stope ubojstava i stope siromaštva", out.width='0.9\\linewidth'}
homicide_gdp <-
  data.frame(x = chicago_pc$Below.Poverty.Level,
             y = chicago_pc$Assault..Homicide.)
plot(homicide_gdp,
     main = "Homicide rate in relation to GDP per capita",
     xlab = "Below poverty level",
     ylab = "Homicide rate (n/100 000)",
     col = "cornflowerblue",
     pch = 16)
```
Pokušajmo sada dodati varijablu u regresijski model.

```{r}
below_poverty <- as.matrix(chicago_pc[5])
regression_2 <-
  lm(homicide ~ sqrt(unemployment) + below_poverty, chicago_pc)
regression_2_sqrt <-
  lm(homicide ~ sqrt(unemployment) + sqrt(below_poverty), chicago_pc)
regression_2_log <-
  lm(homicide ~ sqrt(unemployment) + log(below_poverty), chicago_pc)

summary(regression_2)
summary(regression_2_sqrt)
summary(regression_2_log)
```
Usporedbom 3 mogućnosti kombiniranja varijabli, primijećujemo da zadnja opcija (kombinacija sqrt i log transformacije) daje najbolje rezultate.  
Pokušajmo sada dodati i treću varijablu po koreliranosti, Dependency.

```{r}
dependency <- as.matrix(chicago_pc[7])
regression_3 <-
  lm(homicide ~ sqrt(unemployment) + log(below_poverty) + dependency,
     chicago_pc)

summary(regression_3)
```
Kao i kod analize varijable Firearm related, za varijablu Homicide primijećujemo slično ponašanje - dodavanje varijable Dependency uzrokuje minimalna poboljšanja u modelu te ima vrlo malu signifikantnost.  
Nakon ispitivanja nekoliko modela pokazalo se da varijable per_capita_income i crowded housing također imaju vrlo malu signifikanost te ne poboljšavaju razine signifikantnosti drugih varijabli.
No dodamo li No high school diploma varijablu u model, ona ne da sama sebe ima veliku signifikantnost nego i značajno poboljšava signifikantnost varijabli dependency i below poverty level.  
Sljedeći model pokazao je najbolje performanse:

```{r}
no_diploma <- as.matrix(chicago_pc[8])

regression_4 <-
  lm(
    homicide ~ sqrt(unemployment) + log(below_poverty) + dependency + no_diploma,
    chicago_pc
  )

summary(regression_4)
```

```{r fig.cap = "Prikaz standardiziranih reziduala", out.width='0.9\\linewidth'}
plot(regression_4, 3, col="cornflowerblue", pch=16)
```
Breusch-Paganovim testom možemo provjeriti homoskedastičnost reziduala.

```{r}
bptest(regression_4)
```
Obzirom na činjenicu da linija u grafu iznad nije horizontalna i bilo je za očekivati da varijance nisu homogene.\
Primijenimo sqrt transformaciju na homicide varijablu.

```{r}
regression_4_sqrt <-
  lm(
    sqrt(homicide) ~ sqrt(unemployment) + log(below_poverty) + dependency + no_diploma,
    chicago_pc
  )

summary(regression_4_sqrt)
```
Promotrimo sada graf sa standardiziranim rezidualima:
```{r fig.cap = "Prikaz standardiziranih reziduala", out.width='0.9\\linewidth'}
plot(regression_4_sqrt, 3, col = "cornflowerblue", pch =16)
```
Vidimo da se pravac znatno izravnao u usporedbi s modelom. \
Pokažimo to i testom.

```{r}
bptest(regression_4_sqrt)
```
Test sada prolazi na razini značajnosti 0.05.

Durbin-Watsonovim testom možemo provjeriti postoji li koreliranost između reziduala.
```{r}
dwtest(regression_4_sqrt)
```
Test nam potvrđuje da nema koreliranosti između reziduala što je poželjan ishod kod regresijskih modela. \
Promotrimo sada QQ plot. Vidimo da reziduali poprilično dobro slijede pravac što bi sugeriralo njihovu normalnost.

```{r fig.cap = "Q-Q Plot", out.width='0.9\\linewidth'}
plot(regression_4_sqrt, 2, col = "cornflowerblue", pch=16)
```
Prethodnu pretpostavku možemo provjeriti Kolmogorov-Smirnovljevim testom.

```{r}
ks.test(rstandard(regression_4_sqrt), 'pnorm')
```
Test nam govori da su reziduali normalni. \
Kao dodatnu provjeru možemo provesti i Lillieforseovu inačicu testa.
```{r}
require(nortest)
lillie.test(rstandard(regression_4_sqrt))
```
Vidimo da Lillieforseov test ne prolazi, to se događa jer je on nešto stroži i precizniji od KS testa. \
Promotrimo li histogram reziduala vidimo da oni donekle oponašaju izgled normalne razdiobe.
```{r fig.cap = "Histogram reziduala", out.width='0.9\\linewidth'}
hist((regression_4_sqrt$residuals),
     col = 'cornflowerblue')
```

### Zaključak

Možemo zaključiti da na stopu ubojstava najviše utječu nezaposlenost i stopa ljudi koji žive ispod razine siromaštva, također znatan utjecaj pokazuju manjak više edukacije i radna sposobnost. \
Rezultati su očekivani i stvaraju jasnu poveznicu između negativnih karakteristika kvarta i stope ubojstava.