#------------------------------
#   Zagat's
#------------------------------

zagat = read.csv("http://faculty.chicagobooth.edu/nicholas.polson/teaching/41000/zagat.csv",header = TRUE)
attach(zagat)

# Data summary

summary(zagat)
cor(zagat)

# Plot price versus each of the quality variables
setEPS()
postscript("zagat.eps",width=7,height=3)
par(mfrow=c(1,3))
plot(price~food,main="Price versus Food",pch=20,col=20,bty='n')
plot(price~decor,main="Price versus Decor",pch=20,col=20,bty='n')
plot(price~service,main="Price versus Service",pch=20,col=20,bty='n')
dev.off()
# Multiple regression

zaga = lm(price ~ food + decor + service, zagat)

# Extract coefficients and model summary

coef(zaga)
summary(zaga)

# Extract fitted value for each case

fitted(zaga)

# 4-in-1 residual diagnostics

layout(matrix(c(1,2,3,4),2,2))   
plot(zaga,pch=20)

# Prediction for the case (food = 27, service = 25, decor = 20)

newdata = data.frame(27,25,20)
colnames(newdata) = c("food","service","decor")
predict(zaga,newdata)
# sum(coef(zaga)*c(1,27,20,25))

# Plug-in prediction interval (95%)

# Residual standard error 

se = summary(zaga)$sigma

# Upper bound
sum(coef(zaga)*c(1,27,20,25))+qnorm(0.975)*se

# Lower bound
sum(coef(zaga)*c(1,27,20,25))-qnorm(0.975)*se

detach(zagat)