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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-robust
--
olivier.ren...@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
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))
}
_______________________________________________
R-SIG-Robust@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-robust
