library(MPV) #load the textbook data sets library

Delivery<-p8.3

x1<-Delivery$x1
x2<-Delivery$x2
y<-Delivery$y
par(mfrow=c(1,2))
plot(x1,y,pch=19)
plot(x2,y,pch=19)
par(mfrow=c(1,1))
pairs(cbind(y,x1,x2),pch=19)
dev.print(device=postscript,file='DeliveryData.ps')

cor(cbind(y,x1,x2))

library(car)
scatter3d(y~x1+x2+x1*x2,data=Delivery)

#library(Rcmdr)
#scatter3d(x1,x2,y)

########################################

fit.Del12<-lm(y~x1+x2,data=Delivery)
summary(fit.Del12)
library(scatterplot3d) 
s3d <-scatterplot3d(x1,x2,y, pch=16, grid=FALSE,
	main="Plane of best fit",angle=20)
s3d$plane3d(fit.Del12,col='black')
dev.print(device=postscript,file='DeliveryData-Plane.ps')

s3d <- scatterplot3d(cbind(x1,x2,residuals(fit.Del12)),lab.z=2, 
               pch=16,scale.y=0.7,highlight.3d=TRUE, main="Residuals")
s3d$plane3d(0,0,0, "solid", col="grey")

anova(fit.Del12)

par(mfrow=c(2,2),mar=c(2, 4, 4, 2) + 0.1)
plot(x1,residuals(fit.Del12),pch=19,ylim=range(-10,10),ylab='Residuals')
abline(h=0,lty=2);title('Residuals vs x1')
plot(x2,residuals(fit.Del12),pch=19,ylim=range(-10,10),ylab='Residuals')
abline(h=0,lty=2);title('Residuals vs x2')

plot(fitted(fit.Del12),residuals(fit.Del12),pch=19,ylim=range(-10,10),ylab='Residuals',xlab=expression(hat(y)[i]))
abline(h=0,lty=2);title('Residuals vs fitted values')


dev.print(device=postscript,file='DeliveryData-Residuals.ps',height=9,width=11)


######################################





fit.Del21<-lm(y~x2+x1,data=Delivery)
summary(fit.Del21) #gives same regression results
anova(fit.Del21)  #gives different ANOVA table !

fit.Del1<-lm(y~x1,data=Delivery)
summary(fit.Del1)
anova(fit.Del1)
fit.Del2<-lm(y~x2,data=Delivery)
summary(fit.Del2)
anova(fit.Del2)



par(mfrow=c(1,1))



fit.Del12i<-lm(y~x1+x2+x1:x2,data=Delivery)
summary(fit.Del12i)