#############
## Cviceni 1
#############

delta.t <- c(1:6)*10^4
prodlouzeni <- c(0.18, 0.35, 0.48, 0.65, 0.84, 0.97)
regr.primka.pocatek <- lm(prodlouzeni~delta.t-1)
summary(regr.primka.pocatek)
confint(regr.primka.pocatek)

plot(prodlouzeni~delta.t, main="Tepelna roztaznost")
lines(fitted(regr.primka.pocatek)~delta.t)

#############
## Cviceni 2
#############

matka <- c(40, 64, 34, 15, 57, 45)
dite <- c(33, 46, 23, 12, 56, 40)
regr.primka <- lm(dite~matka)
summary(regr.primka)
confint(regr.primka)

konf.intervaly <- function(x){
  a <- cbind(1, x)
  sd <- qt(0.975, df=4)*sqrt(diag(a%*%vcov(regr.primka)%*%t(a)))
  odhad <- a%*%coef(regr.primka)
  dolni <- odhad - sd
  horni <- odhad + sd
  return(data.frame(dolni, odhad, horni))
}

konf.pasy <- function(x){
  a <- cbind(1, x)
  sd <- sqrt(2*qf(0.95, df1=2, df2=4)*diag(a%*%vcov(regr.primka)%*%t(a)))
  odhad <- a%*%coef(regr.primka)
  dolni <- odhad - sd
  horni <- odhad + sd
  return(data.frame(dolni, odhad, horni))
}

konfidencni.intervaly <- konf.intervaly(xx<-seq(10, 70, length=101))
konfidencni.pasy <- konf.pasy(xx)

plot(dite~matka, main="Koncentrace kyseliny mlecne v krvi", ylim=c(0, 60), asp=1)
lines(fitted(regr.primka)~matka)
lines(konfidencni.intervaly$odhad~xx, col=4)
lines(konfidencni.intervaly$horni~xx, col=4, lty=2)
lines(konfidencni.intervaly$dolni~xx, col=4, lty=2)
lines(konfidencni.pasy$horni~xx, col=6, lty=2)
lines(konfidencni.pasy$dolni~xx, col=6, lty=2)

#############
## Cviceni 3
#############

rychlost <- c(4:10)
palivo <- c(6.1, 5.8, 6.0, 6.5, 6.8, 8.1, 10.0)
regr.parabola <- lm(palivo~rychlost+I(rychlost^2))
summary(regr.parabola)
confint(regr.primka)

regr.parabola.sub1 <- lm(palivo~1)
anova(regr.parabola.sub1, regr.parabola)

plot(palivo~rychlost, main="Trabant", ylab="spotreba [l/100 km]", xlab="rychlost [km/h]")
lines(fitted(regr.parabola)~rychlost)

regr.parabola.sub2 <- lm(palivo~I(rychlost^2))
summary(regr.parabola.sub2)
plot(palivo~rychlost, main="Trabant", ylab="spotreba [l/100 km]", xlab="rychlost [km/h]", ylim=c(0, 10), xlim=c(0, 10))
lines(fitted(regr.parabola)~rychlost)
lines(fitted(regr.parabola.sub2)~rychlost, col=2)
lines(xx<-seq(0, 10, length=501), 4.6+0.045*xx^2, col=2, lty=2)

#############
## Cviceni 4
#############

insulgas <- read.table("insulgas.txt", header=TRUE)
dim(insulgas)
summary(insulgas)

#############
## Cviceni 5
#############

# (a)
# ---
plot(insulgas$Gas~insulgas$Temp, main="Gas consumption", xlab="Outside temperature", ylab="Gas consumption")
points(insulgas$Gas[insulgas$Insulate=="After"]~insulgas$Temp[insulgas$Insulate=="After"], col="springgreen", pch=19)
lines(lowess(insulgas$Gas[insulgas$Insulate=="After"]~insulgas$Temp[insulgas$Insulate=="After"]), col="springgreen")
lines(lowess(insulgas$Gas[insulgas$Insulate=="Before"]~insulgas$Temp[insulgas$Insulate=="Before"]))

# (b)
# ---
model.basic <- lm(Gas~Temp+Insulate, data=insulgas)
summary(model.basic)

est <- coefficients(model.basic)
plot(insulgas$Gas~insulgas$Temp, main="Gas consumption", xlab="Outside temperature", ylab="Gas consumption")
points(insulgas$Gas[insulgas$Insulate=="After"]~insulgas$Temp[insulgas$Insulate=="After"], col="springgreen", pch=19)
abline(est[1], est[2], col="springgreen")
abline(est[1]+est[3], est[2])

plot(model.basic, which=c(1:6))

# (c)
# ---
model.inter <- lm(Gas~Temp*Insulate, data=insulgas)
summary(model.inter)

est <- coefficients(model.inter)
plot(insulgas$Gas~insulgas$Temp, main="Gas consumption", xlab="Outside temperature", ylab="Gas consumption")
points(insulgas$Gas[insulgas$Insulate=="After"]~insulgas$Temp[insulgas$Insulate=="After"], col="springgreen", pch=19)
abline(est[1], est[2], col="springgreen")
abline(est[1]+est[3], est[2]+est[4])

plot(model.inter, which=c(1:6))

# (d)
# ---
anova(model.basic, model.inter)