# 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)
 
Z=read.table("hw3.dat",header=T)[,-1]
Zmat=as.matrix(Z)
bmi=Zmat[,2]/((Zmat[,1]/100)^2); Zb=cbind(bmi,Zmat[,-(1:2)]) 
# Zb, to macierz danych z dodana kolumna bmi - "body mass index"
# oraz usunietymi kolumnami "height" i "weight". Postepuje jak Tracy Bergemann
# wierzac, ze bmi kompresuje informacje zawarta w ht i wt.

# 1. Diagnostyka obserwacji odstajacych
X=Zb; p1=ncol(X)
m1=lm(X[,p1]~X[,-p1])
x11(); par(mfrow=c(3,2))
plot(m1,which=1:6)
# print(summary(influence.measures(m1)))    # obs wplywowe

# diagnostyka obserwacji odstajacych - wersja uproszczona
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[,"fat"],type="n"); abline(0,1,col=2)
text(m1$fit,X[,"fat"],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*86),86-9-1) = 
# -3.375614. Nastepnie sprawdzamy, ze jest jedno residuum przekraczajace te 
# wartosc: res[74] =-3.540783. Zatem odrzucamy H0 i usuwamy obs 74. Ogolnie: 
X=X[ -which( abs(res) > -qt(.1/(2*86),76) ) , ] 

# 2. Transformacje zmiennych
x11(); par(mfrow=c(4,5))
for(i in 1:9)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:9)plot(density(log(X[,i])),main=paste("log-",colnames(X)[i]))
# widac, ze transformacja log X nie zbliza zmiennych do  "normalnosci"

library(MASS); x11(); par(mfrow=c(1,2))
logtrans(X[,p1]~X[,-p1],alpha=seq(-5,40,len=100))
boxcox(X[,p1]~X[,-p1],lambda=seq(-2,3,len=100))
# boxcox wybiera identycznosc, czyli nic nie przeksztalcamy!

# 3. Selekcja cech
print(round(cor(X),2)) # nie ma cech bardzo (|cor|>0.9) ze soba skorelowanych

m2=lm(fat~.,data=data.frame(X)) 
print(summary(m2)) 
# Widac, ze istotne cechy to: tri, bmi, wrist, scap. Do takiego samego modelu
# prowadzi krokowa optymalizacja bayesowskiego kryterium informacyjnego (BIC):
# wykonuje stepAIC z paramentrem k= log(n), gdzie n - liczba obs.
m3=stepAIC(m2,scope=list(upper=.~.,lower=.~1),trace=T,k=log(nrow(X)))
print("**************************")
# Nastepnie sprawdzam istotnosc pozostalych cech, usuwajac po jednej:
m4=update(m3,.~.-wrist); print(dropterm(m4,test="F",k=log(nrow(X))))  # model Tracy Bergemann
m5=update(m3,.~.-scap);  print(dropterm(m5,test="F",k=log(nrow(X))))
m6=update(m5,.~.-wrist); print(dropterm(m6,test="F",k=log(nrow(X))))

print("**************************")
print(summary(m4))
print(round(cor(X[,"fat"],m4$fit),2)) # jakosc modelu Tracy

# Sprawdzilem definicje R-Squared  i   Adjusted R-squared
# R2 - fraction of variance of y which is explained by the model
# R2=cor(X[,"fat"],m3$fit)^2
# R2adj = function(R2,n,p) 1 - (n - 1) / (n - p) * (1 - R2)
# gdzie n - liczba obs, p - liczba parametrow modelu, czyli length(m3$coef)