#------------------------------# #--- Robust Regression in R ---# #------------------------"------# # The robustbase package contains a demonstration dataset called hbk # install.packages("robustbase") library(robustbase) data(hbk) plot(hbk) # using OLS: model.lm = lm(Y~., data = hbk) summary(model.lm) confint(model.lm, level=0.95) hbk.x <- data.matrix(hbk[,1:3]) library(MASS) cat("Measurements of location:") cat("mean(x1) = ", mean(hbk[,1])) cat("median(x1) = ", median(hbk[,1])) cat("20% trimmed mean(x1) = ", mean(hbk[,1], trim = 0.2)) #20% off both ends cat("covMcd location(x1) = ", covMcd(hbk.x)$center[1]) cat("cov.rob location(x1) = ", cov.rob(hbk.x)$center[1]) cat("Measurements of scale:") cat("StdDev(x1) = ", sd(hbk[,1])) cat("MAD(x1) = ", mad(hbk[,1])) # The default constant = 1.4826 cat("0.7413*IQR(x1) = ", 0.7413*IQR(hbk[,1])) cat("covMcd scale(x1) = ", sqrt(covMcd(hbk.x)$cov[1,1])) cat("cov.rob scale(x1) = ", sqrt(cov.rob(hbk.x)$cov[1,1])) bodyfat <- read.csv('bodyfat.csv') model.lm = lm(BodyFat~Abdomen, data = bodyfat) summary(model.lm) model.robust1 = rlm(BodyFat~Abdomen, data = bodyfat) summary(model.robust1) bodyfat2 <- read.csv('bodyfat-reduced.csv') model.lm2 = lm(BodyFat~Abdomen, data = bodyfat2) summary(model.lm2) # robust regression library(MASS) model.robust1 = rlm(Y~., data = hbk) summary(model.robust1) # Fitting is done by iterated re-weighted least squares # method = "M" gives M-estimation (default), method = "MM" gives M-estimation with Tukey's biweight initialized by a specific S-estimator. # psi = psi.huber (default), psi.hampel, or psi.bisquare. # scale.est = method of scale estimation; default = re-scaled MAD of the residuals library(robustbase) model.robust2 = lmrob(Y~., data = hbk, method = "MM") summary(model.robust2) # default method = "MM" gives M-estimation with Tukey's biweight initialized by a specific S-estimator. # install.packages("robust") library(robust) model.robust3 = lmRob(Y~., data = hbk) summary(model.robust3)