# This program produces Figure 6.5 on Page 241
# of our Text.

   t<-matrix(scan("CH06FI05.DAT"),ncol=3,byrow=T)

# The data has 3 columns and read by row

dwaine <-data.frame(TARGTPOP=t[,1],DISPOINC=t[,2],SALES=t[,3])
TXD <- t[,1]*t[,2]
attach(dwaine)

opar <- par(no.readonly = TRUE) # put old par setting
par(mfrow=c(2,2));
plot(TARGTPOP,SALES, main="TARGTPOP vs. SALES")
plot(DISPOINC,SALES, main="DISPOINC vs. SALES")
plot(TARGTPOP,DISPOINC, main="TARGTPOP vs. DISPOINC")
plot(TXD,SALES, main="TXD vs. SALES")
par(opar);
pairs(dwaine,main="Dwaine DATA") # Prelim graphical displays in pairs.

g<- lm(SALES~TARGTPOP*DISPOINC)
g2<- lm(SALES~TARGTPOP+DISPOINC)
gt<- lm(SALES~TARGTPOP)
gd<- lm(SALES~DISPOINC)

# Summary and ANOVA
g2s<-summary(g2);g2s;anova(g2);g2s$cov.unscaled
par(mfrow=c(2,2));
plot.lm(g2)     # Checking residuals

summary(gd);anova(gd)
summary(gt);anova(gt)

names(g)

# Load Mass package to construct confidence intervals for
# beta's

confint(g)
confint(g2)
confint(gt)
confint(gd)

# Estimation of mean response, its confidence intervals
# and prediction intervals
summary(dwaine)
new<- data.frame(TARGTPOP=seq(40,90,5))
predict(lm(SALES~TARGTPOP), new, se.fit = TRUE)
pred.w.plim <- predict(lm(SALES~TARGTPOP), new, interval="prediction")
pred.w.clim <- predict(lm(SALES~TARGTPOP), new, interval="confidence")

pred.w.plim # prediction intervals for TARGTPOP at seq(40,90,5)
pred.w.clim # confidence intervals for TARGTPOP at seq(40,90,5)

# Overlay CI and PI at the given TARGTPOP values
matplot(new$TARGTPOP,cbind(pred.w.clim, pred.w.plim[,-1]),
        lty=c(1,2,2,3,3), type="l", ylab="predicted y")