sap_projekt.Rmd

---
title: "SAP - projekt - Milijarderi"
knit: (function(input_file, encoding) {
  out_dir <- 'docs';
  rmarkdown::render(input_file,
 encoding=encoding,
 output_file=file.path(dirname(input_file), out_dir, 'projekt_sap_dokumentacija.pdf'))})
author: "Dora Bezuk, Marcela Matas, Josip Arelic, Domagoj Marinello"
date: "13.11.2022."
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r error=FALSE, message=FALSE, warning=FALSE, include=FALSE}
library("readxl")
library(tidyverse)
library(dplyr)
library(fastDummies)

```

# Uvod

Pitanja:

1. Ima li neki kontinent statistički značajno više milijarda?

2. Jesu li milijarderi koji su naslijedili bogatstvo statistički značajno bogatiji od onih koji nisu?

3. Možete li iz danih varijabli predvidjeti njihovo bogatstvo?

4. Kada biste birali karijeru isključivo prema kriteriju da se obogatite, koju biste industriju izabrali?

Dodatno pitanje:

5. Jesu li muškarci milijarderi statistički značajno bogatiji od žena milijardera?

# Deskriptivna analiza

Potrebno je učitati podatke.
```{r include=FALSE}
# Učitavanje podataka iz excel datoteke
# Promijeniti path u put do datoteke s podacima
bill_data <- read_excel("billionaires.xlsx")
```

```{r include=FALSE}
dim(bill_data) # dimenzije: 2614 redaka i 22 stupaca
names(bill_data) # imena stupaca
view(bill_data)

# klase pojedinih stupaca
sapply(bill_data, class)

# klasa tablice
class(bill_data) 
# zaključak: bill_data podaci su dobro učitani
```

```{r}
# Pomoćna funkcija za izbacivanje stršećih vrijednosti
remove_outliers <- function(data, data_column) {
  quartiles <- quantile(data_column, probs=c(.25, .75), na.rm = FALSE)
  IQR <- IQR(data_column)
  Lower <- quartiles[1] - 1.5*IQR
  Upper <- quartiles[2] + 1.5*IQR 
  
  return(subset(data, data_column >= Lower & data_column <= Upper))
}

cat('\n Dimenzija podataka: ', dim(bill_data))
```

Svaki milijarder(2614) u danim podatcima sadrži 21 atribut koji ga opisuje. Neki od njih su: broj godina, spol,
državljanstvo, porijeklo bogatstva, struka, vrijednost imovine, itd.

```{r}
for (col_name in names(bill_data)){
  if (sum(is.na(bill_data[,col_name])) > 0){
    cat('Ukupno nedostajućih vrijednosti za varijablu'
        ,col_name, ': ', sum(is.na(bill_data[,col_name])),'\n')
  }
}

```

Naš dataset sastoji se od character i numeric varijabli.
Prvo promotrimo numeričku varijablu wealth.
```{r}

hist(bill_data$`wealth.worth in billions` ,main='wealth worth in billions', xlab='wealth', ylab='Frequency', col="pink")


summary(bill_data$`wealth.worth in billions`)

```



Ovaj histogram nam prikazuje distribuciju bogatstva. Možemo vidjeti da se "wealth" ne ravna po normalnoj distribuciji i da je jako malo onih s velikom količinom bogatstva. To nam pokazuju i rezultati koje smo dobili - srednja vrijednost i medijan.



```{r}

barplot(table(bill_data$wealth.type),las=2,cex.names=.9,main='Wealth type',col="pink")


```



Na ovom grafu vidimo koji tip bogatstva je najzastupljeniji. Najviše ima onih koji su naslijedili bogatstvo.



```{r}
barplot(table(bill_data$wealth.how.industry),las=2,cex.names=.7,main='Industry',col="pink")


```

Na ovom grafu vidimo koje industrije su najzastupljenije kod milijardera. Najviše je onih u potrošačkoj industriji (consumer), a slijede ih oni koji se bave nekretninama, maloprodajom i restoranima. 




```{r}

table(bill_data$demographics.gender)

```


Ovdje vidimo da muškarci značajno prevladavaju u broju milijardera naspram žena i supružnika.

```{r}
median(bill_data$demographics.age)
mean(bill_data$demographics.age)



```

Ako idemo proučiti milijardere po godinama, možemo vidjeti da je izračunata srednja vrijednost za varijablu starost 53 godine. Medijan iznosi 59 godina, odnosno kada bi ih poredali od najmlađeg do najstarijeg, vrijednost u sredini bila bi 59 godina.

# Pitanja

## 1. Ima li neki kontinent statistički značajno više milijardi?

Za početak želimo vidjeti je li svim milijarderima u našem datasetu dodijeljen kontinent. S obzirom na to da kontinent kao varijabla ne postoji, koristit ćemo regiju (location.region). Sada želimo izlistati sve regije koje postoje u datasetu.
```{r}
levels(factor(bill_data$location.region))
```



Ima li nedostajućih vrijednosti? 
```{r}
# is.na ce nam vratiti logical vektor koji ima TRUE na mjestima gdje ima NA:
sum(is.na(bill_data$location.region)) 
 

```

Nema nedostajućih vrijednosti
```{r}
table(bill_data$location.region)



```
S obzirom na to da imamo regiju Middle East/North Africa trebamo ih rastaviti na kontinente kojima pripadaju (Azija i Afrika). Prvo želimo vidjeti koje sve države postoje u toj regiji u našem datasetu pomoću državljanstva.
Također imamo jednu državu čija regija ima vrijednost 0, a država je Bermuda. Dakle, nju ćemo kasnije u kodu svrstati pod kontinent Sjeverna Amerika.


```{r, results='hide'}


bill_data$location.citizenship[bill_data$location.region == "Middle East/North Africa"]
bill_data$location.citizenship[bill_data$location.region == "0"]

```

Sada možemo združiti podatke ovisno o kontinentu.

Kopirajmo najprije podatke u novi data.frame (bill_data_copy).
```{r}
bill_data_copy = data.frame(bill_data)
tracemem(bill_data)==tracemem(bill_data_copy)
untracemem(bill_data_copy)
untracemem(bill_data_copy)
```



```{r,results='hide'}
# Združimo Europu
for (column_name in c("Europe")){
  bill_data_copy$location.region[bill_data_copy$location.region == column_name] = "Europe";
}

# Združimo Afriku
for (column_name in c("Lebanon","Egypt","Morocco","Algeria")){
  bill_data_copy$location.region[bill_data_copy$location.citizenship == column_name] = "Africa";
}

for (column_name in c("Sub-Saharan Africa")){
  bill_data_copy$location.region[bill_data_copy$location.region == column_name] = "Africa";
}

# združimo Sjevernu Ameriku
for (column_name in c("North America")){
  bill_data_copy$location.region[bill_data_copy$location.region == column_name] = "North America";
  
  for (column_name in c("Bermuda")){
  bill_data_copy$location.region[bill_data_copy$location.citizenship == column_name] = "North America";
}
}

# Združimo Južnu Ameriku
for (column_name in c("Latin America")){
  bill_data_copy$location.region[bill_data_copy$location.region == column_name] = "South America";
}

# Združimo Aziju
for (column_name in c("East Asia","South Asia")){
  bill_data_copy$location.region[bill_data_copy$location.region == column_name] = "Asia";
}
for (column_name in c("Saudi Arabia","Kuwait","United Arab Emirates","Israel","Turkey","Oman","Bahrain")){
  bill_data_copy$location.region[bill_data_copy$location.citizenship == column_name] = "Asia";
}

#Združimo Australiju
for (column_name in c("Australia")){
  bill_data_copy$location.region[bill_data_copy$location.citizenship == column_name] = "Australia";
}


bill_data_copy

```
```{r}
tbl = table(bill_data_copy$location.region)
print(tbl)

```





Sada kad smo završili s pripremom podataka za ovaj zadatak, možemo započeti sa statističkim testovima.
ANOVA je parametarski test kojim se uspoređuju srednje vrijednosti više uzoraka te se na temelju F-testa donosi zaključak o postojanju značajnih razlika između tih srednjih vrijednosti. Na taj se način analizira utjecaj jedne ili više nezavisnih varijabli na jednu numeričku kontinuiranu (zavisnu) varijablu.
U ovom slučaju razmatrat ćemo location.region kao varijablu koja određuje grupe (populacije) i wealth kao zavisnu varijablu.
Pretpostavke jednofaktorske ANOVA-e su:

- nezavisnost pojedinih podataka u uzorcima,
- normalna razdioba podataka,
- homogenost varijanci među populacijama. 

Nezavisnost uzoraka je zadovoljena s obzirom na to da jedna osoba ne potječe s više kontinenata.

Kad su veličine grupa podjednake, ANOVA je relativno robusna metoda na blaga odstupanja od pretpostavke normalnosti i homogenosti varijanci. Ipak, dobro je provjeriti koliko su ta odstupanja velika. 

Provjeru normalnosti za svaku pojedinu grupu napravit ćemo Lillieforsovom inačicom KS testa. 

Pretpostavke Lillieforsevog testa:

$$ \begin{aligned}
H_0&: \ podaci\  se \ ravnaju\ po \ normalnoj\ distribuciji\\
H_1&: \ podaci\  se \ ne\  ravnaju\ po \ normalnoj\ distribuciji\\
\end{aligned} $$


```{r, fig.height = 4, fig.width = 5}

# logaritmirali smo wealth kako bi dobili ljepšu distribuciju na grafovima

wealth <- log(bill_data_copy$wealth.worth.in.billions, 2)

require(nortest)
lillie.test(wealth)

lillie.test(wealth[bill_data_copy$location.region=='Africa'])
lillie.test(wealth[bill_data_copy$location.region=='Europe'])
lillie.test(wealth[bill_data_copy$location.region=='South America'])
lillie.test(wealth[bill_data_copy$location.region=='North America'])
lillie.test(wealth[bill_data_copy$location.region=='Asia'])
lillie.test(wealth[bill_data_copy$location.region=='Australia'])


hist(wealth[bill_data_copy$location.region=='Africa'], main = "Histogram of wealth in Africa", xlab="Wealth", col= "pink")
hist(wealth[bill_data_copy$location.region=='Europe'], main = "Histogram of wealth in Europe", xlab="Wealth", col= "pink")
hist(wealth[bill_data_copy$location.region=='South America'], main = "Histogram of wealth in South America", xlab="Wealth", col= "pink")
hist(wealth[bill_data_copy$location.region=='North America'], main = "Histogram of wealth in North America", xlab="Wealth", col= "pink")
hist(wealth[bill_data_copy$location.region=='Asia'], main = "Histogram of wealth in Asia", xlab="Wealth", col= "pink")
hist(wealth[bill_data_copy$location.region=='Australia'], main = "Histogram of wealth in Australia", xlab="Wealth", col= "pink")

```





Po rezultatima testa za normalnost (p vrijednosti manje od 0.05) te dobivenim histogramima vidimo da nam normalnost i nije zadovoljena. Nastavit ćemo s provjerom homogenosti varijanci Bartlettovim testom. Njegove pretpostavke su:
$$ \begin{aligned}
H_0&: \ varijance\ su\ jednake\\
H_1&: \ varijance\ se \ razlikuju \\
\end{aligned} $$

```{r test homogenosti}

# Testiranje homogenosti varijance uzoraka Bartlettovim testom

bartlett.test(wealth ~ bill_data_copy$location.region )


var((wealth[bill_data_copy$location.region=='Africa']))
var((wealth[bill_data_copy$location.region=='Asia']))
var((wealth[bill_data_copy$location.region=='Europe']))
var((wealth[bill_data_copy$location.region=='North America']))
var((wealth[bill_data_copy$location.region=='South America']))
var((wealth[bill_data_copy$location.region=='Australia']))

```


S obzirom na to da na temelju ovih rezultata(p=0.0009 < 0.05) odbacujemo nultu hipotezu da su nam varijance jednake, ne možemo koristiti ANOVA test. Prvo ćemo prikazati na boxplot grafu ovisnost varijable bogatstva o kontinentu kako bi grafički mogli interpretirati rezultat, a zatim ćemo utvrditi sigurnost rezultata koristiti neparametarski Kruskal - Wallis test kao alternativu jednofaktorskoj ANOVA-i.


 

```{r}

# Graficki prikaz podataka
boxplot(wealth ~ bill_data_copy$location.region, xlab="Continent", ylab="Wealth", col= "pink")

```
Izračunajmo sada srednje vrijednosti bogatstva za svaki kontinent.


```{r}

mean_all= mean(bill_data_copy$wealth.worth.in.billions)
mean_all
mean_by_continent <- bill_data_copy %>% 
  group_by(location.region) %>% 
  summarize(mean_continent = mean(wealth.worth.in.billions))
mean_by_continent



```
Razmatrajući graf i numeričke rezultate koje smo dobili utvrdili smo da ima manje razlike između srednjih vrijednosti varijable wealth podijeljene po kontinentima.Sada idemo vidjeti je li ta razlika statistički značajna. Ovaj test služi za testiranje jednakosti srednjih vrijednosti u jednofaktorskoj analizi varijance.
Pretpostavke Kruskal -Wallis testa :


$$ \begin{aligned}
H_0&: \ nema\ razlike\ među \ populacijama\\
H_1&: \ postoji \ razlika \ među\ populacijama\ (barem\ dvije\ \ se \ razlikuju)\\
\end{aligned} $$

```{r}

# Alternativa ANOVI - Kruskal - Wallis test

kruskal.test(bill_data_copy$wealth.worth.in.billions ~ bill_data_copy$location.region)


```

Kako je p-vrijednost manja od nivoa značajnosti od 0.05, možemo zaključiti da ima statistički značajne razlike između milijardi dijeljenim po kontinentima (među barem 2 kontinenta). Dakle, pomoću našeg izračuna srednjih vrijednosti i grafa možemo utvrditi da Sjeverna Amerika prednjači u količini bogatstva, a odmah iza nje nalazi se Europa.



## 2. Jesu li milijarderi koji su naslijedili bogatstvo statistički značajno bogatiji od onih koji nisu?

Potrebno je pripremiti podatke za obradu, razdvojiti podatke iz tablice po polju
how.inherited u dva slučaja: inherited (oni koju su naslijedili bogatstvo) i 
non_inherited (oni koji nisu naslijedili bogatstvo).

```{r, results='hide'}
inherited = bill_data[bill_data$wealth.how.inherited!="not inherited",]
non_inherited = bill_data[bill_data$wealth.how.inherited=="not inherited",]
```

Zatim je potrebno izračunati srednju vrijednost (mean) posebno za svaki slučaj
uzimajući u obzir polje worth.in billions.

```{r}
inherited_mean = mean(inherited$`wealth.worth in billions`)
print(inherited_mean)

non_inherited_mean = mean(non_inherited$`wealth.worth in billions`)
print(non_inherited_mean)

```

Na temelju male razlike u srednjim vrijednostima, ne postoje indikacije da su 
milijarderi koji su naslijedili bogatstvo statistički značajno bogatiji 
od onih koji nisu. No, navedeno je potrebno provjeriti.


Kako bi bolje vizualizirali podatke crtamo histogram i box plot za svaki od
slučaja:

```{r}
hist(inherited$`wealth.worth in billions`, breaks = 20, main = "Histogram of wealth that is inherited", xlab="inherited wealth", col= "pink")
boxplot(inherited$`wealth.worth in billions`, col = "pink")

hist(non_inherited$`wealth.worth in billions`, breaks = 20, main = "Histogram of wealth that is not inherited", xlab="non inherited wealth", col= "pink")
boxplot(non_inherited$`wealth.worth in billions`, col ="pink")
```

Iz prikazane vizualizacije uočavamo kako se podaci ne ravnaju po
normalnoj distribuciji.

Što se može bolje vidjeti sa sljedećih prikaza:
```{r}
qqnorm(inherited$`wealth.worth in billions`, pch = 1, frame = FALSE,main='Inherited')
qqline(inherited$`wealth.worth in billions`, col = "blue", lwd = 2)

qqnorm(non_inherited$`wealth.worth in billions`, pch = 1, frame = FALSE,main='Non inherited')
qqline(non_inherited$`wealth.worth in billions`, col = "red", lwd = 2)

```

Ipak, uočeno je potrebno dodatno ispitati koristeći Kolmogorov–Smirnov test 
kojim se utvrđuje ravna li se distribucija po normalnoj razdiobi.

```{r, warning=FALSE}
ks.test(inherited$`wealth.worth in billions`, y="pnorm")
ks.test(non_inherited$`wealth.worth in billions`, y="pnorm")
```

Iz dobivenih p vrijednosti u oba slučaja odbacujemo mogućnost da se distribucije
ravnaju po normalnoj razdiobi. 

Time je potvrđena pretpostavka da se podaci ne ravnaju po normalnoj distribuciji.

Potrebno je koristiti neparametarski test Mann–Whitney U test, koji se koristi 
kada se podaci se ravnaju po istim distribucijama (obje distribucije su nakošene
u desno) i uzorci su nezavisni iz jedne i druge populacije (jedna osoba ne može
naslijediti i nenaslijediti bogatstvo).


Hipoteze glase:
$$ \begin{aligned}
H_0&: \mu_1 = \mu_2 \\
H_1&: \mu_1 > \mu_2 \quad \quad 
\end{aligned} $$


```{r}
 
wilcox.test(inherited_mean, non_inherited_mean, alt = "greater")

```

Zbog p-vrijednost jednake 0.5, na temelju značajnosti od 50% ne možemo odbaciti 
$H_0$ hipotezu o jednakosti prosječnih vrijednosti bogatstva u korist $H_1$, 
odnosno možemo reći da milijarderi koji su naslijedili bogatstvo nisu statistički
značajno bogatiji od onih koji nisu.



## 3. Možete li iz danih varijabli predvidjeti njihovo bogatstvo?

Cilj ovog pitanja je provjeriti postoji li statistički značajna veza između više ulaznih varijabli (regresora) i izlazne varijable (reakcije,  `wealth.worth in billions`). Korištenjem modela linearne regresije provjerit ćemo koji regresori najviše utječu na izlaznu varijablu.

Pretpostavke modela:

- linearnost veze X i Y
- pogreške nezavisne, homogene i normalno distribuirane s $\epsilon \sim N (0, \sigma^2 )$

Za predobradu podataka radimo sljedeće stvari:

- Izbacujemo nepotrebne regresore:
  - name
  - company.name
  - rank
  - location.gdp: više od pola vrijednosti su 0 (netočan podatak)
  - location.country.code i location.citizenship: koristimo location.region koji je veće granulacije
  - wealth.how.from emerging, wealth.how.was founder, wealth.how.was political: konstantne varijable
  - company.sector: jer ima previše različitih vrijednosti, koje kad bi one hot encodali bi dali previše stupaca
  - wealth.type_inherited, već sadržan u `inherited`
- ignoriramo uzorke s netočnim podacima (kriva dob)
- povećavamo granulaciju varijable relationship (slične/iste vrijednosti svodimu na jednu)
- izbacujemo manji broj uzoraka koji sadrži null vrijednosti

Sve kategorijske varijable obrađujemo tako da ih pretvorimo u dummy varijable. Svaka kategorijska varijabla predstavljena je svojom novom vlastitom varijablom koja poprima vrijednost 1 u slučaju da originalna kategorijska varijabla odgovara novoj varijabli, inače je vrijednost 0.

Za filtrirani podatkovni skup sve iznose `wealth.worth in billions` množimo s 1 + (kupovna moć dolara svedena na godinu 2014).

```{r}
exclude_cols = c("name", "company.name", "rank", "location.gdp", "location.country code", "location.citizenship", "wealth.how.from emerging", "wealth.how.was founder", "wealth.how.was political", "company.sector")

bill_data_clean <- bill_data %>% select(-one_of(exclude_cols)) %>% arrange(year)

bill_data_clean[["company.relationship"]] <- tolower(bill_data_clean[["company.relationship"]] )

bill_data_clean <- bill_data_clean %>% filter(demographics.age > 0)
bill_data_clean <- bill_data_clean %>% filter(!location.region == "0")

# inflation rate $1.00 (1996) -> $1.51 (2014), +50.9%
# inflation rate $1.00 (2001) -> $1.34 (2014), +33.7%
bill_data_clean[bill_data_clean$year == "1996", "wealth.worth in billions"] <- bill_data_clean[bill_data_clean$year == "1996", "wealth.worth in billions"] * 1.509
bill_data_clean[bill_data_clean$year == "2001", "wealth.worth in billions"] <- bill_data_clean[bill_data_clean$year == "2001", "wealth.worth in billions"] * 1.337

# Iskoristili smo godinu da ažuriramo cijene (inflacija), sad ju odbacujemo
bill_data_clean <- bill_data_clean %>% select(., -year)

bill_data_clean$company.relationship <- gsub(".*\b(owner)\b.*", "owner", bill_data_clean$company.relationship)
bill_data_clean$company.relationship <- gsub(".*(ceo|chief executive officeor|chief executive officer|chief executive|exectuitve).*", "ceo", bill_data_clean$company.relationship)
bill_data_clean$company.relationship <- gsub(".*(founder).*", "founder", bill_data_clean$company.relationship)
bill_data_clean$company.relationship <- gsub(".*(chair|chari).*", "chairman", bill_data_clean$company.relationship)
bill_data_clean$company.relationship <- gsub(".*(director).*", "director", bill_data_clean$company.relationship)
bill_data_clean$company.relationship <- gsub(".*(head).*", "head", bill_data_clean$company.relationship)
bill_data_clean$company.relationship <- gsub(".*(president).*", "president", bill_data_clean$company.relationship)

bill_data_clean <- bill_data_clean %>% drop_na()

bill_categorical <- bill_data_clean %>% select(where(is_character))
bill_numeric <- bill_data_clean %>% select(where(is.numeric))
bill_categorical_onehot = dummy_cols(bill_categorical, remove_first_dummy = TRUE, remove_selected_columns = TRUE)
bill_categorical_onehot <- bill_categorical_onehot[, colSums(bill_categorical_onehot) > 5]
bill_data_clean <- bind_cols(bill_numeric, bill_categorical_onehot)
```


Bitna pretpostavka multivarijatne linearne regresije je da ne postoji snažna linearna korelacija regresora modela. U ovom koraku provjerit ćemo postoje li parovi takvih regresora i otkloniti ih ako postoje. Odbacit ćemo sve regresore za koje postoji neki drugi regresor čiji je apsolutna vrijednost Pearsonovog koeficijenta korelacije veća od 0.9. 

```{r}
correlation_threshold = 0.9
tmp <- corr_table <- cor(bill_data_clean)
tmp[upper.tri(tmp)] <- 0
diag(tmp) <- 0  # clean diagonal which is always 1
bill_data_clean <- bill_data_clean[, apply(tmp,2,function(x) all(x<= correlation_threshold))]

bill_data_clean <- remove_outliers(bill_data_clean, bill_data_clean$`wealth.worth in billions`)
wealth <- bill_data_clean$`wealth.worth in billions`
```


Prije stvaranja linearnog modela pogledajmo kojih 5 varijabli najviše linearno korelira sa `wealth.worth in billions`. Rezultate koje dobijemo ne možemo direktno koristiti za statističko zaključivanje, ali možemo kasnije usporediti linearne korelacije s rezultatima i zaključcima koje ćemo dobiti nakon stvaranja linearnog modela.

```{r}
w <- corr_table[, "wealth.worth in billions"]
w <- abs(w)

corr_wealth_vars <- w[order(w, decreasing = TRUE)]
cat("")
corr_wealth_vars[2:6]
```

```{r}
# x setup, y = wealth
normalized<-function(y) {
  x<-y[!is.na(y)]
  x<-(x - min(x)) / (max(x) - min(x))
  y[!is.na(y)]<-x
  return(y)
}

# `wealth.how.industry_Retail, Restaurant` casues fitting issues
exclude_cols = c("wealth.worth in billions", "wealth.how.industry_Retail, Restaurant", "wealth.type_inherited")
x <- bill_data_clean %>% select(-one_of(exclude_cols))
x[, c("company.founded", "demographics.age")] <- apply(x[, c("company.founded", "demographics.age")] , 2 , normalized) # minmax scaling
x <- x[,order(colnames(x))]

```


Prvi linearni model stvorili smo naivno tako da smo iskoristili sve moguće regresore. Provjerom vrijednosti `adj.r.squared` saznat ćemo koliki postotak varijance u podacima opisuje stvoreni linearni model. Također, provjerit ćemo koji regresori objašnjavaju najveći postotak varijance tako da ih poredamo po p vrijednostima.

```{r}
cat("Ukupan broj regresora:", length(colnames(x)), "\n")
p_value_column <- 4
model_all_vars <- lm(wealth ~ . , x)
sa <- summary(model_all_vars)
cat("Postotak varijance objašnjen linearnim modelom", sa$adj.r.squared * 100, "%\n")
coef <- sa$coefficients
coef_sorted <- coef[order(coef[,p_value_column]),]
cat("Prvih 5 regresora sortiranih uzlazno po p vrijednosti:\n")
coef_sorted[1:5,]
```

U ovom slučaju, najveći postotak varijance u podacima za izlaznu varijablu `wealth.worth in billions` objašnjava regresor `demographics.age`. Trenutan model objašnjava svega 7.7% varijance u podacima (Adjusted $R=0.07766$) za reakciju `wealth.worth in billions`. Na žalost, ovaj model ne objašnjava veliki dio varijance u podacima.

Sada ćemo pokušati pronaći najbolje prediktore na sljedeći način: stvarat ćemo model linearne regresije za svaki regresor pojedinačno, i očitavati koliko oni statistički značajno objašnjavaju varijancu u podacima (očitavamo p vrijednosti svakog modela nakon fittanja).


```{r}

n = 10
filtered_col_names = c()
r_squares = c()
ps = c()
col_names=colnames(x)

for(i in 1:ncol(x)){
  
  col_name=col_names[i]
  model=lm(wealth ~ x[[col_name]]) # napravi lienarni model s jednim regresorom
  summary_model = summary(model)
  
  filtered_col_names <- append(filtered_col_names, col_name)
  r_squares <- append(r_squares, summary_model$r.squared)
  ps <- append(ps,  summary_model$coefficients[,4][2])
}

df_g_squares=data.frame(filtered_col_names, r_squares, ps)
df_top_predictors = df_g_squares[order(df_g_squares$ps), ]
df_top_predictors[1:n, ]
top_n_predictors_one_var_lin <- df_top_predictors[1:n, "filtered_col_names"]
```


Možemo zaključiti da kad bismo morali napraviti linearni model koji najbolje predviđa reaktor `wealth.worth in billions`, odabrali bismo upravo regresor `wealth.how.inherited_not inherited`. Međutim ako želimo napraviti multivarijatni linearni model, nije nužno istina da će najbolji model biti onaj za koji uzmemo prvih n regresora iz trenutne liste. Problem koji se može pojaviti takvim pristupom odabira regresora je da postoje regresori koji su međusobno zavisni (iako smo već prethodno otklonili jako zavisne regresore).

Najbolje regresore također možemo pronaći ANOVA-om. Kada dodamo ili izbrišemo prediktivnu varijablu iz linearne regresije, želimo znati je li ta promjena poboljšala model ili nije. ANOVA uspoređuje dva regresijska modela i javlja jesu li značajno različiti. Spojit ćemo najbolje regresore koje smo dobili ANOVA-om i najbolje regresore dobivene u prethodnom koraku da stvorimo novi linearni model s manje regresora.

```{r}
a <- anova(model_all_vars)
ps_a <- a$`Pr(>F)`
ps_a <- head(ps_a, -1) # anova returns NA for last element

ps_a_ord <- order(ps_a)
sorted_cols <- colnames(x)[order(colnames(x))]
top_predictors_anova <- sorted_cols[ps_a_ord][1:n]

top_predictors = c(top_predictors_anova, top_n_predictors_one_var_lin)
top_predictors <- top_predictors[!duplicated(top_predictors)]

cat("Broj regersora", length(top_predictors))
model_top_preds <- lm(wealth ~ . , x[, top_predictors])
summary(model_top_preds)
```

Pokazali smo da smanjenjem broja regresora na 14 i dalje objašnjavamo usporedivo veliki dio varijance (6.1%, originalno 7.7%). Ovisno o namjeni i potrebama možemo se opredijeliti za složeniji ili jednostavniji model. Jednostavniji model je preferiraniji ako je relativno dobar kao neki alternativni složeniji model.


Za kraj, provjerit ćemo pretpostavku linearnog modela (normalnost reziduala) grafički i korištenjem Kolmogorov-Smirnov i Lilliefors testa.

```{r, collapse=TRUE}
require(nortest)
hist(rstandard(model_top_preds), col="pink")
qqnorm(rstandard(model_top_preds))
qqline(rstandard(model_top_preds), col = "blue", lwd = 2)
ks.test(rstandard(model_top_preds),'pnorm')
lillie.test(rstandard(model_top_preds))
```

Iz histogram se može naslutiti da se distribucija reziduala ne ravna po normalnoj distribuciji. Vrijednosti nisu centrirane oko nule i uočavamo debele desne repove. Iz qq grafa jasno uočavamo problem teškog desnog repa i manje problematičnog laganog lijevog repa. Ovaj graf dodatno potvrđuje da se reziduali ne ravnaju po normalnoj distribuciji. Konačno, oba testa za normalnost ukazuju da se reziduali ne ravnaju po normalnoj distribuciji jer je p vrijednost je manja od 0.05.

S obzirom na to da je pretpostavka linearnog modela prekršena i da linearan model objašnjava svega 6.2% varijance za reaktor `wealth.worth in billions` odbacujemo mogućnost da linearnim modelom previđamo bogatstvo koristeći preostale varijable iz skupa podataka.



## 4. Kada biste birali karijeru isključivo prema kriteriju da se obogatite, koju biste industriju izabrali?

Pretpostavljamo da karijerom u određenoj industriji, a ne nasljedstvom zarađujemo novac. Zbog toga gledamo samo milijardere koji nisu naslijedili svoje bogatstvo.
Također, zanimaju nas samo najnoviji milijarderi odnosno oni s popisa iz 2014. godine, jer su vremenski najrelevantniji. 
Uz taj skup milijardera, uzet ćemo i skup milijardera iz 2014. koji nisu bili na popisu 2001., odnosno novonastale milijardere. Na tom skupu vidjet će se u kojim industrijama je nastalo najviše milijardera.
Za kraj ćemo uzeti u razmatranje i skup milijardera koji su bili na popisu 2001. godine, ali zbog određenog razloga više nisu. Tu ćemo vidjeti koje su industrije u tom razdoblju izgubile najviše milijardera.



```{r, results='hide'}
non_inherited_2014 <- non_inherited[non_inherited$year == 2014,]
non_inherited_2001 <- non_inherited[non_inherited$year == 2001,]
non_inherited_2014_new = bill_data[FALSE,]
non_inherited_2001_old = bill_data[FALSE,]

# selekcija novonastalih milijardera iz 2014. koji nisu bili na prethodnoj listi iz 2001.
for(i in 1:nrow(non_inherited_2014)) {
  r <- non_inherited_2014[i,]
  if(sum(str_detect(non_inherited_2001$name, r[[1]])) == 0) {
    non_inherited_2014_new <- rbind(non_inherited_2014_new, non_inherited_2014[i,])
  }
}

# selekcija milijardera iz 2001. koji nisu na listi iz 2014.
for(i in 1:nrow(non_inherited_2001)) {
  r <- non_inherited_2001[i,]
  if(sum(str_detect(non_inherited_2014$name, r[[1]])) == 0) {
    non_inherited_2001_old <- rbind(non_inherited_2001_old, non_inherited_2001[i,])
  }
}
```




```{r}

par(mar=c(10,7,1,1))
barplot(sort(table(subset(non_inherited_2014$wealth.how.industry, non_inherited_2014$wealth.how.industry != "0")), decreasing = TRUE),
        main = "Distribution of billionaires by industry in 2014",
        las = 2, col="pink")
```


Iz stupčastog grafa je vidljivo da su tri najzastupljenije industrije maloprodaja (trgovački lanci, lanci restorana), trgovina nekretninama i računalna tehnologija.



```{r}
par(mar=c(10,5,1,1))
barplot(sort(table(subset(non_inherited_2014_new$wealth.how.industry, non_inherited_2014_new$wealth.how.industry != "0")), decreasing = TRUE),
        main = "Distribution of newly made billionaires by industry in 2014",
        las = 2, col ="pink")
```


Ako usporedimo ovaj graf s prethodnim može se vidjeti da su vrlo slični, jedina razlika je u poretku medijske industrije. To nam govori da se broj milijardera po industrijama mijenja otprilike istom brzinom, odnosno da industrije s najviše milijardera dobivaju najveći broj novih milijardera (i obrnuto).



```{r}
par(mar=c(10,5,1,1))
barplot(sort(table(subset(non_inherited_2001_old$wealth.how.industry, non_inherited_2001_old$wealth.how.industry != "0")), decreasing = TRUE),
        main = "Distribution of ex-billionaires by industry in 2001",
        las = 2, col="pink")
```



Na posljednjem grafu možemo potvrditi prethodno uočeno kretanje medijske industrije. Za razliku od ostalih industrija, broj milijardera u medijskoj industriji jedini je toliko značajno pao da je medijska industrija pala u ukupnom poretku. Sukladno tome na ovom grafu vidimo da je medijska industrija doživjela neproporcionalan pad broja milijardera. Međutim, medijska industrija nije među najvećim industrijama tako da ne utječe na zaključak, odnosno tri najzastupljenije industrije nisu se promijenile.

Zaključno, industrije koje se mogu predložiti na temelju ovih grafova su ponajprije maloprodaja, trgovanje nekretninama i računalna tehnologija. Te industrije najbolji su odabir za početak karijere s ciljem postajanja milijarderom ponajprije zbog najveće količine milijardera u tim industrijama (i sukladno tome najvećeg broja novonastalih milijardera).



## 5. Jesu li muškarci milijarderi statistički značajno bogatiji od žena milijardera?

S obzirom na to da muškaraca milijardera ima značajno više nego žena (rezultate smo dobili u deskriptivnoj analizi), čak 10 puta više, zanima nas jesu li onda oni i uspješniji,odnosno bogatiji u prosjeku od žena milijardera.

Postoje nedostajuće vrijednosti koje moramo izbaciti u varijabli demographics.gender(NA).
```{r}
bill_data_gender= bill_data_copy%>% filter(!is.na(demographics.gender))



women = bill_data_gender[bill_data_gender$demographics.gender == "female",]
men = bill_data_gender[bill_data_gender$demographics.gender == "male",]


```

```{r}
women_mean = mean(women$wealth.worth.in.billions)
women_mean
men_mean = mean(men$wealth.worth.in.billions)
men_mean
```
Izračunali smo srednje vrijednosti bogatstva žena i muškaraca. Iz ovih rezultata možemo vidjeti da se rezultati ne razlikuju puno, no treba provjeriti je li ova mala razlika statistički značajna.
Kako bi mogli provesti t-test, moramo najprije provjeriti pretpostavke normalnosti i nezavisnosti uzorka. Obzirom na to da razmatramo dva uzoraka različitih spolova, možemo pretpostaviti njihovu nezavisnost. Sljedeći korak je provjeriti normalnost podataka koju najčešće provjeravamo: histogramom te KS-testom (kojim provjeravamo pripadnost podataka distribuciji).

```{r}
hist(women$wealth.worth.in.billions,
     main='Histogram of womens wealth worth in billions',
     xlab='Wealth',
     col="pink")

hist(men$wealth.worth.in.billions,
     main='Histogram of mens wealth worth in billions',
     xlab='Wealth',
     col="pink")



```


Iz histograma jasno vidimo da se podaci ne ravnaju po normalnoj distribuciji. Ipak, uočeno je potrebno dodatno ispitati koristeći Kolmogorov–Smirnov test 
kojim se utvrđuje ravna li se distribucija po normalnoj razdiobi.

```{r, warning=FALSE}
ks.test(women$wealth.worth.in.billions, y="pnorm")
ks.test(men$wealth.worth.in.billions, y="pnorm")
```


Iz dobivenih p vrijednosti u oba slučaja odbacujemo mogućnost da se distribucije
ravnaju po normalnoj razdiobi. 



Potrebno je koristiti neparametarski test Mann–Whitney U test, koji se koristi 
kada se podaci se ravnaju po istim distribucijama i uzorci su nezavisni iz jedne i druge populacije (u našem datasetu osoba ima dodijeljen spol male ili female).


Hipoteze glase:
$$ \begin{aligned}
H_0&: \mu_1 = \mu_2 \\
H_1&: \mu_1 > \mu_2 \quad \quad 
\end{aligned} $$


```{r}
 
wilcox.test(men_mean, women_mean, alt = "greater")

```


S obzirom na to da je dobivena p-vrijednost > 0.05, nemamo temelja za odbaciti nultu hipotezu, odnosno muškarci i žene se statistički značajno ne razlikuju u količini bogatstva.