#Problemi sa greskama

#1. Uopsteni najmanji kvadrati (GLS)
#ako disperzija reziduala nije konstantna ili su korelisani
#pritom ta disperzija najcesce nije data nego je ocenjujemo iz uzorka
library(faraway)
data(longley)
g<-lm(Employed~GNP+Population, longley)
summary(g,cor=T)
cor(longley$GNP,longley$Pop)
#korelacija izmedju koeficijenata GNP i Population je jaka negativna, a izmedju tih promenljivih jaka pozitivna
#susedne greske bi mogle biti korelisane
#autoregresioni reziduali - reziduali nisu nezavisni i autoregresiono su korelisani
cor(residuals(g)[-1],residuals(g)[-16])  # bez prvog i bez poslednjeg
#vidimo da korelacija nije zanemarljiva => odlucujemo se za GLS metod
#primena formula pesacki (ako znamo da je korelacija 0.31041 i da AR model)
x<-model.matrix(g)
Sigma<-diag(16)    #jedinicna matrica
Sigma<-0.31041^abs(row(Sigma)-col(Sigma)) #sigma_{ij}=ro^|i-j|
Sigi<-solve(Sigma)
xtxi<-solve(t(x) %*% Sigi %*% x)
(beta<-solve(t(x) %*% Sigi %*% x,t(x) %*% Sigi %*% longley$Empl))
res<-longley$Empl - x %*% beta
(sig<-sqrt((t(res) %*% Sigi %*% res)/g$df))
sqrt(diag(xtxi))*sig
#isto kao prethodno samo krace (regresija S^(-1)y na S^(-1)x)
sm<-chol(Sigma)
smi<-solve(t(sm))
sx<-smi %*% x
sy<-smi %*% longley$Empl
summary(lm(sy ~ sx-1))
#korelacija iz novog modela
cor(res[-1],res[-16])
#prethodna procedura treba da se ponavlja do konvergencije ka nekoj vrednosti
#problem je sto mi ne znamo koji je model u pitanju, pa je bolje koristiti ugradjene funkcije

#drugi nacin: ugradjene funkcije
library(nlme)
# corAR1 nam kaze da je u pitanju AR(1) model, gde se podaci redjaju po godinama
g<-gls(Employed~GNP+Population, correlation = corAR1(form= ~Year), data=longley)
summary(g)
#phi je korelacija
intervals(g)
#na osnovu podataka ne mozemo zakljuciti nista, disperzija ocene je velika
#posto nula ulazi u interval, nema dokaza o korelaciji uzastopnih gresaka
#ne postoji jaka korelacija koja ce izazvati drugacije ocene linearnog modela

#2. Tezinski najmanji kvadrati 
#greske su nekorelisane, ali nemaju jednaku disperziju
#ponderisemo promenljive
data(fpe)
fpe
#kandidati A i B su prosli u drugi krug, A2 i B2 su brojevi glasova u drugom krugu
#disperzija je srazmerna broju glasaca, pa je ponder 1/EI
g<-lm(A2~A+B+C+D+E+F+G+H+J+K+N-1, fpe, weights=1/EI)
coef(g)
lm(A2~A+B+C+D+E+F+G+H+J+K+N-1, fpe)$coef #obicni najmanji kvadrati
#na neke koeficijente nisu mnogo uticale tezine dok su se drugi malo vise promenili

#sve tezine mozemo pomnoziti nekom konstantom i dobicemo iste ocene
lm(A2~A+B+C+D+E+F+G+H+J+K+N-1, fpe, weights=53/EI)$coef

#nema smisla da koeficijent bude negativan 
#vidimo da su vrednosti za A priblizno jednake 1, a za B nula
#nama su interesantni oni koeficijenti koji su izmedju 0 i 1
lm(A2~offset(A+G+K)+C+D+E+F+N-1, fpe, weights=1/EI)$coef #A,G i K imaju koeficijent 1 

#3. Testiranje nedostataka modela
#kako proveriti da li se podaci NE slazu sa modelom?
#sto imamo vise podataka, to pretpostavljamo da cemo imati bolji model
#ne mozemo da proveravamo da li podaci koje imamo odgovaraju modelu koji smo napravili od njih (zato sto je model pravljen na osnovu njih)
#treba oceniti sigma iz modela i proveriti da li je ocena modela nepristasna
#tj. treba da ocenimo model na osnovu nekih podataka na osnovu kojih on NIJE pravljen 
#ako imamo ponovljene vrednosti u uzorku, onda mozemo na nekim od tih vrednosti da testiramo model
data(corrosion) #gubitak tezine predmeta u zavisnosti od kolicine gvozdja
plot(loss~Fe, corrosion,xlab="Iron content",ylab="Weight loss") 
#sa slike vidimo da ima ponovljenih vrednosti za x, sa razlicitim vrednostima za y
#u tom slucaju mozemo te ponovljene vrednosti za x smatrati kao jednu, a za y uzeti sredinu odgovarajucih y
g<-lm(loss~Fe, corrosion)
summary(g)
abline(coef(g)) #regresiona linija
#pravimo model u kom je prediktor faktor (umesto 13 podataka, imamo 7 razlicitih)
ga<-lm(loss~factor(Fe), corrosion)
summary(ga)
points(corrosion$Fe,fitted(ga),pch=18) #za svaku mogucu vrednost Fe, racuna se srednja vrednost loss
#poredimo prethodna dva modela
anova(g,ga)
#=> postoje nedostaci u modelu (bolji je onaj sa faktorom)

# pokusavamo plinomijalnu regresiju
gp<-lm(loss~Fe+I(Fe^2)+I(Fe^3)+I(Fe^4)+I(Fe^5)+I(Fe^6),corrosion)
plot(loss~Fe, data=corrosion,ylim=c(60,130))
points(corrosion$Fe,fitted(ga),pch=18)
grid<-seq(0,2,len=50)
lines(grid,predict(gp,data.frame(Fe=grid)))
#dobijamo model koji idealno prolazi kroz sve tacke
#ali ne postoji realni razlog da se on koristi, pogotovo sto je disperzija velika
#jedino sto se njim postize je veliko R-squared
summary(gp)$r.squared

#4. Robustna regresija
#primenjuje se kad imamo rapodele sa ekstremnim vrednostima koje jako uticu na ocene
#robusne ocene nisu osetljive na narusenost pretpostavki
#umesto da minimiziramo sumu kvadrata, mozemo minimizirati na primer sumu apsolutinh vrednosti
#LSM ako su x-evi ograniceni nekom konstantom, a LAD ako nisu
data(gala)
gl<-lm(Species~Area+Elevation+Nearest+Scruz+Adjacent,gala)
summary(gl)
library(MASS)
#rlm = robustni linearni model, podrazumevan Huberov metod
#konstanta c se dobija kao ocena parametra sigma
gr<-rlm(Species~Area+Elevation+Nearest+Scruz+Adjacent,gala)
summary(gr)
#R^2 i F-test nisu dati jer se ne mogu izracunati (bar ne u istom smislu)
#p-vrednost t-testa takodje nije data, iako se moze koristiti priblizna normalnost ocena

#LAD regresije
install.packages("quantreg")
library(quantreg)
attach(gala)
gq<-rq(Species~Area+Elevation+Nearest+Scruz+Adjacent)
summary(gq)
#nije data stdardana greska ni p vrednost, samo koeficijenti i gornje i donje granice
#=> nema znacajnih razlika u koeficijentima sa LSM, pa se moze koristiti LSM zbog jednostavnosti
detach(gala)

#Zaseceni najmanji kvadrati (LTS)
#pravimo ocenu tako sto minimizujemo sumu prvih q reziduala
# ovo je postojana metoda i bas je dobra ukoliko se zbog prirode podataka zna da postoje veliki autlajeri

install.packages("lqs")
library(lqs)
g<-ltsreg(Species~Area+Elevation+Nearest+Scruz+Adjacent,gala)
coef(g)
g<-ltsreg(Species~Area+Elevation+Nearest+Scruz+Adjacent,gala)
coef(g)
#LTS regression = ltsreg
#algoritam nije deterministicki, zato su navedene dve iste naredbe a dobijeni su razliciti koeficijenti
g<-ltsreg(Species~Area+Elevation+Nearest+Scruz+Adjacent, gala,nsamp="exact") #spor metod za velike baze jer uzima koliko god je potrebno uzoraka
coef(g)
#dobijamo samo ocene, ne i stdardne greske ocena, ne mozemo da nadjemo raspodelu ocena
#intervali poverenja se odredjuju nekim simulacionim metodama (bootstrap metodom)

#Simluacione metode
#pretpostavljamo da epsilon ima neku raspodelu, generisemo veliki uzorak iz te raspodele
#i odredimo ocene za beta a onda imamo i dobru aproksimaciju raspodele za beta
#Butstrep
sample(10,rep=T) #pravi 10 brojeva iz {1,2,...,10} sa ponavljanjem
residuals(g)[sample(30,rep=TRUE)] #slucajan uzorak sa ponavljanjem LTS reziduala
x<-model.matrix( ~ Area+Elevation+Nearest+Scruz+Adjacent,gala)[,-1] #matrica modela
bcoef<-matrix(0,1000,6) #matrica za smestanje rezultata
#ponavljamo butstrep 1000 puta
for(i in 1:1000){
  newy<-predict(g)+residuals(g)[sample(30,rep=T)]
  brg<-ltsreg(x,newy,nsamp="best") #nsamp="best" do 5000 uzoraka
  bcoef[i,]<-brg$coef
}
#95% interval poverenja uzimanjem empirijskih kvantila za Area
quantile(bcoef[,2],c(0.025,0.975)) #nula je van ovog intervala
#raspodela dobijena na osnovu butstrep metode
plot(density(bcoef[,2]),xlab="Coefficient of Area",main="")
#interval poverenja
abline(v=quantile(bcoef[,2],c(0.025,0.975)))
#ranije je dobijeno da je koeficijent uz Area oko -0.02, a sad je značajno veći
#postoji neka vrednost koja ima ogroman rezidual, tek kada se odbaci ta vrednost dobija se potpuno drugaciji koeficijent
#to je ostrvo Izabela, kad se odrede Kukova rastojanja
gli<-lm(Species~Area+Elevation+Nearest+Scruz+Adjacent,  gala,subset=(row.names(gala) != "Isabela"))
summary(gli)
#Adjacent je postao vrlo znacajan

#slicni prikaz na drugoj bazi
data(star)
plot(light ~ temp, star)
#metod najmanjeg kvadrata
gs1<-lm(light ~ temp, star)     
abline(coef(gs1))
#Huberov metod
gs2 <- rlm(light ~ temp, star)    
abline(coef(gs2), lty=2) #skoro uopste se ne razlikuje od prave dobijene pomocu LSM
#zaseceni najmanji kvadrati
gs3<-ltsreg(light ~ temp, star, nsamp="exact")    
abline(coef(gs3), lty=5) #skroz drugacije i deluje kao da se mnogo bolje uklapa
#ignorisani su autlajeri, dobijen je smisleniji model

#ako se ocene dobijene pomocu metoda najmanjeg kvadrata ne razlikuju mnogo od ocena dobijenih pomocu robusnih metoda
#to je onda znak da su ocene po metodu najmanjeg kvadrata dobre