################# -1. ECO


library(faraway)
data(eco)
help(eco)
dim(eco)


m1=lm(income~usborn,eco)
summary(m1)
plot(m1)

library(nortest)
res=residuals(m1)
ad.test(res)
lillie.test(res)
sf.test(res)
pearson.test(res)
## ok


plot(res[1:76],res[2:77]) ### czy cos widac???
dwtest(m1) ### raczej nie
bgtest(m1) ###tak samo jak przy dwtest, to byla fantazja  
library(car)
durbinWatsonTest(m1,max.lag=3)
### OK



library(lmtest)
bptest(m1)
gqtest(m1)
# homoskedastycznosc tak se


resettest(m1)
raintest(m1) 
# jest problem z funkcyjna zaleznoscia


## jak to wyglada w rzeczywistosci?
plot(income~usborn,data=eco)


boxcox(m1)
help(boxcox)
boxcox(m1, lambda = seq(-3, 1, len = 30))

m2=lm(I(income^(-1))~usborn,eco)

summary(m2)
plot(m2)

library(nortest)
res=residuals(m2)
ad.test(res)
lillie.test(res)
sf.test(res)
pearson.test(res)
## slicznie

library(lmtest)
bptest(m2)
gqtest(m2)
# homoskedastycznosc ok (drobna poprawa)

plot(res[1:76],res[2:77]) ### czy cos widac???
dwtest(m2) ### nic a nic
bgtest(m2) ###tak samo jak przy dwtest, czysta sprawa  
library(car)
durbinWatsonTest(m2,max.lag=3)
### OK

resettest(m2)
raintest(m2) 
# problem z funkcyjna zaleznoscia jest znacznie mniejszy 




m2b=lm(I(log(income))~usborn,eco)

summary(m2b)
plot(m2b,1)

library(nortest)
res=residuals(m2b)
ad.test(res)
lillie.test(res)
sf.test(res)
pearson.test(res)
## slicznie


plot(res[1:76],res[2:77]) ### czy cos widac???
dwtest(m2b) ### nic a nic
bgtest(m2b) ###tak samo jak przy dwtest, czysta sprawa  
library(car)
durbinWatsonTest(m2b,max.lag=3)
### OK

library(lmtest)
bptest(m2b)
gqtest(m2b)
# homoskedastycznosc tak se (bardzo drobna poprawa)



resettest(m2b)
raintest(m2b) 
# problem z funkcyjna zaleznoscia jest mniejszy



data<-eco
r1=0
r2=0
r2b=0

for( j in (1:dim(data)[1]))
{
  training=data[-j,]
  test=data[j,]
  
  m1=lm(income~usborn,training)
  m2=lm(I(income^(-1))~usborn,training)
  m2b=lm(I(log(income))~usborn,training)
  
  r1=r1+abs(predict(m1,test)- test$income    )
  r2=r2+abs((predict(m2,test) )^(-1)- test$income    )
  r2b=r2b+abs(exp(predict(m2b,test) )- test$income    )
  
  
}


rmse1=r1/length(data)
rmse2=r2/length(data)
rmse2b=r2b/length(data)

rmse1
rmse2
rmse2b

#### ktory model??

m1=lm(income~usborn,eco)
m2=lm(I(income^(-1))~usborn,eco)
m2b=lm(I(log(income))~usborn,eco)

plot(income~usborn,eco)
plot(I(income^(-1))~usborn,eco)
plot(I(log(income))~usborn,eco)

########################

m3=lm(income~usborn+I(usborn^2),eco)
summary(m3)

library(nortest)
res=residuals(m3)
ad.test(res)
lillie.test(res)
sf.test(res)
pearson.test(res)
## ok

plot(res[1:76],res[2:77]) ### czy cos widac???
dwtest(m3) ### raczej nie
bgtest(m3) ###tak samo jak przy dwtest, to byla fantazja  
library(car)
durbinWatsonTest(m3,max.lag=3)
### OK

library(lmtest)
bptest(m3)
gqtest(m3)
# homoskedastycznosc tak se



resettest(m3)
raintest(m3) 
# nie ma problemu z funkcyjna zaleznoscia


data<-eco
r1=0
r2=0

for( j in (1:dim(data)[1]))
{
  training=data[-j,]
  test=data[j,]
  
  m1=lm(income~usborn,training)
  m2=lm(income~usborn+I(usborn^2),training)
  
  r1=r1+abs(predict(m1,test)- test$income    )
  r2=r2+abs((predict(m2,test) )- test$income    )
  
}


rmse1=r1/length(data)
rmse2=r2/length(data)

rmse1
rmse2

## bez szalu, duzy klopot z interpretacja

##############################################################################
##############################################################################
##############################################################################
##############################################################################



############################# 0. Divorse 

library(faraway)
data=divusa
?divusa
summary(divusa)
m1=lm(divorce~.,divusa)
summary(m1)
par(mfrow=c(3,2))
plot(m1,1) 
plot(m1,2) 
plot(m1,3)
plot(m1,4)
plot(m1,5)
plot(m1,6)
par(mfrow=c(1,1))
qqPlot(m1)

#testy na normalność
library(nortest)
res=residuals(m1)
ad.test(res)
lillie.test(res)
sf.test(res)
pearson.test(res)
##wszystkie p_value powyżej 10% więc residua pochodzą z rokładu normalnego, nie podstaw do odrzucenia


##autokorelacja 
par(mfrow=c(1,1))
plot(res[1:76],res[2:77]) ###jest autokorelacja
dwtest(m1) ###p_value małe czyli jest autokorelacja
bgtest(m1) ###tak samo jak przy dwtest  
library(car)
durbinWatsonTest(m1,max.lag=3)


##homoskedastyczność
library(lmtest)
bptest(m1)
gqtest(m1)
plot(m1,1)
##mamy stałość wariancji?



###funkcyjna zależność
resettest(m1) 
raintest(m1) 
##jest problem...
##współzależność zmiennych 
vif(m1)  ###zalezne year i femlab

##########################################################################

###model bez femlab
m2=lm(divorce~.-femlab,divusa)
summary(m2)
par(mfrow=c(3,2))
plot(m2,1) 
plot(m2,2) 
plot(m2,3)
plot(m2,4)
plot(m2,5)
plot(m2,6)
par(mfrow=c(1,1))

#testy na normalność
qqPlot(m2)
library(nortest)
ad.test(residuals(m2))
lillie.test(residuals(m2))
sf.test(residuals(m2))
pearson.test(residuals(m2))
## p_value powyżej 5% więc residua pochodzą z rokładu normalnego, nie podstaw do odrzucenia

##autokorelacja 
res=residuals(m2)
par(mfrow=c(1,1))
plot(res[1:76],res[2:77]) ###jest autokorelacji
dwtest(m2) ###p_value małe czyli jest autokorelacja
bgtest(m2) ###tak samo jak przy dwtest  


##homoskedastyczność
library(lmtest)
bptest(m2)
gqtest(m2)
##nic nie możemy powiedzieć


###funkcyjna zależność
resettest(m2) 
raintest(m2) 

##współzależność zmiennych 
vif(m2)
###wszystko OK

boxcox(m2)


m3=lm(formula = I(divorce^(-1)) ~ . - femlab, data = divusa)

#testy na normalność
qqPlot(m3)
library(nortest)
ad.test(residuals(m3))
lillie.test(residuals(m3))
sf.test(residuals(m3))
pearson.test(residuals(m3))
## ok

##homoskedastyczność
library(lmtest)
bptest(m3)
gqtest(m3)
##super

##autokorelacja 
res=residuals(m3)
par(mfrow=c(1,1))
plot(res[1:76],res[2:77]) ###jest autokorelacji
dwtest(m3) ###p_value małe czyli jest autokorelacja
bgtest(m3) ###tak samo jak przy dwtest  

###funkcyjna zależność
resettest(m3) 
raintest(m3) 
## znaczna poprawa wlasciwie OK


##współzależność zmiennych 
vif(m3)
###wszystko OK


summary(m3)
# wszystko istotne

step(m3)
#nic nie zmienia


#################################################################

###testowanie krzyżowe modeli m1 i m3
data
r1=0
r2=0

for( j in (1:length(data)))
{
  training=data[-j,]
  test=data[j,]
  
  m1=lm(divorce~year+unemployed+marriage+birth+military,training)
  m2=lm( I(divorce^(-1)) ~ (year + unemployed + femlab +marriage + birth + military) - femlab, training)
  
  
  r1=r1+abs((predict(m1,test) )- test$divorce    )
  r2=r2+abs((predict(m2,test) )^(-1)- test$divorce    )
  
  
}


rmse1=r1/length(data)
rmse2=r2/length(data)
rmse1
rmse2
###wygrywa model 2, czyli m





#############################################################################
#############################################################################
#############################################################################


####################### 1. Chicago

data(chredlin)
help(chredlin)
summary(chredlin)


m2=lm(involact~.-income+log(income)-side,chredlin)
summary(m2)


# normalność rezyduów
res=m2$residuals
library(nortest)
ad.test(res)
lillie.test(res)
sf.test(res)
pearson.test(res)
library(car)
qqPlot(res)
#brak problemów z normalnością rezyduów


#autokorelacja
library(lmtest)
durbinWatsonTest(m2,max.lag=3)
bptest(m2)
#brak problemu z korelacją rezyduów


#homoskedastyczność wariancji
bptest(m2)
gqtest(m2)
plot(m2,1)
hmctest(m2)
#uznajemy, że brak problemu - początkowa linia jest tam gdzie involact=0



# poprawność struktury modelu
resettest(m2)
raintest(m2)
harvtest(m2)
library(MASS)
boxcox(m2)
m3=lm(involact+0.000001~.,chredlin) #śliska sprawa
boxcox(m3)
#odrzucamy transformacje zmiennej objaśnianej ze względu na późniejsze kłopoty z interpretacją

#punkty wpływowe
p=dim(chredlin)[2]-1
n=length(chredlin$race)
2*p/n #granica dużej dźwigni
4/n
plot(m2,3:5)
rstandard(m2)[abs(rstandard(m2))>2]
rstudent(m2)[abs(rstudent(m2))>qt(1-.5/(2*n),n-p-1)]
#zostawiamy je w modelu


vif(m2)
kappa(model.matrix(m2))
# przesłanka do usuwania zmiennych

step(m2)
m3=lm(formula = involact ~ race + fire + theft + age, data = chredlin)
summary(m3)
kappa(model.matrix(m3))

library(leaps)
models=regsubsets(involact~.-income+log(income)-side,force.in=1,chredlin)
(modelssum=summary(models))
cps=modelssum$cp
r2a=modelssum$adjr2
plot(2:5,cps)
abline(0,1)
plot(r2a)

#alternatywny model: 
m4=lm(involact~race+fire+log(income),chredlin)
summary(m4)

#wnioski matematyczne: dobrze dobrany model, konieczność eliminacji zmiennych
#wnioski z analizy: rasa ma małe, ale jednak statystycznie istotne znaczenie przy odmowie zawarcie oferty ubezpieczenia



#######################2. GALA   %%%%%%%%%%%%

data(gala)
summary(gala)


m6=lm(Species~.-Endemics,gala)
summary(m6)


# normalność rezyduów
res=m6$residuals
library(nortest)
ad.test(res)
lillie.test(res)
sf.test(res)
pearson.test(res)
library(car)
qqPlot(res)
#mały problem z normalnością rezyduów


#autokorelacja
library(lmtest)
durbinWatsonTest(m6,max.lag=3)
bgtest(m6)
#brak problemu z korelacją rezyduów


#homoskedastyczność wariancji
bptest(m6)
gqtest(m6)
plot(m6,1)
hmctest(m6)

#ewidentny problem a homoskedastycznością wariancji
plot(res~Area,gala)
plot(res~Nearest,gala)
plot(res~Scruz,gala)
plot(res~Elevation,gala)
plot(res~Adjacent,gala) #widoczna zależność

boxcox(m6)

m7=lm((Species)^(1/3)~.-Endemics,gala)
summary(m6)
summary(m7)

bptest(m7)
gqtest(m7)
plot(m7,1)
hmctest(m6)
#rozwiązany problem z homoskedastycznością wariancji
plot(res~Area,gala)
plot(res~Nearest,gala)
plot(res~Scruz,gala)
plot(res~Elevation,gala)
plot(res~Adjacent,gala) #widoczna zależność

# poprawność struktury modelu
resettest(m7)
raintest(m7)
harvtest(m7)

#punkty wpływowe
p=6
n=length(gala$Area)
2*p/n #granica dużej dźwigni
4/n
plot(m7,3:5)
rstandard(m2)[abs(rstandard(m2))>2]
rstudent(m2)[abs(rstudent(m2))>qt(1-.5/(2*n),n-p-1)]
#usuwamy Isabelę - inna kategoria wysp

gala
gala1=gala[-16,]
#gala1=gala1[-12,] zobaczyć co tutaj się stanie....
gala1

m8=lm((Species)^(1/3)~.-Endemics,gala1)
summary(m8)
summary(m7)

#punkty wpływowe
p=6
n=length(gala1$Area)
2*p/n #granica dużej dźwigni
4/n
plot(m8)
rstandard(m8)[abs(rstandard(m8))>2]
rstudent(m8)[abs(rstudent(m8))>qt(1-.5/(2*n),n-p-1)]
#zostawiamy je w modelu





vif(m8)
kappa(model.matrix(m8))
# przesłanka do usuwania zmiennych

step(m8)
m9=lm((Species)^(1/3) ~ Area + Elevation  + Adjacent, data = gala1)
summary(m9)



#wnioski matematyczne: dość dobrze dobrany model, konieczność stosowania transformacji zmiennej objaśnianej
#wnioski z analizy: decydujący wpływ na ilość gatunków mają wysokość i powierzchnia wyspy

help(gala)



######################3. HAPPY   

data(happy)
help(happy)
summary(happy)


m5=lm(happy~.,happy)
summary(m5)


# normalność rezyduów
res=m5$residuals
library(nortest)
ad.test(res)
lillie.test(res)
sf.test(res)
pearson.test(res)
library(car)
qqPlot(res)
#brak problemów z normalnością rezyduów


#autokorelacja
library(lmtest)
durbinWatsonTest(m5,max.lag=3)
bptest(m5)
#problem z korelacją rezyduów


#homoskedastyczność wariancji
bptest(m5)
gqtest(m5)
plot(m5,1)
hmctest(m5)
#uznajemy, że brak problemu 



# poprawność struktury modelu
resettest(m5)
raintest(m5)
harvtest(m5)
library(MASS)
boxcox(m5)
#odrzucamy transformacje zmiennej objaśnianej ze względu na późniejsze kłopoty z interpretacją

#punkty wpływowe
p=5
n=length(happy$money)
2*p/n #granica dużej dźwigni
4/n
plot(m5,3:5)
rstandard(m5)[abs(rstandard(m5))>2]
rstudent(m5)[abs(rstudent(m5))>qt(1-.5/(2*n),n-p-1)]
happy[36,]
#zostawiamy je w modelu


vif(m5)
kappa(model.matrix(m5))
# przesłanka do usuwania zmiennych

step(m5)
m6=lm(formula = happy ~ money + love + work, data = happy)
summary(m6)
kappa(model.matrix(m6))