##### Zadatak 1.
### 1.b)
model <- lm(Employed ~ ., longley)
summary(model)
# Pravimo interval poverenja za b_pop - b_gnp kao na petom casu
l <- c(0, 0, -1, 0, 0, 1, 0)
X <- model.matrix(model)
XtXi <- solve(t(X)%*%X)
n <- nrow(X)
p <- length(model$coefficients) - 1
estim <- t(l) %*% model$coefficients
stdev <- sqrt(sum(model$res^2)/(n-p-1) * t(l) %*% XtXi %*% l)
lwr <- estim - stdev * qt(0.975, n-p-1)
upr <- estim + stdev * qt(0.975, n-p-1)
c(lwr, upr)

# Ovaj interval ukljucuje nulu, sto znaci da ne mozemo da odbacimo hipotezu da
# je b_pop = b_gnp
# Drugi nacin za testiranje H0: b_pop = b_gnp je koriscenjem linear_hypothesis
linear_hypothesis(model, t(c(0, 0, -1, 0, 0, 1, 0)), 0)
# ili anova...
anova(lm(Employed ~ . -GNP - Population + I(GNP + Population), longley), model)

### 1.g) Multikolinearnost
# Mozemo da gledamo vif
library(car)
vif(model)
# Mnoge vrednosti su velike, dakle postoji prisustvo multikolinearnosti

# Promenljive Unemployed i Armed.Forces imaju vidno manji vif od ostalih, tako
# da cemo ih sacuvati.
# Pogledajmo korelacije ostalih promenljivih
cor(longley[, -c(3, 4, 7)])
# Primecujemo izuzetno velike korelacije, mozemo i nacrtati grafike
plot(longley[, -c(3, 4, 7)])
# Velika linearna zavisnost je prisutna, tako da mozemo ih sve zameniti samo
# jednom od promenljivih, npr. Year.
model2 <- lm(Employed ~ Armed.Forces + Unemployed + Year, longley)
summary(model2)
vif(model2)
# U ovom modelu su vif-ovi svi mali (<30) pa nemamo vise problem
# multikolinearnosti, dok je koeficijent determinacije smanjen sa 0.9955 na
# 0.9928, sto je vrlo mala razlika, pa je dobar smanjeni model.

### 1.zh) Heteroskedasticnost
model <- lm(dist ~ speed + I(speed^2), cars)
# Ispitujemo prisustvo heteroskedasticnosti npr. grafikom zavisnosti apsolutnih
# vrednosti standardizovanih reziduala od ocenjenih vrednosti dist.
plot(abs(rstandard(model)) ~ fitted(model))
# Vidi se da disperzija raste sa porastom dist, tako da imamo problem
# heteroskedasticnosti.

# Trazimo odgovarajucu box-cox transformaciju
library(MASS)
bc <- boxcox(model)
lambda <- bc$x[which.max(bc$y)] # = 0.38
# Pravimo funkciju za transformaciju
BC <- function(y) {
  (y^lambda - 1)/lambda
}
# Fitujemo model sa boxcox transformisanom zavisnom promenljivom
model2 <- lm(BC(dist) ~ speed + I(speed^2), cars)
summary(model2)
plot(abs(rstandard(model2)) ~ fitted(model2))
# Sada je heteroskedasticnost vidno ublazena.

##### Zadatak 2.
library(LSMhelp) # instalacija: devtools::install_github("blaza/lsmhelp")
# Pravimo pocetni model sa svim prediktorima
m <- glm(REMISS ~ ., leukemia, family=binomial)
# U summary vidimo vrednosti Valdovog testa (zvezdice)
summary(m)
# Svi koeficijenti su beznacajni

# Testom kolicnika verodostojnosti proveravamo da li postoji uticaj prediktora
# na zavisnu. Poredimo model koji ukljucuje samo slobodan clan i nas model.
anova(glm(REMISS ~ 1, leukemia, family=binomial), m, test = "Chisq")
# p vrednost < 0.05, pa zakljucujemo da postoji znacajan uticaj

# Poredimo model samo sa LI sa nasim modelom
anova(glm(REMISS ~ LI, leukemia, family=binomial), m, test = "Chisq")
# p vrednost je velika, znaci ne postoji znacajna razlika izmedju dva modela

summary(glm(REMISS ~ LI, leukemia, family=binomial))
# AIC za ovaj model je manji nego za prvobitni, dakle bolji je novi model.