################### ESIMERKKI 4.1 #################################
 hooker <- read.table("http://users.jyu.fi/~junyblom/hooker.dat",
 header=TRUE)

 forbes <- read.table("http://users.jyu.fi/~junyblom/forbes.dat",
 header=TRUE)

 f <- cbind(forbes,0)
 names(f)[3] <- "Group"
 h <- cbind(hooker,1)
 names(h)[3] <- "Group"
 fh <- rbind(f,h)
 attach(fh)
 pres <- 24.5*Pressure   # mmHg
 tmp <- 5*(Temp - 32)/9  # Celsius

 lm.add <- lm(pres ~ tmp)

 oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,4,0,0)) 

 plot(tmp, pres, xlab="", ylab="")

 mtext("Kiehumispiste (Celsius)", font=1, side=1, line=2.5) 
 mtext("Ilmanpaine (mmHg)", font=1, side=2, line=2.5)

 par(mar=c(3,4,0,0))

 plot(fitted(lm.add), resid(lm.add), xlab="", ylab="")

 mtext("Sovite", font=1, side=1, line=2.5) 
 mtext("Jäännös", font=1, side=2, line=2.5)

 par(oldpar)

 lm.log <- lm(log(pres) ~ tmp)

 oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,4,0,0)) 

 plot(tmp, log(pres), xlab="", ylab="")

 mtext("Kiehumispiste (Celsius)", font=1, side=1, line=2.5) 
 mtext("Ilmanpaineen logaritmi", font=1, side=2, line=2.5)

 par(mar=c(3,4,0,0))

 plot(fitted(lm.log), resid(lm.log), xlab="", ylab="")

 mtext("Sovite", font=1, side=1, line=2.5) 
 mtext("Jäännös", font=1, side=2, line=2.5)

 par(oldpar)

  K <- 1/(255.37 + 5*Temp/9) # Kelvin

  lm.kel <- lm(log(pres) ~ K)
  round(cbind(fitted(lm.log), fitted(lm.kel)), 2) 


##################### ESIMERKKI 4.2 #################################

 mesquite <- read.table("w:/mesquite.dat", header=TRUE)
 attach(mesquite)

 lm.0 <- lm(leafwt ~ diam1 + diam2 + canht +  totht +  dens + group)

 lm.1 <- lm(log(leafwt) ~ log(diam1) + log(diam2) + log(totht) + 
 log(canht) + log(dens) + group)

 oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,4,0,0)) 

 plot(fitted(lm.0), resid(lm.0), xlab="", ylab="")

 mtext("Sovite", font=1, side=1, line=2.5) 
 mtext("Jäännös", font=1, side=2, line=2.5)

 par(mar=c(3,4,0,0))

 plot(fitted(lm.1), resid(lm.1), xlab="", ylab="")

 mtext("Sovite", font=1, side=1, line=2.5) 
 mtext("Jäännös", font=1, side=2, line=2.5)

 par(oldpar)


##################### ESIMERKKI 4.3 ##################################

 library(alr3)
 data(transact)
 lm.1 <- lm(Time ~T1 + T2, data=transact)
 res <- resid(lm.1)
 qqnorm(res, main=" ", xlab="Teoreettiset kvantiilit", 
 ylab="Järjestetyt jäännökset"); qqline(res)

 yX <- as.matrix(transact[,c("Time","T1","T2")])
 boot.coef <- matrix(0,ncol=3,nrow=9999)
 colnames(boot.coef) <- c("(Intercept)","T1","T2")
 for (i in 1:9999){
 s <- sample(1:nrow(yX), size=nrow(yX), replace=TRUE)  ## Bootstrap-otosideksit
 y.lm <- lm(yX[s,1]~yX[s,2:3])
 boot.coef[i,] <- coef(y.lm)
                 }
 summary(lm.1)$coef[,"Std. Error"] # tavalliset keskivirheet
 sd(boot.coef)                     # Bootstrap-keskivirheet

##################### ESIMERKKI 4.4 ####################################

 hooker <- read.table("http://users.jyu.fi/~junyblom/hooker.dat",
 header=TRUE)

 forbes <- read.table("http://users.jyu.fi/~junyblom/forbes.dat",
 header=TRUE)

 f <- cbind(forbes,0)
 names(f)[3] <- "Group"
 h <- cbind(hooker,1)
 names(h)[3] <- "Group"
 fh <- rbind(f,h)
 attach(fh)
 pres <- 24.5*Pressure   # mmHg
 tmp <- 5*(Temp - 32)/9  # Celsius

 lm.log <- lm(log(pres) ~ tmp)
 which.max(resid(lm.log))

 48*(1-pt(rstudent(lm.log)[12], df=lm.log$df-1))

 lm.no.12 <- lm(log(pres) ~ tmp, subset = (12 != 1:48))
 coef(lm.no.12)
 coef(lm.log)

##################### ESIMERKKI 4.5 ######################################

 mesquite <- read.table("http://users.jyu.fi/~junyblom/mesquite.dat", 
 header=TRUE)
 attach(mesquite)

 lm.1 <- lm(log(leafwt) ~ log(diam1) + log(diam2) + log(totht) + 
 log(canht) + log(dens) + group)

 plot(hatvalues(lm.1), xlab="Tapaus", ylab="Vipuarvo")

 lm.no.3 <- lm(log(leafwt) ~ log(diam1) + log(diam2) + log(totht) + 
 log(canht) + log(dens) + group, subset = (3 != 1:45))


 round(rbind(coef(lm.1),coef(lm.no.3)), 3)

 plot(hatvalues(lm.B), xlab="Tapaus", ylab="Vipuarvo")
 plot(hatvalues(lm.C), xlab="Tapaus", ylab="Vipuarvo")

 plot(cooks.distance(lm.1))

 source("http://users.jyu.fi/~junyblom/VIF.R")
 VIF(lm.1)
 VIF(lm.B)