R/robR2w.R

In musto101/wilcox_R: What the package does (short line)

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