Robust Regression - Weighting Methods

The stackloss data set

Brownlee’s Stack Loss Plant Data contains operational data of a plant for the oxidation of ammonia to nitric acid.

The variables are:

• Air.Flow Flow of cooling air
• Water.Temp Cooling Water Inlet Temperature
• Acid.Conc. Concentration of acid [per 1000, minus 500]
• stack.loss Stack loss

fit.sl = lm(stack.loss ~ ., data = stackloss)

#attach(stackloss)
#plot(Acid.Conc. , stack.loss, pch = 18, col="red")

par(mfrow=c(2,2))
plot(fit.sl)

par(mfrow=c(1,1))

---

Fitting a robust model (rlm()}

The rlm command in the MASS package command implements several versions of robust regression.

library(MASS)

summary(rlm(stack.loss ~ ., data = stackloss))

## 
## Call: rlm(formula = stack.loss ~ ., data = stackloss)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -8.91753 -1.73127  0.06187  1.54306  6.50163 
## 
## Coefficients:
##             Value    Std. Error t value 
## (Intercept) -41.0265   9.8073    -4.1832
## Air.Flow      0.8294   0.1112     7.4597
## Water.Temp    0.9261   0.3034     3.0524
## Acid.Conc.   -0.1278   0.1289    -0.9922
## 
## Residual standard error: 2.441 on 17 degrees of freedom

---

Using Other Psi Operators

• huber
• bisquare
• hampel

Fitting is done by iterated re-weighted least squares (IWLS).

Psi functions are supplied for the Huber, Hampel and Tukey bisquare proposals as psi.huber, psi.hampel and psi.bisquare. Huber’s corresponds to a convex optimization problem and gives a unique solution (up to collinearity). The other two will have multiple local minima, and a good starting point is desirable.

---

Huber

summary(rlm(stack.loss ~ ., data = stackloss))

## 
## Call: rlm(formula = stack.loss ~ ., data = stackloss)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -8.91753 -1.73127  0.06187  1.54306  6.50163 
## 
## Coefficients:
##             Value    Std. Error t value 
## (Intercept) -41.0265   9.8073    -4.1832
## Air.Flow      0.8294   0.1112     7.4597
## Water.Temp    0.9261   0.3034     3.0524
## Acid.Conc.   -0.1278   0.1289    -0.9922
## 
## Residual standard error: 2.441 on 17 degrees of freedom

summary(rlm(stack.loss ~ ., data = stackloss,  psi = psi.huber))

## 
## Call: rlm(formula = stack.loss ~ ., data = stackloss, psi = psi.huber)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -8.91753 -1.73127  0.06187  1.54306  6.50163 
## 
## Coefficients:
##             Value    Std. Error t value 
## (Intercept) -41.0265   9.8073    -4.1832
## Air.Flow      0.8294   0.1112     7.4597
## Water.Temp    0.9261   0.3034     3.0524
## Acid.Conc.   -0.1278   0.1289    -0.9922
## 
## Residual standard error: 2.441 on 17 degrees of freedom

---

Bisquare

 rlm(stack.loss ~ ., stackloss, psi = psi.bisquare)

## Call:
## rlm(formula = stack.loss ~ ., data = stackloss, psi = psi.bisquare)
## Converged in 11 iterations
## 
## Coefficients:
## (Intercept)    Air.Flow  Water.Temp  Acid.Conc. 
## -42.2852537   0.9275471   0.6507322  -0.1123310 
## 
## Degrees of freedom: 21 total; 17 residual
## Scale estimate: 2.28

---

Hampel

 rlm(stack.loss ~ ., stackloss, psi = psi.hampel, init = "lts")

## Call:
## rlm(formula = stack.loss ~ ., data = stackloss, psi = psi.hampel, 
##     init = "lts")
## Converged in 9 iterations
## 
## Coefficients:
## (Intercept)    Air.Flow  Water.Temp  Acid.Conc. 
## -40.4747819   0.7410853   1.2250731  -0.1455245 
## 
## Degrees of freedom: 21 total; 17 residual
## Scale estimate: 3.09