# Kategoricki prediktori

# Prvo cemo se baviti binarnim prediktorima, tj. prediktorima koji uzimaju samo
# 2 razlicite vrednosti. Najcesce se u modelu predstavljaju kao indikatori,
# odnosno promenljive koje su 0 ili 1.

# Pogledajmo primer podataka koji pokazuju efekat upotrebe cigareta u trudnoci
# na tezinu novorodjenceta.
# Podatke preuzeti sa:
# http://www.matf.bg.ac.rs/p/files/94-birthsmokers.txt
birthsmokers <- read.table("birthsmokers.txt", header = TRUE)
# Promenljiva Smoke je ucitana kao factor, sto u R-u oznacava kategorcku
# promenljivu, koja uzima konacan broj vrednosti koje nisu uredjene.
birthsmokers$Smoke

# Nacrtajmo podatke
plot(birthsmokers)
# Wgt (tezina) ima linearnu vezu sa Gest (nedelja trudnoce) i to cemo
# modelovati, ali cemo posmatrati i uticaj pusenja na tezinu.

# Pravimo model sa oba prediktora
model <- lm(Wgt ~ Gest + Smoke, birthsmokers)
summary(model)
# Bitno je tumacenje koeficijenata modela. Koeficijent koji je vezan za Smoke je
# oznacen kao "Smokeyes", jer je faktorska promenljiva i R automatski kreira
# indikatorsku promenljivu za nivo "yes" koja je 1 ukoliko je Smoke 'yes' i 0
# ako nije. Za nivo "no" ne pravi posebnu promenljivu jer je jednoznacno
# odredjena sa indikatorom za "yes".
# Kako je kodirana promenljiva Smoke u modelu mozemo da vidimo iz model matrice
model.matrix(model)
# ili funkcijom contrasts
contrasts(birthsmokers$Smoke)
# Interpretacija ovog koeficijenta je: Ako je imamo dva podatka za koja je
# vrednost Gest ista, dok je kod prvog Smoke "yes" a drugog "no", za prvi ce
# vrednost zavisne promenljive Wgt biti za 244.544 niza nego za drugi.
# U sustini, ovime pravimo dve regresione prave, jednu za Smoke "yes" i jednu za
# "no", koje dele koeficijent pravca (uz Gest), tj. paralelne su.
# Ilustrujmo na slici:
# Crtamo Wgt u odnosu na Gest i bojimo tacke u zavisnosti od vrednosti Smoke
plot(Wgt ~ Gest, birthsmokers,
     col = ifelse(birthsmokers$Smoke == "yes", "red", "blue"))
# Prava koja bi bila bez razmatranja Smoke je
#abline(lm(Wgt ~ Gest, birthsmokers))
# Prava za Smoke == "no":
abline(model$coef[1], model$coef[2], col = "blue")
# Prava za Smoke == "yes":
abline(model$coef[1] + model$coef[3], model$coef[2], col = "red")
# Prave su paralelne, a Smoke=="yes" prava je ispod Smoke=="no" prave. Ovo znaci
# da je u proseku tezina bebe od majke pusaca 244.5 grama manja od tezine bebe
# majke nepusaca.

# Mozemo testirati da li je razlika izmedju pusaca i nepusaca znacajna
# testiranjem hipoteze H0: b_{smokeyes} = 0, sto vidimo da je znacajno, jer ima
# 3 zvezdice koeficijent uz smokeyes. Mogli smo i da proverimo interval
# poverenja za koeficijente
confint(model)
# Interval poverenja za Smokeyes ne sadrzi nulu, pa odbacujemo H0.

# Mozemo odrediti i interval poverenja za prosecnu vrednost tezine bebe majki
# koje su rodile u 38 nedelji trudnoce, a koje su pusaci ili nepusaci
ci <- predict(model,
              data.frame(Gest=c(38,38), Smoke=c("yes", "no")),
              interval = "confidence")
# Vrednost ocene (fit) se razlikuje upravo za 244.544, dok su sirine intervala
# razlicite (zasto?) pa se granice ne razlikuju za istu vrednost
ci[1, ] - ci[2, ]
ci[ , 2] - ci[, 3] # sirine intervala

# Mogli smo da promenimo referentni nivo faktora da ne bude Somke=="no", tj. da
# kodiramo Smoke=="yes" sa 0, a Smoke=="no" sa 1.
birthsmokers$Smoke2 <- relevel(birthsmokers$Smoke, "yes")
# Kodiranje je sad:
contrasts(birthsmokers$Smoke2)
# Kad napravimo model sa referentnim nivoom "yes", dobijamo
model2 <- lm(Wgt ~ Gest + Smoke2, birthsmokers)
summary(model2)
# Sada pise Smoke2no i koeficijent je 244.544, sto je suprotno od Smokeyes u
# prvobitnom modelu, sto je ocekivano.

# Zasto koristiti kategoricku promenljivu umesto dva razlicita modela, kad vec
# dobijamo dve prave?
# Podelimo podatke na pusace i nepusace i napravimo dva modela Wgt ~ Gest
model_smoke <- lm(Wgt ~ Gest, birthsmokers[birthsmokers$Smoke == "yes", ])
model_nosmoke <- lm(Wgt ~ Gest, birthsmokers[birthsmokers$Smoke == "no", ])
summary(model_smoke)
summary(model_nosmoke)
# Nacrtamo li ih...
abline(model_smoke)
abline(model_nosmoke)
# Dobijamo jako slicne prave kao i sa kategorickim, s tim sto su koeficijenti
# pravca za nijansu drugaciji.

# Standardne greske uz Gest za model, model_smoke, model_nosmoke su, redom,
# 9.12, 14.11, 11.97. Dakle, za prvi model je standardna greska najmanja, sto
# utice i na intervale poverenja. Pogledajmo i intervale poverenja za prosecnu
# vrednost Wgt kada je Gest=38...
predict(model,
        data.frame(Gest=c(38,38), Smoke=c("yes", "no")),
        interval = "confidence")
predict(model_smoke, data.frame(Gest=38), interval = "confidence")
predict(model_nosmoke, data.frame(Gest=38), interval = "confidence")
# Interval poverenja za nepusace je skoro isti, dok je za pusace interval
# poverenja uzi kod modela koji obuhvata sve podatke i indikatorsku promenljivu.


# Interakcije medju prediktorima

# Cesto se desava da prediktori nemaju samo aditivni efekat na zavisnu
# promenljivu, vec postoje i medjusobni uticaji medju njima. Na primer, u
# zavisnosti od neke kategoricke promenljive, koeficijent uz neki prediktor bi
# mozda trebalo promeniti. Taj efekat modelujemo dodavanjem interakcija u model,
# koje se svode na dodavanje novog prediktora koji predstavlja proizvod
# postojecih prediktora.

# Videli smo da kada napravimo dva razlicita modela za pusace i nepusace,
# dobijemo dve prave ciji se koeficijenti pravca malo razlikuju. Nas originalni
# model nam je dozvolio da samo promenimo slobodan clan u zavisnosti od
# promenljive Smoke, cime smo dobili samo dve paralelne prave. Uvodjenjem
# interakcija mozemo da dobijemo i razlicite koeficijente pravca, cime dobijamo
# iste prave kao i kod pojedinacnih modela.
# u R-u se interakcija uvodi znakom mnozenja *
model_int <- lm(Wgt ~ Gest * Smoke, birthsmokers)
summary(model_int)
# sa Gest * Smoke dobili smo model Gest + Smokeyes + Gest*Smokeyes
# Da bismo proverili to, dodacemo promenljivu Smokeyes u podatke
birthsmokers$Smokeyes <- birthsmokers$Smoke == "yes"
summary(lm(Wgt ~ Gest + Smokeyes + I(Gest * Smokeyes), birthsmokers))
# model je isti kao i Gest*Smoke

# Kada nacrtamo prave za odgovarajuce faktore
plot(Wgt ~ Gest, birthsmokers,
     col = ifelse(birthsmokers$Smoke == "yes", "red", "blue"))
# Prava za Smoke == "no":
abline(model_int$coef[1], model_int$coef[2], col = "blue")
# Prava za Smoke == "yes":
abline(model_int$coef[1] + model_int$coef[3],
       model_int$coef[2] + model_int$coef[4],
       col = "red")
# Iste su kao pojedinacni modeli
abline(model_smoke)
abline(model_nosmoke)


# Pokazacemo jos jedan primer, sada sa vise od 2 vrednosti kategoricke
# promenljive. Gledacemo efikasnost 3 razlicite terapije za lecenje depresije.
# Podaci: http://www.matf.bg.ac.rs/p/files/94-depression.txt
depression <- read.table("depression.txt", header = TRUE)
# TRT je kategoricka promenljiva i ima 3 vrednosti - A,B,C. To se moze kodirati
# na mnogo nacina, ali standardan nacin na koji se to radi je da se naprave 2
# indikatorske promenljive (za n kategorija koristi se n-1 indikatora), odredi
# se jedan referentni nivo, koji ce biti predstavljen samo nulama, a ostale
# kategorije su vektori koji imaju 1 na jednom mestu. U ovom primeru, kodiranje
# koje je izabrano je da C bude referentni nivo, a A i B su kodirani kao (1,0) i
# (0,1), respektivno. kodiranje vidimo po promenljivim x2 i x3
depression
# Mi koristimo factor promenljivu TRT u modelu, pa cemo staviti C za referentni
# nivo da bismo dobili isto kodiranje kao u x2, x3. To radimo funkcijom relevel.
depression$TRT <- relevel(depression$TRT, "C")

# Crtamo zavisnost y od age
plot(y ~ age, data = depression, col=TRT)
legend("topleft", legend = c("TRT C", "TRT A", "TRT B"),
       col = 1:3, pch=1)

# Pravimo model koji uzima u obzir sve interakcije izmedju age i TRT
dep_model <- lm(y ~ age * TRT, depression)
summary(dep_model)
# Isti model dobijamo i na sledeci nacin, koristeci kodiranje x2, x3
summary(lm(y ~ age + x2 + x3 + I(age*x2) + I(age*x3), depression))
# Jos jedan nacin na koji smo mogli da postignemo ovo je koriscenje ':' kao znak
# za pojedinacnu interakciju (* uzima u obzir sve moguce interakcije). Znak : se
# moze tumaciti kao mnozenje u vecini slucajeva.
summary(lm(y ~ age + x2 + x3 + age:x2 + age:x3, depression))

# Nacrtacemo sve moguce prave dobijene ovim modelom.
coefs <- dep_model$coefficients
abline(coefs[1], coefs[2], col=1)
abline(coefs[1] + coefs[3], coefs[2] + coefs[5], col=2)
abline(coefs[1] + coefs[4], coefs[2] + coefs[6], col=3)

# Razmotriti interpretacije koeficijenata!

# Nacrtajmo i reziduale
plot(dep_model$res ~ dep_model$fitted)
abline(h=0)
# izgledaju dobro

# Normalnost...
qqnorm(dep_model$res)
qqline(dep_model$res)
shapiro.test(dep_model$residuals)
# izgleda dobro.

# Mozemo odgovoriti na nekoliko pitanja.

# Prvo, da li, za fiksni broj godina, postoji razlika u efikasnosti razlicitih
# vidova terapije (A,B i C)?

# Ovo znaci da se 3 regresione prave, za svaku od terapija A,b i C, poklapaju.
# To ce se desiti samo ako su slobodni clanovi i koeficijenti pravca jednaki.
# To je ispunjeno samo ako vazi da su svi koeficijenti sem age i intercept nule
# (proveriti zasto je to tako!)
# Dakle testiramo hipotezu H0: b_TRTA = b_TRTB = b_age:TRTA = b_age:TRTB = 0
# Mozemo koristiti linear_hypothesis funkciju sa ranijih casova, ili anova f-ju
anova(lm(y ~ age, depression), dep_model)
# p-vrednost testa je mala, pa odbacujemo H0, tj. postoji razlika u efikasnosti
linear_hypothesis(dep_model, matrix(c(0,0,1,0,0,0,
                                      0,0,0,1,0,0,
                                      0,0,0,0,1,0,
                                      0,0,0,0,0,1),
                                    byrow = TRUE, nrow = 4),
                  c(0,0,0,0)) # isti rezultat

# Drugo pitanje: Da li se uticaj broja godina na efikasnost terapije menja u
# zavisnosti od primenjene terapije?

# Ovo je ekvivalentno testiranju H0: b_age:TRTA = b_age:TRTB = 0
anova(lm(y ~ age + TRT, depression), dep_model)
# Opet je odgovor potvrdan - uticaj se menja od terapije do terapije
linear_hypothesis(dep_model, matrix(c(0,0,0,0,1,0,
                                      0,0,0,0,0,1),
                                    byrow = TRUE, nrow = 2),
                  c(0,0)) # isti rezultat


# Sta moze da se desi ako zanemarimo interakcije?

# Generisemo vestacke podatke tako da imamo dve grupe podataka, sa drasticno
# razlicitim pravcima - jedni leze na y=x a drugi na y=10-x.
set.seed(1)
x <- runif(70, min = 0, max = 10)
group <- sample(c(0, 1), 70, replace = TRUE)
y <- (x + rnorm(70)) * group + (10-x + rnorm(70)) * (1 - group)
data <- data.frame(y, x, group)
plot(x, y, col = factor(group))

# Prvo cemo napraviti model bez interakcija
model <- lm(y ~ x + group)
summary(model)
# Svi koeficijenti su neznacajni! Cak i F test to povrdjuje.
# Kako izgledaju prave za grupe 0 i 1?
abline(model$coef[1], model$coef[2], col=1)
abline(model$coef[1] + model$coef[3], model$coef[2], col=2)
legend("top", col=1:2, lty=1, legend=c("group 0", "group 1"))
# Obe su vrlo bliske horizontalnoj liniji i imamo jako los model.

# Pogledajmo i reziduale
plot(model$res ~ model$fitted)
abline(h=0)
# Nisu ni simetricni oko nule, ni homoskedasticni
# Dakle, imamo vrlo slab model.
# Mozda bi bilo pametno dodati i interakciju izmedju x i group
# Kada dodamo interakciju...
model_int <- lm(y ~ x + group + x:group)
summary(model_int)
# Svi koeficijenti postaju znacajni, sa p-vrednocu 0
# Pogledajmo sada dobijene prave za odgovarajuce grupe
plot(x, y, col = factor(group))
# (razmisliti zasto su koeficijenti bas ovakvi)
abline(model_int$coef[1], model_int$coef[2], col = 1)
abline(model_int$coef[1] + model_int$coef[3],
       model_int$coef[2] + model_int$coef[4],
       col = 2)
# Sada model vrlo dobro opisuje podatke.
# Kako izgledaju reziduali?
plot(model_int$res ~ model_int$fitted)
abline(h=0)
# Mnogo bolje, disperzija im je skoro ista, i centrirani su oko nule.