# MLR(2)    for APLM                   by Yu-Ling Tseng                                 
# Multiple regression -- Dwaine Studios Example     
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ch6fi05<-matrix(scan("ch06fi05.txt"),ncol=3, byrow=T) ;
x1 <- ch6fi05[,1];
x2 <- ch6fi05[,2];
y <- ch6fi05[,3];
dum <- data.frame("SALES"=y,"TARGTPOP"=x1,"DISPOINC"=x2)
attach(dum)
rm(x1,x2,y)
fm <- lm(SALES~TARGTPOP+DISPOINC, dum)   # regress SALES on those two predictors #
# To get (a) Multiple Regression Output  on p.237 #
summary(fm)             
anova(fm)        

# Note that the ANOVA table you get from R is a little different from that
# on p.237.  Think: How to get the ANOVA table on p.237 from R's output?
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Load Mass package to construct confidence intervals for
# beta's

confint(fm,level=.95)
confint(fm,level=.9)  # c.f. p. 245 joint conf. int's for beta1 and beta2

# Estimation of mean response, its confidence intervals
# and prediction intervals   
summary(dum)
new<- data.frame(TARGTPOP=c(65.4,53.1),DISPOINC=c(17.6, 17.7))
new
predict(lm(SALES~TARGTPOP+DISPOINC), new, se.fit = TRUE)
pred.w.plim <- predict(lm(SALES~TARGTPOP+DISPOINC), new, interval="prediction",level=.95)
pred.w.clim <- predict(lm(SALES~TARGTPOP+DISPOINC), new, interval="confidence",level=.95)
pred.w.plim 
                #c.f. p. 247 joint 90% prediction intervals for City A and B; notice that
                # level is set to be .95 in the command to get this simultaneous .90 prediction
                # intervals.

pred.w.clim     # p. 246
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Note: ANOVA tables in R
# anova(fittedmodel) of R gives the EXTRA SS's, c.f. chapter 7.  
# ex. anova(fm) v.s. anova(fmDT #
fmDT <- lm(SALES~DISPOINC+TARGTPOP, dum)   # regress SALES on those two predictors 
                                           # w/ different input order #
anova(fm)  
anova(fmDT)  
# but, same outputs with summary(fittedmodel) #
summary(fm)             
summary(fmDT)
 

# Quit R now?  #

# If so......    As usual, before you leave R ......    #
detach()
rm(list=ls())