# 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
median(abs(res-median(res)))*1.4826``````
``[1] 1.512416``
``z=res/stdev``

## 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"))``