#  M-L-R (1)                         Course: APLM
#                                   prepared by Yu-Ling Tseng

# PART I: some matrix operations that you might find useful 

x <- c(1:5)/5;   # say, we have these x values  #
x;

var(x);        # the sample variance of those x's #
x^2;
mean(x);       # x_bar #
sum((x-mean(x))^2); # sum of (x_i-x_bar)^2 i=1,2,..10 #
sum((x-mean(x))^2)/(length(x)-1);  # c.f. var(x) #

y<-x+rnorm(5) 
fm<-lm(y~x)

X <- model.matrix(fm)  # this gives the design matrix  in a S-L-R model, [1, x]#
X
t(X);    # transpose of X #
A <- solve(t(X)%*%X);  # (X^tX)^(-1) #
A;
H <- X%*%A%*%t(X);
H;    # the hat matrix #


# PART II: 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)
dum
attach(dum)
rm(x1,x2,y)

fm <- lm(SALES~TARGTPOP+DISPOINC, dum)   # regress SALES on those two predictors #

# To get (b) Basic Data  on p237 of the text #
fi65b<-data.frame("x1"=TARGTPOP,"x2"=DISPOINC,"y"=SALES,"y_hat"=fitted(fm),"residual"=resid(fm))
fi65b;    

# 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?

# To get (X^tX)^(-1) #
X<-model.matrix(fm)  #the design matrix
solve(t(X)%*%X)

# To get Figure 6.4, p.232 #
pairs(dum, main="Dwaine DATA:Scatter Plot Matrix",lwd=3);   # (a) Scatter plot matrix
cor(dum);              # (b) Correlation matrix


# Quit R now?  #

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