# Modelowanie przeprowadzone w kolejnosci proponowanej przez J. Farawaya:
# diagnostyka obs odstajacych, transformacje, selekcja cech
# (str 134 "Practical Regression and ANOVA using R" - ksiazka dostepna w sieci)

Auta=read.table("samochody.txt",skip=9,header=TRUE)
X=as.matrix(Auta)

# 1. Diagnostyka obserwacji odstajacych - wersja uproszczona
m1=lm(X[,1]~X[,-1])
x11(); par(mfrow=c(1,2))
res=rstudent(m1); hii=hatvalues(m1)
plot(hii,res,type="n"); abline(h=0,col=2)
text(hii,res,as.character(1:nrow(X)),cex=.75)
plot(m1$fit,X[,"Paliwo"],type="n"); abline(0,1,col=2)
text(m1$fit,X[,"Paliwo"],as.character(1:nrow(X)),cex=.75)
# Testuje teraz hipoteze H0 == "nie ma obserwacji odstajacych" na poziomie 
# alfa=0.1 metoda Bonferroniego, czyli wykonuje ciag testow  T[i] dla hipotez
# H[i] == "i-ta obs nie odstaje" na poziomie alfa/n, gdzie n - liczba obs.
# Wiemy, ze res[i] ~ t(n-p-1), gdzie p - liczba param modelu.
# Zatem do testowania H[i] mozemy wykorzystac (dwustronny) test t-studenta.
# Obliczamy wartosc krytyczna qt(.1/(2*n),n-p-1) = qt(.1/(2*24),24-7-1) = 
# -3.338551. Nastepnie sprawdzamy, ze nie ma residuum przekraczajacego te wartosc.
# Zatem nie odrzucamy H0. Ogolnie: X=X[ -which( abs(res) > -qt(.1/(2*n),n-p-1) ) , ] 

# 2. Transformacje zmiennych
x11(); par(mfrow=c(4,4))
for(i in 1:7)plot(density(X[,i]),main=colnames(X)[i])
plot(1,1,type="n",axes=F,xlab="",ylab="") # dla symetrii miedzy X[,i] i log(X[,i]
for(i in 1:7)plot(density(log(X[,i])),main=paste("log-",colnames(X)[i]))
# widac, ze transformacja log(X) zbliza "Paliwo", "Pojemn" i "Moc" do  
# "normalnosci". Zatem
X[,2:3]=log(X[,2:3]) # Paliwo zmienie po "konsultacji" z logtrans i boxcox:

library(MASS); x11(); par(mfrow=c(1,2))
logtrans(X[,1]~X[,-1],alpha=seq(-5.4,0,len=100))
# logtrans uzasadnia przeksztalcenie log(X[,1]-5)
boxcox(X[,1]-5~X[,-1],lambda=seq(-2,2,len=100))
# boxcox potwierdza log(X[,1]-5)
X[,1]=log(X[,1]-5)
colnames(X)[1:3]=c("logPaliwo","logPojemn","logMoc")

# 3. Selekcja cech
print(round(cor(X),2)) 
# Sa cechy b ze soba skorelowane, czyli prawie wspolliniowe. 
print("**************************")

m2=lm(logPaliwo~.,data=data.frame(X)) 
print(summary(m2))
print("**************************")
# Najistotniejsza jest logPojemn, ale tylko na "." - chyba z powodu wspolliniowosci.
# Do takiego samego modelu prowadzi krokowa optymalizacja bayesowskiego kryterium 
# informacyjnego (BIC): wykonuje stepAIC z paramentrem k= log(n).
m3=stepAIC(m2,scope=list(upper=.~.,lower=.~1),trace=T,k=log(nrow(X)))
print("**************************")
print(summary(m3))
# Nastepnie sprawdzam istotnosc pozostalych cech, usuwajac po jednej:
m4=update(m3,.~.-Masa); print(dropterm(m4,test="F")); print(summary(m4))  
print("**************************")
m5=update(m3,.~.-logPojemn);  print(dropterm(m5,test="F")); print(summary(m5))
print("**************************")
# Wybieram model m4, bo nie duzo traci do m2 i m3 na dokladnosci i jest prosty
print(round(cor(X[,"logPaliwo"],m4$fit),2)) # jakosc ostatecznego modelu