# 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)