## Slide 3
require(MASS)
plot(phones$year,phones$calls,pch=19,xlab="Year",ylab="Millions of calls")
fit <- lm(calls~year,phones)
abline(fit$coef,lwd=3)
fit <- rlm(calls~year,phones,maxit=50)
abline(fit$coef,lwd=3,col="red")
fit <- rlm(calls~year,phones,psi=psi.bisquare)
abline(fit$coef,lwd=3,col="blue")
legend("topleft",legend=c("Least squares","Huber","Tukey"),lwd=3,col=c("black","red","blue"))

## Slide 6
u <- seq(-6,6,len=101)
plot(u,u^2,lwd=3,type="l",col="black",ylim=c(0,6),xlab="Residual",ylab="Loss")
lines(u,abs(u),lwd=3,col="green")
y <- numeric(length(u))
p <- function(x){psi.huber(x)*x}
for (i in 1:101) y[i] <- integrate(p,u[1],u[i])$value
y <- y-min(y)
lines(u,y,type="l",lwd=3,col="red")
p <- function(x){psi.bisquare(x)*x}
for (i in 1:101) y[i] <- integrate(p,-3.1,u[i])$value
y <- y-min(y)
lines(u,y,type="l",lwd=3,col="blue")
legend("bottomleft",col=c("black","green","red","blue"),legend=c("Least squares","Absolute value","Huber","Tukey"),lty=1,lwd=3)

## Slide 16-17
fit <- lm(time~dist+climb,hills)
plot(fit$fitted,fit$resid,pch=19,xlab="Fitted value",ylab="Residual")
plot(fit$fitted,fit$resid,pch=19,cex=10*influence(fit)$hat,xlab="Fitted value",ylab="Residual")
summary(fit)
source("http://web.as.uky.edu/statistics/users/pbreheny/764-F11/notes/fun.R")
fit <- rlm(time~dist+climb,hills)
summary(fit)
fit <- rlm(time~dist+climb,hills,psi=psi.bisquare)
summary(fit)