# Handout for APLM :
# on Weighted Multiple Regression -- Blood Pressure Example on p427 of the text
#                                                     prepared by Yu-Ling Tseng  
#              (thanks to Dr. J.J. Tsai for the 1th version of this program) 

# data import 
ch11ta01<-matrix(scan("ch11ta01.txt"),ncol=2, byrow=T) ;
x1 <- ch11ta01[,1];
y <- ch11ta01[,2];

# regress BP on AGE #
dum <- data.frame("Age"=x1, "BP"=y)
dum
attach(dum)
fm <- lm(BP~Age, dum)    

# To get Figure 11.1
windows() 
opar <- par(mfrow = c(2,2), oma = c(0, 0, 2.7, 0));
plot(x1, y, xlab="AGE", ylab="BLOOD PRESSURE", main="Scatter Plot")
abline(fm)
plot(x1,resid(fm), xlab="AGE", ylab="Residul", main="Residual plot against X");
abline(0,0);

# Note: the plot exhibits a megaphone shape => noncontant var => try W.L.S.E. here
# To get the weights w_i's 
# Since the plot exhibits a megaphone shape, apply rule 1 on p.425 to 
#                           get estimtes of the variances:
# => Regress absolute residuals on AGE #

dum1 <- data.frame("Age"=x1, "ABSei"=abs(resid(fm)))
dum1
attach(dum1)
fm1 <- lm(ABSei~Age, dum1)   
plot(x1,abs(resid(fm)), xlab="AGE", ylab="Absolute Residul", main="Absolute Residual plot against X");
abline(fm1);

# To get Table11.1 on p427 of the text #
ch11ta01<-data.frame("x1"=Age,"y"=BP,"e_i"=resid(fm),"ABS_ei"=abs(resid(fm)),"hat{s_i}"=fitted(fm1),"w_i"=1/(fitted(fm1)^2) )
ch11ta01;    

# To get W.L.S.E. using R...  regress regress BP on AGE, weighted by w_i#
fm2 <- lm(BP~Age, dum,weights=1/(fitted(fm1)^2))  

# The equation (11.20) on p428 of the text #
fm2
summary(fm2)             
# you can get the estimate of the standard deviation of b_w1 from this table....