# 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(seq(0, 30, by=10),each=6)
visc <- data.frame(V, f, p)
summary(visc)

plot(visc)
library(lattice)
xyplot(V ~ f | as.factor(p))

library(MASS)  # load the MASS library from Venables and Ripley
boxcox(lm(V ~ f + p), plotit=T)

lnvisc <- data.frame(lnV =log(V), f=f, p=p)
plot(lnvisc)
xyplot(log(V) ~ f | as.factor(p))
par(mfrow=c(2,2)); plot(lm(log(V) ~ f + p), ask=F)
logV.lm = lm(log(V) ~ f + p)

summary(lm(log(V) ~ f + p))
par(mfrow=c(2,2)); plot(lm(V ~ p + f + I(p*f) + I(p^2) + I(f^2)),ask=F)
boxcox(lm(V ~ p + f + I(p*f) + I(p^2) + I(f^2)))
summary(lm(log(V) ~ p + f + I(p*f) + I(p^2) + I(f^2)))k