21

I estimated a robust linear model in R with MM weights using the rlm() in the MASS package. `R`` does not provide an $R^2$ value for the model, but I would like to have one if it is a meaningful quantity. I am also interested to know if there is any meaning in having an $R^2$ value that weighs the total and residual variance in the same way that observations were weighted in the robust regression. My general thinking is that, if, for the purposes of the regression, we are essentially with the weights giving some of the estimates less influence because they are outliers in some way, then maybe for the purpose of calculating $r^2$ we should also give those same estimates less influence?

I wrote two simple functions for the $R^2$ and the weighted $R^2$, they are below. I also included the results of running these functions for my model which is called HI9. EDIT: I found web page of Adelle Coster of UNSW that gives a formula for R2 that includes the weights vector in calculating the calculation of both SSe and SSt just as I did, and asked her for a more formal reference: http://web.maths.unsw.edu.au/~adelle/Garvan/Assays/GoodnessOfFit.html (still looking for help from Cross Validated on how to interpret this weighted $r^2$.)

#I used this function to calculate a basic r-squared from the robust linear model
r2 <- function(x){  
+ SSe <- sum((x$resid)^2);  
+ observed <- x$resid+x$fitted;  
+ SSt <- sum((observed-mean(observed))^2);  
+ value <- 1-SSe/SSt;  
+ return(value);  
+ }  
r2(HI9)  
[1] 0.2061147

#I used this function to calculate a weighted r-squared from the robust linear model
> r2ww <- function(x){
+ SSe <- sum((x$w*x$resid)^2); #the residual sum of squares is weighted
+ observed <- x$resid+x$fitted;
+ SSt <- sum((x$w*(observed-mean(observed)))^2); #the total sum of squares is weighted      
+ value <- 1-SSe/SSt;
+ return(value);
+ }
 > r2ww(HI9)
[1] 0.7716264

Thanks to anyone who spends time answering this. Please accept my apologies if there is already some very good reference on this which I missed, or if my code above is hard to read (I am not a code guy).

CraigMilligan
  • 571
  • 1
  • 4
  • 9
  • put the weights inside lm() and take the r-squared from there (why re-invent the wheel?) – user603 Jan 30 '14 at 17:35
  • 1
    thanks for the tip on a way to do what I did more efficiently. can anyone comment on the meaning of the weighted r-squared that I described/proposed? – CraigMilligan Jan 30 '14 at 19:46
  • @user603: How would you actually go about putting the weights inside lm()? – histelheim Oct 21 '15 at 13:02
  • Just for a compliment, the weighted least square fitted in R is by minimizing sum(w * e^2), where e is the residual. So for you computation code, all weight **w** should be taken a square root. – Yuanhao Lai Dec 10 '16 at 06:49
  • I want to emphasize that we haven't to take a weighted mean, at least I believe that because a program I have written give a r-squared close to 1 with : classic r-squared weighted r-squared but NOT with weighted r-squared where the mean is weighted too, I find -6 it is counter-intuitive even for me, but I believe the experience though – pierre May 28 '17 at 11:33

2 Answers2

26

The following answer is based on: (1) my interpretation of Willett and Singer (1988) Another Cautionary Note about R-squared: It's use in weighted least squates regression analysis. The American Statistician. 42(3). pp236-238, and (2) the premise that robust linear regression is essentially weighted least squares regression with the weights estimated by an iterative process.

The formula I gave in the question for r2w needs a small correction to correspond to equation 4 in Willet and Singer (1988) for r2wls: the SSt calculation should also use a weighted mean: 

the correction is SSt <- sum((x$w*observed-mean(x$w*observed))^2)].

What is the meaning of this (corrected) weighted r-squared? Willett and Singer interpret it as: "the coefficient of determination in the transformed [weighted] dataset. It is a measure of the proportion of the variation in weighted Y that can be accounted for by weighted X, and is the quantity that is output as R2 by the major statistical computer packages when a WLS regression is performed".

Is it meaningful as a measure of goodness of fit? This depends on how it is presented and interpreted. Willett and Singer caution that it is typically quite a bit higher than the r-squared obtained in ordinary least squares regression, and the high value encourages prominent display... but this display may be deceptive IF it is interpreted in the conventional sense of r-squared (as the proportion of unweighted variation explained by a model). Willett and Singer propose that a less 'deceptive' alternative is pseudoR2wls (their equation 7), which is equivalent to my function r2 in the original question. In general, Willett and Singer also caution that it is not good to rely on any r2 (even their pseudor2wls) as a sole measure of goodness of fit. Despite these cautions, the whole premise of robust regression is that some cases are judged 'not as good' and don't count as much in the model fitting, and it may be good to reflect this in part of the model assessment process. The weighted r-squared described, can be one good measure of goodness of fit - as long as the correct interpretation is clearly given in the presentation and it is not relied on as the sole assessment of goodness of fit.

CraigMilligan
  • 571
  • 1
  • 4
  • 9
2

@CraigMilligan. Shouldn't:

  • the weight be outside of the squared parenthesis
  • the weighted mean be calculated as sum(x$w*observed)/sum(x$w) for which we can also use weighted.mean(observed,x$w)

Something like this:

r2ww <- function(x){
  SSe <- sum(x$w*(x$resid)^2)
  observed <- x$resid+x$fitted
  SSt <- sum(x$w*(observed-weighted.mean(observed,x$w))^2)
  value <- 1-SSe/SSt;
  return(value);
}
Julien Massardier
  • 121
  • 4