##  1978 Automobile Data - no missings

##  make            Make and Model
##  price           Price
##  mpg             Mileage (mpg - mile per gallon)
##  rep78           Repair Record 1978
##  headroom        Headroom (in.)
##  trunk           Trunk space (cu. ft.)
##  weight          Weight (lbs.)
##  length          Length (in.)
##  turn            Turn Circle (ft.) 
##  displacement    Displacement (cu. in.)
##  gear_ratio      Gear Ratio
##  foreign         Car type

automobile <- read.table("automobile.txt", header=TRUE,sep=",")
automobile$spotreba <-  1/automobile$mpg * 3.78541178/1.609344 * 100
attach(automobile)



#################################################################
# model s obema prediktory aditivne

model1 <- lm(spotreba ~ weight+foreign,model = TRUE, x = TRUE, y = TRUE)
summary(model1)

plot(spotreba ~ weight, 
     xlab = "Hmotnost [lbs]",
     ylab = "Spotreba [l/100km]",
     type="n")
     
points(weight[foreign=="Domestic"],spotreba[foreign=="Domestic"],pch=1)
points(weight[foreign=="Foreign"],spotreba[foreign=="Foreign"],pch=2)

#slozitejsi zpusob, potreba u slozitejsich modelu
newweightDom <- data.frame(weight = seq(min(weight), max(weight), length.out = 100), foreign = rep("Domestic",100))
newweightFor <- data.frame(weight = seq(min(weight), max(weight), length.out = 100), foreign = rep("Foreign",100))

lines(predict(model1, newweightDom) ~ newweightDom$weight)
lines(predict(model1, newweightFor) ~ newweightFor$weight,lty=2)


#################################################################
# model s obema prediktory - interakce

#dva totozne zapisy
model2 <- lm(spotreba ~ weight+foreign+weight:foreign,model = TRUE, x = TRUE, y = TRUE)
model2 <- lm(spotreba ~ weight+foreign+weight:foreign,model = TRUE, x = TRUE, y = TRUE)

summary(model2)

plot(spotreba ~ weight, 
     xlab = "Hmotnost [lbs]",
     ylab = "Spotreba [l/100km]",
     type="n")
     
points(weight[foreign=="Domestic"],spotreba[foreign=="Domestic"],pch=1)
points(weight[foreign=="Foreign"],spotreba[foreign=="Foreign"],pch=2)

#slozitejsi zpusob, potreba u slozitejsich modelu
newweightDom <- data.frame(weight = seq(min(weight), max(weight), length.out = 100), foreign = rep("Domestic",100))
newweightFor <- data.frame(weight = seq(min(weight), max(weight), length.out = 100), foreign = rep("Foreign",100))

lines(predict(model2, newweightDom) ~ newweightDom$weight)
lines(predict(model2, newweightFor) ~ newweightFor$weight,lty=2)

text(weight[25],spotreba[25],make[25],pos=2)
text(weight[33],spotreba[33],make[33],pos=1)
text(weight[11],spotreba[11],make[11],pos=3)

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

#F-test vyznamnosti noveho clenu

anova(model2)
anova(model1,model2)


################
#GRAFICKE OVERENI PREDPOKLADU

#
plot(resid(model1)~weight,
		xlab = "Weight",
    ylab = "Rezidua",
    main = "Nezávisle proměnné"
)
abline(0,0,lty=2)
text(weight[abs(resid(model1))>3],resid(model1)[abs(resid(model1))>3],make[abs(resid(model1))>3],pos=4)
automobile[abs(resid(model1))>3,]




#zakladni checkovaci funkce v eRku
plot(model1,labels.id=make)

#residuals vs. fitted
plot(model1,labels.id=make, which=1)

#QQ plot
plot(model1,labels.id=make, which=2)

#pakove body
plot(hatvalues(model1)~weight,
		xlab = "Weight",
    ylab = "Leverages"
)
text(weight[abs(hatvalues(model1))>0.08],hatvalues(model1)[abs(hatvalues(model1))>0.08],make[abs(hatvalues(model1))>0.08],pos=2)
automobile[abs(hatvalues(model1))>0.08,]

#DFFITS
plot(dffits(model1)~weight,
		xlab = "Weight",
    ylab = "DFFITS"
)
text(weight[abs(dffits(model1))>0.5],dffits(model1)[abs(dffits(model1))>0.5],make[abs(dffits(model1))>0.5],pos=4)
automobile[abs(dffits(model1))>0.5,]

#DFBETAS
matplot(weight,dfbetas(model1),xlab="Weight",ylab="DFBETAS",type="p")
text(weight[abs(dfbetas(model1))[,2]>0.4],dfbetas(model1)[abs(dfbetas(model1))[,2]>0.4,2],make[abs(dfbetas(model1))[,2]>0.4],pos=4)
automobile[abs(dfbetas(model1))[,2]>0.4,]

#Cookovy vzdalenosti
plot(cooks.distance(model1)~weight,
		xlab = "Weight",
    ylab = "Cookova vzdálenost"
)
text(weight[abs(cooks.distance(model1))>0.1],cooks.distance(model1)[abs(cooks.distance(model1))>0.1],make[abs(cooks.distance(model1))>0.1],pos=4)
automobile[abs(cooks.distance(model1))>0.1,]



#co s tim modelem bez odlehlych hodnot

automobile2 <- automobile[abs(cooks.distance(model1))<=0.1,]
automobile2_out <- automobile[abs(cooks.distance(model1))>0.1,]
model3 <- lm(spotreba ~ weight*foreign,data=automobile2)

plot(spotreba ~ weight, 
     xlab = "Hmotnost [lbs]",
     ylab = "Spotreba [l/100km]",
     type="n")
     
points(weight[foreign=="Domestic"],spotreba[foreign=="Domestic"],pch=1)
points(weight[foreign=="Foreign"],spotreba[foreign=="Foreign"],pch=2)

#slozitejsi zpusob, potreba u slozitejsich modelu
newweightDom <- data.frame(weight = seq(min(weight), max(weight), length.out = 100), foreign = rep("Domestic",100))
newweightFor <- data.frame(weight = seq(min(weight), max(weight), length.out = 100), foreign = rep("Foreign",100))

lines(predict(model2, newweightDom) ~ newweightDom$weight)
lines(predict(model2, newweightFor) ~ newweightFor$weight,lty=2)

###

points(automobile2_out$weight[automobile2_out$foreign=="Domestic"],automobile2_out$spotreba[automobile2_out$foreign=="Domestic"],pch=1, col="red")
points(automobile2_out$weight[automobile2_out$foreign=="Foreign"],automobile2_out$spotreba[automobile2_out$foreign=="Foreign"],pch=2, col="red")

newweightDom <- data.frame(weight = seq(min(weight), max(weight), length.out = 100), foreign = rep("Domestic",100))
newweightFor <- data.frame(weight = seq(min(weight), max(weight), length.out = 100), foreign = rep("Foreign",100))

lines(predict(model3, newweightDom) ~ newweightDom$weight, col="red")
lines(predict(model3, newweightFor) ~ newweightFor$weight,lty=2, col="red")

#vyznamnost???

summary(model2)
summary(model3)

# centrovani - tak aby byl interpretovatelny rozdil mezi primkami
model3 <- lm(spotreba ~ I(weight-3000)*foreign,data=automobile2)












