##------------------------------------------------------------------------## ## Script for web appendix on robust regression ## ## An R and S-PLUS Companion to Applied Regression ## ## John Fox ## ## Sage Publications, 2002 ## ##------------------------------------------------------------------------## # preliminaries options(width=65) options(digits=5) # graphs of objective, psi, and weight functions ls.rho <- function(e) e^2 ls.psi <- function(e) 2*e ls.wt <- function(e) rep(1, length(e)) huber.rho <- function(e, k=1.345){ ifelse (abs(e) <= k, 0.5*e^2, k*abs(e) - 0.5*k^2) } huber.psi <- function(e, k=1.345){ ifelse (e > k, k, ifelse (abs(e) <= k, e, -k)) } huber.wt <- function(e, k=1.345){ ifelse (abs(e) <= k, 1, k/abs(e)) } biwt.rho <- function(e, k=4.685){ ifelse (abs(e) <= k, ((k^2)/6)*(1 - (1 - (e/k)^2)^3), (k^2)/6) } biwt.psi <- function(e, k=4.685){ ifelse (abs(e) <= k, e*(1 - (e/k)^2)^2, 0) } biwt.wt <- function(e, k=4.685){ ifelse (abs(e) <= k, (1 - (e/k)^2)^2, 0) } par(mfrow=c(3,3)) e <- seq(-6, 6, by=.1) plot(e, ls.rho(e), type='l', lwd=2, ylab=expression(rho[LS](e))) plot(e, ls.psi(e), type='l', lwd=2, ylab=expression(psi[LS](e)), main='Least Squares') plot(e, ls.wt(e), type='l', lwd=2, ylab=expression(w[LS](e)), ylim=c(0,1)) plot(e, huber.rho(e), type='l', lwd=2, ylab=expression(rho[H](e))) plot(e, huber.psi(e), type='l', lwd=2, ylab=expression(psi[H](e)), main='Huber') plot(e, huber.wt(e), type='l', lwd=2, ylab=expression(w[H](e)), ylim=c(0,1)) plot(e, biwt.rho(e), type='l', lwd=2, ylab=expression(rho[B](e))) plot(e, biwt.psi(e), type='l', lwd=2, ylab=expression(psi[B](e)), main='Bisquare') plot(e, biwt.wt(e), type='l', lwd=2, ylab=expression(w[B](e)), ylim=c(0,1)) # LS estimator library(car) # mostly for Duncan data set data(Duncan) mod.ls <- lm(prestige ~ income + education, data=Duncan) summary(mod.ls) which.names(c('minister', 'conductor'), Duncan) mod.ls.2 <- update(mod.ls, subset=-c(6,16)) summary(mod.ls.2) # M-estimation library(MASS) mod.huber <- rlm(prestige ~ income + education, data=Duncan) summary(mod.huber) plot(mod.huber\$w, ylab="Huber Weight") identify(1:45, mod.huber\$w, rownames(Duncan)) library(lqs) mod.bisq <- rlm(prestige ~ income + education, data=Duncan, method='MM') summary(mod.bisq, cor=F) plot(mod.bisq\$w, ylab="Bisquare Weight") identify(1:45, mod.bisq\$w, rownames(Duncan)) mod.hampel <- rlm(prestige ~ income + education, data=Duncan, psi=psi.hampel) summary(mod.hampel, cor=F) plot(mod.hampel\$w, ylab="Hampel Weight") identify(1:45, mod.hampel\$w, rownames(Duncan)) # Bounded-influence regression mod.lts <- ltsreg(prestige ~ income + education, data=Duncan) mod.lts