R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 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')

 Robust Regression of Systolic After  
> 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)
Call:
rlm(formula = spa ~ spb, data = bp)
Converged in 6 iterations

Coefficients:
(Intercept)         spb 
  3.5014401   0.8693671 

Degrees of freedom: 15 total; 13 residual
Scale estimate: 8.26 
> # Final Weights
> print(rlmout$w)
 [1] 0.7440184 0.7621694 1.0000000 0.7714647 1.0000000 0.9651177
 [7] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[13] 1.0000000 1.0000000 1.0000000
> 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')
Fake example with added outlier
> 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)
 [1] 0.1040508 0.9636992 1.0000000 1.0000000 1.0000000 1.0000000
 [7] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[13] 1.0000000 1.0000000 1.0000000
> # Now do bisquare
> (rlmout<-rlm(bp$spa~bp$spb,method="MM"))
Call:
rlm(formula = bp$spa ~ bp$spb, method = "MM")
Converged in 5 iterations

Coefficients:
(Intercept)      bp$spb 
 18.0687918   0.7825345 

Degrees of freedom: 15 total; 13 residual
Scale estimate: 9.37 
> rlmout$w
 [1] 0.0000000 0.7885590 0.9973501 0.8138413 0.9968162 0.8826195
 [7] 0.9725022 0.9994477 0.9572922 0.9859642 0.9891965 0.9445466
[13] 0.9630782 0.9753845 0.9412298
> 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))
       1 
2.935084 
> #This is not the right model.  Use Fieller.
> mod <- lm(measured ~ actual, data = arsenic)
> (res <- calibrate(mod, y0 = 3, interval = "inversion", level = 0.9) )
estimate    lower    upper 
2.931449 2.603537 3.258658 
> 
<\/pre>