# Table 8.7.1 data


wt <- c(140,155,153,150,180,160,155,145,163,150,140,170,180,175,155,150,160,145,190,
        228,150,180,165,145)

cs <- c(14.5,15.5,15.5,15.0,16.5,16.5,15.5,14.5,15.0,15.0,15.0,15.5,15.5,15.5,15.5,15.5,15.5,
15.0,16.0,16.5,15.0,15.0,15.0,15.0)

ss <- c(9.5, 9.5,10.5,10.5,11.0, 8.5, 8.5, 9.5, 10.0, 9.0, 8.5, 9.5,11.0,11.0,10.5, 8.5,
10.0, 9.0,12.0,13.0, 8.5, 8.5,11.0, 9.0)

WeightShoes <- data.frame(wt,cs,ss)
WeightShoes
rm(wt,cs,ss)                                  # remove unneeded objects

# -------------- R least squares -------------------------------------

summary(lm(wt ~ cs + ss, data = WeightShoes))

# ******************************************************
#Call:
#lm(formula = wt ~ cs + ss, data = WeightShoes)
#Residuals:
#    Min      1Q  Median      3Q     Max 
#-16.396  -5.548  -1.533   3.776  35.726 
#
#Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept) -156.572     67.886  -2.306  0.03138 *  
#cs            15.074      4.862   3.100  0.00542 ** 
#ss             8.793      2.202   3.993  0.00066 ***
#---
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#
#Residual standard error: 11.6 on 21 degrees of freedom
#Multiple R-Squared: 0.6862,     Adjusted R-squared: 0.6563 
#F-statistic: 22.96 on 2 and 21 DF,  p-value: 5.183e-06 
# ******************************************************

# ---------------------- M estimation ------------------------------

library(MASS)

summary(rlm(wt ~ cs + ss, data = WeightShoes, method="M", psi = psi.bisquare))

# ******************************************************
#Call: rlm(formula = wt ~ cs + ss, data = WeightShoes, psi = psi.bisquare, 
#    method = "M")
#Residuals:
#     Min       1Q   Median       3Q      Max 
#-13.3846  -2.9628   0.3464   5.8195  37.0452 
#
#Coefficients:
#            Value     Std. Error t value  
#(Intercept) -127.4560   47.4269    -2.6874
#cs            13.2687    3.3968     3.9062
#ss             8.3978    1.5383     5.4590
#
#Residual standard error: 6.522 on 21 degrees of freedom
# ******************************************************

WtShoesMBiSq <- rlm(wt ~ cs + ss, data = WeightShoes, method="M", psi = psi.bisquare)

# Here we look at the final weights used to obtain the M estimates, 
# to see how observations were downweighted

plot(1:24,WtShoesMBiSq$w)


summary(rlm(wt ~ cs + ss, data = WeightShoes, method="M", psi = psi.huber))

# ******************************************************
#Call: rlm(formula = wt ~ cs + ss, data = WeightShoes, psi = psi.huber, 
#    method = "M")
#Residuals:
#     Min       1Q   Median       3Q      Max 
#-15.1125  -4.0802  -0.3007   5.3522  37.0346 
#
#Coefficients:
#            Value     Std. Error t value  
#(Intercept) -145.5248   50.7043    -2.8701
#cs            14.1060    3.6316     3.8843
#ss             9.0471    1.6447     5.5009
#
#Residual standard error: 6.734 on 21 degrees of freedom
# ******************************************************

# ---------------------- MM estimation ------------------------------

summary(rlm(wt ~ cs + ss, data = WeightShoes, method="MM"))

# ******************************************************
#Call: rlm(formula = wt ~ cs + ss, data = WeightShoes, method = "MM")
#Residuals:
#      Min        1Q    Median        3Q       Max 
#-15.00946  -3.46271  -0.08995   6.39826  38.50924 
#
#Coefficients:
#            Value     Std. Error t value  
#(Intercept) -160.6751   49.9760    -3.2150
#cs            14.7154    3.5794     4.1111
#ss             9.5805    1.6210     5.9101
#
#Residual standard error: 8.188 on 21 degrees of freedom
# ******************************************************


WtShoesMM <- rlm(wt ~ cs + ss, data = WeightShoes, method="MM")

# Here we look at the final weights used to obtain the MM estimates, 
# to see how observations were downweighted

plot(1:24,WtShoesMM$w)



# ---------------------- Rank-based regression ------------------------------
# Must first load R code from: http://www.stat.wmich.edu/mckean/HMC/Rcode/AppendixB/ww.r
# and also the quantreg package must be installed and loaded


source("http://www.stat.wmich.edu/mckean/HMC/Rcode/AppendixB/ww.r")

library(quantreg)

# The first argument to the wwest function is the set of columns that are the independent variables,
# the second argument is the dependent variable, and the third column is the weight function, here 
# Wilcoxon weights are used
# In addition to printing the coefficient estimates, an option is given to produce residual plots


wwest(WeightShoes[,2:3],WeightShoes$wt,"WIL") 

# ******************************************************
#
#Wald Test of H0: BETA1=BETA2=0
#TS: 32.3743 PVAL: 0 
#
#Drop Test of H0: BETA1=BETA2=0
#TS: 19.4925 PVAL: 0 
#
#            EST      SE    TVAL   PVAL
#BETA0 -135.6250 55.9518 -2.4240 0.0245
#BETA1   13.3333  4.0062  3.3282 0.0032
#BETA2    9.1667  1.8143  5.0524 0.0001
#
#Would you like to see residual plots (y/n)? 
# ******************************************************



# This time Huber weights are used, which can adjust for high leverage points as well as 
# for outliers

wwest(WeightShoes[,2:3],WeightShoes$wt,"HBR") 

# ******************************************************
#Wald Test of H0: BETA1=BETA2=0
#TS: 20.6468 PVAL: 0 
#
#            EST      SE    TVAL   PVAL
#BETA0 -139.3182 57.5179 -2.4222 0.0246
#BETA1   13.6364  3.6593  3.7265 0.0012
#BETA2    9.0909  1.8771  4.8430 0.0001
#
#Would you like to see residual plots (y/n)?
# ******************************************************


#from http://www.stat.wmich.edu/cgi-bin/slab/
#Parameter Estimates and Standard Errors	
#  Wilcoxon R                                                   
#tau-hat=9.40186	                                                 
#	Estimate	SE	t-ratio	p-values	
#Intept	-135.942	55.0374	-2.47	0.0221665	
#Beta 1	13.3594		3.94175	3.3892	0.00276732	
#Beta 2	9.16015		1.78512	5.13138	4.39791e-05