#TESTIRANJE HIPOTEZA

#Poredimo modele m i M
#testiranje: H0(M nije znacajno bolji od m), H1(da je M bolji model)
#test statistika: F=[(RSS_m-RSS_M)/(df_m-df_M)]/[RSS_M/(n-df_M)]~F_(p-q),(n-p) 

#Slucaj 1. M sadrzi sve prediktore, m je model bez X-eva (konstantan model)
#nulta hipoteza: beta1=beta2=...=beta_(p-1)=0
#[(TSS - RSS)/(p-1)/[RSS/(n-p)] ~ F_(p-1,n-p)
#kriticna oblast: W={F>C}

library(faraway)
data(savings)
head(savings)
g<-lm(sr~pop15+pop75+dpi+ddpi, savings) #g<-lm(sr ~ .,savings)

#direktno racunanje preko formula
(tss<-sum((savings$sr-mean(savings$sr))^2)) #spoljne zagrade znace da se ispise rezultat
(rss<-deviance(g))
df.residual(g) #broj stepeni slobode
(fstat<-((tss-rss)/4)/(rss/df.residual(g))) #test statistika     
1-pf(fstat,4,df.residual(g)) #p vrednost testa

summary(g) #direktno dobijanje vrednosti test statistike

#mala p vrednost testa => odbacuje se nulta hipoteza => model iz H0 je los


#Slucaj 2. M sadrzi sve prediktore, m sadrzi sve prediktore sem pop15
g2<-lm(sr~pop75+dpi+ddpi, savings)

#preko F testa
(rss2<-deviance(g2))
(fstat<-(deviance(g2)-deviance(g))/(deviance(g)/df.residual(g)))
1-pf(fstat,1,df.residual(g)) #p-vrednost F testa    
sqrt(fstat) #test statistika t-testa
#preko t testa
(A<-summary(g)$coef)
(tstat<-A[2,1]/A[2,2]) #isto sto i t.st<-A[2,3]
(p.vr<-2*(1-pt(abs(tstat),df.residual(g))))

#poredjenje modela preko ANOVA tabele
anova(g2,g)                         

summary(g)$coefficients #vidi se isto sto i u prethodnom 

#mala p vrednost testa => odbacuje se nulta hipoteza => pop15 je znacajan za model

#Slucaj 3. M sadrzi sve prediktore, m sadrzi sve sem pop15 i ddpi 

g3 <- lm(sr ~ pop15+dpi , savings)
anova(g3,g)

#p vrednost testa vrlo bliska 0.05 pa nije lako odluciti da li da se odbaci hipoteza

#Slucaj 4. M sadrzi sve prdiktore, m sadrzi prediktore pop15 i pop75 kao jedan prediktor 

gr<-lm(sr~I(pop15+pop75)+dpi+ddpi,savings) #I znaci da ih posmatramo zajedno
anova(gr,g)
# p vrednost testa velika => prihvata se nulta hipoteza => moze se koristiti transformisan model

#Slucaj 5. testiramo: koeficijent uz ddpi je 0.5
gr<-lm(sr~pop15+pop75+dpi+offset(0.5*ddpi),savings) #offset znaci da se koeficijent ne ocenju nego se uzima navedena vrednost
anova(gr,g)
# p vrednost testa velika => prihvata se nulta hipoteza => moze se koristiti transformisan model

#preko t testa: beta_ddpi=0.5 protiv beta_ddpi!=0.5
(tstat<-(0.409695-0.5)/0.196197) # (ocenjena_vrednost-pretpostavljena_vrednost)/std_devijacija iz modela g
2*pt(tstat,45)
tstat^2 # Fiserova statistika

# test permutacija

g<-lm(sr~pop75+dpi,savings)
summary(g) #citamo F-statistic  
gs<-summary(g)
gs$fstat #dobijanje F statistike i parametara    

# sample generise slucajne permutacije
# pravimo uzorak od 4000 slucajno odabrani permutacija (ukupno ima 50! permutacija)
fstats<-numeric(4000) #vektor u koji se smestaju permutacije
for(i in 1:4000){
  ge<-lm(sample(sr)~pop75+dpi,data=savings)
  fstats[i]<-summary(ge)$fstat[1]
}
(p<-length(fstats[fstats>2.6796])/4000) # 2.6796 je vrednost F-statistike za model g

# p je priblizno jednako p-vrednosti testa iz modela g => model je koristan

# permutovanje prediktora
# H0(nema razlike ako permutujemo pop75), H1(razlikuju se modeli)
summary(g)$coef[2,3]    # t-statistika iz modela
tstats<-numeric(4000)
for(i in 1:4000){
  ge<-lm(sr~sample(pop75)+dpi, savings)     
  tstats[i]<-summary(ge)$coef[2,3]
}
mean(abs(tstats)>1.6783) #1.6783 je vrednost iz modela g

# p vrednost velika => prihvatamo H0 => Y ne zavisi od pop75

# intervali i regioni poverenja 

g<-lm(sr~ ., savings) #ceo model
# pravljenje intervala poverenja preko formula
t<-qt(0.975,45)
c(-1.69-t*1.08,-1.69+t*1.08) #interval poverenja za pop75   
c(0.41-t*0.196,0.41+t*0.196) #interval poverenja za ddpi
confint(g)     #svi 95%-ni intervali poverenja (za sve bete), ako se doda level menja se nivo poverenja

install.packages("ellipse") #paket za crtanje elipsi poverenja
library(ellipse)
# elipsa poverenja za pop15 i pop75
plot(ellipse(g,c(2,3)),type="l",xlim=c(-1,0))
points(0,0) #dodavanje tačke koja odgovara testirnju
#tacka van elipse => odbacujemo H0 da su parametri jednaki 0
points(coef(g)[2],coef(g)[3],pch=18) #tackasta ocena parametra beta (u centru elipse zato sto je elipsa region poverenja)
#dodavanje pojedinacnih intervala poverenja
abline(v=confint(g)[2,],lty=2)
abline(h=confint(g)[3,],lty=2)

cor(savings$pop15,savings$pop75) #korelacija izmedju provemljivih
summary(g,corr=TRUE)$corr #korelacija koeficijenata
#corr nije ukljucen po difoltu u summary pa moramo da stavimo corr=TRUE

# interval predvidjanja (poverenja zavisne promenljive)
g<-lm(Species~Area+Elevation+Nearest+Scruz+Adjacent,gala)
x0<-c(1, 0.08, 93, 6.0, 12.0, 0.34) #konkretne vrednosti prediktora
(y0 <- sum(x0*coef(g))) #predvidjamo broj vrsta na ostvru sa navedenim karakteristikama
# formiranje intervala za srednju vrednost
summary(g)
(t<-qt(0.975,24)) #vrednost t-testa
x<-model.matrix(g) 
xtxi<-solve(t(x)%*%x) 
(bm<-sqrt(x0 %*% xtxi %*% x0)*t*summary(g)$sigma)  
c(y0-bm,y0+bm) #interval
# formiranje intervala za pojedinacne vrednosti
bm<-sqrt(1+x0 %*% xtxi %*% x0)*t*summary(g)$sigma
c(y0-bm,y0+bm) #leva granica je nema smisla jer brojimo vrste, pa uzimamo da je 0

# isto mozemo uraditi pomocu ugradjenih funkcija
x0<-data.frame(Area=0.08,Elevation=93,Nearest=6.0,Scruz=12,Adjacent=0.34)
str(predict(g,x0,se=TRUE))
predict(g,x0,interval="confidence") #interval poverenja za srednju vrednost
predict(g,x0,interval="prediction") #interval predvidjanja 

# ekstrapolacija
grid<-seq(0,100,1) #niz vrednosti od 1 do 100 sa korakom 1
p<-predict(g,data.frame(Area=0.08,Elevation=93,Nearest=grid,Scruz=12,Adjacent=0.34),se=T,interval="confidence")
matplot(grid,p$fit,lty=c(1,2,2),type="l",xlab="Nearest",ylab="Species")
rug(gala$Nearest)   #gde se nalaze pocetni X-evi, iz baze na osnovu kojih je napravljen model