[RsR] Can I calculate the p-value of the robust regression at this way?

[Olivier Renaud]

Hi,
The output of lmRob DOES provide the robust version of the R-squared 
under "Multiple R-Squared". It does not directly provide a test that 
compares the model with the model with no covariates (equivalent of the 
"F test"), but it can be obtained with the anova function as shown here:

 > creat.lmRob = lmRob(original1 ~ 
approprie1+approprie2+creativite1+creativite2, data=creatif)
 > creat00.lmRob = lmRob(original1 ~ 1, data=creatif)
 > anova (creat00.lmRob, creat.lmRob)

This can also be done with lmrob().

Concerning the R-squared however, we have recently published a paper 
that shows that the robust R-squared provided by lmRob is biased, 
sometimes to a large extent. We provide a consistent and robust 
estimator of R-squared (robR2w.WithCorrection in the output below) and a 
version adjusted for the sample size (robR2w.AdjustedWithCorrection in 
the output below). The code and an example are provided below.

Olivier

ref: Renaud, O. & Victoria-Feser, M.-P. (2010). A robust coefficient of 
determination for regression. Journal of Statistical Planning and 
Inference, 140, 1852-1862. http://dx.doi.org/10.1016/j.jspi.2010.01.008

 > library(robust)
 > source("robR2w.r")

 > creat.lmRob = lmRob(original1 ~ 
approprie1+approprie2+creativite1+creativite2, data=creatif)

 > summary(creat.lmRob)

Call: lmRob(formula = original1 ~ approprie1 + approprie2 + creativite1 +
    creativite2, data = creatif)

Residuals:
        Min          1Q      Median          3Q         Max
-1.96149388 -0.34543174 -0.05500626  0.23168813  1.73781067

Coefficients:
            Value         Std. Error    t value       Pr(>|t|)    
(Intercept) -2.2744444543  1.2656121996 -1.7971100903  0.0784825453
approprie1   0.0914187017  0.1187669959  0.7697315322  0.4451541682
approprie2   0.1505246740  0.0934840433  1.6101643524  0.1137861853
creativite1  0.6270578015  0.1648705921  3.8033332300  0.0003963918
creativite2 -0.2384886952  0.1302348886 -1.8312197120  0.0731512215

Residual standard error: 0.459508 on 49 degrees of freedom
Multiple R-Squared: 0.183627

Test for Bias:
            statistic    p-value
M-estimate   6.527018 0.25825825
LS-estimate 10.368150 0.06545114

 > robR2w(creat.lmRob)
$robR2w.NoCorrection
[1] 0.3732078

$robR2w.WithCorrection
[1] 0.3302368

$robR2w.AdjustedWithCorrection
[1] 0.2604698



friendpine wrote:
> The output for the rlm(),lmRob(),lmrob() didn't give the R-squared and the
> p-value of the equation. There seems to be no function that can do these
> things.Can we calculate it as in the OLS(ordinary least squared) regression? The
> method in OLS use the SSE and SSR to calculate the f for F test. The p-value
> generated in the F test is the p-value for the OLS regression equation. Can we
> do it in the robust regression?
>
> _______________________________________________
> R-SIG-Robust using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-robust

-- 
Olivier.Renaud using unige.ch               http://www.unige.ch/fapse/mad/
Methodology & Data Analysis - Psychology Dept - University of Geneva
UniMail, Office 4164  -  40, Bd du Pont d'Arve   -  CH-1211 Geneva 4


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-robust/attachments/20100406/c2374888/attachment.html>

-------------- next part --------------

robR2w = function (rob.obj, correc=1.2076) {
  ## R2 in robust regression, see
  ## Renaud, O. & Victoria-Feser, M.-P. (2010). A robust coefficient of determination for regression. Journal of Statistical Planning and Inference, 140, 1852-–1862.
  ## rob.obj is an lmRob object. correc is the correction for consistancy. Call:
  ##
  ## library(robust)
  ## creat.lmRob = lmRob(original1 ~ approprie1+approprie2+creativite1+creativite2, data=creatif)
  ## summary(creat.lmRob)
  ## robR2w(creat.lmRob)
  
  ## Weights in robust regression
  wt.bisquare = function(u, c = 4.685) {
    U <- abs(u/c)
    w <- ((1. + U) * (1. - U))^2.
    w[U > 1.] <- 0.
    w
  }
  weight.rob=function(rob.obj){
    resid.rob=rob.obj$resid
    scale.rob=(rob.obj$scale)*rob.obj$df.residual/length(resid.rob)
    resid.rob= resid.rob/scale.rob
    weight=wt.bisquare(resid.rob)
  }

  if (attr(rob.obj, "class") !="lmRob")
    stop("This function works only on lmRob objects")
  pred = rob.obj$fitted.values
  resid = rob.obj$resid
  resp = resid+pred
  wgt = weight.rob(rob.obj)
  scale.rob = rob.obj$scale
  resp.mean = sum(wgt*resp)/sum(wgt)
  pred.mean = sum(wgt*pred)/sum(wgt)
  yMy = sum(wgt*(resp-resp.mean)^2)
  rMr = sum(wgt*resid^2)
  r2 = (yMy-rMr) / yMy
  r2correc= (yMy-rMr) / (yMy-rMr +rMr*correc)
  r2adjcor = 1-(1-r2correc) * (length(resid)-1) / (length(resid)-length(rob.obj$coefficients)-1)  
  return(list(robR2w.NoCorrection=r2, robR2w.WithCorrection=r2correc, robR2w.AdjustedWithCorrection=r2adjcor))
}