Skip to content
Prev 417 / 523 Next

[RsR] Pseudo-R2 and AIC of robust regression

Dear Fadong,

Have you seen this:
https://stat.ethz.ch/pipermail/r-sig-robust/2010/000290.html ?
This has been on my/our to-do list for some time now. I did not have
time to look at it properly.

In the post I linked to above they provide some code that works with
lmRob from the robust package. I took the liberty to modify it so that
it works with lmrob objects from the robustbase package. See below.

Best,
Manuel


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
    ##   (default is for bisquare with 95% efficiency)

    ## Example:
    ## require(robustbase)
    ## m <- lmrob(Y ~ ., data=coleman)
    ## robR2w(m)

    if (attr(rob.obj, "class") != "lmrob")
        stop("This function works only on lmrob objects")
    pred <- fitted(rob.obj)
    resid <- resid(rob.obj)
    resp <- resid+pred
    wgt <- weights(rob.obj, "robustness")
    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))
}
On Tue, Nov 26, 2013 at 9:46 AM, Fadong Chen <fadongchen at gmail.com> wrote: