nflcomb <- read.csv("http://www.stat.ufl.edu/~winner/sta4210/mydata/nfl_combine.csv", header=TRUE) attach(nflcomb); names(nflcomb) #### Matrix form (using lm for |e|,y-hat regressions) ########### n <- length(Weight) X0 <- rep(1,n) X <- as.matrix(cbind(X0,Height,ArmLng,HandLng)) Y <- as.matrix(Weight) p <- ncol(X) #### Fit original regression, and regress |e| on Y-hat b.ols <- solve(t(X) %*% X) %*% t(X) %*% Y yhat.ols <- X %*% b.ols e.ols <- Y - yhat.ols abs.e.ols <- abs(e.ols) e.reg.ols <- lm(abs.e.ols ~ yhat.ols) summary(e.reg.ols) s.ols <- predict(e.reg.ols) w.ols <- 1/s.ols^2 b.old <- b.ols wm.old <- as.matrix(diag(w.ols)) b.diff <- 100 while (b.diff > 0.0001) { b.new <- solve(t(X) %*% wm.old %*% X) %*% t(X) %*% wm.old %*% Y yhat.new <- X %*% b.new abs.e.new <- abs(Y - yhat.new) wm.new <- as.matrix(diag(1/predict(lm(abs.e.new~yhat.new))^2)) b.diff <- sum((b.new-b.old)^2) b.old <- b.new wm.old <- wm.new } b.new <- solve(t(X) %*% wm.old %*% X) %*% t(X) %*% wm.old %*% Y b.wls <- b.new wm.wls <- wm.new mse.w <- (t(Y-X%*%b.wls) %*% wm.wls %*% (Y-X%*%b.wls))/(n-p) s2.b.wls <- mse.w[1,1]*solve(t(X) %*% wm.wls %*% X) s.b.wls <- sqrt(diag(s2.b.wls)) t.b.wls <- b.wls/s.b.wls print(round(cbind(b.wls,s.b.wls,t.b.wls),3))