# An example of simple linear regression with an influential point

xslrinf <- c(3.01, 4.1, 2.1, 3.4, 2.6, 4.9, 1.1, 9.)  
yslrinf <- c(2, 1.5, 2.52, 1.75, 2.25, 1., 2.97, 8.)

# ----------------------Least squares ------------------------------

summary(lm(yslrinf ~ xslrinf))

# ******************************************************
#Call:
#lm(formula = yslrinf ~ xslrinf)
#
#Residuals:
#      Min        1Q    Median        3Q       Max 
#-2.476664 -0.930889  0.005382  1.109672  1.954804 
#
#Coefficients:
#            Estimate Std. Error t value Pr(>|t|)  
#(Intercept)   0.3027     1.1604   0.261   0.8029  
#xslrinf       0.6478     0.2638   2.456   0.0494 *
#---
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#
#Residual standard error: 1.684 on 6 degrees of freedom
#Multiple R-Squared: 0.5013,     Adjusted R-squared: 0.4181 
#F-statistic:  6.03 on 1 and 6 DF,  p-value: 0.04941 
# ******************************************************

plot(xslrinf,yslrinf)
abline(lm(yslrinf ~ xslrinf))

# ---------------------- M estimation ------------------------------
library(MASS)

summary(rlm(yslrinf ~ xslrinf, method="M"))
abline(rlm(yslrinf ~ xslrinf,method="M",lty=4)) 

# ******************************************************
#Call: rlm(formula = yslrinf ~ xslrinf, method = "M")
#Residuals:
#     Min       1Q   Median       3Q      Max 
#-2.50582 -0.95321 -0.01297  1.08606  1.94524 
#
#Coefficients:
#            Value  Std. Error t value
#(Intercept) 0.3066 1.3434     0.2282 
#xslrinf     0.6529 0.3054     2.1380 
#
#Residual standard error: 1.724 on 6 degrees of freedom
#
# ******************************************************


summary(rlm(yslrinf ~ xslrinf, method="MM"))
abline(rlm(yslrinf ~ xslrinf,method="MM",lty=5))

# ******************************************************
#Call: rlm(formula = yslrinf ~ xslrinf, method = "MM")
#Residuals:
#      Min        1Q    Median        3Q       Max 
#-0.056303 -0.028135  0.009656  0.043179  9.102320 
#
#Coefficients:
#            Value    Std. Error t value 
#(Intercept)   3.5723   0.0356   100.2711
#xslrinf      -0.5194   0.0081   -64.1359
#
#Residual standard error: 0.07371 on 6 degrees of freedom
# ******************************************************

# ---------------------- Rank-based regression ------------------------------
# Must first load R code from: http://www.stat.wmich.edu/mckean/HMC/Rcode/AppendixB/ww.r

library(quantreg)
wwest(x=xslrinf,y=yslrinf,bij="WIL")

# ******************************************************
#Wald Test of H0: BETA1=0
#TS: 0.7237 PVAL: 0.4276 
#
#Drop Test of H0: BETA1=0
#TS: 0.5793 PVAL: 0.4754 
#
#         EST     SE   TVAL   PVAL
#BETA0 0.3390 3.1081 0.1091 0.9167
#BETA1 0.6367 0.7485 0.8507 0.4276
# ******************************************************

abline(.3390,.6367,lty=2,lwd=2)
wwest(x=xslrinf,y=yslrinf,bij="HBR")

# ******************************************************
#Wald Test of H0: BETA1=0
#TS: 0.3611 PVAL: 0.5699 
#
#          EST     SE    TVAL   PVAL
#BETA0  3.5792 4.2289  0.8464 0.4298
#BETA1 -0.5184 0.8627 -0.6009 0.5699
# ******************************************************

abline(3.5792,-.5184,lty=6,lwd=3)