> n <- length(Y)
> 
> X2C <- X2-mean(X2)
> 
> library(MASS)
> 
> huber.rreg1 <- rlm(Y ~ X2C + I(X2C^2))
> 
> summary(huber.rreg1)

Call: rlm(formula = Y ~ X2C + I(X2C^2))
Residuals:
      Min        1Q    Median        3Q       Max 
-36.08978  -2.84240  -0.03667   3.00545   8.87016 

Coefficients:
            Value    Std. Error t value 
(Intercept) 259.4211   1.0912   237.7351
X2C           1.5646   0.1536    10.1888
I(X2C^2)      0.0802   0.0184     4.3490

Residual standard error: 4.338 on 37 degrees of freedom
> 
> install.packages("car")
Installing package(s) into ‘C:/Users/Larry/Documents/R/win-library/2.15’
(as ‘lib’ is unspecified)
Warning: package ‘car’ is in use and will not be installed
--- Please select a CRAN mirror for use in this session ---
Error in m[, 1L] : incorrect number of dimensions
> 
> library(car)
> 
> scatterplotMatrix(~Y+X1+X2+X3+X4+X5)
> 
> huber.rreg2 <- rlm(Y ~ X2 + X3 + X4)
> 
> summary(huber.rreg2)

Call: rlm(formula = Y ~ X2 + X3 + X4)
Residuals:
     Min       1Q   Median       3Q      Max 
-22.3984  -2.9999   0.5596   2.6767   9.4139 

Coefficients:
            Value    Std. Error t value 
(Intercept) 207.6806  17.6965    11.7357
X2            0.7972   0.1399     5.6982
X3            0.1609   0.2209     0.7282
X4           -1.1692   0.2231    -5.2412

Residual standard error: 4.342 on 36 degrees of freedom
> 
> ############### Begin brute-force IRWLS - Huber Method - Example 1, ppp. 441-5
> 
> w <- rep(1,n) 
> 
> irwls.1 <- lm(Y ~ X2C + I(X2C^2))
> 
> 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.000001) {
+ 
+ num.iter <- num.iter + 1
+ 
+ irwls.2 <- lm(Y ~ X2C + I(X2C^2),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
[1] 7
> 
> summary(irwls.2)

Call:
lm(formula = Y ~ X2C + I(X2C^2), weights = w)

Weighted Residuals:
     Min       1Q   Median       3Q      Max 
-14.5074  -2.8421  -0.0358   3.0051   7.1923 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 259.42103    1.17608 220.582  < 2e-16 ***
X2C           1.56493    0.15634  10.010 4.47e-12 ***
I(X2C^2)      0.08016    0.02193   3.655 0.000793 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 5.12 on 37 degrees of freedom
Multiple R-squared: 0.7306,     Adjusted R-squared: 0.7161 
F-statistic: 50.18 on 2 and 37 DF,  p-value: 2.892e-11 

> 
> ############### Begin brute-force IRWLS - Huber Method - Example 2, ppp. 445-8
> 
> w <- rep(1,n) 
> 
> irwls.1 <- lm(Y ~ X2 + X3 + 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.000001) {
+ 
+ num.iter <- num.iter + 1
+ 
+ irwls.2 <- lm(Y ~  X2 + X3 + 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
[1] 12
> 
> summary(irwls.2)

Call:
lm(formula = Y ~ X2 + X3 + X4, weights = w)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-11.631  -3.006   0.546   2.677   7.529 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 207.8462    17.5887  11.817 6.00e-14 ***
X2            0.7940     0.1408   5.638 2.12e-06 ***
X3            0.1636     0.2204   0.742    0.463    
X4           -1.1696     0.2189  -5.343 5.25e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 4.192 on 36 degrees of freedom
Multiple R-squared: 0.874,      Adjusted R-squared: 0.8635 
F-statistic: 83.27 on 3 and 36 DF,  p-value: 2.92e-16