# http://www.karlin.mff.cuni.cz/~zvara/regrese/0708/prednaska.html

cuzn<-read.table("~/CuZn.dat",header=TRUE)
cuzn

summary(cuzn)
pairs(cuzn)

attach(cuzn)

hist(protein,col="grey50")
library(lattice)
histogram(~protein|copper+zinc)
boxplot(protein)

cor(copper,protein,method="pearson")
cor(cuzn)

summary(m.inter<-lm(protein~copper+zinc+I(copper*zinc)))
summary(m.lin<-lm(protein~copper+zinc))
summary(m.cu<-lm(protein~copper))
summary(m.zn<-lm(protein~zinc))

anova(m.inter)

boxplot(protein~zinc)
plot(protein~zinc)
abline(m.zn,col=2)

anova(m.zn,m.lin,m.inter)

interaction.plot(zinc,copper,protein)

Copper<-factor(copper)
summary(m.factor<-lm(protein~Copper+zinc),cor=T)

library(car)
Anova(lm(protein~copper*zinc))
Anova(lm(protein~copper*zinc),type="III")

summary(m.inter22<-lm(protein~copper+zinc+I(copper^2)+I(zinc^2)+I(zinc*copper)))
summary(m.inter12<-lm(protein~copper+zinc+I(zinc^2)+I(zinc*copper)))
summary(m.inter11<-lm(protein~copper+zinc+I(zinc*copper)))

anova(m.inter11,m.inter12,m.inter22)

summary(m.inter.infl<-influence.measures(m.inter))

2*(1-pt(abs(rstudent(m.inter)[21]),25-3-1))
2*(1-pt(max(abs(rstudent(m.inter11))),25-3-1))*25 # Bonferroni

plot(m.inter)
plot(m.inter,which=1:4) # Cook (default R <= 2.1.x)

vif(m.inter)
summary(lm(copper~zinc))
summary(lm(zinc~copper))

VIF <- function(lmobj) {
if (!is.null(weights(lmobj)))
    stop("requires unweighted model")
if (!(any(names(coefficients(lmobj))=="(Intercept)")))
    stop("requires model with intercept")
X0 <- scale(model.matrix(lmobj))[,-1]
nam <- labels(terms(lmobj))[-1]
y0 <- scale(lmobj$model[,1])
lmobj0 <- lm(y0~X0)
VIF <- diag(solve(cor(X0)))
tol <- 1/VIF; R2 <- 1-tol
b.star <- coef(lmobj0)[-1]
out <- cbind(b.star,VIF,R2,tol)
return(out)
}

VIF(m.inter)

plot(m.inter$resid~m.inter$fitted)
abline(h=0,col=4,lwd=2)

plot((m.inter$resid)^2~m.inter$fitted)

qqnorm(m.inter$resid)
qqline(m.inter$resid,col="red")
hist(m.inter$resid)

library("lmtest")
bptest(m.inter)

shapiro.test(resid(m.inter))

dwtest(m.inter)