# ANOVA

# Pogledacemo primere analize podataka kada su sve promenljive kategoricke, tj.
# kada je cilj odrediti postojanje razlike u prosecnoj vrednosti zavisne
# promenljive izmedju vise grupa i uticaj raznih kategorickih prediktora na
# prosecnu vrednost zavisne promenljive. Ovakva analiza se cesto naziva ANOVA
# (Analysis of variance). Glavna razlika u odnosu na sta smo ranije radili je
# sto su prediktori sada samo kategoricki i nemamo neprekidne, pa, graficki,
# nemamo prave koje provlacimo kroz podatke, vec gledamo nivoe prosecnih
# vrednosti po grupama.

# Koristicemo Salaries skup podataka iz carData paketa.
salaries <- carData::Salaries
# Ovaj skup podataka sadrzi informacije o visini plate zaposlenih na fakultetu
# (ne matfu) a prediktori na koje cemo se fokusirati su zvanje i disciplina.
?carData::Salaries

# Kao osnovni grafik podataka cemo koristiti boxplot. Za tumacenje, npr. videti
# https://stats.stackexchange.com/questions/189079/how-to-read-a-boxplot-in-r
boxplot(salary ~ rank, salaries)
boxplot(salary ~ discipline, salaries)


# Pogledacemo prvo uticaj zvanja profesora na platu
rank_model <- lm(salary ~ rank, salaries)
summary(rank_model)
# Koeficijenti su znacajni, pa zakljucujemo da postoje znacajne razlike u
# prosecnoj plati izmedju Prof i AssocProf, kao i izmedju AssocProf i AsstProf i
# AsstProf i Prof. F test takodje ima malu p vrednost sto nam kaze da je uticaj
# ove grupe prediktora, tj. ovog kategorickog prediktora, znacajan.

# Bitno je i ovde znati interpretaciju koeficijenata.
# 1. Intercept predstavlja prosenu vrednost plate za referentni nivo (AsstProf)
mean(salaries$salary[salaries$rank == "AsstProf"])
rank_model$coef[1] # jednaki su
# 2. Koeficijent uz prvi prediktor oznacava razliku izmedju proseka za AssocProf
# i AsstProf (ref. nivo), tj. koliko treba dodati proseku za AsstProf da bi se
# dobio prosek za AssocProf
mean(salaries$salary[salaries$rank == "AssocProf"])
mean(salaries$salary[salaries$rank == "AsstProf"]) + rank_model$coef[2]
# 3. Slicno, koeficijent za prediktor rankProf oznacava razliku izmedju proseka
# za Prof i AsstProf (ref. nivo).
mean(salaries$salary[salaries$rank == "Prof"])
mean(salaries$salary[salaries$rank == "AsstProf"]) + rank_model$coef[3]

# Ukljucimo sada i prediktor koji oznacava disciplinu (primenjena vs teorijska)
dis_rank_model <- lm(salary ~ discipline + rank, salaries)
summary(dis_rank_model)
# U ovakvom modelu, ocekivano bi bilo da se prosek plate za rank R i discipline
# D dobija zbirom intercept koeficijenta i koeficijenata uz odgovarajuce
# prediktore za rank R i disciplinu D. Npr. za AssocProf i disciplinu B,
# ocekivali bismo da je prosek plata za primenjene matematicare koji su
# AssocProf jednaka intercept + disciplineB + rankAssocProf. Slicno, prosek
# plata za teorijske matematicare (discipline = A) koji su AsstProf bi trebalo
# da bude jednaka samo koeficijentu uz intercept.
# **Medjutim, to se ne desava!**
mean(salaries$salary[salaries$rank == "AsstProf" & salaries$discipline == "A"])
# Vidimo da je stvarni prosek u drugom slucaju 73935.54
dis_rank_model$coefficients[1]
# dok je koeficijent uz intercept 71944.33 =/= 73935.54
# Ista stvar se desava i za prvi pomenuti slucaj
mean(salaries$salary[salaries$rank == "AssocProf" & salaries$discipline == "B"])
# prosecna plata za primenjenog AssocProf je 101276.4, dok je zbir odgovarajucih
# koeficijenata jednak 99466.83.
dis_rank_model$coef["(Intercept)"] +
  dis_rank_model$coef["disciplineB"] +
  dis_rank_model$coef["rankAssocProf"]
# sum(dis_rank_model$coef[c("(Intercept)", "disciplineB", "rankAssocProf")])

# Razlog za ovo je nedovoljna fleksibilnost ovakvog modela (to ne mora nuzno
# biti lose), jer imamo 4 koeficijenta, a 2*3=6 mogucih kombinacija vrednosti
# kategorickih prediktora, pa nemamo dovoljno koeficijenata da tacno ocenimo
# proseke, vec dajemo priblizne ocene.

# Fleksibilniji model, kod koga se koeficijenti poklapaju sa prosecima, se
# dobija uvodjenjem interakcija izmedju prediktora rank i discipline
dis_rank_int_model <- lm(salary ~ discipline * rank, salaries)
summary(dis_rank_int_model)
# Kako se tumace koeficijenti?
# 1. Koeficijent uz intercept je prosek plate za profesore ciji su i zvanje i
# disciplina na referentnim nivoima, tj. kod nas, prosek plate profesora koji je
# AsstProf i sa disciplinom A (teorijska matematika).
mean(salaries$salary[salaries$rank == "AsstProf" & salaries$discipline == "A"])
dis_rank_int_model$coef[1]
# 2. Posto imamo interakcije, ostale koeficijente mozemo rastaviti na 2 slucaja,
# jedan za referentni nivo discipline (A) i drugi za vrednost discipline "B", tj
# primenjene matematicare. Drugim recimo odvajamo model za teorijske (A) i
# primenjene(B) matematicare.

# 2.1. koeficijenti rank* oznacavaju koliko treba dodati na referentni nivo da
# bi se dobio prosek plate za profesora teorijske matematike odgovarajuceg
# zvanja. Na primer, za teoretskog matematicara koji je AssocProf imamo
mean(salaries$salary[salaries$rank == "AssocProf" & salaries$discipline == "A"])
dis_rank_int_model$coef[1] + dis_rank_int_model$coef["rankAssocProf"]
# 2.2. Kada je disciplineB == 1, tj. kada gledamo primenjene matematicare, svi
# koeficijenti su nam potrebni.
# Koeficijent disciplineB nam kaze koliko treba dodati na intercept da bi se
# dobio prosek plate za AsstProf (ref. nivo zvanja) medju primenjenim
# matematicarima.
mean(salaries$salary[salaries$rank == "AsstProf" & salaries$discipline == "B"])
dis_rank_int_model$coef[1] + dis_rank_int_model$coef["disciplineB"]
# Koeficijenti rank* imaju slicnu svrhu kao i pre, s tim sto se sada sabiraju sa
# koeficijentima disciplineB:rank* (koef. interakcije) da bi se dobilo koliko
# treba dodati na prosek plate za primenjene matematicare da bi rezultat bio
# prosek plate za odgovarajuce zvanje, medju primenjenim matematicarima.
# Na primer, za Prof sa primenjene matematike:
mean(salaries$salary[salaries$rank == "Prof" & salaries$discipline == "B"])
dis_rank_int_model$coef[1] +
  dis_rank_int_model$coef["disciplineB"] +
  dis_rank_int_model$coef["rankProf"] +
  dis_rank_int_model$coef["disciplineB:rankProf"]

# Mozemo odgovoriti na nekoliko pitanja...

# 1. Da li postoji znacajan uticaj discipline na prosecnu platu profesora? Da
# odgovorimo na pitanje, testiramo hipotezu da su koeficijenti koji sadrze
# disciplineB nule, tj.
# H0: b_disciplineB = b_disciplineB:rankAssocProf = b_disciplineB:rankProf = 0
# Ovime gledamo zapravo da li mozemo da izbacimo prediktor discipline iz modela
# ovako...
anova(lm(salary ~ rank, salaries), dis_rank_int_model)
# ili ovako...
linear_hypothesis(dis_rank_int_model, matrix(c(0,1,0,0,0,0,
                                      0,0,0,0,1,0,
                                      0,0,0,0,0,1),
                                    byrow = TRUE, nrow = 3),
                  c(0,0,0)) # isti rezultat
# P vrednosti su jako male, pa zakljucujemo da je znacajan uticaj discipline na
# prosecnu platu. Ovaj rezultat je interesantan, jer su pojedinacni koeficijenti
# u modelu bili neznacajni, ali vidimo da ih ne mozemo izbaciti.

# 2. Da li se razlikuje uticaj zvanja na platu medju razlicitim disciplinama?
# Ovo proveravamo testiranjem da su koeficijenti interakcija nule, tj.
# H0: b_disciplineB:rankAssocProf = b_disciplineB:rankProf = 0
# ovako...
anova(lm(salary ~ rank + discipline, salaries), dis_rank_int_model)
# ili ovako...
linear_hypothesis(dis_rank_int_model, matrix(c(0,0,0,0,1,0,
                                               0,0,0,0,0,1),
                                             byrow = TRUE, nrow = 2),
                  c(0,0)) # isti rezultat
# Sad je p vrednost velika, pa ne mozemo odbaciti hipotezu da je razlika u
# uticajima zvanja na platu oizmedju razlicitih disciplina znacajna. Drugim
# recima, prosecna plata priblizno jednako raste od zvanja do zvanja kod
# teoretskih i primenjenih matematicara.

# Za aktiviste: Analizirati uticaj pola na platu :)