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

fev <- read.table("fev.txt", header=TRUE)
summary(fev)
fev$Height <- fev$Height*2.54

fev.boys <- fev[fev$Sex=="Male" & fev$Smoker=="Non", ]
dim(fev.boys)

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

# (a)
# ---
plot(fev.boys$FEV~fev.boys$Height, xlab="Height[cm]", ylab="FEV[l]", main="Dependence of FEV on Height for non-smoking boys")
lines(lowess(fev.boys$FEV~fev.boys$Height))

# (b)
# ---
# postupne pridavani:
model.boys <- lm(FEV~Height, data=fev.boys)
summary(model.boys)
model.boys <- lm(FEV~Height+I(Height^2), data=fev.boys)
summary(model.boys)
model.boys <- lm(FEV~Height+I(Height^2)+I(Height^3), data=fev.boys)
summary(model.boys)
model.boys <- lm(FEV~Height+I(Height^2), data=fev.boys)

plot(fev.boys$FEV~fev.boys$Height, xlab="Height[cm]", ylab="FEV[l]", main="Dependence of FEV on Height for non-smoking boys")
lines(lowess(fev.boys$FEV~fev.boys$Height))
lines(xx <- seq(min(fev.boys$Height), max(fev.boys$Height), length=101), predict(model.boys, newdata=data.frame(Height=xx)), col=4)

# postupne odebirani:
model.boys <- lm(FEV~Height+I(Height^2)+I(Height^3)+I(Height^4)+I(Height^5), data=fev.boys)
summary(model.boys)
model.boys <- lm(FEV~Height+I(Height^2)+I(Height^3)+I(Height^4), data=fev.boys)
summary(model.boys)

lines(xx, predict(model.boys, newdata=data.frame(Height=xx)), col=6)

# (c)
# ---
# postupne pridavani:
model.boys <- lm(FEV~Height, data=fev.boys)
summary(model.boys)
model.boys <- lm(FEV~poly(Height, 2), data=fev.boys)
summary(model.boys)
model.boys <- lm(FEV~poly(Height, 3), data=fev.boys)
summary(model.boys)
model.boys <- lm(FEV~poly(Height, 2), data=fev.boys)

lines(xx, predict(model.boys, newdata=data.frame(Height=xx)), col=2)
lines(xx, predict(model.boys, newdata=data.frame(Height=xx)), col=4)

# postupne odebirani:
model.boys <- lm(FEV~poly(Height, 5), data=fev.boys)
summary(model.boys)
model.boys <- lm(FEV~poly(Height, 4), data=fev.boys)
summary(model.boys)

lines(xx, predict(model.boys, newdata=data.frame(Height=xx)), col=2)
lines(xx, predict(model.boys, newdata=data.frame(Height=xx)), col=6)

# (d)
# ---
model.boys <- lm(FEV~I(Height-mean(Height))+I((Height-mean(Height))^2), data=fev.boys)
summary(model.boys)

model.boys <- lm(FEV~poly(Height, 4), data=fev.boys)
summary(model.boys)

model.boys <- lm(FEV~poly((Height-mean(Height)), 4), data=fev.boys)
summary(model.boys)

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

# (a)
# ---
model.boys.4 <- lm(FEV~poly(Height, 4), data=fev.boys)
plot(model.boys.4, which=c(1:6))

# (b)
# ---
plot(sqrt(abs(rstandard(model.boys.4)))~fev.boys$Height)
lines(lowess(sqrt(abs(rstandard(model.boys.4)))~fev.boys$Height))

# (c)
# ---
wts <- 1/fev.boys$Height 

S <- diag(sqrt(wts))
Y <- S%*%fev.boys$FEV
X <- S%*%model.matrix(model.boys.4)
model.boys.4.wls <- lm(Y ~ -1 + X[, 1] + X[, 2] + X[, 3] + X[, 4] + X[, 5])
summary(model.boys.4.wls)
summary(model.boys.4.w)

model.boys.4.w <- lm(FEV~poly(Height, 4), data=fev.boys, weights=wts)
summary(model.boys.4.wls)
summary(model.boys.4.w)

plot(fev.boys$FEV~fev.boys$Height, xlab="Height[cm]", ylab="FEV[l]", main="Dependence of FEV on Height for non-smoking boys")
lines(lowess(fev.boys$FEV~fev.boys$Height))
lines(xx, predict(model.boys.4, newdata=data.frame(Height=xx)), col=4)
lines(xx, predict(model.boys.4.w, newdata=data.frame(Height=xx)), col=6)

# (d)
# ---
plot(model.boys.4, which=c(1:6))
plot(model.boys.4.w, which=c(1:6))
plot(model.boys.4.wls, which=c(1:6))

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

spills <- read.table("spills.txt", header=TRUE)
dim(spills)
spills
plot(spills$Spills~spills$Year, xlab="Year", ylab="Number of spills")

# (a)
# ---

spills.lm <- lm(Spills~I(Year-1973), data=spills)
summary(spills.lm)
lines(spills$Year, predict(spills.lm, newdata=data.frame(Year=spills$Year)))

plot(residuals(spills.lm)[-13]~residuals(spills.lm)[-1])
cor(residuals(spills.lm)[-13], residuals(spills.lm)[-1])

library(nlme)
spills.glm <- gls(Spills~I(Year-1973), data=spills, correlation = corAR1(form = ~Year))
summary(spills.glm)
summary(spills.lm)

plot(spills$Spills~spills$Year, xlab="Year", ylab="Number of spills")
lines(spills$Year, predict(spills.lm, newdata=data.frame(Year=spills$Year)))
lines(spills$Year, predict(spills.glm, newdata=data.frame(Year=spills$Year)), col=4)

intervals(spills.glm, which="var-cov")