Robust linear regression with bootstrapped standard errors

rlm_with_bootstrap.R

library(car)
library(MASS)

# generate data with slight non-constant variance
x1 <- gl(n = 3, k = 400, labels = c("A","B","C"))
x2 <- gl(n = 2, k = 600, labels = c("1","2"))
set.seed(1)
y <- 1 + 1.2*(x1 == "B") + 1.3*(x1 == "C") -0.5*(x2 == "2") + 
  c(rnorm(600,sd = 0.5),rnorm(600,sd = 0.7))
DF <- data.frame(y, x1, x2)


# Fit 2-way ANOVA, no interaction -----------------------------------------


aov1 <- aov(y ~ x1 + x2, data = DF)
summary(aov1)


# fit 2-way ANOVA via lm, same thing --------------------------------------

mod1 <- lm(y ~ x1 + x2, data = DF)
# ANOVA table
anova(mod1)
# coefficient table
summary(mod1)
# CI of coefficient estimats
confint(mod1)


# robust regression -------------------------------------------------------

# robust regression: downweight influential, outlying observations
# load MASS for rlm()
mod2 <- rlm(y ~ x1 + x2, data = DF)
summary(mod2)

# Now use bootstrap to evaluate precision of estimators and get CI for
# coefficients

# load car for Boot()
betahat.boot <- Boot(mod2, R=999) 

# bootstrapped CI of coefficients
confint(betahat.boot)
# compare to original CI
confint(mod1)

# not much difference in results


# diagnostic plots --------------------------------------------------------

plot(mod1, which = 1, main = "Regular ANOVA model")
plot(mod2, which = 1, main = "Robust method")