# pdf("F:\\sta4210\\CH11TA01.pdf") bpdata <- read.table("http://www.stat.ufl.edu/~winner/sta4210/data/CH11TA01.txt", header=F, col.names=c("age","dbp")) attach(bpdata) #### Fit Ordinary least squares, plot and regressions of dbp, e, |e|, e^2 vs age bp.mod1 <- lm(dbp ~ age) summary(bp.mod1) e <- residuals(bp.mod1) abs.e <- abs(e) e2 <- e^2 abs_e.mod1 <- lm(abs.e ~ age) # Regression of |e| on age summary(abs_e.mod1) e2.mod1 <- lm(e2 ~ age) # Regression of e^2 on age summary(e2.mod1) par(mfrow=c(2,2)) # Place plots on one page plot(age,dbp,main="scatter plot",xlab="age", ylab="dbp") abline(bp.mod1) plot(age,e,main="Residual Plot",xlab="age",ylab="residual") abline(h=0) plot(age,abs.e,main="Absolute Residual Plot",xlab="age",ylab="abs(residual)") abline(abs_e.mod1) plot(age,e2,main="Squared Residual Plot",xlab="age",ylab="residual^2") abline(e2.mod1) s.hat <- predict(abs_e.mod1) # Save the fitted SD's for W w.hat <- 1/s.hat^2 # Obtain the estimated weights in W print(cbind(age,dbp,e,abs.e,s.hat,w.hat)) X <- as.matrix(cbind(rep(1,length(dbp)),age)) # X-matrix Y <- as.matrix(dbp) # Y-vector W <- as.matrix(diag(w.hat)) # W-hat matrix b_w <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% Y # obtain b_w sse_w <- t(Y - (X %*% b_w)) %*% W %*% (Y - (X %*% b_w)) # Obtain SSE_w mse_w <- sse_w/(nrow(X)-ncol(X)) # Obtain MSE_w se_bw <- sqrt(diag(mse_w[1,1] * solve(t(X) %*% W %*% X))) # Obtain SE{b_w} t_bw <- b_w/se_bw # t-statistics p_bw <- 2*(1-pt(abs(t_bw),nrow(X)-ncol(X))) # P-values print(cbind(b_w,se_bw,t_bw,p_bw)) ### Fit weighted least squares using lm function and w.hat as weight bp.modwls <- lm(dbp ~ age, weight=w.hat) summary(bp.modwls) # dev.off()