wheat <- read.csv("http://www.stat.ufl.edu/~winner/sta4210/mydata/ohio_wheat.csv", header=TRUE) attach(wheat); names(wheat) wheat.mod1 <- lm(whtcnd.y ~ tempon.x1 + rains.x2 + rainon.x3 + sunon.x4) summary(wheat.mod1) library(quantreg) wheat.lad <- rq(whtcnd.y ~ tempon.x1 + rains.x2 + rainon.x3 + sunon.x4) summary(wheat.lad) library(MASS) huber.rreg1 <- rlm(whtcnd.y ~ tempon.x1 + rains.x2 + rainon.x3 + sunon.x4) summary(huber.rreg1) ############### Begin brute-force IRWLS - Huber Method n <- length(whtcnd.y) w <- rep(1,n) irwls.1 <- lm(whtcnd.y ~ tempon.x1 + rains.x2 + rainon.x3 + sunon.x4) res1 <- residuals(irwls.1) b.old <- coef(irwls.1) MAD <- (1/0.6745)*median(abs(res1-median(res1))) u <- res1/MAD for(i in 1:n) w[i] <- min(1,1.345/abs(u[i])) delta_b <- 100.0 num.iter=0 while (delta_b > 0.000000001) { num.iter <- num.iter + 1 irwls.2 <- lm(whtcnd.y ~ tempon.x1 + rains.x2 + rainon.x3 + sunon.x4,weights=w) res2 <- residuals(irwls.2) b.new <- coef(irwls.2) MAD <- (1/0.6745)*median(abs(res2-median(res2))) u <- res2/MAD for(i in 1:n) w[i] <- min(1,1.345/abs(u[i])) delta_b <- sum((b.new-b.old)^2) b.old <- b.new } num.iter summary(irwls.2)