Robust Regression

Robust regression in R

simulate 20 observations from a linear model with errors that follow a normal distribution

set.seed=112358
nobs<-20
sdy<-1
x<-seq(0,1,length=nobs)
y<-10+5*x+rnorm(nobs,sd=sdy)

add outlier at high leverage point

y[nobs]<-7

fit linear model

ols<-lm(y~x)

fit robust linear model

library(MASS)
mEst<-rlm(y~x)

plot results

plot(x,y)
abline(ols,lwd=2)

abline(mEst,col="red",lwd=2)
legend("topleft",legend=c("OLS","M-estimation"),lwd=2,col=1:2)

round(mEst$w,3)

 [1] 1.000 1.000 1.000 1.000 1.000 1.000
 [7] 1.000 1.000 0.954 1.000 1.000 1.000
[13] 0.962 1.000 1.000 1.000 1.000 1.000
[19] 1.000 0.192

Implement it yourself

start from ols fit

lmMod=ols

Use robust variance estimator to calculate the z

res=lmMod$res
stdev=mad(res)
median(abs(res-median(res)))*1.4826

[1] 1.512416

z=res/stdev

calculate weights

use psi.huber function

w=psi.huber(z)
plot(x,y)

plot(x,lmMod$res)

plot(x,lmMod$res,cex=w)

perform a weighted regression use lm with weights=w

lmMod=lm(y~x,weights=w)

plot results

plot(x,y)
abline(ols,lwd=2)

abline(mEst,col="red",lwd=2)
abline(lmMod,col="blue",lwd=2)

legend("topleft",legend=c("OLS","M-estimation","Our Impl"),lwd=2,col=c("black","red","blue"))

repeat this many times

lmMod=ols
for (k in 1:10)
{
######repeat this part several times until convergence
#use robust variance estimator to calculate the z
res=lmMod$res
stdev=mad(res)
median(abs(res-median(res)))*1.4826
z=res/stdev
#calculate weights
#use psi.huber function
w=psi.huber(z)
#perform a weighted regression use lm with weights=w
lmMod=lm(y~x,weights=w)
#plot results
plot(x,y)
abline(ols,lwd=2)
abline(mEst,col="red",lwd=2)
abline(lmMod,col="blue",lwd=2)
legend("topleft",legend=c("OLS","M-estimation","Our Impl"),lwd=2,col=c("black","red","blue"))
####################################
}