# Viscosity, as a function of filler (f) and oil levels (p)

V <- c(26,38,50,76,108,157,17,26,37,53,83,124,13,20,27,37,57,87,NA,15,22,27,41,63)
f <- rep(c(0,12,24,36,48,60),4)
p <- rep(c(0,10,20,30),rep(6,4))

visc <- data.frame(V, f=f, p = p)

pdf("visc.pdf")
plot(visc)
dev.off()

pdf("lm1.pdf")
par(mfrow=c(2,2))
plot(lm(V ~ f + p), ask=F)
dev.off()
# add quadtractic terms
pdf("lm2.pdf")
par(mfrow=c(2,2))
plot(lm(V ~ p + f + I(p*f) + I(p^2) + I(f^2)),ask=F)
dev.off()

# explore transformtions of Y
library(MASS)  # load the MASS library from Venables and Ripley
pdf("boxcox.pdf")
boxcox(lm(V ~ f + p), plotit=T)
dev.off()

lnvisc <- data.frame(lnV =log(V), f=f, p=p)
pdf("logdata.pdf")
plot(lnvisc)
dev.off()
pdf("lmlog.pdf")
par(mfrow=c(2,2))
plot(lm(log(V) ~ f + p), ask=F)
dev.off()
lmlog1 <- lm(log(V) ~ f + p)
summary(lmlog1)

lmlog2 <- lm(log(V) ~ f + p + I(f*p)  + I(f^2) + I(p^2), data=visc)
summary(lmlog2)
anova(lmlog1,lmlog2)