library(MASS) attach(phones) par(mfrow=c(1,1)) plot(year,calls) mtext("Belgium Phone Calls",outer=T,side=3) # # LEAST SQUARES # cat("least squares\n") phones.lm <- lm(calls ~ year) par(mfrow=c(2,2)) print(summary(phones.lm, correlation = F)) plot(year, calls) abline(phones.lm, lty = 1) plot(phones.lm,which=c(1:2,4)) # # ROBST REGRESSION # par(mfrow=c(1,1)) plot(year,calls) mtext("Belgium Phone Calls",outer=T,3) abline(phones.lm, lty = 1) # cat("\ntrimmed means\n") # rrtal.lm <- rreg(year, calls, method = wt.talworth) # print(rrtal.lm$coef) # abline(rrtal.lm, lty = 2) cat("\nHuber\n") rrhub.lm <- rlm(calls~year, psi = psi.huber,maxit=50) print(rrhub.lm$coef) abline(rrhub.lm, lty = 2) cat("\nHampel\n") rrham.lm <- rlm(calls~year, psi=psi.hampel) print(rrham.lm$coef) abline(rrham.lm, lty = 3) cat("\nTukey\n") rrtuk.lm <- rlm(calls~year, psi=psi.bisquare) print(rrtuk.lm$coef) abline(rrtuk.lm, lty = 4) # cat("\nleast absolute deviations\n") # l1.lm <- l1fit(year, calls) # print(l1.lm$coef) # abline(l1.lm, lty = 6) cat("\nleast trimmed squares\n") lts.lm <- lqs(calls~year,method="lts") print(lts.lm$coef) abline(lts.lm, lty = 5) cat("\nleast median of squares\n") lms.lm <- lqs(calls~year,method="lms") print(lms.lm$coef) abline(lms.lm,lty=6) legend(50, 200, legend = c("OLS", "Huber", "Hampel", "Tukey", "least trimmed squares","least median squares"), lty = 1: 6)