##------------------------------------------------------------------------##
## 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