pdf("g13.pdf");options(width=64)
#*************************************************/
# M estimation.  R seems to use Huber by default,*/
# while SAS seems to use bisquare.               */
#*************************************************/
cat('\n Robust Regression of Systolic After  \n')
library(MASS)# For rlm
#**************************************************************/
#* Data on blood presure change after captopril from          */
#* http://www.stat.ucla.edu/~rgould/datasets/bloodpressure.dat*/
#* Variables are                                              */
#* before systolic bp, after systolic bp, systolic bp change  */
#* before diastolic bp, after diastolic bp, diast. bp change  */
#* Example E from Cox and Snell, Applied Statistics: Princi-  */
#* ples and examples.  First line has column headings.  Fields*/
#* are separated by tabs.  SAS will give an error with this.  */
#**************************************************************/
bp<-as.data.frame(scan('bloodpressure.dat',skip=1,
  what=list(spb=0,spa=0,spd=0,dpb=0,dpa=0,dpd=0)))
rlmout<-rlm(spa~spb,data=bp)
library(quantreg)# For rq
rqout<-rq(spa~spb,data=bp)
print(rlmout)
# Final Weights
print(rlmout$w)
plot(bp$spb,bp$spa)
abline(rlmout,lty=2)
abline(rqout,lty=3)
abline(lm(spa~spb,data=bp),lty=4)
legend(min(bp$spb),max(bp$spa),
 legend=c("Huber","Median","OLS"),lty=2:4)
#Create a fake data set with an implausible value for the first response.
cat('Fake example with added outlier\n')
bp$spa[1]<-40
rqout<-rq(bp$spa~bp$spb);rlmout<-rlm(bp$spa~bp$spb);plot(bp$spb,bp$spa)
abline(rlmout,lty=2);abline(rqout,lty=3);
abline(lm(spa~spb,data=bp),lty=4)
legend(min(bp$spb),max(bp$spa),
 legend=c("Huber","Median","OLS"),lty=2:4)
# Final Weights
print(rlmout$w)
# Now do bisquare
(rlmout<-rlm(bp$spa~bp$spb,method="MM"))
rlmout$w
library(investr)# For calibrate
data(arsenic)
#Desired model is measured~actual.
#Estimate the actual corresponding to measured=3
#Swap explanatory and response variables.
predict(lm(actual~measured,data=arsenic),data.frame(measured=3))
#This is not the right model.  Use Fieller.
mod <- lm(measured ~ actual, data = arsenic)
(res <- calibrate(mod, y0 = 3, interval = "inversion", level = 0.9) )
<\/pre>